Coincidence of a high-fluence blazar outburst with a PeV-energy neutrino event

The discovery of extraterrestrial very-high-energy neutrinos by the IceCube collaboration has launched a quest for the identification of their astrophysical sources. Gamma-ray blazars have been predicted to yield a cumulative neutrino signal exceeding the atmospheric background above energies of 100 TeV, assuming that both the neutrinos and the gamma-ray photons are produced by accelerated protons in relativistic jets. Since the background spectrum falls steeply with increasing energy, the individual events with the clearest signature of being of an extraterrestrial origin are those at PeV energies. Inside the large positional-uncertainty fields of the first two PeV neutrinos detected by IceCube, the integrated emission of the blazar population has a sufficiently high electromagnetic flux to explain the detected IceCube events, but fluences of individual objects are too low to make an unambiguous source association. Here, we report that a major outburst of the blazar PKS B1424-418 occurred in temporal and positional coincidence with the third PeV-energy neutrino event (IC35) detected by IceCube. Based on an analysis of the full sample of gamma-ray blazars in the IC35 field and assuming a photo-hadronic emission model, we show that the long-term average gamma-ray emission of blazars as a class is in agreement with both the measured all-sky flux of PeV neutrinos and the spectral slope of the IceCube signal. The outburst of PKS B1424-418 has provided an energy output high enough to explain the observed PeV event, indicative of a direct physical association.

are produced by accelerated protons in relativistic jets. Since the background spectrum falls steeply with increasing energy, the individual events with the clearest signature of being of an extraterrestrial origin are those at PeV energies. Inside the large positional-uncertainty fields of the first two PeV neutrinos detected by IceCube, the integrated emission of the blazar population has a sufficiently high electromagnetic flux to explain the detected IceCube events, but fluences of individual objects are too low to make an unambiguous source association.
Due to the Earth's opacity, the IceCube HESE analysis detects events at PeV energies mainly * from the southern sky 3 . Thus, contemporaneous astronomical data to probe the various source hypotheses can best be obtained via southern-hemisphere monitoring programs. TANAMI † is a multiwavelength program 15,16 that monitors the brightest γ-ray-loud active galactic nuclei (AGN) located in the southern sky at declinations below −30 • . It comprises the ideal database to estimate the diffuse neutrino flux due to the integrated emission of AGN in a given large field at a given time, as well as the maximum possible neutrino flux associated with an individual object of the sample.
Blazars are radio-loud AGN with jets oriented close to the line of sight. This substantially increases the apparent brightness of these objects owing to the Doppler boosting of the emission from the relativistically moving emission zones. A direct association of a PeV-neutrino with an * A different analysis of IceCube muon neutrinos finds an excess signal also from the northern sky 14 . † http://pulsar.sternwarte.uni-erlangen.de/tanami individual γ-ray blazar would have the important implication that a sizeable fraction of their observed γ-ray emission must be due to hadronic decays, and that blazar jets are also sources of ultra-high-energy cosmic rays 17 . The X-ray and γ-ray emission of blazars may originate from the photoproduction of pions by accelerated protons 18 . Protons that are accelerated in the jet (e.g., via shock acceleration) could interact with 'seed' photons (e.g., ultraviolet photons from the accretion disk surrounding the central supermassive black hole). The resulting cascades produce charged and neutral pions, which decay and produce neutrinos and high-energy photons. Simple estimates and detailed Monte Carlo simulations show 5,19 that in this scenario F γ F ν , where the X-ray to γ-ray flux is integrated over the high-energy spectral energy distribution (SED). If the seed photons are provided by a blue/UV bump component, as is typical in the blazar subclass of flat spectrum radio quasars (FSRQs), the neutrino spectrum is expected to peak at PeV energies 5 . Attributing the high-energy electromagnetic emission to these photohadronic processes, the maximum possible neutrino PeV emission can be estimated from the measured integrated flux of high-energy photons.
Using TANAMI multiwavelength data, we previously compiled and discussed the multiwavelength properties of the six radio-and γ-ray-brightest blazars located inside the R 50 fields of the two ∼ 1 PeV events IC 14 and IC 20 from the first two years of IceCube data 5 . We found relatively low maximum neutrino fluxes of these six individual blazars owing to their low fluence built-up over two years, but the diffuse flux due to the integrated emission of all blazars in the fields was found to be sufficiently high to expect up to two events. When the contribution of the large number of fainter sources from the blazar population is taken into account 20 , the maximum possible neutrino flux inside a given field is increased further. A high-angular-resolution point-  Along with the very bright γ-ray emission, an increase in X-ray, optical, and radio emission from PKS B1424−418 has also been reported 28-30 . Figure where F γ is the γ-ray energy flux of all blazars located inside Ω integrated between 5 keV and which includes maximum-possible PeV neutrino counts from all γ-ray blazars from the 2LAC catalog plus a maximum-possible contribution of the large population of faint unresolved blazars.

