ANTARES Search for Point Sources of Neutrinos Using Astrophysical Catalogs: A Likelihood Analysis

An extended maximum likelihood method (Barlow 1990) is used to reject the null hypothesis H₀ where only background is assumed to be present, from an alternative hypothesis H₁ of signal plus background. The log-likelihoods for both hypotheses H₀ and H₁ are written respectively as

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{ln}{ \mathcal L }(x| {{\rm{H}}}_{0}) & = & \displaystyle \sum _{i}^{N}\mathrm{ln}\left[{\mu }_{b}B({x}_{i})\right]-{\mu }_{b}\\ \mathrm{ln}{ \mathcal L }(x| {{\rm{H}}}_{1}) & = & \displaystyle \sum _{i}^{N}\mathrm{ln}\left[{\mu }_{s}S({x}_{i})+{\mu }_{b}B({x}_{i})\right]-{\mu }_{s}-{\mu }_{b},\end{array}\end{eqnarray} \tag{ 1 }$$

where N is the total number of observed neutrino candidates, S(x_i ) is the probability density function (PDF) of the signal, and B(x_i ) the PDF of the background. The free parameters are the estimated numbers of signal μ_s and background μ_b events. The test statistic is then defined as

$$\begin{eqnarray}&&\lambda =\mathrm{ln}\left(\displaystyle \frac{\max ({ \mathcal L }(x| {{\rm{H}}}_{1}))}{\max ({ \mathcal L }(x| {{\rm{H}}}_{0}))}\right).\end{eqnarray} \tag{ 2 }$$

3.2.1. Signal

For a given catalog, the signal term is proportional to the expected number of neutrinos that should be detected as track-like events by ANTARES, from the directions of the N_sources objects, each one of them emitting a ${\nu }_{\mu }+{\bar{\nu }}_{\mu }$ neutrino flux ${{\rm{\Phi }}}_{j}(E)={{\rm{\Phi }}}_{j}^{0}{\left(\tfrac{E}{1{\rm{GeV}}}\right)}^{-\gamma }$, the flux normalization ${{\rm{\Phi }}}_{j}^{0}$ (GeV⁻¹ cm⁻² s⁻¹) being different for each source j. The spectral index γ is assumed to be γ = 2.0 in the analysis, and is fixed in the likelihood maximization.

The signal PDF is then written as

$$\begin{eqnarray}&&S({x}_{i})=\displaystyle \frac{1}{\sum {w}_{j}}\displaystyle \sum _{j=1}^{{N}_{{\rm{sources}}}}{w}_{j}\,{s}_{j}({x}_{i}),\end{eqnarray} \tag{ 3 }$$

where s_j (x_i ) is the contribution of the jth object, with weight w_j . The weight of the jth source is defined as

$$\begin{eqnarray}&&{w}_{j}={w}_{j}^{\mathrm{model}}\,{ \mathcal A }({\delta }_{j}),\end{eqnarray} \tag{ 4 }$$

where ${ \mathcal A }(\delta )$ is the declination-dependent acceptance of the neutrino track sample computed for an E⁻² energy spectrum.

The function ${ \mathcal A }(\delta )$ allows the computation of the number of signal events expected for a single point source having a flux normalization Φ₀ simply via the relation

$$\begin{eqnarray}&&n(\delta )={{\rm{\Phi }}}_{0}\,{ \mathcal A }(\delta ).\end{eqnarray} \tag{ 5 }$$

For a given energy spectrum ${\rm{\Phi }}(E)={{\rm{\Phi }}}_{0}{\left(\tfrac{E}{1{\rm{GeV}}}\right)}^{-\gamma }$, the acceptance is computed from the effective area (Albert et al. 2017) of the detector, ${{ \mathcal A }}_{\mathrm{eff}}(E,\delta ,t)$, with the relation

$$\begin{eqnarray}&&{ \mathcal A }(\delta )=\int \int {\rm{d}}t\,{\rm{d}}E\,{{ \mathcal A }}_{{\rm{eff}}}(E,\delta ,t){\left(\displaystyle \frac{E}{1{\rm{GeV}}}\right)}^{-\gamma },\end{eqnarray} \tag{ 6 }$$

