A Multimessenger Picture of the Flaring Blazar TXS 0506+056: Implications for High-energy Neutrino Emission and Cosmic-Ray Acceleration

IceCube-170922A was an EHE neutrino event (GCN/AMON NOTICE 2016) identified and distributed by the IceCube Observatory via AMON and GCN within of its detection at 20:54:30 UT on 2017 September 22 (GCN/AMON NOTICE IceCube-170922A 2017). A refined localization was reported four hours later (Kopper & Blaufuss 2017): R.A. =  degree, decl. =  degree (J2000; 90% containment ellipse). The maximum likelihood neutrino position is R.A. 05°09^m08784, decl. +05°45′1332 (J2000); see Figure 1 for an illustration of the initial and final localizations.

Zoom In Zoom Out Reset image size

EHE neutrino event reports include the neutrino arrival time, direction (R.A. and decl.), angular error (r₅₀ for 50% containment; r₉₀ for 90% containment), revision number, an estimate of the deposited charge, an estimate of the neutrino energy, the parameter signalness, and an estimate of the probability that the event was due to an astrophysical—rather than atmospheric—neutrino (Aartsen et al. 2017a). Real-time identification, localization, and reporting of IceCube HE neutrinos is enabled by software that has been in place at the South Pole since 2016 April (Aartsen et al. 2017a).

IceCube-170922A triggered the Neil Gehrels Swift Observatory in automated fashion via AMON cyberinfrastructure, resulting in rapid-response mosaic-type follow-up observations, covering a roughly circular region of sky centered on the prompt localization in a 19-point tiling that began 3.25 hr after the neutrino detection. This initial epoch of Swift observations spanned 22.5 hr and accumulated ≈800 s exposure per pointing. The mosaic tiling yielded coverage of a region with radius ≈08 centered on R.A. 05^h09^m08784, decl. +05°45′1332 (J2000), amounting to a sky area of 2.1 deg². XRT data were analyzed automatically, as data were received at the University of Leicester, via the reduction routines of Evans et al. (2009, 2014). Nine X-ray sources were detected in the covered region down to a typical achieved depth of erg cm⁻² s⁻¹(0.3–10.0 keV). Figure 1 shows the exposure map for the 19-point tiling pattern, along with the nine detected X-ray sources. All detected sources were identified as counterparts to known and cataloged stars, X-ray sources, or radio sources (Keivani et al. 2017); fluxes of these X-ray sources were consistent with previously measured values.

Notably, Source 2 from these observations (marked as X2 on Figure 1), located 46 from the center of the neutrino localization, was identified by us as the likely X-ray counterpart to QSO J0509+0541, also known as TXS 0506+056. This was the first report to connect TXS 0506+056 to IceCube-170922A (Keivani et al. 2017).

Following the Fermi report that TXS 0506+056 was in a rare GeV-flaring state (Tanaka et al. 2017), we commenced a Swift monitoring campaign on September 27 (Evans et al. 2017). Swift monitored TXS 0506+056 for 36 epochs by November 30 with 53.7 ks total exposure time (Table 1).

Table 1.  Swift-XRT Monitoring of TXS 0506+056

Epoch	Exposure [ks]	Photon Index
(1)	(2)	(3)	(4)	(5)
58019.4700 ± 0.4630	0.8		65.8 ± 10.1
58023.8535 ± 0.0660	4.9	2.43 ± 0.12	121.0 ± 5.3
58026.2273 ± 0.0391	2.0	2.30 ± 0.33	66.2 ± 7.8
58028.6419 ± 0.0114	2.0	2.73 ± 0.20	117.0 ± 8.2
58029.6745 ± 0.1036	1.1		182.0 ± 14.1
58030.7369 ± 0.0385	1.2		186.0 ± 16.1
58031.7944 ± 0.1025	2.3	2.64 ± 0.13	255.0 ± 11.7
58032.8985 ± 0.3428	2.1		90.1 ± 7.1
58034.4478 ± 0.0337	1.9	2.53 ± 0.21	108.0 ± 8.1
58040.9452 ± 0.0010	0.2
58042.7684 ± 0.1686	2.2	2.00 ± 0.27	61.9 ± 5.8
58044.1300 ± 0.0696	1.8	2.10 ± 0.31	52.6 ± 6.0
58047.2130 ± 0.6360	2.3	2.11 ± 0.28	49.6 ± 5.1
58050.7286 ± 0.0446	2.9		46.6 ± 4.4
58053.3162 ± 0.0335	1.0		54.4 ± 8.2
58059.6620 ± 0.0729	3.3	2.22 ± 0.21	60.2 ± 4.7
58065.6390 ± 0.3371	3.1	2.30 ± 0.20	83.4 ± 6.2
58068.3642 ± 0.0742	3.0	2.38 ± 0.24	55.7 ± 4.9
58069.5359 ± 0.5032	1.8	2.23 ± 0.24	83.3 ± 7.5
58071.1551 ± 0.1299	2.9	2.15 ± 0.19	81.1 ± 5.9
58072.0949 ± 0.0039	0.7		48.5 ± 9.1
58073.0911 ± 0.0042	0.7		55.0 ± 9.7
58074.2484 ± 0.0962	3.0		54.4 ± 4.8
58075.0831 ± 0.0046	0.8	2.51 ± 0.45	63.1 ± 9.6
58075.6655 ± 0.0057	1.0	1.94 ± 0.41	60.6 ± 9.4
58076.0793 ± 0.0049	0.8		56.1 ± 9.0
58077.0756 ± 0.0054	0.9		61.4 ± 8.9
58078.0717 ± 0.0056	1.0		40.5 ± 7.1
58079.0675 ± 0.0057	1.0		49.3 ± 8.6
58080.1320 ± 0.0057	1.0		47.4 ± 9.4
58081.1279 ± 0.0057	1.0		45.8 ± 9.7
58082.0551 ± 0.0057	1.0
58083.1195 ± 0.0057	1.0	2.45 ± 0.56	44.4 ± 8.0
58084.0499 ± 0.0058	1.0		39.0 ± 6.9
58086.1089 ± 0.0059	1.0		45.4 ± 9.8
58087.1560 ± 0.0053	0.9		66.9 ± 9.2

