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

A. Keivani, K. Murase, M. Petropoulou, D.B. Fox, S.B. Cenko, S. Chaty, A.

Coleiro, J.J. Delaunay, S. Dimitrakoudis, P.A. Evans, et al.

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

A. Keivani, K. Murase, M. Petropoulou, D.B. Fox, S.B. Cenko, et al., A Multimessenger Picture of the Flaring Blazar TXS 0506+056: implications for High-Energy Neutrino Emission and Cosmic Ray Acceleration. The Astrophysical Journal, 2018, 864 (1), pp.84. 10.3847/1538-4357/aad59a .
hal-01851256

HAL Id: hal-01851256

https://hal.science/hal-01851256v1

Submitted on 19 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

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

A. Keivani<sup>1,2,17</sup>①, K. Murase<sup>1,2,3,4</sup>①, M. Petropoulou<sup>5</sup>①, D. B. Fox<sup>2,3,6</sup>①, S. B. Cenko<sup>7,8</sup>①, S. Chaty<sup>9</sup>①, A. Coleiro<sup>10,11</sup>①, J. J. DeLaunay<sup>1,2</sup>①, S. Dimitrakoudis<sup>12</sup>①, P. A. Evans<sup>13</sup>①, J. A. Kennea<sup>3</sup>①, F. E. Marshall<sup>14</sup>①, A. Mastichiadis<sup>15</sup>, J. P. Osborne<sup>13</sup>①, M. Santander<sup>16</sup>①, A. Tohuvavohu<sup>3</sup>①, and C. F. Turley<sup>1,2</sup>①

<sup>1</sup>Department of Physics, Pennsylvania State University, University Park, PA 16802, USA; keivani@psu.edu

<sup>2</sup>Center for Particle and Gravitational Astrophysics, Institute for Gravitation and the Cosmos, Pennsylvania State University, University Park, PA 16802, USA; murase@psu.edu

<sup>3</sup>Department of Astronomy and Astrophysics, Pennsylvania State University, University Park, PA 16802, USA

<sup>4</sup>Center for Gravitational Physics, Yukawa Institute for Theoretical Physics, Kyoto, Kyoto 606-8502 Japan

<sup>5</sup>Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA; m.petropoulou@astro.princeton.edu

<sup>6</sup>Center for Theoretical and Observational Cosmology, Institute for Gravitation and the Cosmos, Pennsylvania State University, University Park, PA 16802, USA; dfox@psu.edu

<sup>7</sup>Astrophysics Science Division, NASA Goddard Space Flight Center, Mail Code 661, Greenbelt, MD 20771, USA

<sup>8</sup>Joint Space-Science Institute, University of Maryland, College Park, MD 20742, USA

<sup>9</sup>Laboratoire AIM (UMR 7158 CEA/DRF-CNRS-Université Paris Diderot), Irfu/Département d'Astrophysique, CEA-Saclay, F-91191 Gif-sur-Yvette Cedex, France

<sup>10</sup>APC, Univ Paris Diderot, CNRS/IN2P3, CEA/Irfu, Obs de Paris, Sorbonne Paris Cité, France

<sup>11</sup>IFIC-Instituto de Física Corpuscular (CSIC-Universitat de València), Calle Catedrático José Beltrán, 2 E-46980 Paterna, Valencia, Spain
<sup>12</sup>Department of Physics, University of Alberta, Edmonton, Alberta T6G 2E1, Canada

<sup>13</sup>Department of Physics and Astronomy, University of Leicester, Leicester, LEI 7RH, UK

<sup>14</sup>NASA Goddard Space Flight Center, Mail Code 660.1, Greenbelt, MD 20771, USA

<sup>15</sup>Department of Physics, National and Kapodistrian University of Athens, Panepistimiopolis, GR 15783 Zografos, Greece

<sup>16</sup>Department of Physics and Astronomy, University of Alabama, Tuscaloosa, AL 35487, USA

Received 2018 May 17; revised 2018 July 20; accepted 2018 July 20; published 2018 August 31

## Abstract

Detection of the IceCube-170922A neutrino coincident with the flaring blazar TXS 0506+056, the first and only ~3σ high-energy neutrino source association to date, offers a potential breakthrough in our understanding of high-energy cosmic particles and blazar physics. We present a comprehensive analysis of TXS 0506+056 during its flaring state, using newly collected *Swift*, *NuSTAR*, and X-shooter data with *Fermi* observations and numerical models to constrain the blazar's particle acceleration processes and multimessenger (electromagnetic (EM) and high-energy neutrino) emissions. Accounting properly for EM cascades in the emission region, we find a physically consistent picture only within a hybrid leptonic scenario, with γ-rays produced by external inverse-Compton processes and high-energy neutrinos via a radiatively subdominant hadronic component. We derive robust constraints on the blazar's neutrino and cosmic-ray emissions and demonstrate that, because of cascade effects, the 0.1–100 keV emissions of TXS 0506+056 serve as a better probe of its hadronic acceleration and high-energy neutrino production processes than its GeV–TeV emissions. If the IceCube neutrino association holds, physical conditions in the TXS 0506+056 jet must be close to optimal for high-energy neutrino production, and are not favorable for ultrahigh-energy cosmic-ray acceleration. Alternatively, the challenges we identify in generating a significant rate of IceCube neutrino detections from TXS 0506+056 may disfavor single-zone models, in which γ-rays and high-energy neutrinos are produced in a single emission region. In concert with continued operations of the high-energy neutrino observatories, we advocate regular X-ray monitoring of TXS 0506+056 and other blazars in order to test single-zone blazar emission models, clarify the nature and extent of their hadronic acceleration processes, and carry out the most sensitive possible search for additional multimessenger sources.

**Key words:** BL Lacertae objects: general – BL Lacertae objects: individual (TXS 0506+056) – galaxies: active – gamma rays: galaxies – neutrinos – radiation mechanisms: non-thermal

*Supporting material:* machine-readable table

## 1. Introduction

High-energy (HE; $\epsilon_{\nu} \gtrsim 1$ TeV) neutrinos, as cosmic messenger particles, have the potential to reveal the sources of HE cosmic rays and illuminate their underlying particle acceleration processes. Detection of a nearly isotropic flux of HE cosmic neutrinos has been reported by the IceCube Collaboration (Aartsen et al. 2013a, 2013b); the absence of any identified

point sources (Aartsen et al. 2016b) places firm limits on possible contributions from persistently bright neutrino sources (Murase & Waxman 2016). In this context, multimessenger studies have provided important clues to the origins of the diffuse neutrino, γ-ray, and cosmic-ray backgrounds (Murase et al. 2013; Fang & Murase 2018), and offer a powerful approach for identifying transient or highly variable neutrino sources, including blazar flares (Dermer et al. 2014; Kadler et al. 2016; Petropoulou et al. 2016; Turley et al. 2016). Anticipating these and related opportunities, the Astrophysical

<footnote n="17"><sup>17</sup>Now at Department of Physics, Columbia University, New York, NY 10027; ak4279@columbia.edu.</footnote>
Multimessenger Observatory Network (AMON$^{18}$) was founded to link global HE and multimessenger observatories together into a single network and to distribute relevant alerts to the community in real-time (Smith et al. 2013).

Blazars—active galactic nuclei oriented with a relativistic jet pointing toward Earth—dominate the extragalactic γ-ray sky (Ackermann et al. 2015, 2016). Yet, despite a wealth of electromagnetic (EM) data, blazar radiation mechanism(s) remain unclear, with leptonic and (lepto)hadronic scenarios providing viable explanations (e.g., Boettcher et al. 2013). Since blazars and other jetted active galactic nuclei are proposed ultrahigh-energy cosmic-ray (UHECR; $\epsilon_{cr} \gtrsim 3$ EeV) accelerators (e.g., Murase et al. 2012), the question of hadronic acceleration in these sources has important broader implications, and information from HE neutrinos will likely be crucial to resolving these issues (Murase 2017).

The AMON ICECUBE EHE alert 50579430 (hereafter, IceCube-170922A) was identified by IceCube and publicly distributed via AMON and the Gamma-ray Coordinates Network (GCN) within $\delta t \approx 43$ s of its interaction in the Antarctic ice cap at 20:54:30.43 UT on 2017 September 22 (GCN/AMON NOTICE IceCube-170922A 2017). As with previous likely cosmic events, its location was soon targeted by multiple observatories covering a broad energy range. The Swift-XRT (Keivani et al. 2017) and Fermi Large Area Telescope (LAT; Tanaka et al. 2017) reported an association with a blazar, TXS 0506+056, which showed strong activity in LAT data beginning in 2017 April and significant X-ray variability during Swift monitoring observations. Broadband EM observations of this event, the significance of the blazar flare, and the high-energy neutrino coincidence are discussed in Aartsen et al. (2018a).

The present work is organized as follows. In Section 2, we present a comprehensive analysis of Fermi (γ-ray), NuSTAR (hard X-ray), Swift (X-ray, ultraviolet, optical), and X-shooter (ultraviolet, optical, near-infrared) observations of TXS 0506+056, using these data to construct the source spectral energy distribution (SED) over near-infrared to γ-ray energies for two epochs: a 30 day period centered around the time of the neutrino trigger (Ep. 1) and a 30 day period starting 15 days after the neutrino trigger (Ep. 2). In Section 3, we model the SED of TXS 0506+056 in the flaring phase (Ep. 1) by performing detailed radiative transfer calculations, focusing on leptonic and hadronic single-zone emission scenarios. We discuss the implications of our modeling results in Section 4, and conclude in Section 5.

Throughout the manuscript, we use the notation $Q_x$ as a shorthand for the quantity $Q/10^x$ (with $Q$ in cgs units) unless stated otherwise. We note that given its redshift $z = 0.3365$ (Paiano et al. 2018) and a consensus cosmology, the luminosity distance of TXS 0506+056 is $d_L \approx 1750$ Mpc.

## 2. Observations and Analysis

In this section, we review the IceCube detection of the IceCube-170922A neutrino and present observations and data analysis for EM follow-up observations from the Fermi (γ-ray), NuSTAR (hard X-ray), Swift (X-ray, ultraviolet, optical), and X-shooter (ultraviolet, optical, near-infrared) facilities. We note that while very high-energy ($\epsilon_{\gamma} \gtrsim 100$ GeV) γ-ray observations by multiple facilities have

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Swift-XRT follow-up of IceCube-170922A. X-ray exposure map resulting from the adopted 19-point tiling pattern centered on the initial IceCube neutrino localization is shown in grayscale, and the positions of all detected X-ray sources with red points. The red dashed circle shows the initial 90%-containment region. The red solid ellipse shows the updated 90%-containment region (Kopper & Blaufuss 2017). Grayscale levels indicate achieved exposure at each sky position, as shown by the color bar. White streaks are due to dead regions on the XRT detector caused by a micrometeroid impact (Abbey et al. 2006).</figcaption></figure>