where the integral is performed over a sufficiently large energy range (100 GeV−10 PeV), and for the total livetime considered in the analysis (3125.4 days). The acceptance obtained for a pure E⁻² energy spectrum is shown in Figure 5 of Albert et al. (2017).

The weighting scheme aims to set the weights to be proportional to the true neutrino flux. Since this is unknown, two different weighting assumptions are considered:

• 1.  
  ${w}_{j}^{\mathrm{model}}={{\rm{\Phi }}}_{j}^{0}$
• 2.  
  ${w}_{j}^{\mathrm{model}}=1$ (equal-flux assumption),

where ${{\rm{\Phi }}}_{j}^{0}$ is the photon flux observed in a relevant part of the electromagnetic spectrum, to be used as a proxy for the neutrino flux. As the relation between the flux emitted in a given band of the electromagnetic spectrum and the neutrino flux is rather uncertain, the equal-flux assumption (${w}_{j}^{\mathrm{model}}=1$) provides an additional model-independent test, which could be more sensitive if the weights in the likelihood are very different from the true (unknown) ones.

The individual source contribution to the signal for the ith event is

$$\begin{eqnarray}&&{s}_{j}({x}_{i})={ \mathcal F }({\alpha }_{{ij}},{E}_{i},{\beta }_{i}),\end{eqnarray} \tag{ 7 }$$

where ${ \mathcal F }$ is the energy-dependent point-spread function, defined as the probability density for the reconstructed event direction to lie within an angular distance α_ij from the true source direction.

Figure 1 shows that on average, for an E⁻² spectrum, 50% of the track events in the point-source sample are reconstructed within an angle of 0.4° of the true neutrino direction.

The most relevant parameters that determine the shape of the point-spread function are the estimated energy E_i and the angular reconstruction uncertainty β_i coming from the track reconstruction procedure.

The function ${ \mathcal F }$ is evaluated using a collection of 3D histograms built from Monte Carlo simulations, in declination bands of ${\rm{\Delta }}\sin \delta =0.2$. For a given track event with declination δ, energy E, and angular uncertainty β, the histogram corresponding to the correct declination band is used, and the value of the signal is obtained by linear interpolation in the (α, E, β) space.

3.2.2. Background

A conservative approach is used to determine the background term. As the contribution of an astrophysical signal is expected to be small, real data are used to build the background PDF. The limited statistics available in real data as opposed to Monte Carlo do not allow the construction of a three-dimensional histogram, therefore the expression of the PDF is factorized as

$$\begin{eqnarray}&&B({x}_{i})={ \mathcal B }(\sin {\delta }_{i})\,{f}_{{\rm{b}}}({E}_{i},{\beta }_{i}),\end{eqnarray} \tag{ 8 }$$

where ${ \mathcal B }(\sin \delta )$ is the observed distribution of the sine of the declination, and f_b the distribution of energy E and angular uncertainty β derived in data. The distribution of ANTARES track events measured during more than ten years is independent of the right ascension (see Figure 2 top left), because of the Earth's rotation averaging the nonuniformity of the detector's exposure in local coordinates.

Zoom In Zoom Out Reset image size

Figure 2 (bottom left) shows the rate of events as a function of $\sin \delta$, fitted by a polynomial function that is used in the likelihood computation. The value of the function f_b(E, β) is evaluated by linear interpolation in the 2D histogram shown in Figure 2 (right).

The Fermi 3LAC catalog (Ackermann et al. 2015) contains blazars detected in gamma rays in the range 1–300 GeV by the Fermi-LAT in four years of operation. The previous release (Fermi 2LAC) has been used by the IceCube Collaboration in a stacking analysis (Aartsen et al. 2020a), where several subclasses of objects have been tested. A similar approach is followed here: in addition to the full Fermi 3LAC catalog, two additional subsamples are also considered, as described below.

Blazars have a typical spectral energy distribution (SED) of their electromagnetic spectrum that presents two bumps: a low-energy peak between IR and UV–X due to synchrotron emission of electrons, and a high-energy peak in the hard X-ray/gamma-ray band that can be produced by both hadronic and leptonic processes. Blazars are classified into flat-spectrum radio quasars (FSRQs) and BL Lacs, according to the presence (FSRQs) or absence (BL Lacs) of optical broad emission lines on top of the continuum. When this classification is not possible, a blazar is labelled as a blazar candidate of uncertain type (BCU) in the catalog.