Note. and indicate count rate and energy flux, in units of 10⁻³ ct s⁻¹and 10⁻¹² erg cm⁻² s⁻¹, respectively. Uncertainties are quoted at 90% confidence.

Download table as:  ASCIITypeset image

To characterize the X-ray flux and spectral variability of TXS 0506+056, we performed a power-law fit to each individual Swift-XRT observation (Table 1), as well as to the summed spectrum from all listed epochs, using XSPEC (Arnaud 1996). The observation on October 14 is excluded from the spectral analysis due to low exposure time. The summed spectrum is adequately fit with a single power-law spectral model having the Galactic column density cm⁻², resulting in a photon index and mean flux of erg cm⁻² s⁻¹ (0.3–10.0 keV).

We note that this source has been observed on multiple previous occasions with Swift-XRT, with results published in the 1SXPS catalog (Evans et al. 2014). In past observations, TXS 0506+056 exhibited a typical flux of erg cm⁻² s⁻¹, with one observation at erg cm⁻² s⁻¹(0.3–10.0 keV). The source was thus in an active X-ray flaring state by comparison to historical X-ray measurements (Figure 2, upper left). Photon indices and X-ray flux measurements for each epoch are provided in Table 1, and the variations in photon index are shown in Figure 2 (upper right). The Swift-XRT light curve is shown in Figure 2 (upper left).

Zoom In Zoom Out Reset image size

To further characterize the HE emissions of the source, we requested two observations with the NuSTAR hard X-ray (3.0–100 keV) mission (Harrison et al. 2013).

On 2017 September 29 (02:23–17:48 UTC) NuSTAR carried out a Target of Opportunity observation of TXS 0506+056 (Fox et al. 2017). The full science observation was retrieved from the NuSTAR public archive (ObsID 90301618002). Data were processed within the HEAsoft (Arnaud 1996) software environment using the nupipeline tool with the setting SAAMODE=strict. This yielded exposures of 23.9 ks (24.5 ks) and count rates of 21.3 ct ks⁻¹(20.8 ct ks⁻¹) in the A (B) units, respectively. Level 3 data products for the source were then extracted using the nuproducts tool. Within XSPEC, spectral data from both units were fit to a single power-law model with cm⁻², resulting in a photon index of and a flux of erg cm⁻² s⁻¹ (3.0–79.0 keV).

On 2017 October 19 (10:26 to 21:21 UTC) NuSTAR performed a second Target of Opportunity observation of TXS 0506+056 (ObsID 90301618004). Executing the same reduction as for the first observation yielded exposures of 19.7 ks (19.7 ks) and count rates of 20.2 ct ks⁻¹(19.6 ct ks⁻¹) in the A (B) units at this second epoch. Performing the same spectral fit as for the first observation, we obtain a photon index and flux of erg cm⁻² s⁻¹ (3.0–79.0 keV), consistent with results of the first observation from 20 days earlier.

In order to obtain the energy spectrum for a wider X-ray band (0.2–100 keV), we simultaneously fit data from individual XRT observations and NuSTAR for two main epochs: [−15 days, +15 days] (Ep. 1) and [+15 days, +45 days] (Ep. 2) relative to the neutrino detection. The two epochs include eight and seven XRT observations, respectively, and one NuSTAR observation each. Since the source spectrum over the NuSTAR bandpass does not change from Ep. 1 to Ep. 2, we fit all individual XRT observations together with both NuSTAR observations with a sum of two power laws model, including Galactic absorption frozen at cm⁻², and quote the soft component best-fit parameters when reporting XRT results. A Markov Chain Monte Carlo algorithm is then employed for each fit to provide the 90% confidence levels. We generate 1000 realizations of the spectra from each XRT observation and add all Ep. 1 and Ep. 2 realizations together in order to find the 90% confidence intervals on the flux density versus energy. This joint analysis results in best-fit photon indices () of and . The results are displayed in Figure 3 and used for subsequent SED modeling (see Section 3).

Zoom In Zoom Out Reset image size

The Swift-UVOT (Roming et al. 2005) also participated in the rapid-response follow-up observations of the IceCube-170922A and the subsequent monitoring of the flaring blazar TXS 0506+056 as described in Section 2.2. The u filter was used during all 19 pointings used to tile the region around IceCube-170922A. TXS 0506+056 was readily detected during this initial survey. All six UVOT lenticular filters (v, b, u, uvw1, uvm2, uvw2) were used during the subsequent observations monitoring TXS 0506+056.

UVOT data were analyzed using the standard tool uvotsource of HEAsoft (v6.22.1) and the latest updates to the UVOT CALDB files. uvotsource does aperture photometry (Breeveld et al. 2010, 2011) using user-specified source and background regions. Because TXS 0506+056 is a bright source, a 5″ aperture was used for the source region. A nearby source-free region was used as the background region instead of the usual concentric ring centered on the target to avoid contamination from a read-out streak. The read-out streak is produced by photons from a very bright source near the edge of the UVOT image that arrive during the brief time in which the image frame is transferred for read out of the detector’s CCD (Page et al. 2013). The position of the read-out streak changes from observation to observation as the orientation on the sky of the UVOT image changes. Consequently, a different nearby background region was used for the observations centered near MJD 58065 (Swift sequence 00083368018). Data from observations near MJD 58028.6 (sequence 00083368003) and MJD 58031.8 (sequence 00083368006) are not used because the source region is within the read-out streak.

Table 2 reports the times, exposures, and magnitudes for all the observations. The source varies over a range of at least 0.5 mag in all six filters.

