Dissecting the region around IceCube-170922A: the blazar TXS 0506+056 as the first cosmic neutrino source

We present the dissection in space, time, and energy of the region around the IceCube-170922A neutrino alert. This study is motivated by: (1) the first association between a neutrino alert and a blazar in a flaring state, TXS 0506+056; (2) the evidence of a neutrino flaring activity during 2014 - 2015 from the same direction; (3) the lack of an accompanying simultaneous $\gamma$-ray enhancement from the same counterpart; (4) the contrasting flaring activity of a neighbouring bright $\gamma$-ray source, the blazar PKS 0502+049, during 2014 - 2015. Our study makes use of multi-wavelength archival data accessed through Open Universe tools and includes a new analysis of Fermi-LAT data. We find that PKS 0502+049 contaminates the $\gamma$-ray emission region at low energies but TXS 0506+056 dominates the sky above a few GeV. TXS 0506+056, which is a very strong (top percent) radio and $\gamma$-ray source, is in a high $\gamma$-ray state during the neutrino alert but in a low though hard $\gamma$-ray state in coincidence with the neutrino flare. Both states can be reconciled with the energy associated with the neutrino emission and, in particular during the low/hard state, there is evidence that TXS 0506+056 has undergone a hadronic flare with very important implications for blazar modelling. All multi-messenger diagnostics reported here support a single coherent picture in which TXS 0506+056, a very high energy $\gamma$-ray blazar, is the only counterpart of all the neutrino emissions in the region and therefore the most plausible first non-stellar neutrino and, hence, cosmic ray source.