Predicted all-sky count of PeV neutrinos from blazars and comparison to observation
By extrapolating from the fairly representative field Ω IC 35 , we estimate the maximum number of PeV neutrino events from all blazars (both resolved and unresolved) over three years from the full southern sky to be This number of events would be expected I) if only electron neutrinos would be produced, II) if all blazars harbored dense UV photon fields due to the emission of optically thick accretion disks as is typical for FSRQs, and III) if the neutrino spectrum peaked sharply at PeV energies. All three conditions are clearly not fulfilled as only 3 events have been detected, leading to an empirical scaling factor of This can be compared to a theoretical value f th , which accounts for physically motivated realistic deviations from the three ideal conditions. The theoretical scaling factor allows us to predict the The scaling factor is factorized into a flavor factor f I , a factor accounting for the different classes of blazars f II , and a spectrum factor f III : The IceCube data indicate an equal flavor ratio 4 so that the flavor factor would be 1/3 if only electron neutrinos are accounted for when computing the maximum event numbers. When adding the two other flavors, it has to be considered that the number of detected cascade events due to muon and tau neutrinos is lower than for electron neutrinos because of the energy-dependent cross sections and inelasticities for neutral-current (NC) and charged-current interactions. Assuming an underlying neutrino power law with slope −2.  22 . For our basic model of a sharply peaked neutrino spectrum due to photopion production from monoenergetic UV photons, f III would be equal to unity. In a more realistic scenario, a range of Doppler shifts (depending on the location of the seed-photon sources with respect to the relativistic jet base as discussed in our earlier work 5 ) causes broader spectra extending to lower neutrino energies. Considering also broadening due to the different redshifts of sources, an output range of ∼ 30 TeV to ∼ 10 PeV can be expected. In addition, models which consider proton-proton collisions or assume accretion tori with virial temperatures of ∼ 10 9 K rather than optically thick accretion disks 7, 33, 34 also predict softer spectra. Using a spectral index of −2.3 as measured by IceCube 3 and the (30 TeV to 10 PeV) bandwidth of the spectrum reduces the number of PeV output neutrinos by f III = 0.05, so that we estimate (cf. Eq. 4). Our model thus predicts 0.0125 · 336 ∼ 4 events at PeV energies from the full southern sky, which is remarkably close to the observed three PeV events. We conclude that the measured γ-ray emission of the blazars in the IC 35 field allows us to reproduce both the measured all-sky flux of PeV neutrinos and the measured spectral slope of the IceCube signal assuming a simple photo-hadronic emission model of FSRQs.

Predicted number of PeV neutrinos for individual FSRQs
If Ω becomes small, containing only one individual FSRQ, we can set f II = 1. The predicted number of PeV neutrinos for an individual FSRQ is then Apr 30). During these 9 months, the source increased its predicted neutrino-production rate by more than an order of magnitude so that the Poisson probability to detect a neutrino associated with the 9-month high-fluence outburst of PKS B1424−418 is at a considerable level of about 11 %, which is three times higher than the corresponding probability to detect an event from the integrated emission of all other known γ-ray blazars in the field during this 9-month period. Our model thus allows us to associate an individual blazar during a rare major outburst with the highestenergy extraterrestrial neutrino detected by IceCube to date.
Why don't we detect PeV neutrinos from every bright blazar?
If our model is correct, it also has to explain the non-detection of PeV neutrinos in positional agreement with other high-fluence blazars and with the detection statistics of sub-PeV neutrino events.
We note that the positional uncertainties R 50 given by the IceCube team are median values, which means that only half of all events originate inside their measured R 50 regions while the other half are coming from larger offset angles. Above, we have calculated the maximum number of neutrino events that can be explained by individual astrophysical sources within R 50 for a high-confidence event. When asking for the maximum number of IceCube events that might be associated with a given astrophysical source, a larger radius has to be considered. We have used the Fermi/LAT monitored source list light curves § to identify candidate sources for high keV-to-GeV fluence, compiled the average SEDs over the three years of IceCube integration for the top-ten candidate sources from the whole sky and derived their expected neutrino counts (see Table 2). For northern-and southern-hemisphere events, we have used the effective areas for the appropriate minimum energy provided by the IceCube team 2 . We do not extend the list beyond rank 10, because for the tenth-ranked source, the maximum possible neutrino output has already dropped by more than an order of magnitude This is about two orders of magnitude more constraining than the neutrino-velocity limit derived from SN 1987A 37 . However, as discussed in the Supplementary 'Methods', a ∼ 5 % probability for a chance coincidence remains. It also cannot be excluded that the observed PeV neutrino could be associated with an historical (or future) outburst of the source.