Due to the UVOT’s blue response, extending as far as Å, it is necessary to apply an appropriate extinction correction in order to interpret UVOT observations appropriately and derive physical constraints on the SED of TXS 0506+056. The line-of-sight extinction to TXS 0506+056, as provided by the NASA/IPAC Infrared Science Archive¹⁹ is A_V = 0.286 mag according to all-sky dust maps (Schlegel et al. 1998; Schlafly & Finkbeiner 2011). Although this estimate nominally has little uncertainty ( mag), we note that (1) this quoted value has been corrected from the original published value by 14% (Schlafly & Finkbeiner 2011); and (2) a subsequent recalibration with a similar approach, using PanSTARRS-1 rather than Sloan Digital Sky Survey data (Green et al. 2015), leads to a different extinction estimate, A_V = 0.254 mag. We therefore take a conservative approach and adopt a line-of-sight extinction of mag, with Gaussian form, as our prior for UVOT SED analysis. To calculate extinction as a function of wavelength, we use the Fitzpatrick (1999) extinction law with R_V = 3.1.

On 1 occasion during Ep. 1, on 3 occasions during Ep. 2, and on 20 subsequent occasions before November 30, UVOT observations were carried out using all six UV/optical filters in rapid sequence. We use these observations to characterize the UV/optical SED and UVOT spectral index (Figure 2, lower right). Additional single-filter observations are then used (in tandem with the six-filter epochs) to estimate source variability (Figure 2, lower left).

To fit the SED for each six-filter UVOT observation, we perform a χ² minimization of the predicted versus observed count rates in each filter using a model with three parameters: source flux density F_ν at Hz, UVOT spectral index β, and extinction A_V. Extinction values away from our adopted value of A_V = 0.286 mag are penalized according to our Gaussian extinction uncertainty. The source spectrum is integrated across each filter bandpass using the filter transmission function²⁰ and Fitzpatrick (1999) extinction curve. A fit is considered successful if the total χ² across the six filters gives a p-value . Only one of the six-filter observations (at δt = 64.5 days after IceCube-170922A) violates this constraint, and thus receives treatment as a broken power law across the UVOT filters (forcing a single power-law fit gives a spectral index β = 1.66 ± 0.06 at this epoch, where we quote uncertainties from analysis even though the fit is acknowledged not to be satisfactory). We note that curvature of the source spectrum (in particular, steepening/softening toward the UV) is also observed in the X-shooter data (Section 2.7). Acceptable fits yield a best-fit value for the UVOT spectral index and uncertainty via analysis. These values and uncertainties are reported in Table 3 and plotted in Figure 2 (lower right).

Table 3.  Swift-UVOT Photon Indices for TXS 0506+056

Epoch	Photon Index
58023.7632
58050.6452
58053.2192
58059.5782
58065.5562
58068.2802
58069.9502
58070.9472
58072.0112
58073.0082
58074.1632
58074.9992
58075.5822
58075.9962
58076.9922
58077.9882
58078.9842
58080.0482
58081.0442
58081.9712
58083.9662
58086.0252
58083.0362 (UV)
58083.0362 (Vis)
58087.0722

Note. UVOT photon indices from power-law SED fitting for epochs with data in all six filters. Data for the next-to-last epoch are fitted with a power law to the three UV and three visible filters separately; a forced single power-law fit to this epoch yields a photon index of 2.66, see the text for details.

Download table as:  ASCIITypeset image

We determine UVOT SED bands, the range of fluxes allowed at 90% confidence for each photon energy, by drawing 1000 samples according to the χ² probability function, and generating a power-law spectrum across the UVOT bandpass for each. The allowed range at each energy in the SED is defined as the minimum-width range encompassing 90% of these spectra. For Ep. 2, we draw 1000 samples from each of the three six-filter epochs and combine these before finding the 90%-confidence range; results for this epoch thus account for source flux variability. For Ep. 1, we use only the single six-filter observation to characterize the SED; we estimate accounting for flux variability over Ep. 1 would expand this band by 12%; however, we do not make this correction in our analysis.

We determine UVOT fluxes for single-epoch observations using the β measurement from the temporally proximate six-filter observation, adjusting the flux density at 10¹⁵ Hz to achieve agreement with the observed count rate in the relevant filter. Quoted uncertainties for these flux estimates combine the Poisson count rate uncertainty with the uncertainty in flux for the adopted SED model, in quadrature. Flux values and uncertainties are reported in Table 4 and plotted in Figure 2 (lower left).

Table 4.  Swift-UVOT Extinction-corrected Fluxes at ν = 10¹⁵ Hz

Epoch
58019.4709	22.50 ± 3.29
58023.8470	34.80 ± 3.13
58026.2284	12.80 ± 3.13
58029.6756	15.80 ± 3.13
58030.7380	15.00 ± 3.29
58032.5593	19.30 ± 3.29
58033.1319	11.00 ± 3.13
58034.4489	13.80 ± 3.29
58040.9463	16.00 ± 1.97
58042.7696	15.70 ± 1.97
58044.1311	19.60 ± 1.97
58046.5814	10.90 ± 1.97
58047.8132	13.70 ± 1.97
58050.7285	19.80 ± 1.81
58053.3030	22.40 ± 2.14
58059.6613	26.70 ± 2.47
58065.6395	24.70 ± 2.30
58068.3633	18.60 ± 1.81
58070.0338	19.40 ± 1.81
58071.0301	21.30 ± 2.14
58072.0950	22.70 ± 2.30
58073.0911	22.30 ± 2.14
58074.2466	24.70 ± 2.30
58075.0830	25.60 ± 2.47
58075.6651	23.70 ± 2.30
58076.0792	23.20 ± 2.14
58077.0751	23.80 ± 2.30
58078.0713	26.10 ± 2.47
58079.0670	24.70 ± 2.47
58080.1316	25.10 ± 2.47
58081.1274	24.60 ± 2.30
58082.0546	24.30 ± 2.30
58083.1191	21.40 ± 2.14
58084.0496	27.30 ± 2.63
58086.1085	30.10 ± 2.96
58087.1557	21.70 ± 2.80

Note. is the flux at ν = 10¹⁵ Hz in units of 10⁻¹² erg cm⁻² s⁻¹, derived by SED fitting for six-filter epochs, and adjusted to ν = 10¹⁵ Hz using the nearest best-fit SED index for single-filter epochs.