been reported (de Naurois & H.E.S.S. Collaboration 2017; Mukherjee 2017), including a first detection at these energies by MAGIC (Mirzoyan 2017), details are not yet publicly available. These constraints are not included in our analysis. Given our results, consistency of our models with Fermi data out to $\epsilon_{\gamma} \approx 100$ GeV, and current uncertainties in the necessary extragalactic background light corrections for the source at these energies, we do not expect that inclusion of these constraints would alter our conclusions.

### 2.1. IceCube Data

IceCube-170922A was an EHE neutrino event (GCN/ AMON NOTICE 2016) identified and distributed by the IceCube Observatory via AMON and GCN within $\delta t \approx 43$ s 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. = 77.43$^{+1.3}_{-0.8}$ degree, decl. = +5.72$^{+0.7}_{-0.4}$ degree (J2000; 90% containment ellipse). The maximum likelihood neutrino position is R.A. 05°09'm08''784, decl. +05°45'13"32 (J2000); see Figure 1 for an illustration of the initial and final localizations.

EHE neutrino event reports include the neutrino arrival time, direction (R.A. and decl.), angular error ($r_{50}$ for 50% containment; $r_{90}$ 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

<footnote n="18">18 AMON website: https://www.amon.psu.edu/.</footnote>
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).

## 2.2. Swift-XRT Data

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 ≈0.8 centered on R.A. 05h09m08s784, decl. +05°45′ 13″32 (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 3.8 × 10⁻¹³ 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 4.6 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).

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 $N_H = 1.11 \times 10^{21} \text{ cm}^{-2}$, resulting in a photon index $\alpha_{XRT} = 2.37 \pm 0.05$ and mean flux of $2.27 \times 10^{-12} \text{ erg cm}^{-2} \text{ s}^{-1}$ (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 $1 \times 10^{-12} \text{ erg cm}^{-2} \text{ s}^{-1}$, with one observation at ≈2.8 × 10⁻¹² 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).

