A publishing partnership

The following article is Free article

Patterns in the Multiwavelength Behavior of Candidate Neutrino Blazars

, , , , , , , , ,

Published 2020 April 28 © 2020. The American Astronomical Society. All rights reserved.
, , Citation A. Franckowiak et al 2020 ApJ 893 162DOI 10.3847/1538-4357/ab8307

Download Article PDF
DownloadArticle ePub

You need an eReader or compatible software to experience the benefits of the ePub3 file format.

0004-637X/893/2/162

Abstract

Motivated by the identification of the blazar TXS 0506+056 as the first promising high-energy neutrino counterpart candidate, we search for additional neutrino blazar candidates among the Fermi–Large Area Telescope detected blazars. We investigate the multiwavelength behavior from radio to GeV gamma-rays of blazars found to be in spatial coincidence with single high-energy neutrinos and lower-energy neutrino flare candidates. In addition, we compare the average gamma-ray emission of the potential neutrino-emitting sources to the entire sample of gamma-ray blazars. We find that neutrino-emitting blazar candidates are statistically compatible with hypotheses of both a linear correlation and no correlation between neutrino and gamma-ray energy flux.

Export citation and abstractBibTeXRIS

1. Introduction

After the detection of a diffuse flux of high-energy neutrinos (Aartsen et al. 2013), the most pressing challenge is to identify where these neutrinos are produced. Among the prime candidates are active galactic nuclei (AGNs), especially those with a relativistic jet pointing toward us, the so-called blazars (e.g., Stecker et al. 1991; Mannheim et al. 1992; Mannheim 1993, 1995; Szabo & Protheroe 1994; Mastichiadis 1996; Protheroe 1999; Atoyan & Dermer 2001; Dimitrakoudis et al. 2012; Murase 2017). No significant clusters of neutrinos in either space or time have been identified by all-sky searches of IceCube data (Aartsen et al. 2015, 2017a, 2020). Searching for neutrinos from a predefined list of 110 sources revealed a excess at the position of the Seyfert 2 galaxy NGC 1068 (Aartsen et al. 2020). Combining neutrino with multiwavelength data is the key to probing neutrino emission from various source populations and identifying potential electromagnetic counterparts.

High-energy neutrinos are solely produced in the interaction of cosmic-ray nuclei with ambient matter or photon fields. In either case, both charged and neutral pions are produced. The neutral pions decay into two gamma-rays, and the charged pions produce neutrinos in their decay chain. While gamma-rays can also be produced in leptonic processes such as synchrotron emission, bremsstrahlung, and inverse Compton scattering, neutrinos are exclusively produced in hadronic processes. They are therefore considered the smoking-gun signature for the identification of cosmic-ray accelerators. Gamma-rays produced alongside high-energy neutrinos can cascade down to lower energies through interactions within the source or during propagation. Increased neutrino activity might therefore be accompanied by increased electromagnetic emission that could appear in various wavelength bands.

The first likely extragalactic neutrino counterpart is the gamma-ray blazar TXS 0506+056, which was found to be in a flaring state in spatial and temporal coincidence with the arrival of the 290 TeV neutrino event IC-190722A (Aartsen et al. 2018a) at 3σ significance. This finding motivated an archival search for lower (1–10 TeV) neutrinos from the sky position of TXS 0506+056, which resulted in the detection of a 160 day long neutrino flare in 2014/15 with significance (Aartsen et al. 2018b). Surprisingly, this archival neutrino flare was not accompanied by increased activity in gamma-ray, optical, or radio wavelengths (Aartsen et al. 2018b). Note that no dedicated follow-up campaign was performed at the time of the neutrino flare, and most of the available multiwavelength data were collected by survey instruments. Hints for a hardening of the gamma-ray spectrum during the archival neutrino flare were identified by Padovani et al. (2018) but were not found to be statistically significant () by Garrappa et al. (2019).

These two neutrino observations from the same source are difficult to reconcile through a single emission model. That the neutrino luminosity of the archival flare is more than four times larger than the gamma-ray luminosity (Aartsen et al. 2018b) suggests a hidden mechanism of neutrino production, e.g., through the attenuation of hadronic gamma-rays due to cascades initiated by photons from the jet or the broad-line region (BLR; Reimer et al. 2019; Rodrigues et al. 2019). We note that Petropoulou et al. (2020), Rodrigues et al. (2019), and Reimer et al. (2019) do not find a set of model parameters explaining the large neutrino flux from the archival neutrino flare without overshooting the electromagnetic observations. Furthermore, those hidden source scenarios (Murase et al. 2016) are inconsistent with the association of IC-170922A with a strong gamma-ray flare from TXS 0506+056 (e.g., Ansoldi et al. 2018; Cerruti et al. 2019; Gao et al. 2019). There have been, however, attempts to explain both observations in a single model by Zhang et al. (2020) and Liu et al. (2019), which come at the cost of assuming more complex geometries. The first one assumes a neutral beam scenario, while the second relies on hadronuclear interactions between protons in the jet and material in a dense gas cloud in the vicinity of the black hole.

Accordingly, establishing and understanding either of the scenarios is of significant importance.