Download table as:  ASCIITypeset image

Medium-resolution spectroscopy of TXS 0506+056 was obtained with the X-shooter spectrograph (Vernet et al. 2011) mounted on the Very Large Telescope UT2 at ESO Paranal Observatory on 2017 October 1. The three arms of X-shooter (UV: UVB, optical: VIS and near-infrared: NIR) were used with slit widths of 10, 09, and 09, respectively. These data provide simultaneous 300–2480 nm spectral coverage with average spectral resolutions of 4290, 7410, and 5410, respectively, in each arm. Observing conditions were good, with a clear sky, seeing of ∼08, and an airmass ranging from 1.2 to 1.3. Individual exposure times are 72, 139, and 54 s for the UBV, VIS, and NIR arms, respectively, and sum to total integration times of 1152, 2224, and 864 s. Standard ABBA nodding observing mode was used to allow for an effective background subtraction.

Data were reduced using the ESO X-shooter pipeline (Goldoni et al. 2006; Modigliani et al. 2010) (v.2.9.3) in the Reflex environment (Freudling et al. 2013), producing a background-subtracted, wavelength-calibrated spectrum. The extracted 1D spectrum was flux calibrated with the X-shooter pipeline using a response function produced by observing the HST white dwarf standard GD71 (R.A. 05^h52^m2786, decl. +15°53′138, J2000) just after the observation of TXS 0506+056. To correct results for slit losses, the final spectrum was rescaled to match the broadband B, V, and R magnitudes obtained on 2017 October 29 using the 1 m Kapteyn Telescope at La Palma (Keel & Santander 2017). Overall, the flux calibration is expected to be accurate to 10% in the UVB arm and 15% in both the VIS and NIR arms based on the seeing conditions at the observing time. The telluric absorption lines were removed by using the Molecfit software through a fit of synthetic transmission spectra calculated by a radiative transfer code (Kausch et al. 2015; Smette et al. 2015). Finally, we corrected the spectra for Galactic extinction using the extinction law of Fitzpatrick (1999) with a total extinction at the V filter band A_V = 0.286 mag and a selective-to-total extinction ratio equal to the Galactic average value R_V = 3.1. For the 90%-confidence band plotted in Figure 3, we allowed the extinction to vary by δA_V = ± 0.054 mag, according to our adopted uncertainty (Section 2.6).

Two Galactic interstellar absorption features are observed in the reduced spectrum: Ca K and H absorption lines (at 3933.7 Å and 3968. Å respectively), and the Na ID doublet at 5892.5 Å. No other emission or absorption line is observed. Overall, the spectrum is consistent with the spectrum of a nonthermally dominated blazar and confirms the source to exhibit a bluer spectrum than that published by Halpern et al. (2003), as previously mentioned by Steele (2017).

Both the X-shooter and Swift-UVOT data clearly show that the synchrotron peak is below Hz. This indicates that TXS 0506+056 is an intermediate synchrotron peaked (ISP) or low synchrotron peaked blazar.

The nondetection of Lyα absorption in the X-shooter spectrum provides a rough upper limit on the redshift of TXS 0506+056, , which is compatible with the redshift measurement (z = 0.3365 ± 0.0010) of Paiano et al. (2018).

The Fermi-LAT is a pair conversion telescope sensitive to γ rays in the 20 MeV to >300 GeV (Atwood et al. 2009). In this section, we analyze photons detected by the LAT during our defined Ep. 1 (±15 days from the neutrino detection) and Ep. 2 (15 to 45 days after the neutrino detection). Analysis was performed using version v10r0p5 of the Fermi Science Tools.²¹

Photons with energies between 100 MeV and 300 GeV that were detected within a radius of 15^° from the location of TXS 0506+056 were selected for the analysis, while photons with a zenith angle were discarded to reduce contamination from the Earth’s albedo.

The contribution from isotropic and Galactic diffuse backgrounds was modeled using the parameterization provided in the files iso_P8R2_SOURCE_V6_v06.txt and gll_iem_v06.fits, respectively. Sources in the 3FGL catalog within a radius of 15^° from the source position were included in the model, with spectral parameters fixed to their catalog values, while spectral parameters for sources within 3^° were allowed to vary freely during the fit. The TXS 0506+056 spectral fit was performed with a binned likelihood method using the P8R2_SOURCE_V6 instrument response function. A power-law fit to the spectrum gives a photon index of , consistent with the 3FGL value of 2.04 ± 0.03, and a flux normalization of cm⁻² s⁻¹ MeV⁻¹ at an energy of 1.44 GeV, about a factor of four higher than the 3FGL value of in the same units. The spectral fit was repeated in seven independent energy bins with equal logarithmic spacing in the 100 MeV–300 GeV range to be incorporated in the modeling of the SED. Best-fit flux values and 90% uncertainties, shown in Figure 3 and subsequent figures, are reported for spectral bins with a test statistic (TS) value larger than 9, which corresponds to an excess of ∼3σ. Flux upper limits at the 95% confidence level are quoted otherwise.

We assume that protons and electrons are coaccelerated by some mechanism, whose details lie outside the immediate scope of this work and are subsequently injected isotropically in a spherical region containing a tangled magnetic field. The particle interactions with the magnetic field and with secondary particles leads to the development of a system with five stable particle populations in steady state: protons, which lose energy by synchrotron radiation, Bethe–Heitler pair production, and photomeson production processes; electrons and positrons, which lose energy by synchrotron radiation and IC scattering; photons, which gain and lose energy in a variety of ways; neutrons, which can escape almost unimpeded from the source region, with a certain probability of photohadronic interactions; and neutrinos, which escape without any attenuation. The interplay of the processes governing the evolution of the energy distributions of those five populations is formulated with a set of time-dependent kinetic equations. Through them, energy is conserved in a self-consistent manner, since all the energy gained by a particle type has to come from an equal amount of energy lost by another particle type.