INTRODUCTION
The IceCube Neutrino Observatory at the South Pole 1 has recently reported the detection of a number of high-energy astrophysical neutrinos 2 (Aartsen et al. 2013;IceCube Collaboration 2013, 2015a. These include 82 high-energy starting events collected over six years (IceCube Collaboration 2017b), which are inconsistent with a purely atmospheric origin with a significance greater than 6.5 σ. The Ice-to greater than 300 GeV (Atwood et al. 2009) that has been surveying the γ-ray sky for the past almost 10 years. The constant monitoring and archiving of all-sky γ-ray data permits unprecedented investigations of variable sources.
To explore the complexity of the multi-wavelength sky, we make use of innovative tools that are under development within "Open Universe", a new initiative under the auspices of the United Nations Committee On the Peaceful Uses of Outer Space (COPUOS). The goal of Open Universe is to stimulate a large increase of the accessibility and usability of space science data in all sectors of society from the professional scientific community, to universities, schools, museums, and citizens. A web portal of the Open Universe initiative, developed at the Italian Space Agency, is available at openuniverse.asi.it. In this paper we make use of software, such as the VOU-BLAZAR tool, which has been specifically designed to identify blazars based on multifrequency information in large error regions, and spectral energy distribution (SED) animations.
Blazars are active galactic nuclei (AGN; see Padovani et al. 2017, for a recent AGN review) having a relativistic jet that is seen at a small angle with respect to the line of sight. The jet contains charged particles moving in a magnetic field emitting non-thermal radiation over the entire electromagnetic spectrum (Urry & Padovani 1995;Padovani et al. 2017). Since the energy distribution of these particles can significantly differ from object to object, the electromagnetic emission exhibits a wide range of intensity levels and spectral slopes across the spectrum. This results in observational properties that depend strongly on the energy band where blazars are discovered. In a series of papers Giommi et al. (2012aGiommi et al. ( , 2013; ;  proposed a new blazar paradigm (but see Ghisellini et al. 2017, for an alternative scenario) based on dilution by the jet and the host galaxy, minimal assumptions on the physical properties of the non-thermal jet emission, and unified schemes. By means of detailed Monte Carlo simulations, it was shown that this scenario is consistent with the complex observational properties of blazars as we know them in all parts of the electromagnetic spectrum.
The possibility that blazars could be the sources of highenergy neutrinos has been investigated by many authors, even long before the IceCube detections (e.g. Mannheim 1995;Halzen & Zas 1997;Mücke et al. 2003;Padovani & Resconi 2014;Petropoulou et al. 2015;Tavecchio & Ghisellini 2015). Padovani et al. (2016) have correlated the second catalogue of hard Fermi-LAT sources (2FHL, E > 50 GeV, Ackermann et al. 2016) and other catalogues, with the publicly available high-energy neutrino sample detected by Ice-Cube. The chance probability of association of 2FHL highenergy peaked blazars (HBL/HSP, i.e. sources with the peak of the synchrotron emission ν S peak > 10 15 Hz 3 ; Padovani & Giommi 1995) with IceCube events was 0.4 per cent, which becomes 1.4 per cent (2.2 σ) by evaluating the impact of trials (Resconi et al. 2017). This hint appears to be strongly dependent on γ-ray flux. The corresponding fraction of the 3 Blazars can be further divided into low (LBL: ν S peak < 10 14 Hz) and intermediate (IBL: 10 14 Hz< ν S peak < 10 15 Hz) energy peaked sources.
IceCube signal explained by HBL is however only ∼ 10 − 20 per cent, which agrees with the results of Aartsen et al. (2017); IceCube Collaboration (2017c), who by searching for cumulative neutrino emission from blazars in the second Fermi-LAT AGN (2LAC; Ackermann et al. 2011) and other catalogues (including also the 2FHL), have constrained the maximum contribution of known blazars to the observed astrophysical neutrino flux to < 27 per cent.
High-energy astrophysical neutrinos originate in cosmic ray interactions providing a natural link with high-energy and possibly ultrahigh-energy cosmic ray (UHECR) detection. Resconi et al. (2017) have presented a hint of a connection between HBL, IceCube neutrinos, and UHECRs (E ≥ 52 × 10 18 eV) with a probability ∼ 0.18 per cent (2.9 σ) after compensation for all the considered trials. Even in this case, HBL can account only for ≈ 10 per cent of the UHECR signal.
None of the possible neutrino counterparts in Padovani et al. (2016) and Resconi et al. (2017) are tracks, as they are all cascade-like events 4 . This indicates that by using tracks we are still limited in sensitivity to the HBL neutrino signal.
Although tracks trace only about 1/6 of the astrophysical signal for a flavour ratio ν e : ν µ : ν τ = 1 : 1 : 1, standard neutrino cross-sections, and IceCube event selection efficiencies, after a long enough exposure a track IceCube signal from blazars should also start to appear. This is of great interest, because false (random) associations of tracks with a blazar are unlikely due to the better defined position of this eventclass with respect to cascades.
Recently Lucarelli et al. (2017a) have found a transient γ-ray (> 100 MeV) AGILE source positionally coincident with an IceCube track with a post-trial significance ∼ 4 σ and possibly associated with an HBL. However, no other space missions nor ground observatories have reported any detection of transient emission consistent with this event.
The most probable hint of an association (3 − 3.5 σ) reported so far (IceCube Collaboration 2018a) between an IceCube astrophysical neutrino and an extragalactic object is that of the neutrino IceCube-170922A and the radio bright (∼ 1 Jy at 5 GHz) and γ-ray flaring BL Lac object TXS 0506+056 (also known as 5BZB J0509+0541, 2FHL J0509.5+0541, and 3FGL J0509.4+0541). Moreover, IceCube has reported in IceCube Collaboration (2018b) an independently observed 3.5 σ excess of neutrinos from the direction of TXS 0506+056 between October 2014 and February 2015 providing further indication of a high-energy neutrino association.
The presence of several non-thermal objects including variable blazars within 80 arc-minutes of IceCube-170922A (that is within the size of the Fermi-LAT point spread function [PSF 5 ], which is ∼ 2.8 • [95 per cent containment] at E = 1 GeV) makes the γ-ray emission from this area quite complex, with possible source confusion. Different objects 4 The topology of IceCube detections can be broadly classified in two types: (1) cascade-like, characterized by a compact spherical energy deposition, which can only be reconstructed with a spatial resolution ≈ 15 • ; (2) track-like, defined by a dominant linear topology from the induced muon, with positions known typically within one degree or less. 5 The Fermi PSF for this analysis event selection has been determined using the gtpsf Fermi Science Support Center tool. could in fact contribute to the overall γ-ray flux at different levels in a time-dependent manner. For this reason, we report here on what we have called the dissection of the region around IceCube-170922A taking into account the fact that all sources present in the area could be in principle contributors to the neutrino emission observed in 2014-2015 and in 2017. We use innovative software tools that exploit all the publicly available multi-frequency data to study in detail the area around the position of IceCube-170922A at all energies and in the time domain, together with a very careful analysis of the γ-ray emission, providing a wide perspective in space and time.
Section 2 describes the multi-messenger data we used, while Section 3 puts them all together to study the relevant sources in the area, their γ-ray light curves and SEDs. Section 4 gives our results, which are discussed in Section 5. Section 6 summarizes our conclusions. We use a ΛCDM cosmology with Hubble constant H 0 = 70 km s −1 Mpc −1 , matter density Ω m,0 = 0.3, and dark energy density Ω Λ,0 = 0.7. IceCube Collaboration (2018a) have derived the coincidence probability as a measure of the likelihood that a neutrino alert like IceCube-170922A is correlated by chance with a flaring blazar, considering the large number of known γray sources and the modest number of neutrino alerts. This has been done through several hypothesis tests covering a range of assumptions on the spatial and temporal signal distribution and neutrino emission scenarios. The derived probability was trial-corrected by multiplying the p-value by 51, which correspond to the number of alerts issued by IceCube (10) plus the 41 inspected archival events, which would have triggered alerts if the realtime system had been operational. The final post-trial coincidence probability ranges between 2.5 × 10 −4 and 1.3 × 10 −3 (3 − 3.5 σ).

The IceCube neutrino flare in 2014-2015
In contrast with the neutrino alert discussed above, which is a single event identified in real time and satisfying stringent selection criteria, we define here a neutrino flare as a statistically significant (> 3 σ) accumulation of neutrinos coming from a specific direction over a well-defined time period.

Radio and optical monitoring data
The radio (15 GHz) and optical (V mag ) data have been taken from the Owens Valley Radio Observatory (OVRO 7 ) database, the Catalina Real time Transient Survey (CRTS 8 ) and from the All Sky Automatic Survey (ASAS 9 ; Kochanek et al. 2017) online services.

X-ray and optical/UV data
The Neil Gehrels Swift observatory (Geherls et al. 2004) carried out a total of 92 observations within 80 arc-minutes of IceCube-170922A. Of these, 35 were performed before the arrival of the neutrino and were mostly pointed at the nearby flat spectrum radio quasar (FSRQ) PKS 0502+049. The remaining pointings have been carried out either as a Target of Opportunity follow up mapping the error region of IceCube-170922A a few hours after the event (Keivani et al. 2017), or as part of the monitoring program of the blazar TXS 0506+056 triggered by the IceCube neutrino alert.
We have analyzed all 92 X-Ray Telescope (XRT; Burrows et al. 2005) observations using the latest version of the Swift data reduction software (HEADAS 6.22) applying standard procedures. This led to the detection of 251 X-ray sources, which were combined with those of existing X-ray catalogues to build Fig.1 (left). X-ray spectral data were used together with the available multi-frequency data to assemble the SEDs of all interesting sources in the field (see Sect. 3.1).
All optical and UV data of the Swift Ultra-Violet and Optical telescope (UVOT; Roming et al. 2005) for TXS 0506+056 and PKS 0502+049 were analyzed using the SSDC online interactive archive 10 .
The NuSTAR hard X-ray observatory (Harrison et al. 2013) was pointed twice, on September 29 and October 19 2017, at TXS 0506+056 following the detection of IceCube-170922A. A few days after the observations the data were made openly available. We have analyzed these data sets using the online analysis tool of the SSDC archives following the standard procedure. In both observations the spectral shape shows a sharp hardening at about 4 − 5 keV. Figure 1. Left: Radio and X-ray sources within 80 arc-minutes of the position of IceCube-170922A. Symbol diameters are proportional to source intensity. Radio sources appear as red filled circles, X-ray sources as open blue circles. Right: Known and candidate blazars around IceCube-170922A as detected by a tool developed within the framework of the Open Universe initiative as described in the text. Dark blue circles represent LBL type candidates, that is sources with flux ratio in the range observed in the sample of LBL blazars of the latest edition of the BZCAT catalogue (Massaro et al. 2015), cyan symbols are for IBL type candidates, and orange symbols are for HBL candidates. Known blazars are also marked by a red diamond if they are included in the BZCAT catalogue or a star if they are part of the 2WHSP sample (see text for sources n. 2, 4, 6, and 7). The diameters of filled and open circles are proportional to radio flux density and X-ray flux, respectively. The dashed line shows the 90 per cent error contour of IceCube-170922A. The localization of the Fermi sources is such that the error ellipses are smaller than the size of the symbols.

γ-ray data
We used publicly available Fermi-LAT Pass 8 (with the P8R2 SOURCE V6 instrument response functions) data acquired in the period from August 4, 2008 to February 10, 2018 and followed the standard procedures suggested by the Fermi-LAT team. Only the events with a high probability of being photons (evclass = 128, evtype = 3 [FRONT+BACK]) in the energy range of 100 MeV -300 GeV from a region of interest (ROI) defined as a circle of radius 12 • centred at the γ-ray position of TXS 0506+056 (RA, Dec = 77.364, 5.699) were analyzed. We removed a possible contamination from the Earth limb by cutting out all the events with zenith angle > 90 • and only used the time intervals in which the data acquisition of the spacecraft was stable (DATA QUAL>0 && LAT CONFIG==1). Consistently with the event selection we used the standard Galactic (gll_iem_v06 11 ) and isotropic (iso_P8R2_SOURCE_V6_v06 11 ) models to describe the diffuse background emissions.

Test Statistic maps
To evaluate the presence of relevant γ-ray signatures around the arrival direction of IceCube-170922A, we built test statistics (TS) maps of the region (Mattox et al. 1996). The 11 https://fermi.gsfc.nasa.gov/ssc/data/access/lat/ BackgroundModels.html test statistic for all the γ-ray analysis in this paper is defined as where L (source) represents the nested likelihood of the data given a specific source hypothesis and L (nosource) the likelihood of the background model. In our TS maps the signal hypothesis of a γ-ray point source is tested against a background model consisting of a diffuse Galactic and a diffuse isotropic component, as well as all the the Fermi-LAT third source catalogue (3FGL; Acero et al. 2015) sources that lie outside the region of the TS map. While the point source is modeled using a power-law with free normalization and spectral index, the parameters of the background sources remain fixed. According to Wilks' theorem the test-statistic distributions follows a X 2 distribution with two degrees of freedom. Hence TS values of 30 and 8 are equivalent to a 5 and 2 σ significance, respectively. We centred our 80×80 arc-minute maps at IceCube-170922A and used an equally spaced grid with 0.05 • step size in right ascension and declination. For each of the grid points the Fermi Science Tools unbinned likelihood analysis tool gtlike was used to maximize the likelihoods in eq. (1) with respect to the free parameters. We built TS maps for different time windows and energy cuts, resolving the γ-ray activity during the periods of the neutrino detections and for energy thresholds 1 GeV, 2 GeV, and 5 GeV. Our choices are driven by the the need for sufficiently high space resolution to distinguish different γ-ray sources (as the PSF decreases with energy) and by the fact that, given the possible neutrino production scenarios, we are mostly interested in the photons at the highest energies.

Light curve and photon index
γ-ray flux and photon index variations were investigated using light curves generated with a fixed time binning of 55 (28 12 ) days for TXS 0506+056 (PKS 0502+049) and with an adaptive binning method with a constant relative flux uncertainty (Lott et al. 2012). In the fixed time binned light curve, photons in the 2 -100 GeV and 0.1 -100 GeV ranges were used for TXS 0506+056 and PKS 0502+049, respectively (see below). In the adaptive binning light curve analysis for PKS 0502+049 the fluxes are computed above an optimum energy of E min = 214 MeV in order to reach the required constant relative flux uncertainty of 15 per cent. For all the light curves the flux normalization and photon index of the target source are determined by applying the gtlike tool in each time bin. The model file, describing the ROI, contains point sources from the 3FGL catalogue within ROI+5 • from the target, as well as the Galactic and isotropic γ-ray background models. It is generated using the user contributed 12 This difference is due to the fact that the photon statistics for PKS 0502+049 is better than for TXS 0506+056 for the chosen energy thresholds. The bin sizes used have been chosen in order to have enough photon statistics in the majority of the bins but also to avoid any fine tuning.
make3FGLxml.py 13 tool. For the case of TXS 0506+056 the chosen energy threshold of 2 GeV efficiently removes any source confusion (see Fig. 4 and section 3.3), hence only this source is left free for the fit. The resulting light curve is robust against mis-modelling and strong time-variability of the other sources in the ROI. The light curve of PKS 0502+049, on the other hand, is derived by reaching lower energies, since the majority of the photons are below 1 GeV (Acero et al. 2015). In the source model both sources are fitted at the same time. The diffuse background components are fixed to their nine year value, since they are not expected to vary on the time scales of this analysis. Additionally, since we are using short integration times, we model PKS 0502+049 with a power-law.

PUTTING IT ALL TOGETHER
3.1 Relevant astronomical sources in the region of IceCube-170922A We have searched for non-thermal emission in blazar-like sources in the vicinity of IceCube-170922A using a tool that is being developed in the framework of the United Nations Open Universe initiative 14 using Virtual Observatory protocols. This was used to find sources in all available lists of radio and X-ray emitters. Fig. 1 (left) shows a plot of the 637 radio and X-ray sources present in those catalogues or which were detected in dedicated Swift observations within 80 arcminutes of the position of IceCube-170922A. Of these, only 7 emit in both bands and are all blazar-like in their X-ray-toradio flux ratio, as determined from the BZCAT (Massaro et al. 2015) and 2WHSP (Chang et al. 2017) samples. These are shown in Fig. 1 (right). Three known objects, i.e. TXS 0506+056 (an IBL/HBL 15 at z = 0.3365: Paiano et al. 2018), PKS 0502+049 (also known as 5BZQ J0505+0459, an LBL/FSRQ at z = 0.954), and 2WHSP J050833.3+05310 (an HBL), i.e. sources no. 5, 1, and 3 respectively in Fig. 1 (right), and four additional blazar candidates are present in the area. The first two blazars are also bright γ-ray emitters (Sect. 3.2). Visual inspection of the SED of the other sources allowed us to confirm that source no. 4 is a good candidate HBL object, while source 7 is likely a cluster of galaxies (due to its extended X-ray emission), source 6 is a steep radio spectrum object, and source 2 is a nearby (z = 0.03677) elliptical galaxy showing low luminosity X-ray emission (L ∼ 10 41 erg s −1 at 1 keV) that could be due to a jet or even to non-nuclear sources.

γ-ray emission near IceCube-170922A
The γ-ray emission near IceCube-170922A is dominated at various times either by TXS 0506+056 or by PKS 0502+049 but there are also times when the two sources have roughly equal fluxes, as shown in Fig. 2. To investigate which of the   (October 19, 2014-February 6, 2015. In this period the brightest source was PKS 0502+049 at lower energies while TXS 0506+056 dominated at higher energies. See Fig. 2 for more information.
sources dominate the γ-ray emission during the IceCube-170922A event and the neutrino flare period we constructed TS maps for these two periods. We observe the following: (i) Fig. 3 shows the TS maps during the period contemporaneous with and before the IceCube-170922A event. From the maps it appears that TXS 0506+056 dominates the photon flux of the region at energies > 1 GeV; (ii) during the time of the neutrino flare, on the other hand, the situation is different, with PKS 0502+049 being the brightest source at E > 1 GeV and TXS 0506+056 progressively taking over above 2 and 5 GeV, as shown in Fig. 4.
We tried to unveil any evidence of γ-ray emission asso-ciated with the neighbour HBL 2WHSP J050833.3+05310, which is also within the 90 per cent error contour of IceCube-170922A, by conducting a series of eleven unbinned likelihood analyses, covering 9 years of observation with Fermi-LAT. We set energy cuts at E >1 GeV and >5 GeV, integrating over 400 days intervals 16 . We kept all 3FGL sources in the field, setting both normalizations and photon indexes as free parameters. At the 2WHSP J050833.3+05310 position 16 Each bin has 100 days superposed to the previous/next bin, meaning we cover the following time windows: MJD 54700 to 55100, 55000 to 55400, 55300 to 55700... and so on. an additional power-law source was included in the model. The strongest signature for γ-ray emission was found between MJD 55900 to 56300 for the 1 GeV energy cut, reaching TS ∼ 5. This result was confirmed by the residual map of the region.

γ-ray light curves
As discussed in Sect. 2.4.2, we have derived the γ-ray light and spectral index curves of TXS 0506+056 and PKS 0502+049. As inferred from Fig. 4, we have evidence that PKS 0502+049 dominates the γ-ray sky at low energies during the neutrino flare, possibly contaminating TXS 0506+056 below 2 GeV. We build the light curves for TXS 0506+056 integrating photons above 2 GeV to avoid any bias from PKS 0502+049 and also study the highest energies (above 10 GeV) during the specific neutrino flare period. This is the best compromise for the energy threshold. Ideally, one would like to sample as high energies as possible, to profit from the smaller Fermi-LAT PSF and source containment region and therefore reduce contamination from PKS 0502+049. However, the larger the energy, the smaller the photon statistics. In contrast, we build light curves for PKS 0502+049 above 0.1 GeV with fixed time bins and above 214 MeV with the adaptive binning method to study both long and short term structures. Fig. 5 (top) shows the light curve of TXS 0506+056. The blue band denotes the neutrino flare, while the red line indicates the IceCube-170922A event. Significant flux variations are visible, as typical of blazars. The corresponding photon flux and spectral index distributions are also shown on the right. It is interesting to note that TXS 0506+056 was at its hardest in the Fermi-LAT band during the neutrino flare, while being relatively faint (a "low/hard" state), and at its brightest during the IceCube-170922A event, while being softer (a "high/soft" state). Note that, based on the overall distributions, a spectral index as Figure 7. The SED of TXS 0506+056 assembled using all available public data. The blue measurements in the γ-ray band show the variability at 0.1, 1, and 10 GeV using Fermi-LAT data collected between January 1 and December 31, 2017 (see Fig. 10 for the 1 and 10 GeV light curves). Green and brown points at X-ray frequencies are from our analysis of Swift and NuSTAR data respectively. MAGIC data (IceCube Collaboration 2018a) are shown as purple crosses and were de-absorbed to correct for the extragalactic background light following Domínguez et al. (2011). Red points at optical and UV frequencies are from Swift-UVOT. The other non-simultaneous multi-frequency measurements are from catalogues and online archives (grey points).
hard as observed during the neutrino flare is expected with a probability of only ∼ 2 per cent, while a flux as high at that during the IceCube-170922A event has a probability of only ∼ 1 per cent to be detected. The average photon index during the entire duration of the neutrino flare for E > 2 GeV is 1.62 ± 0.20.
Since we want to concentrate on the highest energy photons we zoom in on the period around the neutrino flare to investigate in more detail the source variability looking for the most extreme emission using only events above 10 GeV. The period from MJD 56850 to MJD 56750 was then divided into half-and one-day bins and an unbinned maximum likelihood analysis was performed. Next, the nearby bins with TS > 0 were merged and new light curves were calculated. In order to improve the statistics, the length of the time periods with TS > 0 was then progressively increased by adding 1-hour intervals. As a result, we find two periods with significant emission above 10 GeV (Fig. 5, bottom)  days and we did not consider small γ-ray fluctuations over much shorter time scales). PKS 0502+049 presents no particular states in photon flux and/or spectral index during IceCube-170922A and the neutrino flare period. Fig. 7 shows the SED of TXS 0506+056. The NuSTAR data show a hardening at ∼ 10 18 Hz (∼ 4 − 5 keV) most likely due to the onset of the inverse Compton component (e.g. Padovani et al. 2017). Fig. 8 shows the SED of PKS 0502+049. Two points can be made: (1) during flares PKS 0502+049 is a brighter γ-ray source than TXS 0506+056; (2) the SED of PKS 0502+049 cuts off at E ≈ 10 GeV, while that of TXS 0506+056 reaches E 100 GeV (compatibly with the likely extragalactic background light absorption at its redshift: Paiano et al. 2018).

Spectral energy distributions
To study the time evolution of the SED of TXS 0506+056 we have built an animation of nearly simultaneous data, which include 15 GHz monitoring data from OVRO, the ASAS V mag light curve (Kochanek et al. 2017), the optical/UV data of the Swift-UVOT (Roming et al. 2005) analyzed with the SSDC online tool, and the X-ray . The hybrid photon -neutrino SED of TXS 0506+056. The red points (OVRO at 15 GHz and ASAS V mag ) are simultaneous with neutrinos, grey ones refer to historical data, while the black ones are Fermi data. The red bands for the γ-ray flux show the 1 σ error bounds on the best fit, while upper limits are given at 95 per cent C.L.. Fermi data points were de-absorbed to correct for the extragalactic background light following Domínguez et al. (2011). Left: the MJD 57908 -58018 period (June 4 -September 22, 2017). The neutrino flux has been derived by IceCube Collaboration (2018a) over the 200 TeV -7.5 PeV range (see text for more details); we give here the all-flavour flux. The vertical upper limit is drawn at the most probable neutrino energy. The average Fermi-LAT photon index for E > 2 GeV is 2.16 ± 0.10. Right: the MJD 56949 -57059 period (October 19, 2014-February 6, 2015. The neutrino flux has been derived by IceCube Collaboration (2018b) over the 32 TeV -3.6 PeV range; the error is the combined error on the spectral index and the normalization. The average Fermi-LAT photon index for E > 2 GeV is 1.62 ± 0.20. and γ-ray data described in sect 2.3 and 2.4. The animation is available here: https://youtu.be/lFBciGIT0mE. Fig. 9 shows the hybrid photon -neutrino SED of TXS 0506+056 for the period around the IceCube-170922A event (left) and the neutrino flare (right), based on the concept introduced by Padovani & Resconi (2014). The red points are the electromagnetic emission simultaneous with the neutrinos. The detection of high-energy neutrinos above ∼ 30 TeV implies the existence of protons up to at least 3 × 10 14 − 3 × 10 15 eV, which then collide with other protons (pp collisions) or photons (pγ collisions). High-energy γ-rays with energy and flux about a factor two higher than the neutrinos at the source are then expected as secondary products in both cases (Kelner, Aharonian, & Bugayov 2006;Kelner & Aharonian 2008). Indeed, in both cases the (linearly extrapolated) γ-ray and neutrino fluxes are comparable, consistently with the hypothesis that they are produced by the same physical process. This is especially true for the neutrino flare, when the neutrino flux has a relatively small uncertainty being derived from ∼ 13 events within a well-defined period of 110 +35 −24 days (IceCube Collaboration 2018b). This is different for the IceCube-170922A event, given the large uncertainty on the neutrino flux since we are dealing with a single event over an ill-defined period of time. To estimate the neutrino flux from one neutrino event IceCube Collaboration (2018a) had to assume a spectral emission shape, an emission time τ, and an energy emission range. The corresponding mean number of ν µ events N expected in IceCube is

The hybrid SED of TXS 0506+056
where A eff is the effective area of the IceCube detector and 1 3 is the flavour ratio assumed. For a source described by a single power-law distribution the flux producing one neutrino event is where N = 1 and τ is taken as 0.5 years (of the same order as the duration of the γ-ray flare). Interpreting the observation of one IceCube alert event as an upward Poissonian fluctuation, then the flux value calculated can be understood as an upper limit on the neutrino flux (see also IceCube Collaboration 2018a).

RESULTS
Following up the IceCube-170922A event observed in coincidence with a γ-ray flare of TXS 0506+056 (IceCube Collaboration 2018a), the IceCube collaboration has also detected a neutrino flare in late 2014 -early 2015 from the same direction (IceCube Collaboration 2018b). Given the complexity of the γ-ray sky in this area, both spatially and temporally, we have carefully dissected the region and found the following: (i) out of the 637 radio and X-ray sources within 80 arcminutes of the IceCube-170922A event position, only 7 are both radio and X-ray emitters and therefore likely nonthermal sources. As it turns out, the X-ray-to-radio flux ratios of these 7 sources are blazar-like; (ii) out of these 7 sources, 4 are blazars, two of which are very bright γ-ray sources, namely TXS 0506+056 and PKS 0502+049, competing for dominance; (iii) while TXS 0506+056 dominates in all γ-ray bands during the IceCube-170922A event, the situation is more complex during the neutrino flare, as PKS 0502+049 dominates up to E ∼ 1 − 2 GeV but TXS 0506+056 takes over at E 2 − 5 GeV. The γ-ray spectrum of PKS 0502+049, in fact, cuts off at high energy even during flares, a behaviour typical of LBL blazars; (iv) PKS 0502+049 is flaring right before and right after the neutrino flare (but not in coincidence with it) while TXS 0506+056 was at its hardest in that time period but in a relatively faint state, suggesting a shift to high energies of the γ-ray SED; (v) the hybrid γ-ray -neutrino SED of TXS 0506+056 during the neutrino flare is as expected for lepto-hadronic models since the photon and neutrino fluxes are at the same level (Petropoulou et al. 2015). We note that the hybrid SEDs of Padovani & Resconi (2014) and Padovani et al. (2016) were based on one shower-like IceCube event, which could in principle have been emitted over the full IceCube detection live time, and were therefore affected by a very large uncertainty. In the case of the neutrino flare, instead, a sizable (∼ 13) number of neutrinos has been detected within a well-defined time window and good spatial resolution.
In short, all spatial, timing, and energetic multimessenger diagnostics point to TXS 0506+056 as the first identified non-stellar neutrino (and therefore cosmic ray) source.

Source properties
We now explore in more detail the properties of TXS 0506+056. First, we note that this source is a very strong γ-ray source, having an average flux of 7.1 × 10 −8 ph cm −2 s −1 above 100 MeV, which puts it among the top 4 per cent of the Fermi 3LAC catalogue . Moreover, it also belongs to the 2FHL sample (Ackermann et al. 2016), which includes all sources detected above 50 GeV by Fermi-LAT in 80 months of data. TXS 0506+056 also has a large radio flux density ∼ 1 Jy at 6 cm (Gregory & Condon 1991), and ∼ 537 mJy at 20 cm, which makes it one of the brightest radio sources (in the top 0.3 per cent) of the NRAO VLA Sky Survey, which covers 82 per cent of the sky (Condon et al. 1998). Fig. 7 shows the overall SED of the source in luminosity, based on the redshift of 0.3365 recently reported by Paiano et al. (2018).
The peak luminosities of ∼ 2 × 10 46 erg s −1 in the synchrotron peak, and almost 10 47 erg s −1 at 10 GeV, place this object among the most powerful BL Lacs known, particularly in the high-energy/very high-energy γ-ray band. For comparison, the corresponding maximum luminosities ever observed in MKN 421 (and PKS 2155−304) are ∼ 4 × 10 45 (∼ 2 × 10 46 ) and ∼ 1.5 × 10 45 (10 46 ) erg s −1 , a factor of ∼ 5 (1) and ∼ 50 (10) lower than TXS 0506+056 (Giommi et al., in preparation). What seems to be peculiar in this source is the very large luminosity at ∼ 10 GeV compared to other similar sources. From the overall SED point of view TXS 0506+056 shows a variability range in the γ-ray band (almost a factor 1,000 at 10 GeV: see Fig. 7) much larger than that observed at the peak of the synchrotron emission. Even during the large γ-ray flaring event observed close to the detection of IceCube-170922A the peak of the synchrotron emission (located in the UV band) did not vary by more than a factor of 2, nor did the X-ray flux, at the tail of the synchrotron peak, change by a large factor. This behaviour is consistent with an excess of hard γ-ray radiation possibly associated with hadronic processes.

Physical implications
Our results have several physical implications. Namely: (i) with neutrinos we are exploring an energy range which is inaccessible with photons at this redshift. Therefore, the first evidence for neutrino emission from the direction of TXS 0506+056 opens a new window on blazar physics; (ii) the derived L ν during the neutrino flare is quite large and even larger than L γ (2 GeV -1 TeV), which might imply that a sizeable fraction of the neutrino-related γ-rays have energies above the Fermi-LAT energy band, as expected in the case of pγ collisions (Sect. 3.4.1); (iii) the ratio between the two SED humps (high-energy to low-energy, usually known as Compton dominance [CD]: e.g. Giommi et al. 2012b) for TXS 0506+056 is > 1, while typical HBL have CD ∼ 0.1. This could be due to the electromagnetic emission coming from different blobs e.g. one dominant at optical frequencies and the other dominant at highenergy γ-rays, or to the presence of an additional hadronic component in the γ-ray band; (iv) a discrete cross-correlation analysis using the ztransformed discrete correlation function (ZDCF; Alexander 1997) among the radio (15 GHz), optical, X-ray and γ-ray light curves of TXS 0506+056, was built with the data described in the previous paragraphs (Fig. 10). The ZDCF shows a strong correlation (∼ 10 σ) with a time lag of ≈ 4 months between the radio and optical emission (with the Figure 10. The radio (15 GHz), optical (V mag ), X-ray and γ-ray light curves of TXS 0506+056. The radio data have been taken from the OVRO database, the visual magnitude data are from the Catalina Real Time Transient Survey (CRTS) and from the All Sky Automatic Survey ASAS (Kochanek et al. 2017). The γ-ray light curves have been produced using Fermi-LAT data with the adaptive-bin method (Lott et al. 2012). The blue band denotes the neutrino flare, while the red line indicates the IceCube-170922A event.
optical band leading) but no correlation between optical/Xray/radio and γ-ray bands. This is at variance with typical single-zone leptonic (synchrotron self-Compton or external Compton) models where all energy bands are expected to be well correlated; (v) the SED of at least one blazar now has to be modeled within a lepto-hadronic scenario (e.g. Petropoulou et al. 2015). This is far from trivial and goes well beyond the scope of this paper; (vi) the neutrino sky has been populated so far by only two sources: the Sun (e.g. BOREXINO Collaboration et al. 2014) and SN 1987A (Hirata et al. 1987;Bionta et al. 1987;Alekseev et al. 1987). TXS 0506+056 is now a third plausible candidate, whose neutrino energies are, however, more than six orders of magnitude larger than those of the two stellar sources.

Future searches
Our results suggest two periods of neutrino emission for TXS 0506+056, which appear to be in connection with two very different γ-ray states, namely one high/soft, connected to the IceCube-170922A event, and another one low/hard, related to the neutrino flare. This implies that the search for multi-messenger sources needs to be carried out not only for flaring γ-ray sources but also for relatively hard emitters. These criteria can also be used by neutrino telescopes to look for other neutrino sources.
Why is this the only strong IceCube neutrino source candidate? What makes it special? We believe this is due to a series of factors. First, as discussed in Sect. 5.1, TXS 0506+056 is a very strong γ-ray and radio source. Moreover, its declination and energy range happened to be in the regions of parameter space where IceCube reaches maximum sensitivity and the duration of the flare was also long enough for IceCube to (marginally) detect it (see also IceCube Collaboration 2018b).
Future searches for cosmic neutrino sources should emphasize: (1) extreme blazars of the BL Lac type, as hinted at by Padovani et al. (2016); (2) high/soft -low/hard γ-ray states; (3) regions of parameter space where the neutrino detectors are most sensitive.

CONCLUSIONS
The IceCube-170922A event and the neutrino flare at the end of 2014 have been linked to the same source, TXS 0506+056, a blazar of the BL Lac type at z = 0.3365. This is the most plausible association so far between Ice-Cube neutrinos and an extragalactic object. The area near TXS 0506+056 is quite complex due to the presence of several non-thermal sources, which in principle could all contribute to the overall γ-ray flux. We have therefore carefully dissected this region using a multi-messenger approach, obtaining the following results: (i) TXS 0506+056 was the brightest Fermi source in the region of interest at energies above 1 GeV during the IceCube-170922A event but only above 2 − 5 GeV during the neutrino flare. PKS 0502+049, a nearby blazar of the FSRQ type offset by ∼ 1.2 • , dominated the γ-ray region at lower energies in the latter period.
(ii) We have observed two periods of significant neutrino emission consistent with the position of TXS 0506+056 in connection with two very different γ-ray states, namely one, high/soft, connected to the IceCube-170922A event, and another one, low/hard, related to the neutrino flare. PKS 0502+049 was also flaring right before and right after the neutrino flare but not in coincidence with it.
(iv) Both the lack of a correlation between the γ-ray and radio/optical flux and the SED shape of TXS 0506+056, which is unusual in terms of its Compton dominance, appear not to be consistent with simple leptonic models.
(v) All of the above is fully consistent with the hypothesis that TXS 0506+056 has undergone a hadronic flare during the neutrino detections.
In short, all spatial, timing, and energetic multimessenger diagnostics point to TXS 0506+056 as the only counterpart of all the neutrinos observed in the vicinity of IceCube-170922A and therefore the first non-stellar neutrino (and hence cosmic ray) source. The emergent picture is that extreme blazars, i.e., strong, very high energy γ-ray sources with the peak of the synchrotron emission > 10 14 − 10 15 Hz (Padovani et al. 2016), are the first class of sources with evident contribution to the IceCube diffuse signal (IceCube Collaboration 2017a,b).
Future searches for cosmic neutrino sources, concentrated on similar classes of sources, using additional high-energy track-like events, and based on detailed multimessenger analysis, will likely provide further associations.
Urry C. M., Padovani P., 1995, PASP, 107, 803 This paper has been typeset from a T E X/L A T E X file prepared by the author.