In addition, a second classification can be made using the position of the synchrotron peak to distinguish low synchrotron peak (LSP) from intermediate and high synchrotron peaks (ISP and HSP respectively). This classification was proposed as an extension of the LBL/HBL classification of BL Lacs originally introduced by Padovani & Giommi (1995).

In addition, the objects that are flagged as doubtful in the Fermi catalog are removed from the analysis (see Ackermann et al. 2015, Section 2). The different (sub)samples that could be used in the analysis are:

• 1.  
  Fermi 3LAC: 1420 objects, with 1255 in the field of view (FoV) of ANTARES,
• 2.  
  FSRQ: 414 objects (382 in FoV),
• 3.  
  BL Lacs: 604 objects (509 in FoV),
• 4.  
  LSP: 666 objects (604 in FoV),
• 5.  
  ISP+HSP: 685 objects (590 in FoV).

The overlap between subsamples is indicated in Table 1.

As can be seen, most of the FSRQs are classified as LSP and the majority of BL Lacs are classified as ISP+HSP. In order to reduce the number of trials, the analysis is then restricted to the FSRQ and BL Lac subsamples only, in addition to the full Fermi 3LAC sample.

For the implementation in the likelihood, the measured slope γ of the gamma-ray energy spectrum is taken into account in the weights. This is done by using the energy flux: ${w}^{{\rm{model}}}=\int E\times {{\rm{\Phi }}}_{0}{\left(E/{E}_{0}\right)}^{-\gamma }\,{dE}$ where the integral runs over the energy range 1–100 GeV.

The authors of Maggi et al. (2016) suggest considering a specific class of AGNs where the jet axis points toward the Earth but passes through a large quantity of dust or gas, thus potentially enhancing the neutrino production rate via hadronic p–p interactions.

The sample is built by cross-correlating objects from a catalog of radio-emitting galaxies (van Velzen et al. 2012) with the Fermi 2LAC catalog (Ackermann et al. 2011). The morphologies that are incompatible with a jet viewed face-on are rejected. The final selection criterion uses the observed X-ray flux of the objects, which should be low compared to the radio flux because of the attenuation in the surrounding matter.

The 25% weakest X-ray emitters are then selected, leading to a small catalog of 15 objects, containing three FSRQs, one ultraluminous IR galaxy, and 11 BL Lacs.

To take into account the distance and the intrinsic luminosity, the radio flux density measured at 1.4 GHz is used as a weight in the likelihood: w^model = Φ_Radio.

Radio galaxies are a subclass of active galaxies that could be considered as possible sources of high-energy neutrinos (Hooper 2016). The neutrinos could be produced by the hadronic interaction of accelerated cosmic rays within the jets or in the giant lobes at the ends of the jets.

Radio galaxies have an active nucleus emitting jets of relativistic plasma, but contrary to blazars, the direction of the jet axis is not close to our line of sight. Radio galaxies are then in general less luminous than blazars individually, but they are more abundant in the local universe.

A sample of 65 hard X-ray selected radio galaxies suggested by L. Bassani is considered in this study. The authors of Bassani et al. (2016) looked at extended radio galaxies in the sample of soft gamma-ray (20–100 keV) selected AGNs in either the INTEGRAL or Swift-BAT survey. The selection procedure yields a sample that contains the brightest and most accretion-efficient radio galaxies in the local sky, with very powerful jets producing an extended double-lobed radio morphology.

Among the published list of objects, one source without redshift information is removed, and Centaurus A is excluded because it cannot be considered as a point source (the angular extent of its radio lobes is ∼10°).

The weight w^model = L_γ /d² is used, where L_γ is the luminosity measured in soft gamma rays by Swift-BAT (except for five galaxies where the luminosity measured by INTEGRAL is used) and d is the luminosity distance to the source. The luminosity distance is computed for each galaxy adopting the redshift indicated in Bassani et al. (2016), using a standard flat ΛCDM cosmology with parameters H₀ = 71 km s⁻¹ Mpc⁻¹, Ω_m = 0.27, and Ω_Λ = 0.73 (Hinshaw et al. 2013).