To simultaneously solve the coupled kinetic equations for all particle types, we use the time-dependent code described in Dimitrakoudis et al. (2012). Photomeson production processes are modeled using the results of the Monte Carlo event generator SOPHIA (Mücke et al. 2000), while the Bethe–Heitler pair production is similarly modeled with the Monte Carlo results of Protheroe & Johnson (1996) and Mastichiadis et al. (2005). The only particles that are not modeled with kinetic equations are muons, pions, and kaons (Dimitrakoudis et al. 2014; Petropoulou et al. 2014); their energy losses can be safely ignored for the parameter values relevant to this study (see also Murase et al. 2014 for numerical calculations where the kinetic equations for these particles are explicitly solved).

The parameters that describe the source (i.e., Doppler factor δ, comoving magnetic field strength B′, and comoving blob size R′) as well as those of accelerated (i.e., primary) particle distributions can often be constrained by multiwavelength data (Takahashi et al. 1996; Mastichiadis & Kirk 1997; Romanova & Lovelace 1997; Kirk et al. 1998; Li & Kusunose 2000); a complete list of model parameters is provided in Table 5.

Table 5.  Physical Parameters (Description, Symbol, and Units) Used in Blazar Leptonic and Hadronic Modeling

Parameter	Symbol	Unit [in cgs]
Doppler factor	δ	n/a
Magnetic field strength	B′	G
Blob radius	R′	cm
Electron injection luminosity		erg s−1
Minimum electron Lorentz factor		n/a
Maximum electron Lorentz factor		n/a
Break electron Lorentz factor		n/a
Power-law electron index below the break		n/a
Power-law electron index above the break		n/a
Proton injection luminosity		erg s−1
Minimum proton Lorentz factor		n/a
Maximum proton Lorentz factor		n/a
Power-law proton index	sp	n/a
Energy density of external radiation		erg cm−3
Effective temperature of blackbody external radiation		K
Photon index of power-law external radiation	α	n/a
Minimum photon energy of power-law external radiation		keV
Maximum photon energy of power-law external radiation		keV

Note. Primed quantities are measured in the jet comoving frame. Parameters describing the relativistic particle distributions refer to their properties at injection.

Download table as:  ASCIITypeset image

We search for models that adequately describe the multiwavelength data (i.e., the model curve passes through most of the instrument-specific SED bands in Figure 3). We begin the parameter space search using values that we obtain analytically from expressions that relate observables to model parameters, as described in Murase et al. (2012) and Petropoulou et al. (2015). As we do not perform a statistical fit to the whole multiwavelength data in the strict sense (i.e., by maximizing a likelihood function), no uncertainty ranges for the model parameters can be formally computed. However, thanks to the detailed quasi-simultaneous X-ray data obtained in this work, we can place limits on the HE neutrino flux without depending on details of the model uncertainties (see subsequent sections). Quantitative upper limits on the proton and neutrino luminosities are placed by the requirement that the EM cascade does not overproduce emission in the X-ray regime (0.1–100 keV), where the source SED exhibits a prominent dip.

The resulting upper limits are quite robust, as they depend on the energy flux ratio of the EM and neutrino components—determined by well-known particle interactions—as well as the properties of EM cascades, which reliably yield a flat, broadband component by redistributing energy from high to low energies.

In the leptonic scenario, the blazar’s SED (optical to γ-rays) is explained by synchrotron and IC processes of accelerated (primary) electrons (Maraschi et al. 1992; Dermer & Schlickeiser 1993; Sikora et al. 1994). The radiation produced by relativistic protons in the source, which are necessary for the production of HE neutrinos, may not be directly observed due to the two-photon annihilation process and subsequent EM cascades inside the source. We coin these hybrid scenarios “LMs,” which stand for LMs, in reference to the leptonic origin of the γ rays. Significant intra-source γ-ray attenuation at sufficiently high energies and the associated EM cascade is unavoidable in single-zone models, because target photons responsible for photohadronic interactions hinder HE γ rays from leaving the source. This implies that a source with efficient HE neutrino production can be γ-ray dark and may even be regarded as a hidden cosmic-ray accelerator (Murase et al. 2016).

The photomeson production process also leads to the production of γ-ray photons from neutral pion decay. Moreover, the decay of charged pions leads to the production of secondary electrons and positrons, which also emit HE photons via synchrotron and IC processes. The HE photons can be attenuated by low-energy photons in the source, while enhancing the number of secondary electron–positron pairs. The total absorbed photon luminosity will eventually be redistributed at lower photon energies through the development of an EM cascade (Mannheim et al. 1991; Mannheim 1993).

The IC emission of primary electrons explains the HE peak of the SED, and the emission from the EM cascade should be subdominant. We can therefore set an upper limit on the power of the cosmic-ray proton component by requiring that any proton-induced emission does not fill in the dip (in hard X-rays for ISPs, as here) between the two peaks of the SED. In turn, this translates into an upper limit on the blazar’s neutrino flux.

We first derive the maximum neutrino flux expected in the leptonic scenario by assuming that the proton distribution is a power law with a proton index of s_p = 2, extending from to . From the X-ray and γ-ray light curves, we infer a variability timescale of . Our choice of cm is broadly consistent with the size inferred from the variability, namely cm. We also consider an arbitrary external photon field with a blackbody-like energy distribution that can be described by only two free parameters: its characteristic temperature and energy density , as measured in the comoving frame of the source. We also neglect any angular dependencies of the external radiation field, which is assumed to be isotropic in the rest frame of the supermassive black hole. Such an additional photon field has also been shown to be necessary in the leptonic SED modeling of other ISP blazars (Boettcher et al. 2013). Furthermore, inclusion of external photon fields has been shown to significantly enhance the efficiency of HE neutrino production (Atoyan & Dermer 2001; Dermer et al. 2014; Murase et al. 2014).