## 2.3. NuSTAR Data

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 nupipe-line 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 $N_H = 1.11 \times 10^{21} \text{ cm}^{-2}$, resulting in a photon index of $\alpha_{NuSTAR} = 1.69 \pm 0.12$ and a flux of $4.27_{-0.58}^{+0.50} \times 10^{-12} \text{ erg cm}^{-2} \text{ s}^{-1}$ (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 $\alpha_{NuSTAR} = 1.68 \pm 0.14$ and flux of $3.65_{-0.59}^{+0.54} \times 10^{-12} \text{ erg cm}^{-2} \text{ s}^{-1}$ (3.0–79.0 keV), consistent with results of the first observation from 20 days earlier.

## 2.4. Joint Swift-XRT and NuSTAR Analysis

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 $N_H = 1.11 \times 10^{21} \text{ cm}^{-2}$, 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 ($F_{\epsilon_\gamma} \propto \epsilon^{-\alpha}$) of $\alpha_{XRT} = 2.37 \pm 0.05$ and $\alpha_{NuSTAR} = 1.68 \pm 0.14$. The results are displayed in Figure 3 and used for subsequent SED modeling (see Section 3).

## 2.5. Swift UltraViolet and Optical Telescope (UVOT) Data

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
<table><caption>**Table 1**

Swift-XRT Monitoring of TXS 0506+056</caption><thead><tr><th>Epoch<br>(1)</th><th>Exposure [ks]<br>(2)</th><th>Photon Index<br>(3)</th><th>R<sub>X-3</sub><br>(4)</th><th>F<sub>X-12</sub><br>(5)</th></tr></thead><tbody><tr><td>58019.4700 ± 0.4630</td><td>0.8</td><td>1.83<sup>+0.43</sup><sub>-0.42</sub></td><td>65.8 ± 10.1</td><td>2.33<sup>+1.07</sup><sub>-0.72</sub></td></tr><tr><td>58023.8535 ± 0.0660</td><td>4.9</td><td>2.43 ± 0.12</td><td>121.0 ± 5.3</td><td>3.51<sup>+0.32</sup><sub>-0.29</sub></td></tr><tr><td>58026.2273 ± 0.0391</td><td>2.0</td><td>2.30 ± 0.33</td><td>66.2 ± 7.8</td><td>1.55<sup>+0.43</sup><sub>-0.33</sub></td></tr><tr><td>58028.6419 ± 0.0114</td><td>2.0</td><td>2.73 ± 0.20</td><td>117.0 ± 8.2</td><td>2.92<sup>+0.40</sup><sub>-0.35</sub></td></tr><tr><td>58029.6745 ± 0.1036</td><td>1.1</td><td>2.46<sup>+0.22</sup><sub>-0.21</sub></td><td>182.0 ± 14.1</td><td>4.96<sup>+0.80</sup><sub>-0.70</sub></td></tr><tr><td>58030.7369 ± 0.0385</td><td>1.2</td><td>2.82<sup>+0.26</sup><sub>-0.25</sub></td><td>186.0 ± 16.1</td><td>4.41<sup>+0.74</sup><sub>-0.65</sub></td></tr><tr><td>58031.7944 ± 0.1025</td><td>2.3</td><td>2.64 ± 0.13</td><td>255.0 ± 11.7</td><td>6.47<sup>+0.57</sup><sub>-0.53</sub></td></tr><tr><td>58032.8985 ± 0.3428</td><td>2.1</td><td>2.36<sup>+0.22</sup><sub>-0.21</sub></td><td>90.1 ± 7.1</td><td>2.56<sup>+0.44</sup><sub>-0.37</sub></td></tr><tr><td>58034.4478 ± 0.0337</td><td>1.9</td><td>2.53 ± 0.21</td><td>108.0 ± 8.1</td><td>2.87<sup>+0.44</sup><sub>-0.38</sub></td></tr><tr><td>58040.9452 ± 0.0010</td><td>0.2</td><td>2.22<sup>+1.04</sup><sub>-0.97</sub></td><td>82.9<sup>+29.8</sup><sub>-24.1</sub></td><td>2.15<sup>+2.35</sup><sub>-1.06</sub></td></tr><tr><td>58042.7684 ± 0.1686</td><td>2.2</td><td>2.00 ± 0.27</td><td>61.9 ± 5.8</td><td>2.20<sup>+0.59</sup><sub>-0.44</sub></td></tr><tr><td>58044.1300 ± 0.0696</td><td>1.8</td><td>2.10 ± 0.31</td><td>52.6 ± 6.0</td><td>1.70<sup>+0.49</sup><sub>-0.37</sub></td></tr><tr><td>58047.2130 ± 0.6360</td><td>2.3</td><td>2.11 ± 0.28</td><td>49.6 ± 5.1</td><td>1.61<sup>+0.41</sup><sub>-0.32</sub></td></tr><tr><td>58050.7286 ± 0.0446</td><td>2.9</td><td>2.08<sup>+0.25</sup><sub>-0.24</sub></td><td>46.6 ± 4.4</td><td>1.56<sup>+0.35</sup><sub>-0.28</sub></td></tr><tr><td>58053.3162 ± 0.0335</td><td>1.0</td><td>2.20<sup>+0.42</sup><sub>-0.41</sub></td><td>54.4 ± 8.2</td><td>1.68<sup>+0.66</sup><sub>-0.46</sub></td></tr><tr><td>58059.6620 ± 0.0729</td><td>3.3</td><td>2.22 ± 0.21</td><td>60.2 ± 4.7</td><td>1.84<sup>+0.32</sup><sub>-0.27</sub></td></tr><tr><td>58065.6390 ± 0.3371</td><td>3.1</td><td>2.30 ± 0.20</td><td>83.4 ± 6.2</td><td>2.36<sup>+0.39</sup><sub>-0.33</sub></td></tr><tr><td>58068.3642 ± 0.0742</td><td>3.0</td><td>2.38 ± 0.24</td><td>55.7 ± 4.9</td><td>1.55<sup>+0.30</sup><sub>-0.25</sub></td></tr><tr><td>58069.5359 ± 0.5032</td><td>1.8</td><td>2.23 ± 0.24</td><td>83.3 ± 7.5</td><td>2.52<sup>+0.52</sup><sub>-0.42</sub></td></tr><tr><td>58071.1551 ± 0.1299</td><td>2.9</td><td>2.15 ± 0.19</td><td>81.1 ± 5.9</td><td>2.51<sup>+0.42</sup><sub>-0.35</sub></td></tr><tr><td>58072.0949 ± 0.0039</td><td>0.7</td><td>2.80<sup>+0.60</sup><sub>-0.59</sub></td><td>48.5 ± 9.1</td><td>1.24<sup>+0.55</sup><sub>-0.37</sub></td></tr><tr><td>58073.0911 ± 0.0042</td><td>0.7</td><td>1.98<sup>+0.47</sup><sub>-0.45</sub></td><td>55.0 ± 9.7</td><td>1.84<sup>+0.89</sup><sub>-0.59</sub></td></tr><tr><td>58074.2484 ± 0.0962</td><td>3.0</td><td>2.03<sup>+0.24</sup><sub>-0.23</sub></td><td>54.4 ± 4.8</td><td>1.82<sup>+0.39</sup><sub>-0.32</sub></td></tr><tr><td>58075.0831 ± 0.0046</td><td>0.8</td><td>2.51 ± 0.45</td><td>63.1 ± 9.6</td><td>1.75<sup>+0.66</sup><sub>-0.46</sub></td></tr><tr><td>58075.6655 ± 0.0057</td><td>1.0</td><td>1.94 ± 0.41</td><td>60.6 ± 9.4</td><td>2.13<sup>+0.90</sup><sub>-0.62</sub></td></tr><tr><td>58076.0793 ± 0.0049</td><td>0.8</td><td>1.71<sup>+0.45</sup><sub>-0.44</sub></td><td>56.1 ± 9.0</td><td>2.26<sup>+1.14</sup><sub>-0.74</sub></td></tr><tr><td>58077.0756 ± 0.0054</td><td>0.9</td><td>2.27<sup>+0.41</sup><sub>-0.40</sub></td><td>61.4 ± 8.9</td><td>1.81<sup>+0.66</sup><sub>-0.48</sub></td></tr><tr><td>58078.0717 ± 0.0056</td><td>1.0</td><td>2.42<sup>+0.52</sup><sub>-0.51</sub></td><td>40.5 ± 7.1</td><td>1.15<sup>+0.52</sup><sub>-0.34</sub></td></tr><tr><td>58079.0675 ± 0.0057</td><td>1.0</td><td>2.02<sup>+0.47</sup><sub>-0.46</sub></td><td>49.3 ± 8.6</td><td>1.68<sup>+0.80</sup><sub>-0.53</sub></td></tr><tr><td>58080.1320 ± 0.0057</td><td>1.0</td><td>2.18<sup>+0.55</sup><sub>-0.53</sub></td><td>47.4 ± 9.4</td><td>1.58<sup>+0.86</sup><sub>-0.54</sub></td></tr><tr><td>58081.1279 ± 0.0057</td><td>1.0</td><td>2.51<sup>+0.69</sup><sub>-0.70</sub></td><td>45.8 ± 9.7</td><td>1.33<sup>+0.93</sup><sub>-0.48</sub></td></tr><tr><td>58082.0551 ± 0.0057</td><td>1.0</td><td>3.43<sup>+1.16</sup><sub>-0.93</sub></td><td>30.0<sup>+9.4</sup><sub>-7.8</sub></td><td>0.85<sup>+0.47</sup><sub>-0.34</sub></td></tr><tr><td>58083.1195 ± 0.0057</td><td>1.0</td><td>2.45 ± 0.56</td><td>44.4 ± 8.0</td><td>1.26<sup>+0.65</sup><sub>-0.39</sub></td></tr><tr><td>58084.0499 ± 0.0058</td><td>1.0</td><td>1.58<sup>+0.50</sup><sub>-0.51</sub></td><td>39.0 ± 6.9</td><td>1.87<sup>+1.22</sup><sub>-0.69</sub></td></tr><tr><td>58086.1089 ± 0.0059</td><td>1.0</td><td>2.49<sup>+0.67</sup><sub>-0.63</sub></td><td>45.4 ± 9.8</td><td>1.44<sup>+0.84</sup><sub>-0.52</sub></td></tr><tr><td>58087.1560 ± 0.0053</td><td>0.9</td><td>2.53<sup>+0.39</sup><sub>-0.38</sub></td><td>66.9 ± 9.2</td><td>1.79<sup>+0.55</sup><sub>-0.42</sub></td></tr></tbody></table>

Note. $R_{X-3}$ and $F_{X-12}$ indicate count rate and energy flux, in units of $10^{-3}$ ct s$^{-1}$ and $10^{-12}$ erg cm$^{-2}$ s$^{-1}$, respectively. Uncertainties are quoted at 90% confidence.

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 UVOT-SOURCE 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.

## 2.6. Swift-UVOT Analysis

Due to the UVOT's blue response, extending as far as $\lambda \approx 1600$ Å, 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
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 2.** (Left top) Swift-XRT light curve. Each bin corresponds to one observation in the 0.3–10 keV energy range. The horizontal bands show the XRT historical data (four observations) of TXS 0506+056: the mean historical quiescent flux from combining three data points, and one showing the rate from the outburst observation. (Left bottom) Swift-UVOT light curve for all 36 observations performed on TXS 0506+056. The dashed line shows the IceCube-170922A arrival time. (Right top) Swift-XRT photon index variation during the XRT monitoring campaign of TXS 0506+056. The solid horizontal line shows the photon index of the stacked X-ray spectrum over the 2 epochs, while the dashed lines represent the uncertainties. (Right bottom) Swift-UVOT photon index variations obtained from a power-law fit to the energy flux spectrum (εγ vs. Fγν). In all plots, Ep. 1 and Ep. 2 are, respectively, defined as [−15 days, +15 days] and [+15 days, +45 days] time windows with respect to the IceCube-170922A arrival time.</figcaption></figure>

the NASA/IPAC Infrared Science Archive<sup>19</sup> 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 ($\delta A_V = 0.008$ 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 $A_V = 0.286 \pm 0.032$ 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 $\chi^2$ minimization of the predicted versus observed count rates in each filter using a model with three parameters: source flux density $F_\nu$ at $\nu = \epsilon_\gamma/h = 10^{15}$ Hz, UVOT spectral index $\beta$, and extinction $A_\nu$. 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<sup>20</sup> and Fitzpatrick (1999) extinction curve. A fit is considered successful if the total $\chi^2$ across the six filters gives a p-value $p(\chi^2 > \chi^2_{\text{obs}}) > 5\%$.

Only one of the six-filter observations (at $\delta 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 $\beta = 1.66 \pm 0.06$ at this epoch, where we quote uncertainties from $\delta\chi^2$ 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 $\Delta\chi^2$ analysis. These values and uncertainties are reported in Table 3 and plotted in Figure 2 (lower right).

We determine UVOT SED bands, the range of fluxes allowed at 90% confidence for each photon energy, by drawing 1000 samples according to the $\chi^2$ 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 $\beta$ measurement from the temporally proximate six-filter observation, adjusting the flux density at $10^{15}$ 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).

<footnote n="19"><sup>19</sup> NASA/IPAC Infrared Science Archive: https://irsa.ipac.caltech.edu/.</footnote>

<footnote n="20"><sup>20</sup> UVOT Calibrations Database (CALDB): https://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/data/swift/uvota/index.html.</footnote>
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 3.** Observations and spectral energy distribution (SED) for TXS 0506+056 in its high-flux state. Left panel: timeline of observations by *Fermi*, *Swift*, *NuSTAR*, and X-shooter. Observations are divided into two 30 day epochs each for analysis and discussion purposes; the vertical dashed line shows the IceCube-170922A detection time. Right panel: multiwavelength SED for TXS 0506+056; data with the 90%-confidence bands on source emission are shown separately for the two epochs for each facility. The SEDs for Ep. 1 and Ep. 2 are broadly similar, with the source fading somewhat at optical through X-ray energies, and the ultraviolet/optical spectrum softening.</figcaption></figure>

<table><caption>**Table 2**

Swift-UVOT Monitoring of TXS 0506+056</caption><thead><tr><th>Epoch<sup>a</sup> (MJD)<br>(1)</th><th>Exposure (s)<br>(2)</th><th>Filter<br>(3)</th><th>Magnitude<sup>b</sup><br>(4)</th></tr></thead><tbody><tr><td>58019.4699 ± 0.4632</td><td>780</td><td>u</td><td>14.31 ± 0.03</td></tr><tr><td>58023.8555 ± 0.0641</td><td>809</td><td>uvw2</td><td>14.58 ± 0.03</td></tr><tr><td>58023.8257 ± 0.0306</td><td>157</td><td>v</td><td>14.62 ± 0.04</td></tr><tr><td>58023.8364 ± 0.0403</td><td>2930</td><td>uvw1</td><td>14.50 ± 0.03</td></tr><tr><td>58023.8515 ± 0.0639</td><td>472</td><td>u</td><td>14.36 ± 0.04</td></tr><tr><td>58023.8529 ± 0.0635</td><td>236</td><td>b</td><td>14.27 ± 0.04</td></tr><tr><td>58023.8539 ± 0.0635</td><td>236</td><td>b</td><td>15.08 ± 0.03</td></tr><tr><td>58026.2274 ± 0.0391</td><td>1954</td><td>uvw2</td><td>14.81 ± 0.04</td></tr></tbody></table>

**Notes.**

<sup>a</sup> MJD at the middle of the observation.

<sup>b</sup> Errors at $1\sigma$ uncertainty.

(This table is available in its entirety in machine-readable form.)

## 2.7. X-shooter Data

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 1′′0, 0′′9, and 0′′9, respectively. These data provide simultaneous 300–2480 nm spectral coverage with average spectral resolutions $\lambda/\Delta\lambda$ of 4290, 7410, and 5410, respectively, in each arm. Observing conditions were good, with a clear sky, seeing of ~0′′8, 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

<table><caption>**Table 3**

Swift-UVOT Photon Indices for TXS 0506+056</caption><thead><tr><th>Epoch</th><th>Photon Index</th></tr></thead><tbody><tr><td>58023.7632</td><td>2.16<sup>+0.11</sup><sub>-0.09</sub></td></tr><tr><td>58050.6452</td><td>2.46<sup>+0.09</sup><sub>-0.11</sub></td></tr><tr><td>58053.2192</td><td>2.35<sup>+0.11</sup><sub>-0.13</sub></td></tr><tr><td>58059.5782</td><td>2.32<sup>+0.08</sup><sub>-0.11</sub></td></tr><tr><td>58065.5562</td><td>2.38<sup>+0.11</sup><sub>-0.09</sub></td></tr><tr><td>58068.2802</td><td>2.41<sup>+0.11</sup><sub>-0.10</sub></td></tr><tr><td>58069.9502</td><td>2.32<sup>+0.14</sup><sub>-0.12</sub></td></tr><tr><td>58070.9472</td><td>2.35<sup>+0.11</sup><sub>-0.12</sub></td></tr><tr><td>58072.0112</td><td>2.27<sup>+0.13</sup><sub>-0.12</sub></td></tr><tr><td>58073.0082</td><td>2.41<sup>+0.13</sup><sub>-0.12</sub></td></tr><tr><td>58074.1632</td><td>2.27<sup>+0.11</sup><sub>-0.09</sub></td></tr><tr><td>58074.9992</td><td>2.32<sup>+0.10</sup><sub>-0.13</sub></td></tr><tr><td>58075.5822</td><td>2.38<sup>+0.12</sup><sub>-0.10</sub></td></tr><tr><td>58075.9962</td><td>2.18<sup>+0.10</sup><sub>-0.13</sub></td></tr><tr><td>58076.9922</td><td>2.44<sup>+0.12</sup><sub>-0.11</sub></td></tr><tr><td>58077.9882</td><td>2.30<sup>+0.13</sup><sub>-0.10</sub></td></tr><tr><td>58078.9842</td><td>2.35<sup>+0.12</sup><sub>-0.11</sub></td></tr><tr><td>58080.0482</td><td>2.38<sup>+0.12</sup><sub>-0.10</sub></td></tr><tr><td>58081.0442</td><td>2.38<sup>+0.10</sup><sub>-0.13</sub></td></tr><tr><td>58081.9712</td><td>2.30<sup>+0.11</sup><sub>-0.12</sub></td></tr><tr><td>58083.9662</td><td>2.27<sup>+0.13</sup><sub>-0.10</sub></td></tr><tr><td>58086.0252</td><td>2.13<sup>+0.10</sup><sub>-0.13</sub></td></tr><tr><td>58083.0362 (UV)</td><td>3.22<sup>+0.37</sup><sub>-0.34</sub></td></tr><tr><td>58083.0362 (Vis)</td><td>2.35<sup>+0.25</sup><sub>-0.24</sub></td></tr><tr><td>58087.0722</td><td>2.10<sup>+0.12</sup><sub>-0.10</sub></td></tr></tbody></table>

**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.

HST white dwarf standard GD71 (R.A. 05h52m27s86, decl. +15°53′13″8, 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
<table><caption>**Table 4**
Swift-UVOT Extinction-corrected Fluxes at $\nu = 10^{15}$ Hz</caption><thead><tr><th>Epoch</th><th>F<sub>UV,-12</sub></th></tr></thead><tbody><tr><td>58019.4709</td><td>22.50 ± 3.29</td></tr><tr><td>58023.8470</td><td>34.80 ± 3.13</td></tr><tr><td>58026.2284</td><td>12.80 ± 3.13</td></tr><tr><td>58029.6756</td><td>15.80 ± 3.13</td></tr><tr><td>58030.7380</td><td>15.00 ± 3.29</td></tr><tr><td>58032.5593</td><td>19.30 ± 3.29</td></tr><tr><td>58033.1319</td><td>11.00 ± 3.13</td></tr><tr><td>58034.4489</td><td>13.80 ± 3.29</td></tr><tr><td>58040.9463</td><td>16.00 ± 1.97</td></tr><tr><td>58042.7696</td><td>15.70 ± 1.97</td></tr><tr><td>58044.1311</td><td>19.60 ± 1.97</td></tr><tr><td>58046.5814</td><td>10.90 ± 1.97</td></tr><tr><td>58047.8132</td><td>13.70 ± 1.97</td></tr><tr><td>58050.7285</td><td>19.80 ± 1.81</td></tr><tr><td>58053.3030</td><td>22.40 ± 2.14</td></tr><tr><td>58059.6613</td><td>26.70 ± 2.47</td></tr><tr><td>58065.6395</td><td>24.70 ± 2.30</td></tr><tr><td>58068.3633</td><td>18.60 ± 1.81</td></tr><tr><td>58070.0338</td><td>19.40 ± 1.81</td></tr><tr><td>58071.0301</td><td>21.30 ± 2.14</td></tr><tr><td>58072.0950</td><td>22.70 ± 2.30</td></tr><tr><td>58073.0911</td><td>22.30 ± 2.14</td></tr><tr><td>58074.2466</td><td>24.70 ± 2.30</td></tr><tr><td>58075.0830</td><td>25.60 ± 2.47</td></tr><tr><td>58075.6651</td><td>23.70 ± 2.30</td></tr><tr><td>58076.0792</td><td>23.20 ± 2.14</td></tr><tr><td>58077.0751</td><td>23.80 ± 2.30</td></tr><tr><td>58078.0713</td><td>26.10 ± 2.47</td></tr><tr><td>58079.0670</td><td>24.70 ± 2.47</td></tr><tr><td>58080.1316</td><td>25.10 ± 2.47</td></tr><tr><td>58081.1274</td><td>24.60 ± 2.30</td></tr><tr><td>58082.0546</td><td>24.30 ± 2.30</td></tr><tr><td>58083.1191</td><td>21.40 ± 2.14</td></tr><tr><td>58084.0496</td><td>27.30 ± 2.63</td></tr><tr><td>58086.1085</td><td>30.10 ± 2.96</td></tr><tr><td>58087.1557</td><td>21.70 ± 2.80</td></tr></tbody></table>

**Note.** $F_{UV,-12}$ is the $\nu F_{\nu}(=\varepsilon\gamma F_{\gamma})$ flux at $\nu = 10^{15}$ Hz in units of $10^{-12}$ erg cm$^{-2}$ s$^{-1}$, derived by SED fitting for six-filter epochs, and adjusted to $\nu = 10^{15}$ Hz using the nearest best-fit SED index for single-filter epochs.

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 $\delta A_V = \pm 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 $3 \times 10^{14}$ Hz. This indicates that TXS 0506+056 is an intermediate synchrotron peaked (ISP) or low synchrotron peaked blazar.

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

## 2.8. Fermi Data

The Fermi-LAT is a pair conversion telescope sensitive to $\gamma$ 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 ($\pm 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.<fnref n="21" />

Photons with energies between 100 MeV and 300 GeV that were detected within a radius of $15^{\circ}$ from the location of TXS 0506+056 were selected for the analysis, while photons with a zenith angle $> 90^{\circ}$ 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^{\circ}$ from the source position were included in the model, with spectral parameters fixed to their catalog values, while spectral parameters for sources within $3^{\circ}$ 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 $\alpha_{LAT} = 2.05 \pm 0.05$, consistent with the 3FGL value of $2.04 \pm 0.03$, and a flux normalization of $(1.42 \pm 0.11) \times 10^{-11} \text{ cm}^{-2} \text{ s}^{-1} \text{ MeV}^{-1}$ at an energy of 1.44 GeV, about a factor of four higher than the 3FGL value of $(3.24 \pm 0.10) \times 10^{-12}$ 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 $\sim 3\sigma$. Flux upper limits at the 95% confidence level are quoted otherwise.

## 3. Multimessenger Modeling

Traditionally, blazar SEDs are interpreted in two different ways. In the leptonic scenario, the $\gamma$-ray component is interpreted as synchrotron self-Compton (SSC) emission or external inverse-Compton (EIC) emission (e.g., Maraschi et al. 1992; Dermer & Schlickeiser 1993; Sikora et al. 1994). In the SSC model, the seed photons for Compton scattering are produced internally in the blazar jet; in particular, these are the synchrotron photons produced by nonthermal electron–positron pairs accelerated in the jet. In the EIC model, the seeds for Compton scattering are provided by external radiation fields, such as scattered accretion disk radiation, broadline/dust

<footnote n="21">²¹ Fermi Science Tools can be downloaded from https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/.</footnote>
emission, and soft radiation from the sheath region of a structured jet.

It is natural that protons and nuclei are also accelerated in the jet, leading to the so-called leptohadronic scenario<fnref n="22" /> where the γ-ray emission is explained by processes related to relativistic protons: proton-induced EM cascades (Mannheim 1993; Mücke et al. 2003), proton synchrotron emission (Aharonian 2000; Mücke et al. 2003), or intergalactic magnetic cascades induced by UHECRs (Essey et al. 2011; Murase et al. 2012). In the presence of relativistic protons, theory has predicted that PeV–EeV neutrinos can be produced via the photomeson production process between cosmic-ray protons and target photons provided by the intra-jet and/or external radiation fields (see Murase 2017 for a recent review on AGN neutrinos and references therein). For example, a neutrino with $\epsilon_\nu \approx 0.1–1$ PeV implies a parent proton with energy $\epsilon_p \approx 2.0–20$ PeV, for which photomeson production mainly occurs with target photons with UV or greater energies.

HE neutrinos generated by photohadronic interactions must be accompanied by EM emission of secondary electron–positron pairs and pionic γ rays. EM cascades redistribute energy from high energies (e.g., PeV) to lower energies (e.g., keV–MeV) and exhibit $F_\gamma \sim F_\nu$. These cascade effects are included in our detailed numerical calculations, as presented in the following sections.

### 3.1. Model Description

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).

<caption>**Table 5**</caption>

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

<table><thead><tr><th>Parameter</th><th>Symbol</th><th>Unit [in cgs]</th></tr></thead><tbody><tr><td>Doppler factor</td><td>δ</td><td>n/a</td></tr><tr><td>Magnetic field strength</td><td>B'</td><td>G</td></tr><tr><td>Blob radius</td><td>R'</td><td>cm</td></tr><tr><td>Electron injection luminosity</td><td>L'<sub>e</sub></td><td>erg s<sup>-1</sup></td></tr><tr><td>Minimum electron Lorentz factor</td><td>γ'<sub>e,min</sub></td><td>n/a</td></tr><tr><td>Maximum electron Lorentz factor</td><td>γ'<sub>e,max</sub></td><td>n/a</td></tr><tr><td>Break electron Lorentz factor</td><td>γ'<sub>e,br</sub></td><td>n/a</td></tr><tr><td>Power-law electron index below the break</td><td>s<sub>e,1</sub></td><td>n/a</td></tr><tr><td>Power-law electron index above the break</td><td>s<sub>e,2</sub></td><td>n/a</td></tr><tr><td>Proton injection luminosity</td><td>L'<sub>p</sub></td><td>erg s<sup>-1</sup></td></tr><tr><td>Minimum proton Lorentz factor</td><td>γ'<sub>p,min</sub></td><td>n/a</td></tr><tr><td>Maximum proton Lorentz factor</td><td>γ'<sub>p,max</sub></td><td>n/a</td></tr><tr><td>Power-law proton index</td><td>s<sub>p</sub></td><td>n/a</td></tr><tr><td>Energy density of external radiation</td><td>u'<sub>ext</sub></td><td>erg cm<sup>-3</sup></td></tr><tr><td>Effective temperature of blackbody external radiation</td><td>T'</td><td>K</td></tr><tr><td>Photon index of power-law external radiation</td><td>α</td><td>n/a</td></tr><tr><td>Minimum photon energy of power-law external radiation</td><td>ε'<sub>min</sub></td><td>keV</td></tr><tr><td>Maximum photon energy of power-law external radiation</td><td>ε'<sub>max</sub></td><td>keV</td></tr></tbody></table>

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

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.

We search for models that adequately describe the multi-wavelength 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.

<footnote n="22">²² We will refer to this scenario simply as hadronic for simplicity.</footnote>
## 3.2. Leptonic Models (LMs)

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 $\gamma'_{p,min} = 1$ to $\gamma'_{p,max} = 1.6 \times 10^7$. From the X-ray and γ-ray light curves, we infer a variability timescale of $t_{var} \lesssim 10^5$ s. Our choice of $R' = 10^{17}$ cm is broadly consistent with the size inferred from the variability, namely $R' \approx \delta ct_{var}/(1+z) \simeq 0.56 \times 10^{17}(\delta/25)(t_{var}/10^5)$ s 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 $T'$ and energy density $u'_{ext}$, 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

<figure><img src="image_9.png" /><figcaption>**Figure 4.** Leptonic model (LMBB2b) for the TXS 0506+056 flare (Ep. 1). Two SED cases (gray lines) are plotted against the observations (colored points, showing allowed ranges at 90% confidence), one with a hadronic component set to the maximum allowed proton luminosity $L_p^{(max)} \approx 2 \times 10^{50}$ erg s$^{-1}$ (solid gray), and the other set to twice this maximal value (dashed gray line). Corresponding all-flavor neutrino fluxes for the maximal (solid red) and “twice maximal” (dashed line) cases are also shown. Photon attenuation at $\epsilon_\gamma \gtrsim 3 \times 10^{11}$ eV due to interactions with the extragalactic background light is not included here. All-flavor neutrino-flux upper limits of producing an event similar to the IceCube-170922A are shown in blue (from Figure 4 of Aartsen et al. 2018a) for 0.5 (solid blue line) and 7.5 years (dashed blue line).</figcaption></figure>

<table><caption>**Table 6**

Parameter Values Common to All Leptonic Models (LMs) for TXS 0506+056</caption><tr><td>B'<sup>*</sup> [G]</td><td>0.4</td></tr><tr><td>R'<sup>*</sup> [in cm]</td><td>10<sup>17</sup></td></tr><tr><td>&delta;</td><td>24.2</td></tr><tr><td>L'<sub>e</sub><sup>*</sup> [in erg s<sup>-1</sup>]</td><td>2.2 &times; 10<sup>42</sup></td></tr><tr><td>s<sub>e,1</sub></td><td>1.9</td></tr><tr><td>s<sub>e,2</sub></td><td>3.6</td></tr><tr><td>&gamma;'<sub>e,min</sub></td><td>1</td></tr><tr><td>&gamma;'<sub>e,b</sub></td><td>5 &times; 10<sup>3</sup></td></tr><tr><td>&gamma;'<sub>e,max</sub></td><td>8 &times; 10<sup>4</sup></td></tr></table>

**Note.** The isotropic-equivalent electron luminosity is $L_e = \delta^4 L'_e$. Parameter definitions are provided in Table 5.

explained by the IC scattering of the fiducial external photon field by the same electron population. The steepening of the γ-ray spectrum at $\gtrsim 10$ GeV 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 ($R' \ll 10^{17}$ 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.

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.
<table><caption>**Table 7**
Model-specific Parameter Values for Leptonic Models (LMs) for TXS 0506+056 Discussed in the Text</caption><thead><tr><th></th><th>LMBB1a</th><th>LMBB1b</th><th>LMBB1c</th><th>LMBB2a</th><th>LMBB2b</th><th>LMBB2c</th><th>LMPL1a</th><th>LMPL1b</th><th>LMPL2a</th><th>LMPL2b</th></tr></thead><tbody><tr><td>L'<sub>p</sub><sup>(max)</sup> [10<sup>44</sup> erg s<sup>-1</sup>]</td><td>0.54</td><td>0.27</td><td>0.34</td><td>1</td><td>5.4</td><td>10</td><td>0.54</td><td>0.54</td><td>10</td><td>10</td></tr><tr><td>s<sub>p</sub></td><td>2</td><td>2.5</td><td>3</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td></tr><tr><td>γ'<sub>p,min</sub></td><td>1</td><td>3 × 10<sup>6</sup></td><td>3 × 10<sup>6</sup></td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>γ'<sub>p,max</sub> [10<sup>8</sup>]</td><td>30</td><td>30</td><td>30</td><td>1.6</td><td>0.16</td><td>0.016</td><td>30</td><td>30</td><td>0.016</td><td>0.016</td></tr><tr><td>u'<sub>ext</sub> [erg cm<sup>-3</sup>]</td><td colspan="6">0.033</td><td>0.033</td><td>0.067</td><td>0.04</td><td>0.08</td></tr><tr><td>T' [K]</td><td colspan="6">3 × 10<sup>5</sup></td><td colspan="4">n/a</td></tr><tr><td>α</td><td colspan="6">n/a</td><td>3</td><td>2</td><td>3</td><td>2</td></tr><tr><td>ε'<sub>min</sub> [keV]</td><td colspan="6">n/a</td><td colspan="4">0.05</td></tr><tr><td>ε'<sub>max</sub> [keV]</td><td colspan="6">n/a</td><td colspan="4">5</td></tr></tbody></table>

**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 $T'$, while in LMPL models, it is a power law between comoving energies $\epsilon'_{min}$ and $\epsilon'_{max}$, with photon index $\alpha$. In all cases, $u'_{ext}$ is the comoving energy density of the external photon field. Note that the isotropic-equivalent cosmic-ray proton luminosity is $L_p = \delta^4 L'_p$.

The photon spectrum obtained with $L_p = 2 \times L_p^{(max)}$ 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 $\epsilon_\nu F_{\epsilon_\nu}^{(max)} \sim (2-3) \times 10^{-12} \text{ erg cm}^{-2} \text{ s}^{-1}$.

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 $3 \times 10^{18} \text{ 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 $\gamma'_{p,max} = 1.6 \times 10^8$, $1.6 \times 10^9$, and $3 \times 10^9$—see Table 7. Our results on the neutrino fluxes are presented in Figure 5.

Neutrino spectra in the LMBB1x models are more extended in energy compared to the default case (LMBB2b). They peak around $10 \text{ PeV}$ (100 PeV) for $\gamma'_{p,max} = 1.6 \times 10^7$ ($1.6 \times 10^8$) 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., $\gamma'_{p,max} = 1.6 \times 10^6$, should not lead to higher upper limits on the neutrino flux. The reason is that protons with $\gamma'_{p,max} \sim 10^6$ 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 ($\nu \sim 10^{18} \text{ 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 $s_p = 2.5$ and $s_p = 3$. For particle distributions with soft spectra (i.e., $s_p > 2$), the total energy in protons is determined by the low-energy cutoff ($\gamma'_{p,min}$) of the distribution. These low-energy protons, however, are of no interest for HE neutrino production. In an

<figure><img src="image_10.png" /><figcaption>**Figure 5.** Upper limits on the all-flavor neutrino ($\nu + \bar{\nu}$) fluxes predicted for our modeling of the SED in the leptonic (LMx) and hadronic (Hmx) models. The dashed blue line shows the value of producing an event similar to the IceCube-170922A in 7.5 years (Aartsen et al. 2018a) of data.</figcaption></figure>

attempt to minimize the energy budget, while retaining the HE neutrino fluxes for $s_p > 2$, one has to assume $\gamma'_{p,min} \gg 1$—see Table 7. The large $\gamma'_{p,min}$ 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., $\gamma'_{p,min} \rightarrow \gamma'_{p,br}$). 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 $\gamma$-ray range (Tavecchio et al. 2014; Tavecchio & Ghisellini 2015). From Figure 5, we see that the results for LMPL1x models with $\gamma'_{p,max} = 3 \times 10^9$ 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 $\gamma'_{p,max} = 1.6 \times 10^7$ 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 $\varepsilon_\nu F_{\varepsilon_\nu} \sim 10^{-12} \text{ erg cm}^{-2} \text{ s}^{-1}$, similar to that for LMBB2b. We thus conclude that the neutrino-flux upper limit is $\varepsilon_\nu F_{\varepsilon_\nu}^{(\text{max})} \sim a \times 10^{-12} \text{ erg cm}^{-2} \text{ s}^{-1}$ 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 $\gamma$ 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 $\gamma$-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 $\gamma$ rays. The maximal all-flavor neutrino flux over $0.1 \text{ PeV} < \varepsilon_\nu < 10 \text{ PeV}$ is $\varepsilon_\nu F_{\varepsilon_\nu}^{(\text{max})} \approx 3.6 \times 10^{-12} \text{ erg cm}^{-2} \text{ s}^{-1}$, implying a Poisson probability to detect one event with IceCube over the six-month duration of the TXS 0506+056 $\gamma$-ray flare of at most $\sim 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 $L_p^{(\text{max})} \approx 2 \times 10^{50} \text{ erg s}^{-1}$. Cases with proton luminosities exceeding $L_p^{(\text{max})}$ 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.

### 3.3. Hadronic Models (HMs)

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 $\gamma$-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

<figure><img src="image_11.png" /><figcaption>**Figure 6.** Hadronic model (HM3) for the SED of TXS 0506+056 flare (Ep. 1), as computed for different values of the Doppler factor (gray curves), together with resulting all-flavor neutrino fluxes (red curves) and electromagnetic observations (colored points, showing allowed ranges at 90% confidence). Photon attenuation at $\varepsilon_\gamma \gtrsim 3 \times 10^{11} \text{ eV}$ due to interactions with the extragalactic background light is not included here. All-flavor neutrino-flux upper limits of producing an event similar to the IceCube-170922A are shown in blue (from Figure 4 of Aartsen et al. 2018a) for 0.5 (solid blue line) and 7.5 years (dashed blue line).</figcaption></figure>

<table><caption>**Table 8**

Parameter Values for Hadronic Models (HMs) for TXS 0506+056 Discussed in the Text and Presented in Figure 6</caption><thead><tr><th></th><th>HM1</th><th>HM2</th><th>HM3</th></tr></thead><tbody><tr><td>B' [G]</td><td></td><td>85</td><td></td></tr><tr><td>R' [in 10<sup>16</sup>cm]</td><td>2</td><td>3</td><td>4.5</td></tr><tr><td>δ</td><td>5.2</td><td>10</td><td>15</td></tr><tr><td>L'<sub>e</sub> [in 10<sup>43</sup> erg s<sup>-1</sup>] 9.3</td><td>0.6</td><td>0.06</td><td></td></tr><tr><td>s<sub>e,1</sub></td><td></td><td>1.8</td><td></td></tr><tr><td>s<sub>e,2</sub></td><td>4.2</td><td>3.6</td><td>3.6</td></tr><tr><td>γ'<sub>e,min</sub> [in 10<sup>2</sup>]</td><td>6.3</td><td>1</td><td>1</td></tr><tr><td>γ'<sub>e,br</sub> [in 10<sup>2</sup>]</td><td>7.9</td><td>6.3</td><td>5</td></tr><tr><td>γ'<sub>e,max</sub></td><td></td><td>10<sup>4</sup></td><td></td></tr><tr><td>L'<sub>p</sub> [in 10<sup>46</sup> erg s<sup>-1</sup>]</td><td>2.7</td><td>0.1</td><td>0.01</td></tr><tr><td>s<sub>p</sub></td><td></td><td>2.1</td><td></td></tr><tr><td>γ'<sub>p,min</sub></td><td></td><td>1</td><td></td></tr><tr><td>γ'<sub>p,max</sub></td><td></td><td>2 × 10<sup>9</sup></td><td></td></tr></tbody></table>

**Note.** Parameter definitions are provided in Table 5.

et al. 2013; Cerruti et al. 2015). As a starting point, we search for combinations of $\delta$ 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: $\delta_{eq} \sim 5$, $B'_{eq} \sim 80 \text{ G}$, $R'_{eq} \sim 10^{16} \text{ cm}$, and $\varepsilon'_{p,max} \sim 10^9 \text{ 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
$L_j \approx 2\pi c R'^2 (\delta/2)^2 (u'_p + u'_e + u'_B) \sim 4 \times 10^{47}$ erg s$^{-1}$, with $u'_p \approx 2u'_B \sim 500$ erg cm$^{-3}$, where $u'_p, u'_e, u'_B$ 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 $\delta = \delta_{eq}$ (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.

The EM cascade emission can be suppressed if the source becomes less opaque to the intra-source $\gamma\gamma$ absorption for HE photons. This can be achieved for larger values of the Doppler factor since $\tau_{\gamma\gamma} \propto \delta^{-4}$ (see also Murase et al. 2016; Petropoulou et al. 2017, for analytical expressions). The photon and neutrino spectra for $\delta = 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 $\delta \ge 3\delta_{eq}$ (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 $\gamma$ rays in the Fermi-LAT band (Murase et al. 2016). Furthermore, this model unavoidably leads to a higher jet power, i.e., $L_j \sim 6 \times 10^{48}$ erg s$^{-1}$ and $u'_B \gg u'_p$ (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 $\epsilon_p \gtrsim 3$ 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 $\epsilon_\gamma F_{\epsilon_\gamma} \sim \epsilon_\nu F_{\epsilon_\nu}$ unavoidably overshoots the observed X-ray flux, giving $\epsilon_\gamma F_{\epsilon_\gamma} \approx 8 \times 10^{-12}$ erg cm$^{-2}$ s$^{-1}$, which is strongly excluded. Alternatively, the proton synchrotron model can explain the TXS 0506+056 $\gamma$-ray emission, but gives a maximal neutrino flux $\epsilon_\nu F_{\epsilon_\nu}^{(\text{max})} \approx 2 \times 10^{-15}$ erg cm$^{-2}$ s$^{-1}$, which implies a very low probability for IceCube neutrino detection, $p_{\text{IC}} < 10^{-5}$. 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.

### 3.4. Implications of IceCube-170922A

Relativistic protons of energy $\epsilon_p$ can interact with photons in the source and produce neutrinos with energy $\sim \epsilon_p/20$.<fnref n="23" /> 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 $\Delta$ resonance is usually valid (Murase & Nagataki 2006; Murase et al. 2014), and the characteristic proton Lorentz factor interacting with photons of the frequency $\nu_{\text{syn}} = \epsilon_{\text{syn}}/h$ is given by (Murase 2017)

$$ \begin{aligned} \gamma'_{p,b} &\approx 0.5 \delta \bar{\epsilon}_{\Delta} (\epsilon_{\text{syn}})^{-1} (1+z)^{-1} \\ &\approx 1.3 \times 10^9 \delta_1 \nu_{\text{syn},14.5}^{-1} (1+z)^{-1}, \end{aligned} \quad (1) $$

where $\bar{\epsilon}_{\Delta} \sim 0.3$ GeV is the resonance energy. The respective neutrino energy, in the observer's frame, is then given by $\epsilon_{\nu,b} \approx 0.05 \gamma'_{p,b} \delta m_p c^2/(1+z)$. This is also an estimate of the peak energy of the HE neutrino spectrum:

$$ \epsilon_{\nu,b} \approx 0.05 \epsilon_p^b \approx 600 \text{ PeV} \delta_1^2 \nu_{\text{syn},14.5}^{-1} (1+z)^{-2}. \quad (2) $$

Neutrinos with lower energies than $\epsilon_{\nu,b}$ 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 $\epsilon_{\nu,b}$ due to the decreasing number density of target photons, unless the proton distribution is a super-soft power law (i.e., $s_p \gg 2$). 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 $\gamma'_{p} \sim \gamma'_{p,\text{max}}$ 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 $\gamma$-ray data is primarily explained by the proton-induced cascade emission itself, is strongly ruled out. Although this model naturally predicts $\epsilon_\nu F_{\epsilon_\nu} \sim \epsilon_\gamma F_{\epsilon_\gamma} \sim 10^{-10}$ erg cm$^{-2}$ s$^{-1}$, 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, $\epsilon_\nu F_{\epsilon_\nu} \sim 10^{-15}$ erg cm$^{-2}$ s$^{-1}$.

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, $\epsilon_\nu F_{\epsilon_\nu} \sim (2-3) \times 10^{-12}$ erg cm$^{-2}$ s$^{-1}$, the corresponding muon neutrino fluence is estimated to be:

$$ \epsilon_{\nu}^{2} \phi_{\nu_{\mu}}^{(\text{max})} \approx 1.6 \times 10^{-5} \text{ erg cm}^{-2} \\ \times \left( \frac{\epsilon_{\nu} F_{\epsilon_{\nu}}^{(\text{max})}}{3 \times 10^{-12} \text{ erg cm}^{-2} \text{ s}^{-1}} \right) \left( \frac{\Delta T}{0.5 \text{ year}} \right), \quad (3) $$

<footnote n="23">²³ Hadronuclear reactions such as proton–proton collisions in blazar jets are expected to be too inefficient compared to the photomeson production process, for typical values of the jets’ plasma density (Atoyan & Dermer 2003; Murase et al. 2014). Also, even ad hoc high-density environments are similarly (~order of magnitude) constrained by the cascade bound.</footnote>
<caption>Table 9</caption>

The Upper Limit on the All-flavor Neutrino Flux $F_{\nu}^{(\text{max})}$ for the Different Models that Satisfactorily Explain the Flaring SED of TXS 0506+056—for Details, See Sections 3.2 and 3.3

<table><thead><tr><th></th><th colspan="2">∫F<sub>εν</sub><sup>(max)</sup> dε<sub>ν</sub> [erg cm<sup>-2</sup> s<sup>-1</sup>]</th><th>N<sub>ν</sub></th></tr></thead><tbody><tr><td>LMBB1a</td><td>100 TeV-1 PeV</td><td>100 TeV-10 PeV</td><td>≤10 PeV</td></tr><tr><td>LMBB1b</td><td>5.2 × 10<sup>-14</sup></td><td>4.5 × 10<sup>-13</sup></td><td>1 × 10<sup>-3</sup></td></tr><tr><td>LMBB1c</td><td>9.1 × 10<sup>-14</sup></td><td>2.7 × 10<sup>-12</sup></td><td>4 × 10<sup>-3</sup></td></tr><tr><td>LMBB2a</td><td>4.5 × 10<sup>-14</sup></td><td>1.1 × 10<sup>-12</sup></td><td>6 × 10<sup>-3</sup></td></tr><tr><td>LMBB2b</td><td>1.8 × 10<sup>-13</sup></td><td>3.6 × 10<sup>-12</sup></td><td>3 × 10<sup>-3</sup></td></tr><tr><td>LMBB2c</td><td>2.5 × 10<sup>-14</sup></td><td>7.3 × 10<sup>-14</sup></td><td>8 × 10<sup>-3</sup></td></tr><tr><td>LMPL1a</td><td>3.1 × 10<sup>-14</sup></td><td>5.2 × 10<sup>-13</sup></td><td>2 × 10<sup>-4</sup></td></tr><tr><td>LMPL1b</td><td>9 × 10<sup>-14</sup></td><td>6.3 × 10<sup>-13</sup></td><td>1 × 10<sup>-3</sup></td></tr><tr><td>LMPL2a</td><td>2.5 × 10<sup>-13</sup></td><td>5.2 × 10<sup>-13</sup></td><td>5 × 10<sup>-3</sup></td></tr><tr><td>LMPL2b</td><td>1.2 × 10<sup>-12</sup></td><td>2 × 10<sup>-12</sup></td><td>1 × 10<sup>-2</sup></td></tr><tr><td>HM3</td><td>1.6 × 10<sup>-16</sup></td><td>2 × 10<sup>-15</sup></td><td>4 × 10<sup>-6</sup></td></tr></tbody></table>

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 $\Delta T = 0.5$ year.

where $\Delta T \sim 0.5$ year is the flare duration and the flavor ratio is assumed to be $\nu_e : \nu_\mu : \nu_\tau \approx 1 : 1 : 1$. Then, using the effective area for EHE real-time alerts, $A_{\text{eff}} \sim 10^6 \text{ cm}^2$ in the PeV range (Aartsen et al. 2017a), the expected number of muon neutrinos is estimated to be:

$$ \mathcal{N}_{\nu} \sim (\varepsilon_{\nu} \phi_{\nu_{\mu}}^{(\text{max})}) \Delta \varepsilon_{\nu} A_{\text{eff}} \approx 0.02 \times \left( \frac{\varepsilon_{\nu} F_{\varepsilon_{\nu}}^{(\text{max})}}{3 \times 10^{-12} \text{ erg cm}^{-2} \text{ s}^{-1}} \right) \left( \frac{\Delta T}{0.5 \text{ year}} \right) \left( \frac{A_{\text{eff}}}{10^6 \text{ cm}^2} \right), \quad (4) $$

where the width of the neutrino spectrum is assumed to be $\Delta\varepsilon_\nu \sim \ln(10)$. 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 $\Delta 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 $p_{\text{IC}} \sim 1\%$ 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 $2L_p^{(\text{max})}$ as the most conservative limit, which gives $p_{\text{IC}} < 2\%$.

## 4. Discussion

The contribution of blazar jets to the diffuse neutrino flux has been calculated based on both leptonic and hadronic scenarios for the observed γ-rays, under the common assumption of a constant neutrino flux (Murase et al. 2014; Tavecchio et al. 2014; Padovani et al. 2015; Tavecchio & Ghisellini 2015; Rodrigues et al. 2018). Most of the most optimistic scenarios for the diffuse neutrino emission from blazars (for a review, see Murase 2017) can now be constrained by IceCube (e.g., Aartsen et al. 2016a, 2016b, 2017b). In addition, all of the model-independent analyses (stacking, multiplet, and auto-correlation analyses) have disfavored the blazar population as the dominant (~100%) origin of IceCube’s neutrinos, implying that their contribution is less than ~3%–30% of the diffuse neutrino intensity in the 0.1–1 PeV range (Aartsen et al. 2015, 2017b, 2017c; Murase & Waxman 2016). Importantly, even if the blazar contribution to the diffuse neutrino flux is subdominant, the flaring blazar-associated neutrinos are detectable with the current detector since the atmospheric backgrounds can be reduced by temporal and spatial coincidence. Although blazars do not exhibit the extreme variability of nonrepeating transients like gamma-ray bursts and supernovae, they are highly variable on different timescales and across the EM spectrum (Kataoka et al. 2001; Abdo et al. 2010; González-Martín & Vaughan 2012; Sobolewska et al. 2014), and have long been one of the most promising possibilities for HE neutrino production (Atoyan & Dermer 2001; Dermer et al. 2012, 2014; Halzen & Kheirandish 2016; Kadler et al. 2016; Petropoulou et al. 2016; Turley et al. 2016; Gao et al. 2017).

We presented the time-averaged data for Ep. 1 and Ep. 2 and used the former for multimessenger modeling of the source SED. Although we have argued that this is a reasonable approach, it has some limitations. In particular, we cannot exclude the possibility that physical conditions of the source change drastically on a short timescales. Although there is currently no evidence of such rapid variability, this could affect the results of our time-average SED modeling. Nevertheless, the EM cascade effects are inevitable and the resulting X-ray component must appear. We thus expect that our conclusions are robust, since the upper limits on neutrino fluence are basically set by the observed X-ray fluence.

We modeled the SED for Ep. 1 with a 30 day window. The X-ray flux varies by a factor of 6 within 60 days, but the time-average SED for Ep. 2 does not differ much from that for Ep. 1. This indicates that the physical conditions are approximately similar during Ep. 1 and Ep. 2, which justifies the setup of our calculations. The flare duration that is relevant for estimates of the signal neutrino fluence is at least $\Delta T \sim 60$ days and the γ-ray data suggest that a duration of $\Delta T \sim 0.5–1$ year is possible.

As noted above, it is the neutrino fluence (i.e., the product of the duration of neutrino emission and the neutrino flux) that matters in the calculation of the expected number of events. It is therefore likely that the expected number of events from a nonflaring blazar integrated over the IceCube lifetime is larger than the one expected from a flare. In our optimal case, if the SED shape and X-ray flux in the steady state remain similar to those in the flaring state, the expectation value of the number of muons that can be found in the eight-year point-source analysis on upgoing muons is ~1 event. Thus, time-averaged X-ray fluxes had to be higher in the past to obtain ~13
neutrino events (Aartsen et al. 2018b). To properly address this question, one would need to have a good description of the nonflaring SED of the source, especially in the X-ray range. We plan to compare the flaring emission with the nonflaring emission of TXS 0506+056 in a dedicated future work (M. Petropoulou et al. 2018 in preparation).

Although our results on the neutrino-flux upper limit in the LMs are insensitive to details of the external photon field, we briefly discuss possible origins of the external photons. The typical photon energy of the external radiation field in the black hole rest frame is $\epsilon_{\gamma,ext} = 3k_B T'/\Gamma = 2.5 T'_{5.5}/\Gamma_{1.5}$ eV and its energy density is $u_{ext} \approx u'_{ext}/\Gamma^2 \simeq 3.3 \times 10^{-5} \Gamma_{1.5}^{-2} \text{ erg cm}^{-3}$ (see Table 7). The putative external photon field is compatible with scattered disk emission or soft emission from the sheath region of the blazar jet (Dermer et al. 2014; Tavecchio & Ghisellini 2015). Any additional external component should not exceed the observed fluxes in the optical, UV, and X-ray data. Thus, its luminosity should be $L_{ext} \lesssim 10^{45}-10^{46} \text{ erg s}^{-1}$, depending on the spectral shape.

In the scattered disk emission case, if there is a scattering region with the Thomson optical depth $\tau_T$ at radius $R$, the energy density of the scattered emission can be $u_{ext} \lesssim (3-30) \times 10^{-5}(\tau_T/0.1)(R/3 \times 10^{18} \text{ cm})^{-2} \text{ erg cm}^{-3}$. Alternatively, dissipation in the sheath region of the jet may lead to electron acceleration and associated synchrotron emission with a peak energy of $\epsilon_{\gamma,ext} \sim 20 (\Gamma_s/2)^2 \Gamma_{1.5}^{-1} \text{ eV}$ and luminosity of $L_{ext} \sim 7 \times 10^{45} (R/3 \times 10^{18} \text{ cm})^2 \Gamma_{1.5}^{-2} (\Gamma_s/2)^4 \text{ erg s}^{-1}$, where $\Gamma_s$ is the Lorentz factor of the sheath region.

Third, external photons can be provided by the possible broadline region, and the energy density of the broadline region can be written as $u_{BLR} \approx 0.26 f_{cov} \text{ erg cm}^{-3}$, where $f_{cov}$ is the covering fraction (Ghisellini & Tavecchio 2008). However, there are two drawbacks. The lack of broadline signatures in the optical spectrum of TXS 0506+056 and other BL Lac objects suggests that such line emissions are weak.

Also, such emissions will only be important when the blob is located in the broadline region. If we follow the treatment in Murase et al. (2014), the observed $\gamma$-ray luminosity indicates that the broadline region is located at $R \lesssim 10^{16} \text{ cm}$, so that the typical emission radius would be larger than the radius of the possible broadline region.

Even though blazars like TXS 0506+056 could make a significant contribution to the diffuse neutrino flux, it is premature to extrapolate our findings from the modeling of the TXS 0506+056 flare to other blazars. A more dedicated study is left for future work, after a representative set of flares have been modeled individually (see also Padovani et al. 2015). So far, the only other flare that was modeled with the same numerical code is the 13-day (2010) flare of the high-frequency peaked blazar Mrk 421 (Petropoulou et al. 2016). In that case, the available EM data could be explained with a higher neutrino-to-$\gamma$-ray luminosity ratio than we find here. If we were to take the neutrino-to-$\gamma$-ray luminosity ratio for the TXS 0506+056 flare as representative of flaring blazars, their contribution to the diffuse neutrino flux would be smaller than previous estimates.

Finally, as we emphasized, our theoretical interpretation is based on the single-zone model. In this assumption, the broadband EM and HE neutrino emissions are produced in the same localized region in the blazar jet. Efficient HE neutrino production requires large optical depths for photomeson production. This, in turn, implies that the emitting region is

optically thick for HE $\gamma$-rays (Murase et al. 2016). This is especially the case for blazars that have soft photon spectra, in which the EM energy of the attenuated very HE photons will reappear at lower energies through an EM cascade. The tight constraints stemming from the large optical depth of a localized region could be alleviated, if the photon and neutrino emissions originate from different regions in the jet, as in multizone models (Dermer et al. 2012, 2014; Murase et al. 2014). Although a detailed discussion must be deferred to future work (Murase et al. 2018), the constraints from the EM cascade are unavoidable even for such multizone models, so that further ad hoc adjustments of the source parameters seem necessary.

## 5. Conclusions

We have used the best available multiwavelength data to construct the broadband SED of TXS 0506+056, over 10 orders of magnitude in photon energy from $\approx 1$ eV to $>10$ GeV, proximate to its likely ($\sim 3\sigma$; Aartsen et al. 2018a) emission of the high-energy neutrino IceCube-170922A. Working with this SED and the likely neutrino association, we have explored multimessenger models for TXS 0506+056 to evaluate whether this neutrino association is physically reasonable, and if so, under what conditions of the blazar jet and jet environment.

We find that a leptonic scenario with a radiatively subdominant hadronic component provides the only physically consistent single-zone picture for this source's multimessenger (EM and neutrino) emissions. If IceCube-170922A is associated with this flaring blazar, then physical conditions were close to optimal for neutrino production during its flare. We find a maximal all-flavor neutrino flux over $0.1 \text{ PeV} < \epsilon_{\nu} < 10 \text{ PeV}$ of $\epsilon_{\nu} F_{\epsilon_{\nu}}^{(\text{max})} \approx (2-4) \times 10^{-12} \text{ erg cm}^{-2} \text{ s}^{-1}$. The inferred ratio of proton to electron luminosities is large, $L_p/L_e \sim 250-500$, with smaller values prohibited because increased optical depth to $\gamma\gamma \rightarrow e^{+}e^{-}$ would suppress the observed $\gamma$-rays.

Under these conditions, we find a probability of $p_{\text{IC}} \approx 1\%$ to 2% for IceCube to detect an HE muon neutrino in real time at some point during the blazar's six-month flare.

Since the blazar's $>\text{GeV}$ emissions are dominated by leptonic processes, and since EM cascades efficiently reilibute hadronic EM emissions across the spectrum, we find that the SED exhibits its greatest sensitivity to hadronic acceleration processes across its $0.1-100 \text{ keV}$ “dip.” Flux variations over this energy range are more likely to reflect the source’s high-energy neutrino emissions than its GeV–TeV flux state. We thus find that, going forward, regular X-ray monitoring of TXS 0506+056 and related blazars, in conjunction with continued monitoring by high-energy neutrino observatories, will provide a critical test of single-zone blazar models. Moreover, careful selection of temporal acceptance windows via X-ray observations, as in Turley et al. (2016), will likely yield the most sensitive search for further multimessenger sources.

Finally, we find that under the observed flaring conditions, assuming the IceCube-170922A association holds, TXS 0506+056 was not a significant UHECR accelerator. This is because a proton spectrum extending to $\epsilon_p \gtrsim 3 \text{ EeV}$ would yield a neutrino spectrum peaking above $100 \text{ PeV}$ (Murase et al. 2014); with the neutrino peak flux bounded by X-ray observations via cascade effects, this would strongly suppress the $0.1-10 \text{ PeV}$ neutrino flux (Figure 5).
It is possible that multizone models, which more readily decouple blazar EM and neutrino emissions, may ultimately be required to explain multimessenger observations of TXS 0506 +056 and other blazars. Independent of whether this particular source association holds, our results demonstrate that the detection of even one or two coincident neutrinos can grant us deep insight into a source, and should energize future searches for further multimessenger sources.

We note that the neutrino upper flux data points in Figures 4–6 have been added after the completion of the modeling using data points from Aartsen et al. (2018a).

The authors thank the referee for valuable comments that helped to improve the manuscript. The authors thank the IceCube Collaboration for publicly distributing HE neutrinos in real-time, the *Swift* and *NuSTAR* teams for their rapid responses to our target of opportunity requests, and the *Fermi* collaboration for their publicly available data and analysis software. The authors acknowledge helpful discussions with L. Hagen, J. Charlton, and M. Eracleous. A.K., D.B.F., J.J.D., and C.F.T. acknowledge support from the National Science Foundation under grant PHY-1412633; A.K. acknowledges support from the National Aeronautics and Space Administration *Swift* Guest Investigator Program under grant NNX17AI95G. The work of K.M. is supported by NSF Grant No. PHY-1620777 and the Alfred P. Sloan Foundation. A.K. and K.M. gratefully acknowledge support from the Institute for Gravitation and the Cosmos at the Pennsylvania State University. J.K. and A.T. acknowledge support from NASA contract NAS5-00136. P.A.E. and J.P.O. acknowledge support from the UK Space Agency. M.P. acknowledges support from the L. Jr. Spitzer Postdoctoral Fellowship. S.C. thanks the Centre National d'Etudes Spatiales (CNES) for support and funding. Based on observations collected at the European Organization for Astronomical Research in the Southern Hemisphere under ESO program 099.D-0640(A).

Correspondence regarding this work should be sent to Kohta Murase (murase@psu.edu) or Maria Petropoulou (m.petropoulou@astro.princeton.edu) for theoretical interpretations, and Azadeh Keivani (ak4279@columbia.edu) or Derek Fox (dfox@psu.edu) for observations.

Software: XSPEC (Arnaud 1996), HEASoft (v6.22.1), X-shooter pipeline (v.2.9.3; Goldoni et al. 2006; Modigliani et al. 2010), Reflex (Freudling et al. 2013), Molecfit (Kausch et al. 2015; Smette et al. 2015), Fermi Science Tools (https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/), SOPHIA (Mücke et al. 2000).

## ORCID iDs

A. Keivani https://orcid.org/0000-0001-7197-2788

K. Murase https://orcid.org/0000-0002-5358-5642

M. Petropoulou https://orcid.org/0000-0001-6640-0179

D. B. Fox https://orcid.org/0000-0002-3714-672X

S. B. Cenko https://orcid.org/0000-0003-1673-970X

S. Chaty https://orcid.org/0000-0002-5769-8601

A. Coleiro https://orcid.org/0000-0003-0860-440X

J. J. DeLaunay https://orcid.org/0000-0001-5229-1995

S. Dimitrakoudis https://orcid.org/0000-0002-3368-3739

P. A. Evans https://orcid.org/0000-0002-8465-3353

J. A. Kennea https://orcid.org/0000-0002-6745-4790

F. E. Marshall https://orcid.org/0000-0002-3559-6305

J. P. Osborne https://orcid.org/0000-0002-1041-7542

M. Santander https://orcid.org/0000-0001-7297-8217

A. Tohuvavohu https://orcid.org/0000-0002-2810-8764

C. F. Turley https://orcid.org/0000-0002-9689-8075

## References

Aartsen, M. G., Abbasi, R., Abdou, Y., et al. 2013a, *Sci*, 342, 1242856

Aartsen, M. G., Abbasi, R., Abdou, Y., et al. 2013b, *PhRvL*, 111, 021103

Aartsen, M. G., Abraham, K., Ackermann, M., et al. 2016a, *PhRvL*, 117, 241101

Aartsen, M. G., Abraham, K., Ackermann, M., et al. 2016b, *ApJ*, 833, 3

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2015, *ApH*, 66, 39

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2017a, *ApH*, 92, 30

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2017b, *arXiv*:1710.01179

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2017c, *arXiv*:1710.01191

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2018a, *Sci*, 361, eaat1378

Aartsen, M. G., Ackermann, M., Adams, J., et al. 2018b, *Sci*, 361, eaat2890

Abbey, T., Carpenter, J., Read, A., Wells, A. & Xmm Science Centre, & Swift Mission Operations Center 2006, in ESA Special Publication 604, The X-ray Universe 2005, ed. A. Wilson (Paris: ESA), 943

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, *ApJ*, 722, 520

Ackermann, M., Ajello, M., Albert, A., et al. 2015, *ApJ*, 799, 86

Ackermann, M., Ajello, M., Albert, A., et al. 2016, *PhRvL*, 116, 151105

Aharonian, F. A. 2000, *NewA*, 5, 377

Arnaud, K. A. 1996, in ASP Conf. Ser. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes (San Francisco, CA: ASP), 17

Atoyan, A., & Dermer, C. D. 2001, *PhRvL*, 87, 221102

Atoyan, A. M., & Dermer, C. D. 2003, *ApJ*, 586, 79

Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, *ApJ*, 697, 1071

Boettcher, M., Reimer, A., Sweeney, K., & Prakash, A. 2013, *ApJ*, 768, 54

Böttcher, M., Reimer, A., Sweeney, K., & Prakash, A. 2013, *ApJ*, 768, 54

Breeveld, A. A., Curran, P. A., Hoversten, E. A., et al. 2010, *MNRAS*, 406, 1687

Breeveld, A. A., Landsman, W., Holland, S. T., et al. 2011, in AIP Conf. Ser. 1358, Gamma Ray Bursts 2010, ed. J. E. McEnery, J. L. Racusin, & N. Gehrels (Melville, NY: AIP), 373

Cerruti, M., Zech, A., Boisson, C., & Inoue, S. 2015, *MNRAS*, 448, 910

de Neurois, M. & H.E.S.S. Collaboration 2017, *ATel*, 10787, 1

Dermer, C. D., Murase, K., & Inoue, Y. 2014, *JHEAp*, 3, 29

Dermer, C. D., Murase, K., & Takami, H. 2012, *ApJ*, 755, 147

Dermer, C. D., & Schlickeiser, R. 1993, *ApJ*, 416, 458

Dimitrakoudis, S., Mastichiadis, A., Protheroe, R. J., & Reimer, A. 2012, *A&A*, 546, A120

Dimitrakoudis, S., Petropoulou, M., & Mastichiadis, A. 2014, *ApH*, 54, 61

Essey, W., Kalashev, O., Kusenko, A., & Beacom, J. F. 2011, *ApJ*, 731, 51

Evans, P. A., Beardmore, A. P., Page, K. L., et al. 2009, *MNRAS*, 397, 1177

Evans, P. A., Keivani, A., Kennea, J. A., et al. 2017, *ATel*, 10792, 1

Evans, P. A., Osborne, J. P., Beardmore, A. P., et al. 2014, *ApJS*, 210, 8

Fang, K., & Murase, K. 2018, *NatPh*, 14, 396

Fitzpatrick, E. L. 1999, *PASP*, 111, 63

Fox, D. B., DeLaunay, J. J., Keivani, A., et al. 2017, *ATel*, 10845, 1

Freudling, W., Romaniello, M., Bramich, D. M., et al. 2013, *A&A*, 559, A96

Gao, S., Pohl, M., & Winter, W. 2017, *ApJ*, 843, 109

GCN/AMON NOTICE IceCube-170922A, 2017, IceCube-170922A, https://gcn.gsfc.nasa.gov/notices_amon/50579430_130033.amon

GCN/AMON NOTICE, 2016, Webpage, https://gcn.gsfc.nasa.gov/amon.html

Ghisellini, G., & Tavecchio, F. 2008, *MNRAS*, 387, 1669

Goldoni, P., Royer, F., François, P., et al. 2006, *Proc. SPIE*, 6269, 62692K

González-Martín, O., & Vaughan, S. 2012, *A&A*, 544, A80

Green, G. M., Schlaflay, E. F., Finkbeiner, D. P., et al. 2015, *ApJ*, 810, 25

Halpern, J. P., Eracleous, M., & Mattox, J. R. 2003, *AJ*, 125, 572

Halzen, F., & Kheirandish, A. 2016, *ApJ*, 831, 12

Harrison, F. A., Craig, W. W., Christensen, F. E., et al. 2013, *ApJ*, 770, 103

Kadler, M., Krauß, F., Mannheim, K., et al. 2016, *NatPh*, 12, 807

Kataoka, J., Takahashi, T., Wagner, S. J., et al. 2001, *ApJ*, 560, 659

Kausch, W., Noll, S., Smette, A., et al. 2015, *A&A*, 576, A78

Keel, W., & Santander, M. 2017, *ATel*, 10831, 1

Keivani, A., Evans, P. A., Kennea, J. A., et al. 2017, GCN Circular Service, 21930, 1

Kirk, J. G., Rieger, F. M., & Mastichiadis, A. 1998, *A&A*, 333, 452

Kopper, C., & Blaufuss, E. 2017, GCN Circular Service, 21916, 1

Li, H., & Kusunose, M. 2000, *ApJ*, 536, 729

Mannheim, K. 1993, *A&A*, 269, 67

Mannheim, K., Biermann, P. L., & Kruells, W. M. 1991, *A&A*, 251, 723
Maraschi, L., Ghisellini, G., & Celotti, A. 1992, ApJL, 397, L5

Mastichiadis, A., & Kirk, J. G. 1997, A&A, 320, 19

Mastichiadis, A., Protheroe, R. J., & Kirk, J. G. 2005, A&A, 433, 765

Mirzoyan, R. 2017, ATel, 10817, 1

Modigliani, A., Goldoni, P., Royer, F., et al. 2010, Proc. SPIE, 7737, 773728

Mücke, A., Engel, R., Rachen, J. P., Protheroe, R. J., & Stanev, T. 2000, CoPhC, 124, 290

Mücke, A., Protheroe, R. J., Engel, R., Rachen, J. P., & Stanev, T. 2003, APh, 18, 593

Mukherjee, R. 2017, ATel, 10833, 1

Murase, K. 2017, arXiv:1511.01590

Murase, K., Ahlers, M., & Lacki, B. C. 2013, PhRvD, 88, 121301

Murase, K., Dermer, C. D., Takami, H., & Migliori, G. 2012, ApJ, 749, 63

Murase, K., Guetta, D., & Ahlers, M. 2016, PhRvL, 116, 071101

Murase, K., Inoue, Y., & Dermer, C. D. 2014, PhRvD, 90, 023007

Murase, K., & Nagataki, S. 2006, PhRvD, 73, 063002

Murase, K., Oikonomou, F., & Petropoulou, M. 2018, ApJ, in press (arXiv:1807.04748)

Murase, K., & Waxman, E. 2016, PhRvD, 94, 103006

Padovani, P., Petropoulou, M., Giommi, P., & Resconi, E. 2015, MNRAS, 452, 1877

Page, M. J., Kuin, N. P. M., Breeveld, A. A., et al. 2013, MNRAS, 436, 1684

Paiano, S., Falomo, R., Treves, A., & Scarpa, R. 2018, ApJL, 854, L32

Petropoulou, M., Coenders, S., & Dimitrakoudis, S. 2016, APh, 80, 115

Petropoulou, M., Coenders, S., Vasilopoulos, G., Kamble, A., & Sironi, L. 2017, MNRAS, 470, 1881

Petropoulou, M., & Dermer, C. D. 2016, ApJL, 825, L11

Petropoulou, M., Dimitrakoudis, S., Padovani, P., Mastichiadis, A., & Resconi, E. 2015, MNRAS, 448, 2412

Petropoulou, M., Giannios, D., & Dimitrakoudis, S. 2014, MNRAS, 445, 570

Petropoulou, M., & Mastichiadis, A. 2015, MNRAS, 447, 36

Petropoulou, M., Nalewajko, K., Hayashida, M., & Mastichiadis, A. 2017, MNRAS, 467, L16

Protheroe, R. J., & Johnson, P. A. 1996, APh, 4, 253

Rodrigues, X., Fedynitch, A., Gao, S., Boncioli, D., & Winter, W. 2018, ApJ, 854, 54

Romanova, M. M., & Lovelace, R. V. E. 1997, ApJ, 475, 97

Roming, P. W. A., Kennedy, T. E., Mason, K. O., et al. 2005, SSRv, 120, 95

Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103

Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525

Sikora, M., Begelman, M. C., & Rees, M. J. 1994, ApJ, 421, 153

Smette, A., Sana, H., Noll, S., et al. 2015, A&A, 576, A77

Smith, M. W. E., Fox, D. B., Dowen, D. F., et al. 2013, APh, 45, 56

Sobolewska, M. A., Siemiginowska, A., Kelly, B. C., & Nalewajko, K. 2014, ApJ, 786, 143

Steele, I. A. 2017, ATel, 10799, 1

Takahashi, T., Tashiro, M., Madejski, G., et al. 1996, ApJL, 470, L89

Tanaka, Y. T., Buson, S., & Kocevski, D. 2017, ATel, 10791, 1

Tavecchio, F., & Ghisellini, G. 2015, MNRAS, 451, 1502

Tavecchio, F., Ghisellini, G., & Guetta, D. 2014, ApJL, 793, L18

Turley, C. F., Fox, D. B., Murase, K., et al. 2016, ApJ, 833, 117

Vernet, J., Dekker, H., D'Odorico, S., et al. 2011, A&A, 536, A105