A catalog of 64 star-forming galaxies (SFGs) published by the Fermi-LAT Collaboration (Ackermann et al. 2012) was used recently by the Pierre Auger Collaboration, reporting a 3.9σ correlation with ultrahigh-energy cosmic rays with energy above 39 EeV (Aab et al. 2018). Each SFG does not necessarily have an active nucleus, but has a strong star formation rate traced by the intense infrared radiation. The dense interstellar gas can act as target material for the cosmic rays that are trapped in the galactic magnetic field, thus potentially producing neutrinos.

The dominant galaxies in terms of emitted infrared flux are M82, NGC 253, and NGC 4945. Due to the selection procedure (detection in gamma rays by Fermi-LAT), most of the objects present in the catalog are very nearby: five objects are located within a radius of 5 Mpc, 14 galaxies are within 10 Mpc, and 51 within 100 Mpc. A caveat of the use of this catalog is then its incompleteness due to the gamma-ray selection.

The following weights are used in the likelihood: w^model = L_IR/d² where L_IR is the total IR luminosity (8–1000 μm) and d is the luminosity distance to the source.

The high-energy tracks measured by IceCube are an interesting sample for a stacking analysis: they are likely to be of astrophysical origin and have for the majority a median angular resolution of ∼1°. Their directions are then in general a sufficiently good proxy for the location of neutrino sources.

The sample consists of 56 events, 55 of which are in the FoV of ANTARES: 35 tracks from the IceCube 8 yr sample of up-going muons with deposited energy above 200 TeV (Haack & Wiebusch 2017) and 21 additional tracks from the 6 yr HESE sample (Kopper 2017). The majority of those events are located in the Northern Hemisphere, and have been reconstructed with angular errors varying from ∼0.5° up to ∼5°. A complementary study of the correlation between IceCube high-energy tracks and ANTARES point-source data is reported in Illuminati et al. (2019), where each track is considered separately.

The IceCube candidate neutrinos have individual angular errors that need to be accounted for when searching for point-source excesses. The likelihood method is then modified as follows: an independent fit is performed for each IceCube track in the sample, where the location in (R.A., δ) is let free to vary in the fit within a range of ±2σ separately for each coordinate, where σ is the angular error reported by IceCube (for the 8 yr up-going muon sample, the 50% containment value reported in the article is multiplied by a factor 1.5).

The test statistic is then defined as the sum over all the individual test statistics:

$$\begin{eqnarray}&&\lambda =\displaystyle \frac{1}{\displaystyle {\sum }_{j}{w}_{j}}\displaystyle \sum _{j=1}^{{N}_{{\rm{sources}}}}{w}_{j}{\lambda }_{j},\end{eqnarray} \tag{ 9 }$$

where λ_j is the local test statistic obtained for the jth IceCube track, carrying a weight ${w}_{j}={ \mathcal A }({\delta }_{j})$ equal to the ANTARES acceptance at the corresponding declination.

The global significance is evaluated via the total test statistic λ; nevertheless, the relative contribution of each IceCube track can be examined by looking at the distribution of individual test statistics λ_j .

To search for a possible dominant contribution of a small number of sources, an individual likelihood fit has been performed for each object in the radio galaxy and Fermi 3LAC catalogs. The resulting p-values associated with each source are shown in Figure 4. Apart from the bulk of the distribution, a single source is present in the tail of the distribution of p-values for both the radio galaxies and the Fermi blazars. Their properties are thus investigated in more details in the following sections.

Zoom In Zoom Out Reset image size

Radio galaxies. The radio galaxy with the smallest p-value is 3C 403, which is found with a pre-trial p-value of p = 2.3 ⨯ 10⁻⁴, equivalent to 3.7σ. The probability of finding by chance a similar excess among the N = 56 radio galaxies in the field of view of ANTARES has been computed using P(n ≥ 1) = 1 − (1 − p)^N , leading to a value P = 1.3 ⨯ 10⁻², equivalent to 2.5σ.

The arrival directions of ANTARES events around this source are shown in Figure 5(left), where one can see that the excess comes from the presence of two events lying at less than 0.5° from the source. Those track events have an estimated angular uncertainty of 0.2°, and reconstructed energies of ∼5 TeV and ∼10 TeV, respectively.

Zoom In Zoom Out Reset image size