The respective photon spectrum and the maximum predicted neutrino flux for this parameter set (LMBB2b model) are presented in Figure 4 (solid curves) and the parameter values are summarized in Tables 6 and 7. We find that the X-ray flux in the NuSTAR energy band is dominated by the SSC emission of the accelerated electrons, whereas the γ-ray emission is explained by the IC scattering of the fiducial external photon field by the same electron population. The steepening of the γ-ray spectrum at is due to the Klein–Nishina cross section. Intriguingly, because of the steep Swift-XRT spectrum and the low synchrotron peak-frequency revealed by our X-shooter data, the HE peak of the SED cannot be explained by the SSC emission alone. In addition, any attempt to describe the emission from a more compact region ( cm) fails because of the emergence of the SSC component, which has a different photon index than the observed one in the NuSTAR band. This also demonstrates the importance of the detailed X-ray data provided by this work.

Zoom In Zoom Out Reset image size

Table 6.  Parameter Values Common to All Leptonic Models (LMs) for TXS 0506+056

B′ [G]	0.4
R′ [in cm]	1017
δ	24.2
[in erg s−1]
	1.9
	3.6
	1

Note. The isotropic-equivalent electron luminosity is . Parameter definitions are provided in Table 5.

Download table as:  ASCIITypeset image

As noted in the previous section, HE photons produced via photohadronic interactions are attenuated in the source and induce an EM cascade whose emission should emerge in the Swift-XRT and NuSTAR bands. As a result, the neutrino and proton luminosities are strongly constrained by the X-ray data. The photon spectrum obtained with already violates the observed X-ray data points. In Figure 4, the upper limit on the all-flavor neutrino flux at the neutrino peak energy is .

In what follows, we show that our neutrino-flux limits are fairly insensitive to the exact parameter values that may affect the photomeson production optical depth.

Proton maximum energy—Motivated by the hypothesis that blazars are UHECR accelerators, i.e., at energies above eV (Murase et al. 2012), we explore the effect of the proton maximum energy on the neutrino-flux upper limits. We thus explore cases with , and —see Table 7. Our results on the neutrino fluxes are presented in Figure 5.

Zoom In Zoom Out Reset image size

Table 7.  Model-specific Parameter Values for Leptonic Models (LMs) for TXS 0506+056 Discussed in the Text

	LMBB1a	LMBB1b	LMBB1c	LMBB2a	LMBB2b	LMBB2c	LMPL1a	LMPL1b	LMPL2a	LMPL2b
[1044 erg s−1]	0.54	0.27	0.34	1	5.4	10	0.54	0.54	10	10
sp	2	2.5	3	2	2	2	2	2	2	2
	1	3 × 106	3 × 106	1	1	1	1	1	1	1
[108]	30	30	30	1.6	0.16	0.016	30	30	0.016	0.016
[erg cm−3]	0.033						0.033	0.067	0.04	0.08
[K]							n/a
α	n/a						3	2	3	2
[keV]	n/a						0.05
[keV]	n/a						5

Note. See Table 5 for parameter definitions and Table 6 for parameter values common to all LMs. In LMBB models, the external photon field is blackbody-like with comoving temperature , while in LMPL models, it is a power law between comoving energies and , with photon index α. In all cases, is the comoving energy density of the external photon field. Note that the isotropic-equivalent cosmic-ray proton luminosity is .

Download table as:  ASCIITypeset image

Neutrino spectra in the LMBB1x models are more extended in energy compared to the default case (LMBB2b). They peak around 10 PeV (100 PeV) for () for LMBB2b (LMBB2a), respectively. The number density of target photons decreases fast with increasing energy, while the photomeson production efficiency increases with energy (Murase 2017). However, the upper limits imposed on the proton luminosity and the peak neutrino flux are comparable in the LMBB2a and LMBB2b models. This is because the peak neutrino flux is bounded by the X-ray data points through EM cascades, even though the photomeson production optical depths are quite different. As such, even lower maximum proton energies, e.g., , should not lead to higher upper limits on the neutrino flux. The reason is that protons with will produce electron–positron pairs (via the Bethe–Heitler process) on the synchrotron photons from the peak of the spectrum. Meanwhile, the photomeson interactions of the same protons on the X-ray photons ( Hz) are less efficient (Dimitrakoudis et al. 2012). The proton luminosity cannot be arbitrarily large in this regime because the synchrotron emission from the Bethe–Heitler pairs will overshoot the X-ray data.

Proton spectral index—The slope of the power-law proton distribution is hardly constrained from the SED fitting. Here, we investigate its effects on the neutrino spectrum by considering two additional cases with and s_p = 3. For particle distributions with soft spectra (i.e., ), the total energy in protons is determined by the low-energy cutoff () of the distribution. These low-energy protons, however, are of no interest for HE neutrino production. In an attempt to minimize the energy budget, while retaining the HE neutrino fluxes for , one has to assume —see Table 7. The large can also be justified if the proton distribution has a broken power law and the lower-energy segment has s_p < 2 below the break (i.e., ). Our results on the neutrino flux are presented in Figure 5 and compared to those obtained for s_p = 2. The neutrino spectra become more sharply peaked as the proton distribution becomes softer, while the constraints on the 0.1–10 PeV neutrino flux approach those of our fiducial model (LMBB2b).

External radiation spectrum—Importantly, our results on the neutrino-flux upper limit are insensitive to details of the unknown photon spectrum of external radiation fields. In addition to the external blackbody spectrum, we also consider a power-law spectrum. Such a broadband spectrum might be produced, for example, by electrons accelerated with a hard spectrum in the sheath region of a structured jet, with the associated synchrotron photons—with a low synchrotron peak—serving as seeds for the EIC emission in the γ-ray range (Tavecchio et al. 2014; Tavecchio & Ghisellini 2015). From Figure 5, we see that the results for LMPL1x models with do not differ much from those for LMBB1x models. This is because the relativistic protons at ultrahigh energies mainly interact with target photons around the synchrotron peak. On the other hand, LMPL2x models with give more optimistic neutrino fluxes than LMBB2c, because the photomeson interaction rate is enhanced compared to the photopair production rate (Petropoulou & Mastichiadis 2015). However, the neutrino-flux upper limits are still saturated at erg cm⁻² s⁻¹, similar to that for LMBB2b. We thus conclude that the neutrino-flux upper limit is erg cm⁻² s⁻¹ whether the unknown target spectrum of the external radiation is described by a broadband power law or a narrower blackbody. In this work, we use the blackbody spectrum as a fiducial case, which is conservative in the sense that it introduces fewer free parameters.