The detection of the archival neutrino flare from the direction of TXS 0506+056 motivated a follow-up analysis (O'Sullivan & Finley 2019) searching for similar neutrino flares from the position of all sources in the third catalog of AGNs detected with the Large Area Telescope (LAT) on board the Fermi Gamma-ray Space Telescope (3LAC; Ackermann et al. 2015). The most significant neutrino flare candidate was derived for each source without accounting for electromagnetic observations of the sources. With neutrino data alone, no significant excess of flares was found above the expected atmospheric background.

A similar situation occurred in the case of TXS 0506+056, where IC-170922A with a signalness of 56% by itself was not significant, and the archival neutrino flare in 2014/15 was not found significant in an all-sky neutrino-only search. Only with the added information through multiwavelength data was it possible to identify TXS 0506+056 as the first promising candidate high-energy neutrino source.

The goal of this paper is to better understand blazars as possible source candidates for cosmic neutrinos and their emission mechanisms through the study of their electromagnetic activity. We search for coincidences of single well-reconstructed high-energy (100 TeV) neutrino events with blazars detected by the Fermi-LAT. Furthermore, we investigate the multiwavelength behavior for the most significant sources reported by O'Sullivan & Finley (2019) and sources found spatially consistent with single high-energy neutrinos.

We search for gamma-ray, X-ray, optical, and radio activity correlated with the neutrino emission.

Finally, we study the ensemble of candidate sources by testing for the generic properties expected for neutrino-emitting source populations.

We describe the sample of potential neutrino source candidates in Section 2. Section 3 describes the multiwavelength data used to compile light curves and spectral energy distributions (SEDs) for this study, and Section 4 discusses the statistical methods applied. We present our results in Section 5 and conclude in Section 6.

2. Source Sample

The neutrino sample used in this paper comes from the IceCube Neutrino Observatory, a cubic kilometer–scale Cerenkov detector located at the geographic South Pole. A complete description of the IceCube detector is provided in Aartsen et al. (2017b).

2.1. Neutrino Flare Candidates

O'Sullivan & Finley (2019) used a sample of well-reconstructed muon tracks from atmospheric and astrophysical neutrinos in the time period from 2012 April 26 to 2017 May 11. The sample covers the northern sky at declinations above −5°. The positions of 1023 sources from the 3LAC catalog are searched for a time-dependent neutrino signal in an unbinned maximum-likelihood analysis. The 11 most significant neutrino flares are reported in O'Sullivan & Finley (2019). Two of the sources, B2 1126+37 and MG2 J112910+3702, are the two possible counterparts of 3FGL J1129.0+3705, which corresponds to 4FGL J1129.1+3703 in the fourth catalog of AGNs detected by Fermi-LAT (4LAC; Fermi-LAT collaboration 2019). In 4LAC, the source is associated only with CRATES J112916+370317, which has the same coordinates as MG2 J112910+3702. We therefore only keep CRATES J112916+370317 in our sample; i.e., we study the remaining 10 sources reported by O'Sullivan & Finley (2019).

The temporal profile of each neutrino flare candidate is described by a Gaussian, and for each flare, the best-fit central value T0 and width TW of the Gaussian are reported. The latter is defined as twice the standard deviation of the Gaussian (see Table 3). O'Sullivan & Finley (2019) reported pretrial p-values for the neutrino flare candidates ranging from to , but after trial correction, none are significant. They performed a binomial test to assess the statistical significance of the ensemble yielding a p-value of 11%, which increases to 24% if TXS 0506+056 is removed from the sample, compatible with expectations from the background.

2.2. Single High-energy Neutrinos

The IceCube real-time program selects high-energy (≳100 TeV) starting and throughgoing muon track events (Aartsen et al. 2017c). A sample of real-time and archival events that would have qualified as a real-time alert but was recorded before the real-time system was operational was searched for blazar–neutrino coincidence (Aartsen et al. 2018a; Garrappa et al. 2019). Both studies focused on well-reconstructed events with a 90% containment radius of less than 5 deg2. In addition to the coincidence of IC-170922A with TXS 0506+056, the ∼100 TeV neutrino IC-141209A was identified in spatial coincidence with the BL Lac object GB6 J1040+0617. A detailed description of the multiwavelength behavior of the two sources can be found in Aartsen et al. (2018a) and Garrappa et al. (2019), respectively. Here we study the IceCube real-time alerts (see Table 1) and archival neutrino events that would have passed the same selection criteria (see Table 2). Between 2016 April and 2019 May, the IceCube collaboration operated two high-energy neutrino alert streams: the extremely high-energy (EHE) stream and the high-energy starting-track stream (HESE). In 2019 June, the alert streams were unified to the GOLD and BRONZE streams, defined by a purity of 50% and 30%, respectively. Similarly to what was done in Garrappa et al. (2019), we exclude sources with a 90% angular uncertainty larger than 5 deg2 to remove events for which no significant association would be possible given their poor localization. Such events will typically be coincident with many blazars, resulting in a poor association probability.

Table 1.  Real-time Neutrino Alerts

IceCube Alert Name Signalness Alert Type R.A.(deg) Decl.(deg) Coincident 4LAC Source/Comments GCN Circular
IC-191001A 58.9% GOLD Large angular uncertainty 25913
IC-190922B 50.5% GOLD 25806
IC-190922A 20.2% GOLD Large angular uncertainty 25802
IC-190819A 29.2% BRONZE Large angular uncertainty 25402
IC-190730A 67.2% GOLD 4FGL J1504.4+1029 25225
IC-190712A 30.3% BRONZE Large angular uncertainty 25057
IC-190704A 48.6% BRONZE Large angular uncertainty 24981
IC-190629A 33.9% BRONZE 27.22 Decl. value too close to pole for accurate error on R.A. 24910
IC-190619A 54.5% GOLD Large angular uncertainty 24854
IC-190529A 53% HESE Retracted 24674
IC-190504A 63% HESE 65.77 −37.44 No detailed angular uncertainty provided (IceCube Collaboration 2019b) 24392
IC-190503A 36.6% EHE 24378
IC-190331A 57% HESE 24028
IC-190221A 37% HESE 4FGL J1750.4–1721, 4FGL J1758.7–1621 23918
IC-190205A 84% HESE Retracted 23876
IC-190124A 91% HESE 23785
IC-190104A 35% HESE 23605
IC-181031A 87% HESE Retracted 23398
IC-181023A 28.0% EHE Large angular uncertainty 23375
IC-181014A 10% HESE Large angular uncertainty 23338
IC-180908A 34.4% EHE 23214
IC-180423A 34% HESE Retracted 22669
IC-171106A 74.6% EHE 22105
IC-171028A 30% HESE Retracted 22065
IC-171015A 51% HESE Large angular uncertainty 22016
IC-170922A 56.5% EHE 4FGL J0509.4+0542 21916
IC-170506A 35% HESE Consistent with atmospheric muon background 21075
IC-170321A 28.0% EHE 20929
IC-170312A 78% HESE Consistent with atmospheric muon background 20857
IC-161210A 49.0% EHE 20247
IC-161103A 30% HESE 4FGL J0244.7+1316 20119
IC-160814A 12% HESE Large angular uncertainty
IC-160806A 28.0% EHE 19787
IC-160731A 84.9% EHE/HESE
IC-160427A 92% HESE 19363

Note. Coordinates are reported in J2000 epoch with 90% uncertainties. Alerts with a 90% angular error larger than 5 deg2 are excluded from the analysis. For alerts shown in bold, a 4LAC source was identified located within the 90% uncertainty region. The signalness is added for completion but not used for further analysis. Alerts are taken from https://gcn.gsfc.nasa.gov/amon_hese_events.html, https://gcn.gsfc.nasa.gov/amon_ehe_events.html, and https://gcn.gsfc.nasa.gov/amon_icecube_gold_bronze_events.html.

Download table as:  ASCIITypeset image

Table 2.  Archival Neutrino Alerts

IceCube Event Name Alert Type R.A.(deg) Decl.(deg) Coincident 4LAC Source/Comments
IC-160510A EHE
IC-160128A EHE
IC-151207A HESE Bad angular resolution would have been retracted
IC-151122A EHE
IC-150926A EHE 4FGL J1258.7–0452
IC-150923A EHE
IC-150911A HESE Large angular uncertainty
IC-150831A EHE
IC-150812A EHE
IC-150428A HESE
IC-141209A HESE 4FGL J1040.5+0617
IC-141109A HESE No coincident sources
IC-140923A EHE
IC-140611A EHE
IC-140420A HESE Large angular uncertainty
IC-140203A EHE Large angular uncertainty
IC-140122A HESE Large angular uncertainty
IC-140109A EHE
IC-140108A EHE
IC-131204A EHE
IC-131202A HESE Large angular uncertainty
IC-131023A EHE
IC-130907A EHE
IC-130627A HESE No coincident sources
IC-130408A HESE Large angular uncertainty
IC-121011A EHE
IC-120922A EHE Large angular uncertainty
IC-120523A EHE
IC-120501A HESE Bad angular resolution would have been retracted
IC-120301A EHE
IC-111228A HESE Bad angular resolution would have been retracted
IC-110930A EHE
IC-110714A HESE
IC-110304A EHE
IC-110216A HESE Bad angular resolution would have been retracted
IC-110128A EHE
IC-101112A HESE
IC-101028A EHE
IC-101009A EHE
IC-100912A HESE Bad angular resolution would have been retracted

Note. Coordinates are reported in J2000 epoch with 90% uncertainties. Alerts with a 90% angular uncertainty larger than 5 deg2 are excluded from the analysis. For alerts shown in bold, a 4LAC source was identified located within the 90% uncertainty region. The signalness is added for completion but not used for further analysis. Archival events are taken from https://icecube.wisc.edu/science/data/TXS0506_alerts.

Download table as:  ASCIITypeset image

Between 2016 April and 2019 October, a total of 35 alerts were issued, and 16 survive our selection criteria (see Table 1). Forty archival events have been identified between 2010 September and 2016 May,10 and 28 pass our selection (see Table 2). Four additional coincidences are identified. The neutrino event IC-190730A was reported to be in spatial coincidence with the bright gamma-ray blazar PKS 1502+106 (IceCube Collaboration 2019a). Another coincidence was found with the HESE event IceCube-190221A and two 4FGL sources, 4FGL J1758.7–1621 and 4FGL J1750.4–1721. The first one is associated with a counterpart named AT20G J175841–161703 and classified as a blazar of uncertain type (BCU), while the second one is unassociated. The neutrino best-fit position is located just 4° from the Galactic plane, where the source density is high and the large diffuse emission complicates the detection and association. Since our work focuses on AGNs, we only consider the BCU.

In this work, we identify two additional new coincidences with archival neutrino events using 4LAC compared to the ones reported in Garrappa et al. (2019), where 3LAC was used to search for coincidences. They are IC-150926A, spatially coincident with 4FGL J1258.7–0452, and IC-161103A, spatially coincident with 4FGL J0244.7+1316. Both sources are included in 4LAC but not in 3LAC.

The ANTARES collaboration (Aublin 2019) searched for an excess of neutrinos from the positions of 3LAC sources and reported a hot spot of ANTARES neutrinos from the direction of MG3 J225517+2409. The object MG3 J225517+2409 is also spatially coincident with the 340 TeV neutrino IC-100608A with 65% signalness (event number three in Aartsen et al. 2016) and in flaring state during the IceCube neutrino arrival time (Aublin 2019). We note that IC-100608A would not have passed our selection criteria outlined above due to its large 90% angular uncertainty of roughly 30 deg2 (assuming an elliptical shape). Therefore, the source was excluded from our source sample test presented in Section 4.2. However, the temporal coincidence with IC-100608A and spatial coincidence of the ANTARES hot spot make this source interesting as a potential counterpart.

Kun et al. (2017) reported a spatial coincidence between the HESE event IC-101112A and flat-spectrum radio quasar (FSRQ) PKS 0723–008. We do not consider this source here because it lies outside the reported 90% uncertainty region,11 which is smaller than the originally published uncertainty radius in Aartsen et al. (2014).

3. Multiwavelength Data

In the following, we motivate why different wavelengths may provide relevant information connected to high-energy neutrino emission.

High-energy neutrinos are produced together with high-energy photons of similar energy in hadronic processes. The TeV to PeV photons are quickly absorbed within the source or in interactions with the extragalactic background light through photon–photon annihilation and cascade down to lower energies. Hence, GeV gamma-rays detected by Fermi-LAT provide the all-sky data set closest in energy to the neutrinos of interest. However, if the source environment is optically thick to GeV gamma-rays due to high densities of photons in the keV range, then those gamma-rays will cascade down to even lower energies, which then become an important tracer of the source activity as well. Furthermore, TeV instruments relying on the imaging atmospheric Cerenkov technique have a limited field of view and therefore even with archival data do not provide all-sky coverage. All-sky TeV instruments such as HAWC (Abeysekara et al. 2013) have limited sensitivity to extragalactic sources due to extragalactic background light (EBL) absorption.

It is possible that X-rays would be a good tracer for hadronic interactions in sources where the GeV emission is dominated by leptonic processes (Keivani et al. 2018; Gao et al. 2019).

Increased radio emission was found from TXS 0506+056 at the arrival time of IC-170922A and PKS B1424–418 in coincidence with the arrival time of a PeV neutrino (Kadler et al. 2016). We note that the chance coincidence of the neutrino association with the latter was relatively large (5%). Britzen et al. (2019) used radio data to suggest a possible collision of two jets in TXS 0506+056. However, Ros et al. (2020) excluded the presence of a secondary jet core with higher-resolution radio data but found signs of a spine-sheath structure of the jet, which could be relevant for neutrino production (see also Ghisellini et al. 2005; Tavecchio et al. 2014; Ansoldi et al. 2018).

Gamma-ray and X-ray polarization data could be used to pinpoint the leptonic and/or hadronic blazar radiation mechanisms in the high-energy bands and infer the magnetic field strength in the emission region (Zhang et al. 2019), but they are not available for the sources in our sample.

Finally, archival optical data are available for all sources in our sample. In combination with gamma-ray data, optical data can be useful to identify high-energy flares without low-energy counterparts, which could be due to hadronic interaction (Krawczynski et al. 2004).

3.1. Fermi-LAT Data

The Fermi-LAT is a pair-conversion telescope sensitive to gamma-rays with energies from MeV to greater than GeV (Atwood et al. 2009). It has a field of view >2 sr and scans the entire sky every 3 hr during standard operations. We use almost 11 yr of Pass 8 data collected between 2008 August 4 and 2019 May 30 (MJD 54,682–58,633) with an exception for the source PKS 1502+106, for which we use data up to 2019 July 31 (MJD 58,695) in order to include the arrival time of IC-190730A. We select photons from the event class developed for point-source analyses12 in the energy range from 100 MeV to 800 GeV binned into 10 logarithmically spaced energy intervals per decade. We select a region of interest (ROI) of 15° × 15° centered on the gamma-ray source position, binned in size pixels. The binning is applied in celestial coordinates using a Hammer–Aitoff projection. We perform a maximum-likelihood analysis using the standard Fermi-LAT ScienceTools package version v11r04p00 available from the Fermi Science Support Center13 (FSSC) and the P8R3_SOURCE_V2 instrument response functions, together with the fermipy package v0.17.4 (Wood et al. 2017).

We use standard data-quality cuts to select events observed when the detector was in a normal operation mode. In order to obtain a sample of events for each analysis with reduced contamination from gamma-rays produced in the Earth's upper atmosphere, we apply an additional instrument zenith angle cut of . We also remove time periods coinciding with bright solar flares and gamma-ray bursts detected by the LAT. The input model for the ROI includes all known gamma-ray sources from the 4FGL catalog in a region of 20° × 20°, slightly larger than the ROI, and the isotropic and Galactic diffuse gamma-ray emission models provided by the standard templates iso_P8R3_SOURCE_V2_v01.txt and gll_iem_v07.fits.14 The effect of energy dispersion is included in the fits performed with the Fermi-LAT ScienceTools for all point sources and the Galactic diffuse gamma-ray emission model. We use an iterative source-finding algorithm to scan the ROI and include in the model sources that are significantly () detected over the full data set time range but not over the 8 yr data that produced the 4FGL catalog. New putative point sources are modeled with a single power-law spectrum with the index fixed to 2 and the normalization free to vary in the fit. The search procedure is iterated until no further significant excess is found. The new point sources significantly detected in the longer-integration time data set are accounted for by the final ROI model.

The definition of test statistics (TS) from Mattox et al. (1996) is used to measure the detection level of each source. The minimum separation allowed between two independent point-source detections is set to . We compute the light curve for each source using the adaptive binning algorithm from Lott et al. (2012) with the prescriptions outlined in Garrappa et al. (2019), in order to better resolve the flaring activities of the target sources. Statistically significant variations in the light curve's behavior are detected in this work with the Bayesian Blocks algorithm (Scargle et al. 2013), for which we use its Astropy implementation.15 We adopt a prior that makes the algorithm sensitive to variations that are significant at a 95% confidence level.

All reported gamma-ray fluxes are in the analysis energy range from 100 MeV to 800 GeV.

3.2. Neil Gehrels Swift Observatory Data

While no sensitive all-sky X-ray monitor exists, we can take advantage of the pointed observations of the Swift X-Ray Telescope (XRT) collected in target-of-opportunity and monitoring operations. We first reprocessed the Swift-XRT data to calibrate and clean the event files using the task xrtpipeline.

The pipeline xrtgrblc was adopted to extract the source and background spectra and ancillary response files used for the light-curve generation. This tool automatically adjusts the source and background region sizes based on the source count rate.16 Due to the low photon statistics of the individual observation IDs, we fit a simple absorbed power-law model in XSPEC (Arnaud 1996) while taking the Galactic neutral hydrogen column density along the line of sight from Kalberla et al. (2005).

For the broadband SEDs, the event files were combined with xselect. Exposure maps and ancillary response files were extracted with the tasks ximage and xrtmkarf. The source region was chosen as a circle of radius centered at the target, whereas the background region has an annular shape with inner and outer radii of and , respectively, centered at the source of interest. We tested both an unbroken and a broken power law, taking into account the Galactic neutral hydrogen column density along the line of sight (Kalberla et al. 2005), and report the spectral parameters for the model that represents the data better. Depending on the source brightness, the source spectra are rebinned to have at least one or 20 counts per bin–1. The spectral analysis is performed in XSPEC.

Snapshot observations from the UltraViolet and Optical Telescope (UVOT) on board the Swift satellite during each pointing to the target source are first combined using the tool uvotimsum. To derive the source instrumental magnitude using uvotsource, we adopt a circular source region of radius centered at the object position, and a nearby source-free region of radius is considered to derive the background contamination. The computed magnitudes are converted to energy flux units using the zero-points and calibrations of Breeveld et al. (2011) corrected for the Galactic reddening following Schlafly & Finkbeiner (2011).

3.3. ASAS-SN and CSS Optical Data

Optical data in the V and g bands from the All-Sky Automated Survey for Supernovae (ASAS-SN; Shappee et al. 2014; Kochanek et al. 2017) are processed by the fully automatic ASAS-SN pipeline using the ISIS image subtraction package (Alard & Lupton 1998; Alard 2000). We then perform aperture photometry on the subtracted science image using the IRAF apphot package, adding back in the flux from the reference image. The photometry is calibrated using the AAVSO Photometric All-Sky Survey (APASS; Henden et al. 2015).

Additional V-band data from the Catalina Sky Survey  (CSS; Drake et al. 2009) are available from the public database and based on aperture photometry. To mitigate color-dependent differences between the UVOT, CSS, and ASAS-SN V-band filters, we add an offset to the ASAS-SN and UVOT data to match the CSS data in regions with overlapping exposure. A similar offset was applied to the ASAS-SN g-band observations to line them up with the V-band data. The applied shift is a constant in flux space and indicated in the legend of the corresponding light-curve figures. Since our study only relies on the shape of the light curve rather than the absolute optical flux level, and given that none of the neutrino flares occurred in the transition region between ASAS-SN and CSS data, this shift is not critical for our results.

3.4. Radio Data

Owens Valley Radio Observatory (OVRO) 15 GHz radio monitoring data (Richards et al. 2011) are available for nine sources of the sample (one of them is TXS 0506+056, which was already presented in Aartsen et al. 2018a).

3.5. Other Data

We collect archival spectral observations with the Space Science Data Center SED builder tool17 to supplement the data analyzed in this work. This allows us to cover the broadband SED of the target objects as well as possible, admittedly using nonsimultaneous data sets. However, considering that these observations represent the "average" activity of the sources, we can use them to compare the existing data acquired contemporaneously to the reported neutrino events.

4. Methods

4.1. Quantifying the Gamma-Ray Activity during the Neutrino Arrival Time

For each individual source, we calculate the chance probability

to find the neutrino in a period of gamma-ray activity larger than the gamma-ray energy flux in the time bin t overlapping with the neutrino arrival. Here is a Gaussian function with mean and standard deviation evaluated at Fx; i.e., we assume that the flux uncertainty is normally distributed. The index i runs over all time bins of the source of interest, and the pγ for all sources is reported in Table 3. Low values of pγ indicate that the source was in a high gamma-ray flux state during the neutrino arrival time compared to the other time bins in the 11 yr light curve, while high values indicate that the source did not show an excess in gamma-rays in temporal coincidence with the neutrino emission. We use the adaptive bins that were used to compile the gamma-ray light curves. Due to noncontinuous exposure and gaps in the data, we do not perform a similar analysis for optical and X-ray data. We note that the optical data show in general a similar temporal behavior to the gamma-ray data, as was found in previous studies (see, e.g., Cohen et al. 2014).

Table 3.  Neutrino Source Candidates

Source Name 4FGL Name Class Redshift T0 (MJD) Tw (days) pγ (MJD) Lγ (erg s−1)
Single High-energy Neutrinos
MG3 J225517+2409 J2255.2+2411 BL Lac 1.37a 55,355.49 0.04 [55,346.73, 55,403.54]
GB6 J1040+0617 J1040.5+0617 BL Lac 0.73b 57,000.14311 0.17 [56,997.67, 57,055.08]
1RXS J125847.7–044746 J1258.7–0452 BL Lac 0.586c 57,291.90119
GB6 J0244+1320 J0244.7+1316 BCUd 57,695.38
TXS 0506+056 J0509.4+0542 BL Lace 0.336f 58,018.87 0.009 [58,016.57, 58,019.94]
AT20G J175841–161703 J1758.7–1621 BCU 58,535.35 0.39 [58,304.43, 58,633.01]
PKS 1502+106 J1504.4+1029 FSRQ 1.839 58,694.8685 0.75 [58,603.54, 58,695.14]
Neutrino Flare Candidates
4C +20.25 J1125.9+2005 FSRQ 0.133 56,464.1 5.2 0.64 [56,369.45, 57,248.31]
CRATES J112916+370317 J1129.1+3703 BL Lac 0.445 56,501.385 0.45 [56,404.68, 57,066.59]
MG2 J112758+3620 J1127.8+3618 FSRQ 0.884 56,501.385 0.24 [56,482.90, 56,555.93]
TXS 0506+056 J0509.4+0542 BL Lace 0.336 57,000 120 0.92 [56,965.28, 57,089.28]
1H 0323+342 J0324.8+3412 NLSY1g 0.061 57,326.2938 0.08 [57,326.10, 57,333.17]
RBS 1467 J1508.8+2708 BL Lac 0.27 57,440 170 0.53 [56,474.88, 58,633.01]
S4 1716+68 J1716.1+6836 FSRQ 0.777 57,469.17919 0.48 [57,378.18, 57,510.76]
M87 J1230.8+1223 Radio galaxy 0.00428 57,730.0307 0.55 [57,724.77, 57,847.51]
GB6 J0929+5013 J0929.3+5014 BL Lac 0.37h 57,758.0 1.2 0.44 [57,647.78, 57,759.66]
1ES 0927+500 J0930.5+4951 BL Lac 0.187 57,758.0 1.2 0.49 [57,031.36, 58,633.01]

Notes.  The first seven sources were found in coincidence with single high-energy neutrinos, while the remaining sources were found coincident with neutrino flares by O'Sullivan & Finley (2019). The source classes and redshifts (if not noted otherwise) are reported in the 4LAC catalog (Fermi-LAT Collaboration 2019). Here T0 is the central time of the reported neutrino flare, Tw is twice the standard deviation of the Gaussian flare, pγ is the probability that the neutrino flare center is coincident with a gamma-flare of the found or larger flux, and is the time window used to calculate the gamma-ray flux in which the neutrino arrived. The gamma-ray luminosity is the 8 yr average calculated from the 4LAC values. Luminosity is only calculated when a redshift measurement is available. The objects 1RXS J125847.7–044746 and GB6 J0244+1320 are too dim in gamma-rays to study the variability.

aRedshift from 4LAC, which is taken from SDSS, where it is flagged as "chi-squared of best fit is too close to that of second best ( in reduced chi-squared)." Paiano et al. (2019) found that the redshift is . bRedshift from Ahn et al. (2012). cRedshift from Bauer et al. (2000). dBlazar of uncertain type. eNote that TXS 0506+056 was reclassified by Padovani et al. (2019) as "masquerading BL Lac," i.e., intrinsically an FSRQ with hidden broad lines and a standard accretion disk. fRedshift from Paiano et al. (2018). gNarrow-line Seyfert 1. hRedshift from Abazajian et al. (2009).

Download table as:  ASCIITypeset image

4.2. Comparison of Neutrino Blazar Candidates to the Gamma-Ray Blazar Sample

In addition to studying the multiwavelength behavior of individual sources, we study the average gamma-ray properties of the sources identified as potential neutrino emitters and compare them to the entire gamma-ray blazar population. Figure 1 shows the time-integrated gamma-ray energy flux in the energy range from 100 MeV to 100 GeV as a function of redshift for all blazars in 4LAC (including BCUs). All values are taken from 4LAC. We have added the redshift of four sources (see Table 3). We apply a Kolmogorov–Smirnov (K-S) test to determine how compatible the gamma-ray energy flux distribution of the candidate neutrino blazars is with the expected distribution of gamma-ray blazars under a given hypothesis. To verify that the K-S test p-value is not biased, we performed a sanity check with randomized data. We generate a background K-S p-value distribution by randomly selecting N blazars from the entire blazar sample and calculating the K-S p-value for those. Here N is the number of identified neutrino blazar candidates. A calibrated p-value for the measurement is then calculated as the ratio of background p-values smaller than the measured K-S p-value. We note that no significant bias was found, and the calibrated p-value is similar to the one obtained directly from the K-S test method.

Figure 1. Refer to the following caption and surrounding text.

Figure 1. Comparison of candidate neutrino blazars with all blazars in the 4LAC AGN sample (shown in gray). The gamma-ray energy flux is shown as a function of redshift. Sources identified in the neutrino flare search are displayed with blue circles. The object 1H 0323+342 is highlighted in cyan. Sources associated with single high-energy neutrinos are marked by colored stars. The side panels show projections of the distributions. The dashed lines in the projection panels are the values of individual blazars associated with single high-energy neutrinos, and the blue distribution shows the histogram of the neutrino flare candidate sources.

Standard image High-resolution image

We compare the observed gamma-ray energy flux of candidate neutrino blazars to the expectation for three separate scenarios. First, we test the uncorrelated case, in which all neutrino blazar coincidences occur by chance. In that case, we expect the gamma-ray flux of the candidate neutrino blazars to follow the distribution of the gamma-ray blazar population as a whole. Second, we test the hypothesis of a linear correlation of the neutrino flux with the gamma-ray energy flux of blazars. In that case, we expect the candidate neutrino blazars to have preferentially higher gamma-ray energy flux. In the third case, we test whether the neutrino flux is proportional to the square of the gamma-ray energy flux, as has been suggested in Oikonomou et al. (2019). Here we expect the candidate neutrino blazar distribution to be skewed toward even higher gamma-ray energy fluxes.

For internal consistency, the single high-energy neutrino blazar candidates are compared to the 4LAC blazar population, while the neutrino flare blazar candidates are compared to 3LAC, because only 3LAC source positions were searched for neutrino flares in O'Sullivan & Finley (2019). A large K-S p-value implies that the data are well described by a given hypothesis, while a small one indicates that that hypothesis is disfavored.

The same K-S test is also applied to the candidate neutrino flare sources. As pointed out in O'Sullivan & Finley (2019), one pair and one triplet of sources are correlated. After associating 4FGL J1129.1+3703 with CRATES J112916+370317, the triplet becomes a pair, because the second possible counterpart can be discarded. The position of CRATES J112916+370317 is correlated with MG2 J112758+3620, which is associated with the gamma-ray source 4FGL J1127.8+3618. The second correlated source positions are GB6 J0929+5013 and 1ES 0927+500 (associated with the gamma-ray sources 4FGL J0929.3+5014 and 4FGL J0930.5+4951, respectively). We recalculate the K-S test using only one of the correlated source positions and quote a range of p-values, which brackets the outcome of removing a different set of sources from the test.

Because MG3 J225517+2409 did not fulfill our angular uncertainty criteria, it was excluded from the K-S test.

The results of the K-S test are presented in Table 4 and split into BL Lacs, FSRQs, and all blazars (including BCUs) combined.

Table 4.  K-S Test p-values

  BL Lacs FSRQs All Blazars
  Uncorrelated Correlated Uncorrelated Correlated Uncorrelated Correlated
Single neutrinos 0.32 0.45 (0.0013) 0.10 0.36 (0.28) 0.126 0.64 (0.00032)
Neutrino flares 0.37–0.98 0.027–0.533 0.01–0.36 0.0075–0.023 0.39–0.98 0.0039–0.021

Note. The range of p-values for the neutrino flare case comes from removing different combinations of the correlated source positions. Different columns represent the uncorrelated and linearly correlated hypothesis; values in parentheses represent the quadratically correlated case. Note that neutrino flare candidate blazars are compared to the 3LAC population and single high-energy neutrino candidate blazars with the 4LAC population.

Download table as:  ASCIITypeset image

5. Results

5.1. Individual Sources

The collected multiwavelength light curves are presented in multipanel Figures in the Appendix. We do not show the light curves of TXS 0506+056 and GB6 J1040+0617, because they were already discussed in detail in Garrappa et al. (2019).

All sources are detected in GeV gamma-rays, which is expected, since they are selected from the 3LAC or 4FGL catalog. However, some of them are too faint to resolve temporal structure. We present both the flux and spectral index variations assuming a power-law spectrum for the source in each bin.

Most sources have good coverage in the optical during the neutrino arrival times. Radio data from the OVRO monitoring program are available for nine out of 14 sources. The X-ray data are sparse and only available for eight sources. Only 1H 0323+342 has good coverage in X-rays during the neutrino flare.

In the following, we discuss the three most interesting sources. We discuss the brightest source in gamma-rays, PKS 1502+106, and the two sources, 1H 0323+342 and MG3 J225517+2409, that show gamma-ray flares during the neutrino arrival time, reflected by small pγ of 8% and 4%, respectively, while the other sources showed p-values ranging from 17% to 92%. However, given that we have performed this calculation for 15 sources, these findings are well compatible with the background expectations.

5.1.1. 1H 0323+342

The radio-loud narrow-line Seyfert 1 galaxy 1H 0323+342 at z = 0.061 (Zhou et al. 2007; Abdo et al. 2009) shows increased gamma-ray activity during the reported neutrino flare time (see Figure 2). The gamma-ray counts map integrated over 11 yr of data is shown in Figure 3. The neutrino arrived during a mild excess in gamma-rays of ph cm−2 s−1 and roughly 1 month after a flare in the X-ray, UV, and optical (see Figure 2). The chance probability to find the neutrino in a period of increased gamma-ray activity at the level of or higher is . The neutrino flare arrives in the time bin just next to the peak. We note that the source shows even stronger flares at earlier times, which are not found connected to neutrino emission.

Figure 2. Refer to the following caption and surrounding text.

Figure 2. Multiwavelength light curve of 1H 0323+342. The duration of the neutrino flare is short (Tw = 147 s), and its arrival time is shown as an orange line. An excess in gamma-rays is found coincident with the neutrino arrival time, and an excess in X-ray emission is visible roughly 1 month before the neutrino arrival time. The Fermi-LAT gamma-ray light curve covers the energy range from 100 MeV to 800 GeV, the Swift X-ray light curve is from 0.3 to 10 keV, and the OVRO radio data are at 15 GHz.

Standard image High-resolution image
Figure 3. Refer to the following caption and surrounding text.

Figure 3. Gamma-ray count maps of 1H 0323+342, MG3 J225517+2409, and PKS 1502+106 integrated over 11 yr of Fermi-LAT data. The green cross and green line show the best-fit neutrino position and 90% uncertainty, respectively. White plus signs are 4FGL sources included in the background model. The count maps cover the energy range of 100 MeV to 800 GeV, except for 1H 0323+342, where we start at 1 GeV to suppress the significant Galactic diffuse emission at the source's latitude of . The 1H 0323+342 count map is not overlaid with a neutrino contour, since it was identified in the neutrino flare search from 3LAC sources; i.e., the neutrino flare candidate is by definition located at the position of 1H 0323+342.

Standard image High-resolution image

Figure 4 (upper left panel) shows the broadband SED of 1H 0323+342.

Figure 4. Refer to the following caption and surrounding text.

Figure 4. SEDs of 1H 0323+342, MG3 J225517+2409, and PKS 1502+106. Archival data are shown in orange. The 11 yr Fermi-LAT SEDs are overlaid as blue circles. The black data points refer to the observations contemporaneous to the epoch of neutrino detection. The lower right panel shows the SEDs of PKS 1502+106 during four selected epochs (see Table 5 for details).

Standard image High-resolution image

The X-ray spectrum reveals a break at ∼3 keV. The spectrum before the break energy is soft (; see Table 5), possibly due to coronal emission (Abdo et al. 2009; Paliya et al. 2014, 2019). Note that various spectral features are detected in the X-ray spectrum of this source, which includes a soft X-ray excess below 2 keV, an Fe K-alpha emission line at ∼6 keV, and a possible Compton hump at higher frequencies (see, e.g., Paliya et al. 2014, 2019; Ghosh et al. 2018; Kynoch et al. 2018 for details). Covering these aspects is beyond the scope of this work. Furthermore, the broadband SED modeling of this object suggests that the gamma-ray emission region lies well within the BLR, i.e., close to the central black hole (Abdo et al. 2009; Paliya et al. 2014; Kynoch et al. 2018). If so, the X-ray photons from the corona could constitute a target photon field for photohadronic interactions producing high-energy neutrinos. In particular, the interaction of the protons with the thermal continuum with a characteristic temperature (T*) can produce neutrinos with energy TeV ( K)−1 (see, e.g., Rodrigues et al. 2019). The X-ray coronal photons would also absorb the gamma-rays via pair production, leading to the steepening of the gamma-ray spectrum that is observed (Figure 4; Ghisellini & Tavecchio 2009; Paliya et al. 2014; Rodrigues et al. 2019). Another observational signature for this process is the detection of a bright X-ray emission with a soft spectral shape (see Ghisellini & Tavecchio 2009 for details), which is reflected in the X-ray spectrum of 1H 0323+342 (Figure 4). A quantitative discussion will be the subject of a separate publication.

Table 5.  Summary of SED Analysis

      Fermi-LAT      
Name Time Window Flux Power-law Index log-parabola Indices TS
  MJD 10−7 ph cm−2 s−1   α β  
11 yr Averaged
1H 0323+342 54,682−58,633 0.45 ± 0.02 2.77 ± 0.04 0.09 ± 0.04 1027
PKS 1502+106 54,682−58,695 2.97 ± 0.02 2.12 ± 0.01 0.10 ± 0.01 95412
MG3 J225517+2409 54,682−58,633 0.11 ± 0.01 2.03 ± 0.04 955
Contemporaneous
1H 0323+342 57,263−57,392 1.24 ± 0.21 3.25 ± 0.47 0.24 ± 0.25 90
PKS 1502+106 58,664−58,724 0.86 ± 0.22 2.31 ± 0.16 0.01 ± 0.07 97
MG3 J225517+2409 55,346−55,501 0.41 ± 0.08 2.02 ± 0.09 181
Interesting Multiwavelength Features of PKS 1502+106
PKS 1502+106 54,682−54,692 19.48 ± 0.81 1.87 ± 0.04 0.12 ± 0.02 4205
  55,266−57,022 0.88 ± 0.03 2.27 ± 0.02 0.07 ± 0.01 56318
  57,210−57,219 14.15 ± 0.69 1.74 ± 0.05 0.09 ± 0.02 4830
  58,107−58,125 4.80 ± 0.35 2.16 ± 0.06 0.06 ± 0.04 25676
Swift-XRT
Name Exposure Flux Normalization /dof
  ks     10−12 erg cm−2 s−1 10−4 ph cm−2 s−1 keV−1  
Contemporaneous
1H 0323+342 32.81 366.91/285
PKS 1502+106 5.61   56.59/66
MG3 J225517+2409 1.65   12.57/12
Interesting Multiwavelength Features of PKS 1502+106
PKS 1502+106 36.67   63.62/69
  13.50   12.51/7
  11.06   20.60/12
  14.69   42.60/29
Swift-UVOT
Name V B U W1 M2 W2
Contemporaneous
1H 0323+342 19.58 ± 0.23 19.85 ± 0.19 22.82 ± 0.26 21.58 ± 0.30 24.94 ± 0.38 23.47 ± 0.31
PKS 1502+106     0.90 ± 0.06   0.67 ± 0.06  
MG3 J225517+2409 2.91 ± 0.32 1.90 ± 0.13 2.26 ± 0.12 1.92 ± 0.12 2.74 ± 0.25 1.54 ± 0.10
Interesting Multiwavelength Features of PKS 1502+106
PKS 1502+106 4.71 ± 0.11 4.45 ± 0.08 3.97 ± 0.08 2.59 ± 0.06 2.67 ± 0.07 2.21 ± 0.06
  0.55 ± 0.06 0.67 ± 0.05 0.61 ± 0.06 0.55 ± 0.05 0.51 ± 0.04 0.43 ± 0.03
  6.38 ± 0.24 6.24 ± 0.18 5.32 ± 0.16 3.71 ± 0.14 3.57 ± 0.11 2.98 ± 0.11
  8.06 ± 0.20 7.20 ± 0.15 6.42 ± 0.15 5.03 ± 0.16 4.68 ± 0.13 4.34 ± 0.10

Note. The first block shows the Fermi-LAT SED results performed in a given time window. The quoted gamma-ray flux is integrated in the 0.1−800 GeV energy range. The SED is modeled with a power law, unless a log-parabola description with spectral parameters α and β results in a significantly better fit to the data. The second block shows the results of the Swift-XRT spectral analysis. Here is the photon index of a power-law model or photon index before the break energy in a broken power-law model, while is the photon index after the break energy in a broken power-law model. The flux is integrated in the 0.3−10 keV energy range, and the normalization is defined at 1 keV in units of 10−4 ph cm−2 s−1 keV−1. Absorption by Galactic neutral hydrogen is taken into account using the following column densities along the line of sight (Kalberla et al. 2005): (1H 0323+342), (PKS 1502+106), and cm−2 (MG3 J225517+2409). The third block shows the results of the Swift-UVOT analysis and gives the average flux in the Swift V, B, U, W1, M2, and W2 bands in units of 10−12 erg cm−2 s−1. The Swift analyses were performed in the same time periods specified for the Fermi-LAT results.

Download table as:  ASCIITypeset image

5.1.2. MG3 J225517+2409

The distant BL Lac object MG3 J225517+2409 shows a major flare coincident with the neutrino arrival time (see Figure A1). A redshift of 1.37 (Fermi-LAT Collaboration 2019), which is taken from SDSS, is reported by 4LAC. However, the extracted redshift is flagged as "chi-squared of best fit is too close to that of second best ( in reduced chi-squared)." Paiano et al. (2019) found that the redshift is . The gamma-ray flare reaches a flux level of ph cm−2 s−1 and lasts roughly 140 days (see Figure A1). The chance probability to find the neutrino in a period of increased gamma-ray activity at this level or higher is .

Figure 4 (upper right panel) shows the broadband SED of MG3 J225517+2409, and the best-fit spectral values for the gamma-ray, X-ray, and UV bands are provided in Table 5.

5.1.3. PKS 1502+106

The FSRQ PKS 1502+106 was found to be located within the 50% uncertainty region of IC-190730A. The neutrino was reported with a signalness of 67% and an energy of 300 TeV (IceCube Collaboration 2019a). It is the 15th brightest out of 2863 sources in the 4LAC catalog in terms of gamma-ray energy flux at despite its large redshift of 1.84 (Hewett & Wild 2010), suggesting an extremely high intrinsic luminosity.

It was found to be in a low-activity state during the arrival time of the high-energy neutrino (see Figures 4 and A5). However, the OVRO radio light curve of PKS 1502+106 shows a long-term outburst starting in 2014 and reaching the highest flux density ever reported from this source (since the beginning of the OVRO measurements in 2008) during the arrival of the 300 TeV neutrino IC-190730A (Kiehlmann et al. 2019). The object TXS 0506+056 showed a similar increase in the radio emission observed by OVRO in coincidence with IC-170922A (Aartsen et al. 2018a; Kiehlmann et al. 2019). A strong increase in radio emission was also determined in very long baseline interferometry (VLBI) data for another blazar, PKS B1424–418, which was found coincident with a high-energy but poorly reconstructed neutrino event (Kadler et al. 2016). Plavin et al. (2020) found a correlation of IceCube neutrinos with radio-bright AGNs with a 0.2% p-value. Quantifying the chance coincidence of a radio flare with the arrival time of a neutrino is outside the scope of this paper.

Figure 4 (lower left panel) shows the broadband SED of PKS 1502+106. The 11 yr averaged gamma-ray spectrum of this source reveals a significant curvature/break that could be reflecting the shape of the particle spectrum or due to extrinsic absorption by the BLR photons. Interestingly, the EBL absorption is not significant below 50 GeV at z = 1.84 (Fermi-LAT Collaboration et al. 2018), and the gamma-ray emission from PKS 1502+106 has been explained by the interaction of the jet electrons with the BLR photons (e.g., Abdo et al. 2010). Therefore, the observed spectral curvature could be due to gamma-ray absorption by the BLR photons via the pair-production process and/or a transition from the Thomson to Klein–Nishina regime. If so, the same BLR photon field could also act as a target photon field for neutrino production by interacting with the hadrons present in the jet (see, e.g., Rodrigues et al. 2019).

The gamma-ray spectral index shows a variation in time (see Figure A5, second panel). The spectrum tends to harden when the gamma-ray flux increases. The hard spectral regions indicate an increase in high-energy emission and are therefore promising targets for follow-up searches of (TeV) neutrinos. Since PKS 1502+106 is the most interesting source of our sample (due to its high gamma-ray energy flux), we study the spectral behavior during the multiwavelength flares in more detail to give guidance for future neutrino searches. We split the 11 yr light curve into four regions of interest, where we obtain the gamma-ray spectral shape (see Table 5). We select one period from MJD 55,266–57,022 to cover the quiet state in gamma-rays and three short periods of roughly 10 days in length chosen to cover the three bright X-ray flares, which are also accompanied by optical flares (see Figure A5). We find that during the gamma-ray quiet state, the flux values in each wavelength reach a minimum flux level (shown in green). Interestingly, the highest flare in the X-ray and optical (cyan) does not correspond to the highest flare in gamma-rays, while the highest gamma-ray activity is also accompanied by a significant increase in the optical and X-ray. The different flaring behavior indicates different conditions of the emission region in the source. Detailed time-dependent modeling, which is outside the scope of this work, could give deeper insight into the variable nature of the source.

5.2. Source Population

As a result of the K-S test (see Table 4), we find that the blazars found in coincidence with single high-energy neutrinos are well described by both the gamma-ray energy flux distribution expected in the case of a linear correlation between neutrino and gamma-ray energy flux and the hypothesis of no correlation between the two fluxes. If all sources are combined, the single neutrino source candidates are compatible with the no-correlation hypothesis with a p-value of and consistent with the expectation in the case of a linear correlation between neutrino and gamma-ray energy flux (p-value of 64%). They show a mismatch with the hypothesis of a quadratical correlation (p-value of ). The gamma-ray energy flux distribution is illustrated in Figure 5.

Figure 5. Refer to the following caption and surrounding text.

Figure 5. Expected gamma-ray energy flux (in the 100 MeV–100 GeV energy range) distribution of all blazars in 3LAC (left) and 4LAC (right) according to the no-correlation hypothesis shown in blue. The expected distribution according to the hypothesis of a linear (quadratical) correlation between neutrino and gamma-ray energy flux is shown in green (cyan). Left: observed gamma-ray energy flux of neutrino flare candidate blazars from 3LAC overlaid in red. The red distribution matches the blue one but not the green or cyan. Right: observed gamma-ray energy flux of 4LAC blazars found in spatial coincidence with single high-energy neutrinos (red).

Standard image High-resolution image

The object MG3 J225517+2409 would have failed our angular uncertainty requirement selection and was not included in the K-S test.

At the same time, the candidate neutrino flare sources show a good match with the random distribution (p-value of 39%–98%) but are less well described by the energy flux-weighted distribution (p-value of 0.4%–2.1%).

6. Summary and Conclusions

6.1. Individual Sources

In summary, the available data do not show any significant temporal correlation with the neutrino arrival time considering all blazars studied here, which would allow us to identify one of the sources as a potential cosmic-ray acceleration site. This is consistent with findings by Righi et al. (2019), who studied the gamma-ray light curves of seven BL Lac objects and did not find a clear pattern in common among the sources.

Most sources are well observed in GeV gamma-rays and optical wavelengths, where no significant temporal correlation with the neutrino emission is found. For the sources monitored by OVRO, no short-term features related to the neutrino arrival time are observed. Three out of five sources coincident with single high-energy neutrinos are monitored by OVRO in radio, and two (TXS 0506+056 and PKS 1502+106) show a long-term increase of the radio flux density, which peaks during the neutrino arrival time. The third one (MG3 J225517+2409) is only covered by OVRO observations 70 days after the neutrino arrival time but might be compatible with a radio flux increase assuming a smooth extrapolation of the temporal behavior to earlier times. The five neutrino flare source candidates, which are monitored in radio, show no correlation with long-term radio activity. Radio monitoring of future neutrino blazar coincidences could reveal whether there is indeed a connection between the radio and single high-energy neutrino emission.

6.2. Source Population

We find that the single high-energy neutrino coincidences with blazars are consistent with a p-value of with being due to random chance. At the same time, they are consistent with the hypothesis that the single high-energy neutrino emission is correlated linearly with the gamma-ray brightness of the blazars. This interpretation is consistent with the association of IC-170922A and the bright gamma-ray flare of TXS 0506+056.

We note that the contribution of gamma-ray blazars to the diffuse neutrino flux is constrained by blazar stacking analyses (Aartsen et al. 2017d), which constrain the blazar contribution to less than 27% under the assumption that the neutrino spectrum follows an unbroken power law. Assuming a steeper power law with an index of 2, which is compatible with the diffuse flux fit above ∼200 TeV, weakens the constraints on the diffuse flux contribution to 40%–80% (Aartsen et al. 2018b) and leaves room for a significant contribution from gamma-ray blazars. Krauß et al. (2018) studied the gamma-ray and X-ray emission of 3LAC blazars in the vicinity of high-energy starting events detected by IceCube and found no direct correlation between the Fermi-LAT gamma-ray flux and the IceCube neutrino flux.

Assuming that the observed linear correlation of single high-energy neutrinos and the average gamma-ray energy flux is genuine, our results would have many broad implications. They would allow us to utilize neutrino blazar coincidences for the study of cosmic-ray acceleration. Furthermore, they would allow for an effective search for more coincidences to further characterize the population of sources of high-energy neutrinos.

At the same time, the candidate neutrino flare–emitting blazars are compatible with the background hypothesis, and the data do not support the hypothesis that neutrino emission is correlated to the average gamma-ray energy flux. This could indicate that either neutrino flares are not accompanied by strong gamma-ray emission or these coincidences are of a random nature. The first case could be realized in sources where neutrinos are produced in regions optically thick to gamma-rays, where gamma-ray emission is absorbed (so-called hidden sources) and cascades to the X-ray band, where we do not have good observational constraints from archival data (see also Keivani et al. 2018; Gao et al. 2019).

This work was supported by the Initiative and Networking Fund of the Helmholtz Association. We thank Xavier Rodrigues, Shan Gao, and Walter Winter for fruitful discussions. B.J.S. is supported by NSF grants AST-1908952, AST-1920392, and AST1911074. S.K. acknowledges support from the European Research Council under the European Union's Horizon 2020 research and innovation program, under grant agreement No. 771282. W.M. acknowledges support from CONICYT project Basal AFB-170002. T.H. was supported by Academy of Finland projects 317383 and 320085.

Fermi-LAT: The Fermi-LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT, as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States; the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France; the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy; the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK), and Japan Aerospace Exploration Agency (JAXA) in Japan; and the K. A. Wallenberg Foundation, the Swedish Research Council, and the Swedish National Space Board in Sweden.

Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France. This work was performed in part under DOE contract DE-AC02-76SF00515.