Summary and outlook
Tentative associations of high-energy neutrinos with flaring blazars have been suggested before 39,41 but it remained questionable whether a high-enough neutrino flux could be produced in the candidate flares 40 . Here, we have identified for the first time a single source that has emitted a sufficiently high fluence during a major outburst to explain an observed coinciding PeV neutrino event. There Competing Interests The authors declare that they have no competing financial interests.
Correspondence Correspondence and requests for materials should be addressed to M.K. (email: matthias. kadler@astro.uni-wuerzburg.de).

VLBI Imaging
The TANAMI VLBI images of PKS B1424−418 (Supplementary Data Fig. 1 Table 3 gives an overview of the array configuration and image parameters. The noise level for

Fermi/LAT γ-ray data analysis
For the analysis of Fermi/LAT γ-ray data, we used the Fermi Science Tools (v9r32p5) with the reprocessed Pass 7 data and the P7REP SOURCE V15 instrument response functions 35  has been calculated using 14 day bins. We applied a Bayesian-blocks analysis 36 to the light curve.
We used the change points to determine the onset of the outburst as 2012 Jul 16, which is marked as a shaded blue region in Fig. 1.
All γ-ray sources detected by Fermi /LAT with flares that exceeded the pre-defined threshold of 10 −6 cm −2 s −1 (> 100 MeV) are included in the LAT Monitored Source List * * . By definition, this sample represents the brightest Fermi /LAT sources. We used this list to identify the sources with the highest γ-ray fluences and analyzed their γ-ray light curves in detail as described above. We

Broadband SED
The broadband SED of PKS B1424−418 (Fig. 2) has been built from quasi-simultaneous data from Swift/XRT data were reduced with standard methods, using the most recent software packages (HEASOFT 6.15.1) and calibration databases. Spectra were grouped to a minimum signalto-noise ratio of 5. Spectral fitting was performed with ISIS 1.6.2 32 . We fitted the (0.5-10) keV energy band with an absorbed power law model, which yielded good results for all time ranges.
The source showed no evidence of intrinsic X-ray absorption in excess of the Galactic value 40 . Xray data were deabsorbed using the best available abundances 41 and cross-sections 42 . Swift/UVOT data were extracted following standard methods. Optical, infrared, and ultraviolet data were dereddened using the same absorbing columns 43  in peak flux in order to still describe the slope at X-ray energies and reduces the overall fit quality.
The increase in peak flux leads to greater numbers of neutrinos (see below), with the exception of the short flare, so that our method to ignore the upper limits in the fit is more conservative for our purposes.

Neutrino event prediction
To derive maximum neutrino events to be expected from the measured SEDs, we applied the techniques that we developed for our analysis of the AGN in the fields of the first two PeV IceCube events. We used 2FGL GeV γ-ray spectra 47 and extracted X-ray spectra taken during the same period by Swift 48 . For sources not observed by Swift, we used X-ray data or upper limits from the ROSAT all-sky survey 49 . For PKS B1424−418, we calculated in addition Fermi/LAT γ-ray spectra for the various time ranges of interest according to the methods described above. We fit a single log parabola to the high-energy data. From the integrated 1 keV to 5 GeV SEDs, we can derive a maximum-possible neutrino flux for each source. From this we can derive the maximum-possible number of counts for three years of IceCube integration (see Supplementary Data Table 1) using the appropriate effective area for northern/southern-hemisphere sources for electron neutrinos of the IceCube 998-day analysis 3 . Formal uncertainties of neutrino counts can be derived from the uncertainties of the integrated electromagnetic energy output but are dominated by systematic uncertainties due to the simplicity of the model.

Calculation of chance coincidences of major blazar outbursts and PeV neutrino events
We estimate the number of chance coincidences of major blazar outbursts and PeV neutrino events (2 periods), and PKS B1510−089 (3 periods). Hence, the rate of such outbursts occurring within Ω ν is expected to beṄ outburst ∼ 0.1 yr −1 sr −1 × Ω ν = 0.016 yr −1 . Using these numbers we end up with N coinc ∼ 0.05 as the mean number of chance coincidences. The Poisson probability to observe one or more such coincidences is about 5%. However, because of the lack of a pre-defined statistical test, we cannot formally use this value to test the hypothesis of a chance coincidence.

Coincidences of sub-PeV neutrino events with high-fluence blazars
All high-fluence blazars from     Table 3 for the individual observations. All contours start at 3σ RMS and increase logarithmically by factors of 2 .