In summary, in the LMs (Figure 4), the γ rays are explained by the EIC emission, while there is a small contribution of the SSC component to the hard X-ray band. We note that the SSC component alone cannot explain the γ-ray component of the SED, mainly because of (i) the separation of the low- and high-energy humps of the SED, (ii) the steep Swift-XRT spectrum with the low synchrotron peak inferred by the X-shooter data, and (iii) the flat broad Fermi-LAT spectrum. Accelerated protons, generating HE neutrinos by photohadronic processes, are also present in this scenario, but with an associated EM component that is subdominant in γ rays. The maximal all-flavor neutrino flux over is erg cm⁻² s⁻¹, implying a Poisson probability to detect one event with IceCube over the six-month duration of the TXS 0506+056 γ-ray flare of at most ∼1% under our assumed conditions, which are subject to model and observational constraints but otherwise optimal for HE neutrino production. See Table 9 and Section 3.4 below for estimates of the expectation number of HE muon neutrinos for different model cases. The maximum proton isotropic-equivalent luminosity consistent with the SED is erg s⁻¹. Cases with proton luminosities exceeding lead to higher neutrino fluxes, but they are bounded by the observed X-ray data due to EM cascade effects—as shown in the inset plot, the “twice maximal” case already violates these constraints. We studied different parameters to investigate the parameter dependence, and considered both blackbody-like and power-law spectra for the external target radiation field. As a result, we find that the LM can provide at most a few percent expectation of an associated HE neutrino detection by IceCube.

In hadronic scenarios, while the low-energy peak in the blazar’s SED is explained by synchrotron radiation from relativistic primary electrons, the HE peak is explained by EM cascades induced by pions and muons as decay products of the photomeson production (Mannheim 1993; Mücke et al. 2003), or synchrotron radiation from relativistic protons in the ultrahigh-energy range (Aharonian 2000; Mücke et al. 2003). We coin this scenario “HM,” which stands for Hadronic Model, in reference to the hadronic origin of the γ-rays. The synchrotron and IC emission of secondary pairs may provide an important contribution to the bolometric radiation of the source. In contrast to the leptonic scenario (Section 3.2), the parameters describing the proton distribution can be directly constrained from the NuSTAR and Fermi-LAT data. For the TXS 0506+056 flare, in the hadronic scenario, the SED can be fully explained without invoking external radiation fields.

There are different combinations of parameters that can successfully explain the SED in the HM scenario (Böttcher et al. 2013; Cerruti et al. 2015). As a starting point, we search for combinations of δ and B′ that lead to rough energy equipartition between the magnetic field and protons, since the primary electron energy density is negligible in this scenario. With analytical calculations, we derive rough estimates of the parameter values for equipartition: , G, cm, and GeV (Petropoulou & Dermer 2016).

The parameter values obtained by numerically modeling the SED (see Figure 6) are summarized in Table 8 and are similar to the estimates provided above. The jet power computed for this parameter set (HM1) is close to the minimum value expected in the hadronic scenarios. More specifically, the absolute power of a two-sided jet inferred for these parameters is erg s⁻¹, with erg cm⁻³, where , , are comoving energy densities of relativistic protons, electrons, and magnetic fields, respectively. As demonstrated in Figure 6, the emission from the EM cascade forms a “bridge” between the low-energy and high-energy peaks of the SED for (gray dotted line). Despite minimizing the power of the jet, the adopted set of parameters for HM1 cannot explain the SED due to the associated significant EM cascade component.

Zoom In Zoom Out Reset image size

Table 8.  Parameter Values for Hadronic Models (HMs) for TXS 0506+056 Discussed in the Text and Presented in Figure 6

	HM1	HM2	HM3
B′ [G]	85
R′ [in 1016cm]	2	3	4.5
δ	5.2	10	15
[in 1043 erg s−1] 9.3	0.6	0.06
	1.8
	4.2	3.6	3.6
[in 102]	6.3	1	1
[in 102]	7.9	6.3	5
	104
[in 1046 erg s−1]	2.7	0.1	0.01
sp	2.1
	1
	2 × 109

Note. Parameter definitions are provided in Table 5.

Download table as:  ASCIITypeset image

The EM cascade emission can be suppressed if the source becomes less opaque to the intra-source absorption for HE photons. This can be achieved for larger values of the Doppler factor since (see also Murase et al. 2016; Petropoulou et al. 2017, for analytical expressions). The photon and neutrino spectra for δ = 10 and 15 are also shown in Figure 6, while the respective parameter sets (HM2 and HM3) are listed in Table 8. The SED is compatible with (gray solid line). However, the photomeson production optical depth becomes lower as the two-photon annihilation optical depth decreases. In fact, the sub-PeV neutrino production efficiency is related to the opaqueness for γ rays in the Fermi-LAT band (Murase et al. 2016). Furthermore, this model unavoidably leads to a higher jet power, i.e., erg s⁻¹ and (Petropoulou et al. 2017). Moreover, as the Doppler factor increases, the peak of the neutrino energy spectrum is pushed to the EeV energy range (Dimitrakoudis et al. 2014), while the neutrino flux in the 100 TeV–10 PeV range decreases due to the low efficiency of photomeson interactions.

HM3 demonstrates that the SED of the TXS 0506+056 flare can nicely be explained by the proton synchrotron model, but with the consequence that the HE neutrino production inside the source is very inefficient because of the X-ray constraints on EM cascade emission. The acceleration of UHECRs with EeV is promising in this model, but cannot be reconciled with an IceCube-170922A association, since the predicted neutrino flux is too low in the 0.1–10 PeV energy range.