ASAS-SN: ASAS-SN is supported by the Gordon and Betty Moore Foundation through grant GBMF5490 to the Ohio State University and NSF grants AST-1515927 and AST-1908570. Development of ASAS-SN has been supported by NSF grant AST-0908816, the Mt. Cuba Astronomical Foundation, the Center for Cosmology and AstroParticle Physics at the Ohio State University, the Chinese Academy of Sciences South America Center for Astronomy (CAS-SACA), the Villum Foundation, and George Skestos.

Others: This research has made use of data from the OVRO 40 m monitoring program (Richards et al. 2011), which is supported in part by NASA grants NNX08AW31G, NNX11A043G, and NNX14AQ89G and NSF grants AST-0808050 and AST-1109911.

The CSS survey is funded by the National Aeronautics and Space Administration under grant No. NNG05GF22G issued through the Science Mission Directorate Near-Earth Objects Observations Program. The CRTS survey is supported by the U.S. National Science Foundation under grants AST-0909182 and AST-1313422.

Part of this work is based on archival data, software, or online services provided by the Space Science Data Center (SSDC).

Software: Fermi-LAT ScienceTools (v11r04p00), fermipy (v0.17.4; Wood et al. 2017), astropy (Astropy Collaboration et al. 2013, 2018), IRAF (Tody 1986, 1993), ISIS (Alard & Lupton 1998; Alard 2000), XSPEC (Arnaud 1996), HEASoft (v6.26.1).