The radio galaxy is located at a declination of +2.5° and its distance is evaluated to be ∼260 Mpc. The galaxy 3C 403 is of type FR-II and presents a peculiar radio morphology, called "X-shape" or "Winged" (see Figure 1 of Dennett-Thorpe et al. 2002) due to a recent change in the orientation of the jet axis, which could be the signature of a binary black hole merger or fusion with a small galaxy. The interaction of the newly formed jet with the radiation and matter environment around the black hole could lead to high-energy neutrino emission, as suggested in Kun et al. (2018).

Blazars. The most significant blazar from the Fermi 3LAC catalog is the BL Lac MG3 J225517+2409 (Fermi source 3FGL J2255.1+2411, or 4FGL J2255.2+2411 in the latest version of the catalog). A recent compilation of multiwavelength observations for this object is reported in Figure 8 of Giommi et al. (2020), showing a synchrotron peak around 10¹⁵ Hz, indicating that the blazar can be classified in the ISP category, similar to TXS 0506+056. The redshift of MG3 J225517+2409 is currently unknown, but a limit of z > 0.86 has been recently reported in Paiano et al. (2019).

The gamma-ray flux of MG3 J225517+2409 measured by Fermi-LAT in eight years of observation (Abdollahi et al. 2020) is well described by a simple power law without cutoff, with a spectral index γ = 2.10 ± 0.05 and an integral flux of Φ = 1.0 ⨯ 10⁻⁹ photons cm⁻² s⁻¹ for energies between 1 and 100 GeV.

It is noteworthy that another candidate blazar, 3FGL J2258.8+2437, is located only ∼1° away from MG3 J225517+2409, and ∼1.2° away from the position of the IceCube high-energy track event IC-100608A. However, this source does not have an identified counterpart and is not present in the 10 yr Fermi-LAT catalogs, suggesting that 3FGL J2258.8+2437 might have been a spurious detection. For the present study, this source is thus not considered in the calculation of chance coincidence.

The result of the individual likelihood fit for this source gives a pre-trial p-value of p = 1.4 ⨯ 10⁻⁴, equivalent to 3.8σ. The probability of finding by chance a similar excess among the N = 1255 blazars in the field of view of ANTARES is evaluated by generating pseudo-experiments, giving a post-trial value P = 0.15 (equivalent to 1.4σ). The simple binomial evaluation gives a very similar value, P(n ≥ 1) = 1 − (1 − p)^N ≃ 0.16, because the overlap between sources is very small: there are only 58 pairs of 3LAC objects that are separated by an angular distance less than 1°.

The distribution of ANTARES events around the blazar is shown in Figure 5(right). There are five ANTARES tracks lying at less than 1° from the source, with estimated energies ranging from ∼3 to ∼40 TeV and with estimated angular uncertainty between 0.2° and 0.5°. This particular region of the sky has been previously reported by the ANTARES Collaboration (Illuminati et al. 2019) to contain the hottest spot of the full-sky point-source blind search; the maximum excess lies at <0.7° from the position of the blazar MG3 J225517+2409.

Under the (restrictive) assumption that a neutrino flux is produced steadily by the source with an E⁻² spectrum, an estimation of the flux normalization at 1 GeV of MG3 J225517+2409 with the ANTARES data would be Φ₀ ∼ (0.6–3.0) ⨯ 10⁻⁸ GeV⁻¹ cm⁻² s⁻¹, depending on the number of neutrino candidates that are considered as signal (between 1 and 5) from this source.

In addition to the ANTARES measurement, a high-energy muon track IC-100608A from the IceCube 8 yr up-going muon sample (Haack & Wiebusch 2017) is found to lie at an angular distance of 1.1° from the source. This event has a deposited energy in the IceCube detector of about ∼300 TeV, but it has a poor angular resolution: the 50% CL containment region has a radius ∼1°, but the estimated 90% CL region has a radius larger than ∼3°. This poor angular reconstruction explains why the spatial correlation with MG3 J225517+2409 has not been reported before, the IceCube muon event IC-100608A being excluded from the selection in previous blazar correlation searches (Padovani et al. 2016; Giommi et al. 2020).

The possible association between ANTARES and IceCube neutrinos and the blazar MG3 J225517+2409 is investigated in more detail in this section, by further looking at the time evolution of the gamma-ray emission of the source. The latest data from the Fermi 4FGL catalog (Abdollahi et al. 2020) are used to compare the evolution of the gamma-ray flux with the arrival times of ANTARES and IceCube events.