In summary, we find that no reconciliation of the EM and neutrino observations is possible in hadronic models (HMs; Figure 6). The proton-induced cascade model that predicts ∼ unavoidably overshoots the observed X-ray flux, giving erg cm⁻² s⁻¹, which is strongly excluded. Alternatively, the proton synchrotron model can explain the TXS 0506+056 γ-ray emission, but gives a maximal neutrino flux erg cm⁻² s⁻¹, which implies a very low probability for IceCube neutrino detection, . If the Doppler factor is sufficiently large, the proton-induced cascade emission is suppressed and can avoid overproduction of X-rays (see the main and inset plots over 0.3–100 keV), but at the price of a reduced neutrino flux; hence, only the low neutrino-flux case (red solid curve) is viable. Such low neutrino-flux cases, leading to negligible HE muon neutrino detection probabilities, cannot accommodate production of IceCube-170922A.

Relativistic protons of energy can interact with photons in the source and produce neutrinos with energy .²³ The targets for photohadronic interactions can be synchrotron and IC photons emitted by primary and secondary electrons as well as external photons to the source, if present. For a typical synchrotron spectrum around the SED peak, the rectilinear approximation around the Δ resonance is usually valid (Murase & Nagataki 2006; Murase et al. 2014), and the characteristic proton Lorentz factor interacting with photons of the frequency is given by (Murase 2017)

where is the resonance energy. The respective neutrino energy, in the observer’s frame, is then given by . This is also an estimate of the peak energy of the HE neutrino spectrum:

Neutrinos with lower energies than can still be produced by interactions of lower-energy protons with higher-energy photons. However, the neutrino flux at such lower energies is expected to be lower than the flux at due to the decreasing number density of target photons, unless the proton distribution is a super-soft power law (i.e., ). Similar estimates can be derived in the presence of external radiation fields, as demonstrated in Murase et al. (2014). Note that the main target photons for the photomeson production process in the LMBB2 and LMPL2 models are photons with energies above the synchrotron peak, since even protons with do not satisfy the photomeson production threshold for the peak synchrotron photons—see Equation (1).

As described above, EM cascade emission induced by cosmic-ray protons readily fills the dip between the two peaks of the SED (keV to MeV energies). Thanks to the optical depth correspondence between photomeson production and two-photon annihilation, efficient production of HE neutrinos can only be achieved for conditions that lead to a stronger EM cascade emission inside the source. In particular, the proton-induced cascade model, where the Fermi-LAT γ-ray data is primarily explained by the proton-induced cascade emission itself, is strongly ruled out. Although this model naturally predicts erg cm⁻² s⁻¹, which is consistent with the observation of IceCube-170922A (Aartsen et al. 2018a), the EM cascade component typically has a broadband energy spectrum that extends through the X-ray range. Thus, this model does not provide a proper description of the broadband SED. On the other hand, the proton synchrotron model can explain the blazar’s SED, but the 0.1–10 PeV neutrino flux is predicted to be very low, erg cm⁻² s⁻¹.

Figure 5 presents the upper limits on the all-flavor neutrino fluxes obtained in the leptonic and hadronic models for a wide set of parameters. For the maximum neutrino flux displayed in the figure, , the corresponding muon neutrino fluence is estimated to be:

where is the flare duration and the flavor ratio is assumed to be . Then, using the effective area for EHE real-time alerts, in the PeV range (Aartsen et al. 2017a), the expected number of muon neutrinos is estimated to be:

where the width of the neutrino spectrum is assumed to be . We numerically confirm this analytical estimate for the different models presented in Figure 5.

Table 9 summarizes our results on the upper limits obtained for the integrated all-flavor neutrino fluxes by modeling the flaring SED of TXS 0506+056. Assuming a flare duration of ΔT = 0.5 year, we also evaluate the expected number of muon neutrinos by using the effective area for EHE real time alerts (Aartsen et al. 2017a), taking into account a correction due to the Earth attenuation toward the direction to TXS 0506+056. The exact effective area may be slightly different for IceCube-170922A (Aartsen et al. 2018a), but the results do not change within a factor of two. In the LMBB2b model or LMPL2b model, which give the most optimistic neutrino fluxes among our parameter sets, for example, the probability for IceCube to observe one event that physically originates from the TXS 0506+056 flare is evaluated to be % based on Poisson statistics. This value is already achieved for the near-optimal case, but we allow a factor of two uncertainty for several possible reasons. First, the duration of the HE neutrino flare may be a bit longer (although the confidence of association would be reduced for longer durations). Second, the EHE effective area for IceCube-170922A may be slightly different. Then, taking into account these variations as well as model uncertainties, we may regard the case for the proton luminosity of as the most conservative limit, which gives %.

Table 9.  The Upper Limit on the All-flavor Neutrino Flux for the Different Models that Satisfactorily Explain the Flaring SED of TXS 0506+056—for Details, See Sections 3.2 and 3.3

	[erg cm−2 s−1]
	100 TeV–1 PeV	100 TeV–10 PeV	≤10 PeV
LMBB1a	1.6 × 10−14	4.5 × 10−13	1 × 10−3
LMBB1b	5.2 × 10−14	1.7 × 10−12	4 × 10−3
LMBB1c	9.1 × 10−14	2.7 × 10−12	6 × 10−3
LMBB2a	4.5 × 10−14	1.1 × 10−12	3 × 10−3
LMBB2b	1.8 × 10−13	3.6 × 10−12	8 × 10−3
LMBB2c	2.5 × 10−14	7.3 × 10−14	2 × 10−4
LMPL1a	3.1 × 10−14	5.2 × 10−13	1 × 10−3
LMPL1b	9 × 10−14	6.3 × 10−13	1 × 10−3
LMPL2a	2.5 × 10−13	5.2 × 10−13	5 × 10−3
LMPL2b	1.2 × 10−12	2 × 10−12	1 × 10−2
HM3	1.6 × 10−16	2 × 10−15	4 × 10−6

Note. The reported values correspond to the neutrino fluxes integrated over 0.1–1 PeV and 0.1–10 PeV ranges. The last column shows the expected number of muon neutrinos below 10 PeV for a flare duration ΔT = 0.5 year.

Download table as:  ASCIITypeset image