Appendix: Multiwavelength Light-curve Plots

Here, we present multiwavelength light curves (Figures A1A14) for all sources in our sample.

Figure A1. Refer to the following caption and surrounding text.

Figure A1. Multiwavelength light curve of MG3 J225517+2409. The orange line indicates the arrival time of the high-energy neutrino IC-100608A.

Standard image High-resolution image
Figure A2. Refer to the following caption and surrounding text.

Figure A2. Multiwavelength light curve of 1RXS J125847.7–044746. The orange line indicates the arrival time of the high-energy neutrino IC-150926A. The source is too dim in gamma-rays to resolve the temporal variability.

Standard image High-resolution image
Figure A3. Refer to the following caption and surrounding text.

Figure A3. Multiwavelength light curve of GB6 J0244+1320. The orange line indicates the arrival time of the high-energy neutrino IC-161103A. The source is too dim in gamma-rays to resolve the temporal variability.

Standard image High-resolution image
Figure A4. Refer to the following caption and surrounding text.

Figure A4. Multiwavelength light curve of AT20G J175841–161703. The orange line indicates the arrival time of the high-energy neutrino IC-190221A.

Standard image High-resolution image
Figure A5. Refer to the following caption and surrounding text.

Figure A5. Multiwavelength light curve of PKS 1502+106. The orange line indicates the arrival time of the high-energy neutrino IC-190730A. The green region in the first panel marks the quiescent state, and vertical red, cyan, and black lines mark three flaring states, which are accompanied by X-ray and optical flares, selected for a dedicated gamma-ray spectral analysis (see Table 5 and Figure 4, which uses the same color code).