The Fermi 4FGL catalog provides measurements with twice as much exposure as the Fermi 3FGL, with 20% larger acceptance, improved angular resolution above 3 GeV, a better model of the background galactic diffuse emission, and a refined likelihood analysis. The 4FGL catalog was not available at the time of the stacking analysis; it is used only in this section for the time analysis.

The time evolution of the gamma-ray flux measured by Fermi-LAT for the blazar 4FGL J2255.2+2411 above 100 MeV is shown in Figure 6, along with the arrival times of spatially coincident IceCube and ANTARES neutrino events. The data are presented in 48 bins, each two months wide, from 2008 August to 2016 August. A Bayesian block algorithm (Scargle et al. 2013) has been applied to the flux data in order to pinpoint possible flaring states. A significant increase of the flux by a factor of ∼4 with respect to the average value is evident in a time window of ∼4 months, centered at about MJD = 55,400.

Zoom In Zoom Out Reset image size

5.2.1. Antares Events

The arrival times of the five ANTARES events around the blazar MG3 J225517+2409 are clustered at earlier observation times, in a large time window encompassing the flare but with none of them actually occurring during the period of high gamma-ray activity (the two latest events were observed ∼2 and ∼3 months after the end of the gamma-ray flare). Therefore, including a time-dependent term in the likelihood analysis does not strengthen the potential association between ANTARES candidate neutrino events and the gamma-ray emission of the blazar.

Considering only the time information, a very simple and model-independent test is to compare the average time to the flare center (t_F )

$$\begin{eqnarray}&&\tau =\displaystyle \frac{1}{N}\,\displaystyle \sum _{i=1}^{N=5}| {t}_{i}-{t}_{F}| ,\end{eqnarray} \tag{ 11 }$$

observed for the five ANTARES events to the distribution that is expected from random pseudo-experiments (following the true time distribution of ANTARES events).

The value obtained for the data is τ ≃ 400 days, and the fraction of pseudo-experiments having a lower or equal value is found to be ≃2%, thus leading to an estimated p-value of p = 0.02 (2.3σ).

5.2.2. IceCube Event

To evaluate the degree of space and time correlation between the gamma-ray flux of MG3 J225517+2409 and the neutrino candidate of IceCube, an additional likelihood analysis has been performed for this specific event.

The likelihood of Equation (1) is modified as follows:

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal L }(x| {{\rm{H}}}_{1})=\mathrm{ln}\left[{\mu }_{s}S({x}_{i}){f}_{S}({t}_{i})+{\mu }_{b}B({x}_{i}){f}_{B}({t}_{i})\right]-{\mu }_{s}-{\mu }_{b},\end{eqnarray} \tag{ 12 }$$

where the new terms f_S (t_i ) and f_B (t_i ) represent the signal and background time PDFs respectively.

The signal f_S (t) is chosen to be proportional to the Bayesian block estimation (Scargle et al. 2013) of the gamma-ray flux of MG3 J225517+2409, shown as a solid red curve in Figure 6. The background f_B (t) is proportional to the average number of events detected per day by IceCube, using the point-source selection cuts (Aartsen et al. 2019). The value of f_B (t) is assumed to be constant during the different data-taking periods, from the IC59 sample (2009–2010) where f_B (t) ≃ 60 events/day, up to the IC86 sample where f_B (t) ≃ 210 events/day.

The spatial PDF S(x) is taken as a two-dimensional Gaussian function of width σ = 1.8°, corresponding to the 90% CL statistical angular error of ∼3° reported in Haack & Wiebusch (2017). The background distribution B(x) is assumed to be uniform in right ascension, and proportional to the value of a Bayesian block fit of the distribution in declination of the 56 IceCube high-energy tracks considered in the stacking analysis, as shown in Figure 7.

Zoom In Zoom Out Reset image size

The result of the likelihood fit gives a fitted number of signal events n_s = 0.68 and a test statistic equal to λ = 0.39. The associated p-value is computed by generating pseudo-experiments, leading to a value p = 5.3 ⨯ 10⁻², equivalent to 1.9σ. As the considered source has been specified a priori using the ANTARES result, there are no trial factors to account for.