Standard image High-resolution image
Figure A6. Refer to the following caption and surrounding text.

Figure A6. Multiwavelength light curve of 4C +20.25. The duration of the neutrino flare is short (Tw = 5.2 days), and its arrival time is shown as an orange line.

Standard image High-resolution image
Figure A7. Refer to the following caption and surrounding text.

Figure A7. Multiwavelength light curve of CRATES J112916+370317. The duration of the neutrino flare is short (Tw = 1.4 hr), and its arrival time is shown as an orange line.

Standard image High-resolution image
Figure A8. Refer to the following caption and surrounding text.

Figure A8. Multiwavelength light curve of MG2 J112758+3620. The duration of the neutrino flare is short (Tw = 1.4 hr), and its arrival time is shown as an orange line.

Standard image High-resolution image
Figure A9. Refer to the following caption and surrounding text.

Figure A9. Multiwavelength light curve of 1H 0323+342. The duration of the neutrino flare is short (Tw = 147 s), and its arrival time is shown as an orange line.

Standard image High-resolution image
Figure A10. Refer to the following caption and surrounding text.

Figure A10. Multiwavelength light curve of RBS 1467. The orange shaded region is centered on the mean of the Gaussian describing the neutrino emission with a width corresponding to twice the standard deviation.

Standard image High-resolution image
Figure A11. Refer to the following caption and surrounding text.

Figure A11. Multiwavelength light curve of S4 1716+68. The duration of the neutrino flare is short (Tw = 4.7 s), and its arrival time is shown as an orange line.

Standard image High-resolution image
Figure A12. Refer to the following caption and surrounding text.

Figure A12. Multiwavelength light curve of M87. The duration of the neutrino flare is short (Tw = 3.9 minutes), and its arrival time is shown as an orange line.

Standard image High-resolution image
Figure A13. Refer to the following caption and surrounding text.

Figure A13. Multiwavelength light curve of GB6 J0929+5013. The duration of the neutrino flare is short (Tw = 1.2 days), and its arrival time is shown as an orange line.

Standard image High-resolution image
Figure A14. Refer to the following caption and surrounding text.

Figure A14. Multiwavelength light curve of 1ES 0927+500. The duration of the neutrino flare is short (Tw = 1.2 days), and its arrival time is shown as an orange line.

Standard image High-resolution image

Footnotes

10.3847/1538-4357/ab8307
undefined