A Search for Neutrino Point-source Populations in 7 yr of IceCube Data with Neutrino-count Statistics

IceCube is a cubic-kilometer Cerenkov detector, which is composed of 5160 digital optical modules (DOMs) embedded in the Antarctic ice at the South Pole Achterberg et al. (2006). These DOMs are attached to 86 strings of cable at depths between 1450 and 2450 m beneath the surface of the ice. Most of the DOMs have a vertical spacing of 17 m along the strings and the average distance between neighboring strings is ∼125 m. Each module consists of a photomultiplier tube, onboard digitization board, and a separate board with LEDs for calibration (Abbasi et al. 2009, 2010).

Construction of IceCube started in 2004 and was completed in 2010 December. Before the full detector was completed, data were being taken in partial configurations with fewer than 86 strings. In the present work, we make use of 7 yr of IceCube data. The first three years were taken during the 40-string (IC-40), 59-string (IC-59), and 79-string (IC-79) configurations, as described in Aartsen et al. (2013c) and Abbasi et al. (2011). The subsequent four years of data exploited the full 86-string (IC-86) configuration, as outlined in Aartsen et al. (2014b, 2017a).

Neutrinos are notoriously difficult to detect, and just because they point back to their origin does not mean we can necessarily extract that direction. As we cannot detect the neutrinos directly, the challenge is in inferring the direction from visible products left behind from a neutral or charged-current interaction the neutrino undergoes within or in the vicinity of IceCube. In order to enhance our sensitivity, we choose to focus on events where the flight direction of the neutrino can be accurately reconstructed, as explained below.

In detail, there are three different event topologies within IceCube produced by neutrino interactions: tracks, cascades, and double-bangs. Tracks result from muons traversing the detector, while cascades result from the charged-current interactions of electron or tau neutrinos or neutral-current interactions of any neutrino. The interactions in cascade events produce an almost spherical light emission making directional reconstruction difficult. Tracks from muons of ∼TeV or greater energies can travel several kilometers while constantly emitting Cerenkov light, making them ideal candidates for an accurate directional reconstruction. Double-bangs result from charged-current interactions of tau neutrinos at very high energies where the tau lepton decays to hadrons far enough from the initial interaction to create two distinct cascades. The first candidates for such events have now been identified (Stachurska 2018). Since tracks are the optimal topology for directional accuracy, in this paper, we will only use tracks.

The only neutrino interactions that can create tracks at TeV energies are charged-current interactions of muon neutrinos (and muon antineutrinos). Track events can be further divided into two subclasses: starting tracks occur when a muon neutrino has its charged-current interaction inside of the detector volume, while through-going tracks occur when that interaction occurs outside of the detector volume. Through-going tracks can result from any high-energy muon, including muons produced by cosmic-ray interactions in the atmosphere, but they also have a much higher effective area for muon neutrinos than starting tracks since the charged-current interaction can occur in a much larger volume than the detector volume. In this paper, we will consider only events reconstructed as through-going tracks, in order to take advantage of the high number of events resulting from the much larger effective area; however, the lower purity of astrophysical neutrinos in the sample creates a background we will need to account for in the NPTF analysis.

With these motivations, the specific data set used in this analysis was the through-going tracks in IceCube's 7 yr point-source sample (Aartsen et al. (2017a), and we refer to Section 2.2 of that reference for a detailed account of the data set. The 7 yr point-source sample has been accrued through several different event selections. Events from the year of IC-40 data were selected using fixed selection criteria on several parameters (Abbasi et al. 2011), while events from the remaining years (Aartsen et al. 2013c, 2014b, 2017a) were selected using multivariate boosted decision trees (BDTs) to classify events as signal or background. The BDTs were trained with separate background and signal data sets. In the background case, the BDT was trained on the data themselves, a procedure that is justified as the data are known to be strongly dominated by the background. For the signal, the BDT was instead trained on muon neutrino Monte Carlo simulations. These same Monte Carlo–simulated muon neutrinos are also used to calculate the effective area of the detector and the PSF. More information on the Monte Carlo simulation can be found in Aartsen et al. (2016b). Through these selection processes, the sample is divided into two regions, up-going (with decl. δ > −5° ) and down-going (δ < −5°). The down-going region is dominated by atmospheric muons, while the up-going region is shielded from these muons by the Earth. The up-going region is dominated by atmospheric neutrinos. Despite this distinction, in our analysis, we will not distinguish between up- and down-going events, instead performing a full sky analysis. As in Aartsen et al. (2017a), the total livetime is 2431 days, with 422,791 up-going events, and 289,078 down going. The directions of the events in this sample were reconstructed using a likelihood-based method, which uses information on how the photons scatter and get absorbed with the ice (Ahrens et al. 2004). The reconstructions of events in the IC-86 data used a more advanced description of the ice (Aartsen et al. 2013e) and a better parameterization of how the photons interact with it Aartsen et al. (2014b). The muon energy is estimated by approximating its energy loss along its reconstructed track Aartsen et al. (2014c). In order to fine tune this sample for our analysis, we make a cut on the reconstructed muon energy. The expected energy spectrum for an astrophysical neutrino flux is harder than the energy spectrum for the atmospheric neutrino flux. It should also be noted that the reconstructed muon energy is only an estimate of the muon energy at a point where it enters the detector and, thus, can only provide an estimated lower-bound of the primary neutrino energy, as much of the energy of the event can be deposited outside the detector. Nevertheless, a cut on the reconstructed muon energy can still effectively increase the sample's purity with respect to the harder-spectrum astrophysical flux. We determine the optimal energy cut to be 10^0.5 TeV ≈3 TeV. This value is determined by maximizing the sensitivity for the full sky analysis presented in this work. We note that as the energy distribution of the background in the northern and southern hemispheres is quite different, if we performed an analysis restricting to either of these, the optimal energy cut would vary. We also make a spatial cut around the poles of the detector ($| \delta | \gt 85^\circ$), as the scrambling procedure is less effective, and the reconstructions can be poorly behaved, in these regions. These cuts lead us to our final sample of 309,134 events.

Since the NPTF is a binned analysis, we must spatially bin the data. For this purpose, we use HEALPix Gorski et al. (2005) skymaps to bin the data into pixels of equal solid angle. There is still freedom in terms of how large to choose these bins, controlled by the nside parameter; however we find that a value of 64 maximizes our full-sky sensitivity, which corresponds to 49,152 bins, each of size approximately 0.84 square degrees.

In this section, we provide a brief, quantitative review of the NPTF likelihood framework, particularly emphasizing aspects that will be relevant to an application at IceCube. A more comprehensive description of the method and a derivation of all quoted expressions can be found in Mishra-Sharma et al. (2017).

Our ultimate goal is to write down a likelihood for a set of model parameters ${\boldsymbol{\theta }}$, given the data d described in Section 2—i.e., we want a function ${ \mathcal L }(d| {\boldsymbol{\theta }})$. Let us start with a description of the model parameters. In the case where we only have neutrinos originating from a single point-source population, then the model parameters specify the source-count distribution dN/dF. In principle, it is possible to keep the form of dN/dF very general, but a particularly simple analytic expression for the source-count function that is likely a good approximation to many realistic neutrino source classes is a broken power law⁶⁵

$$\begin{eqnarray}\displaystyle \frac{{{dN}}_{p}}{{dF}}(F;{\boldsymbol{\theta }})=A\,{T}_{p}\left\{\begin{array}{lc}{\left(F/{F}_{b}\right)}^{-{n}_{1}} & F\geqslant {F}_{b},\\ {\left(F/{F}_{b}\right)}^{-{n}_{2}} & F\lt {F}_{b}.\end{array}\right.\end{eqnarray} \tag{ 1 }$$

The use of a common functional form to describe the flux of both galactic and extragalactic sources is motivated by the fact that this ansatz appears to be a reasonable description of both populations in gamma-rays observed by Fermi (see, for example, Lee et al. 2016). In Equation (1), we have added the term T_p to the source-count function. T_p is a template, or pixelated spatial map describing the spatial distribution of the sources on the sky. It is the only term with an explicit pixel-by-pixel variation. As an example, for an isotropic or extragalactic distribution of sources, we have T_p = 1 for every value of p. This template could be considered as a model parameter and fitted for, but in our analysis, we will take T_p to be fixed before performing any likelihood analysis. With the template fixed, there are four model parameters for a singly broken power law: the normalization A, location of the break F_b, and power-law indices above and below the break n₁ and n₂. Thus, formally ${\boldsymbol{\theta }}=\{A,{F}_{b},{n}_{1},{n}_{2}\}$.

To provide some intuition for these parameters, note that we can calculate the total expected number of point sources and also the total expected flux from the population in each pixel from direct integration of the source-count function over all possible fluxes as follows

$$\begin{eqnarray}\begin{array}{rcl}{N}_{p}^{\mathrm{PS}} & = & {\int }_{0}^{\infty }{dF}\,\displaystyle \frac{{{dN}}_{p}}{{dF}}=\displaystyle \frac{A\,{T}_{p}\,{F}_{b}({n}_{1}-{n}_{2})}{({n}_{1}-1)(1-{n}_{2})},\\ {F}_{p}^{\mathrm{PS}} & = & {\int }_{0}^{\infty }{dF}\,F\displaystyle \frac{{{dN}}_{p}}{{dF}}=\displaystyle \frac{A\,{T}_{p}\,{F}_{b}^{2}({n}_{1}-{n}_{2})}{({n}_{1}-2)(2-{n}_{2})}.\end{array}\end{eqnarray} \tag{ 2 }$$

In performing these integrals, we require n₁ > 1 and n₂ < 1 to obtain a finite ${N}_{p}^{\mathrm{PS}}$, whereas for a finite ${F}_{p}^{\mathrm{PS}}$ we need ${n}_{1}\gt 2$ and n₂ < 2. There is an important distinction between the number of sources and flux as given in Equation (2). The total expected flux per pixel, ${F}_{p}^{\mathrm{PS}}$, is related to the expected number of neutrinos observed through the IceCube effective area, and in this sense is tied to an observable in the data. Yet the total expected number of sources, ${N}_{p}^{\mathrm{PS}}$, is not tied to an observable derivable from a map of neutrino counts. In this sense, a real map could have a best-fit value n₂∈ [1, 2], which corresponds to an infinite number of sources but a finite flux. While in practice, the total number of sources should be finite, the effective number of sources, for realistic source populations such as star-forming galaxies (Tamborra et al. 2014), can appear infinite. This effect occurs because the cutoff in the source-count distribution, dN/dS, that makes the number of sources finite appears at flux values well below the level where the sources contribute, on average, more than one photon or neutrino (Lisanti et al. 2016).

When presenting results, it will be helpful to use a different set of variables instead of $\{A,{F}_{b},{n}_{1},{n}_{2}\}$. Specifically, we will replace A and F_b with the expected number of point sources across the whole sky, N^PS, and the expected flux per source, ${\bar{F}}^{\mathrm{PS}}$. The change of variables can be implemented as follows:

$$\begin{eqnarray}\begin{array}{rcl}{N}^{\mathrm{PS}} & = & \displaystyle \sum _{p}{N}_{p}^{\mathrm{PS}}=\displaystyle \frac{{{AF}}_{b}({n}_{1}-{n}_{2})}{({n}_{1}-1)(1-{n}_{2})}\displaystyle \sum _{p}{T}_{p},\\ {\bar{F}}^{\mathrm{PS}} & = & \displaystyle \frac{\sum _{p}{F}_{p}^{\mathrm{PS}}}{\sum _{p}{N}_{p}^{\mathrm{PS}}}={F}_{b}\displaystyle \frac{({n}_{1}-1)(1-{n}_{2})}{({n}_{1}-2)(2-{n}_{2})}.\end{array}\end{eqnarray} \tag{ 3 }$$

Given a template and values of n₁ and n₂, we can then change variables to these more intuitive quantities, which we will exploit when presenting our results.

At this point, we comment on an important assumption that will be used throughout our analysis. We have repeatedly discussed the flux of sources, F, eschewing the question of what energy this flux is being measured at. To resolve this, we will assume that our point-source population follows the canonical astrophysical expectation of E⁻². In detail the flux from an individual source is given by

$$\begin{eqnarray}&&\displaystyle \frac{d{\rm{\Phi }}(E)}{{dE}}=F{\left(\displaystyle \frac{E}{1\mathrm{TeV}}\right)}^{-2},\end{eqnarray} \tag{ 4 }$$

which serves to contextualize the quantity F discussed throughout this section. It should be thought of as the normalization constant for the energy dependent flux, and we will assume an identical E⁻² scaling for all sources. If the actual spectrum is softer than this, the events will be shifted to lower energies where the backgrounds are higher, and generically, we would expect a reduction of the sensitivities shown here.

Returning to our derivation of the likelihood, the source-count function parameterizes our model prediction for the population of sources, but our goal is to embed this into a likelihood that can be fit to neutrino-count data. As a first step toward this, we need to address the fact that the discussion so far has been couched in the language of fluxes per neutrino source, Φ, commonly quoted with units of [neutrinos cm⁻² s⁻¹], whereas what is observed in the instrument is an integer number of neutrino events. The conversion between these two variables is provided by the effective area matrix, which accounts for the fact that two-point sources with equal flux at different locations on the sky will contribute a different number of detected neutrinos within IceCube. The effective area is the amalgamation of the detection efficiency for neutrinos incident on the IceCube detector from different directions, as well as an accounting for the fact the detector has a fixed location at the South Pole.

The conversion from flux, F, in units of [neutrinos cm⁻² s⁻¹] to counts, S, is achieved with the combination of the effective area and collection time, usually called an exposure map, which we denote by ${{ \mathcal E }}_{p}$—a pixel dependent quantity due to the spatial dependence in the effective area matrix. We defer the discussion of how the appropriate ${{ \mathcal E }}_{p}$ map for our analysis was derived until the following subsection. Assuming for now that the map is known, then using this, we can then convert to a source-count function in terms of counts rather than flux as follows:

$$\begin{eqnarray}&&\displaystyle \frac{{{dN}}_{p}}{{dS}}(S;{\boldsymbol{\theta }})=\displaystyle \frac{1}{{{ \mathcal E }}_{p}}\displaystyle \frac{{{dN}}_{p}}{{dF}}(F=S{/{ \mathcal E }}_{p};{\boldsymbol{\theta }}).\end{eqnarray} \tag{ 5 }$$

As ${{ \mathcal E }}_{p}$ can take on a different value in every single pixel, the conversion from flux to counts should in truth be performed in every pixel. Nevertheless, this is in practice often unnecessary. It is usually sufficient to divide the full map into a number of subregions, which each have comparable ${{ \mathcal E }}_{p}$ values, take the mean ${{ \mathcal E }}_{p}$ in this region, and perform the conversion once per region. Within the NPTFit framework, the number of subregions is controlled by the keyword nexp, and results are commonly convergent for small values of this parameter. As the appropriate ${{ \mathcal E }}_{p}$ for our data set varied significantly, albeit smoothly, as a function of decl., we chose to use 50 exposure regions in order to ensure the transformation from flux to counts was accurately performed.

From dN_p/dS, we can then move toward the NPTF likelihood by deriving the following useful quantity: the expected number of sources that will contribute m neutrinos in a pixel p, ${x}_{p,m}({\boldsymbol{\theta }})$. To do so, note that dN_p/dS evaluated at a particular S provides the expected number of sources that contribute an expected number of counts S, where of course S does not need to be an integer. The probability that one such source provides m neutrinos is then determined by the Poisson distribution, specifically S^me^−S/m!. From here, ${x}_{p,m}({\boldsymbol{\theta }})$ is given by weighting this factor by the source-count distribution and integrating over all S, as each value could Poisson fluctuate to m. In detail,

$$\begin{eqnarray}&&{x}_{p,m}({\boldsymbol{\theta }})={\int }_{0}^{\infty }{dS}\,\displaystyle \frac{{{dN}}_{p}}{{dS}}(S;{\boldsymbol{\theta }})\displaystyle \frac{{S}^{m}{e}^{-S}}{m!}.\end{eqnarray} \tag{ 6 }$$

This expression, while intuitive, has an inherent assumption that will be invalidated in most real applications: if a point source is located in pixel p, then it deposits all of its observed neutrinos in that same pixel. This neglects the finite PSF at IceCube, but we will hold off on a discussion of how to address this until the next subsection.

Our final goal of this subsection is to move from x_{p, m} to the probability of observing k neutrinos in a pixel p, ${p}_{p,k}$. Combining ${p}_{p,k}$ with the observed number of neutrinos in the data, d_p, and then taking the product over all pixels p, exactly gives us the likelihood through which we can constrain dN/dF, or more specifically ${\boldsymbol{\theta }}=\{A,{F}_{b},{n}_{1},{n}_{2}\}$. Being fully explicit, we have:

$$\begin{eqnarray}&&{ \mathcal L }\left(d| {\boldsymbol{\theta }}\right)=\prod _{p}{p}_{p,k={d}_{p}}({\boldsymbol{\theta }}),\end{eqnarray} \tag{ 7 }$$

where the product is taken over all pixels in the data set analyzed. In order to derive an expression for ${p}_{p,k}$, we use the concept of probability generating functions.⁶⁶ If we have a discrete probability distribution described by a set $\left\{{p}_{p,k}\right\}$ known for all k, then the generating function in a given pixel is defined as

$$\begin{eqnarray}&&{P}_{p}(t)\equiv \sum _{k=0}^{\infty }{p}_{p,k}{t}^{k},\end{eqnarray} \tag{ 8 }$$

where t is an auxiliary variable. The probabilities can be recovered from P_p(t) via

$$\begin{eqnarray}&&{p}_{p,k}={\left.\displaystyle \frac{1}{k!}\displaystyle \frac{{d}^{k}{P}_{p}(t)}{{{dt}}^{k}}\right|}_{t=0}.\end{eqnarray} \tag{ 9 }$$

In the case of models described by the Poisson distribution with expected counts μ_p, substituting the Poisson distribution into Equation (8) reveals the associated generating function to be ${P}_{p}(t)=\exp [{\mu }_{p}(t-1)]$.

Next, the aim is to construct the associated non-Poissonian generating function, starting with the expected number of m-neutrino sources, ${x}_{p,m}({\boldsymbol{\theta }})$, as given in Equation (6). Of course, m can take on any integer value, but for the moment, let us take it to be fixed, and we will determine the generating function for m-neutrino sources, denoted ${P}_{p}^{(m)}(t)$. From the definition in Equation (8), we need to know the probability of seeing k neutrinos in the pixel p, given by ${p}_{p,k}$, a value that will depend on how many m-neutrino sources there are. Specifically, k must be some integer n_m multiple of m, where n_m drawn from a Poisson distribution with mean ${x}_{p,m}({\boldsymbol{\theta }})$,

$$\begin{eqnarray}&&{p}_{p,{n}_{m}}=\displaystyle \frac{{x}_{p,m}^{{n}_{m}}{e}^{-{x}_{p,m}}}{{n}_{m}!},\end{eqnarray} \tag{ 10 }$$

where we have left the dependence on the parameters ${\boldsymbol{\theta }}$ implicit. In terms of this, the probability for seeing k neutrinos in pixel p is simply ${p}_{p,{n}_{m}}$ for the case that k = n_m × m, or zero otherwise as we are still keeping m fixed. Substituting this information into Equation (8), we obtain

$$\begin{eqnarray}\begin{array}{rcl}{P}_{p}^{(m)}(t) & = & \displaystyle \sum _{{n}_{m}}{t}^{m{n}_{m}}\,\displaystyle \frac{{x}_{p,m}^{{n}_{m}}{e}^{-{x}_{p,m}}}{{n}_{m}!}\\ & = & \exp \left[{x}_{p,m}({t}^{m}-1)\right],\end{array}\end{eqnarray} \tag{ 11 }$$

where we only included the nonzero values in the sum. Now, this result was obtained for a fixed m, but in order to obtain the full non-Poissonian generating function, we need to account for all possible m. As each source is independent, each value of m is independent also. We can then make use of the fact that the generating function of a sum of independent random variables is given by the product of each variable's generating function. Accordingly, the full non-Poissonian generating function is given by

$$\begin{eqnarray}&&{P}_{p}(t)=\exp \left[\sum _{m=1}^{\infty }{x}_{p,m}\left({t}^{m}-1\right)\right].\end{eqnarray} \tag{ 12 }$$

From here, using the inversion formula in Equation (9), the probability of observing k neutrinos in a pixel p is

$$\begin{eqnarray}&&{p}_{p,k}={\left.\displaystyle \frac{1}{k!}\displaystyle \frac{{d}^{k}}{{{dt}}^{k}}\exp \left[\sum _{m=1}^{\infty }{x}_{p,m}\left({t}^{m}-1\right)\right]\right|}_{t=0},\end{eqnarray} \tag{ 13 }$$

or combining this with our earlier results

$$\begin{eqnarray}\begin{array}{rcl}{p}_{p,k}({\boldsymbol{\theta }}) & = & \displaystyle \frac{1}{k!}\displaystyle \frac{{d}^{k}}{{{dt}}^{k}}\exp \left[\displaystyle \sum _{m=1}^{\infty }\left({t}^{m}-1\right)\right.\\ & & \times {\left.\left.{\int }_{0}^{\infty }{dS}\displaystyle \frac{1}{{{ \mathcal E }}_{p}}\displaystyle \frac{{{dN}}_{p}}{{dF}}\left(\displaystyle \frac{S}{{{ \mathcal E }}_{p}};{\boldsymbol{\theta }}\right)\displaystyle \frac{{S}^{m}{e}^{-S}}{m!}\right]\right|}_{t=0}.\end{array}\end{eqnarray} \tag{ 14 }$$

This expression, combined with Equation (7), gives us our full NPTF likelihood as a function of the source-count function d N_p/dF, which is exactly what we want to constrain. Although this expression contains an unevaluated integral, in the case of a multiply broken power law, it can be efficiently implemented as the integral can be calculated analytically, and furthermore, the p_p,k can be evaluated recursively in k. All of these details are described in Mishra-Sharma et al. (2017).

Below, we will outline how this likelihood can be extended to account for the presence of additional contributions to the data beyond just point sources. But before doing so, we turn to the practicalities of implementing the NPTF in the specific case of IceCube.

There are two immediate obstacles to implementing the NPTF at IceCube, beyond the basic outline described in the previous subsection. First, we need to calculate an appropriate effective area matrix in order to determine ${{ \mathcal E }}_{p}$, which appears in the conversion between the ${dN}/{dF}$ and ${dN}/{dS}$, as per Equation (5). Second, we must incorporate the real IceCube PSF and thereby remove the assumption hidden in Equation (6), that a source deposits all of its neutrinos in the pixel it is located. For some of the neutrinos recorded at IceCube, there is a small but nonzero probability that their true incident direction is separated from the reconstructed value by a significant (angular) distance, so this effect must be accounted for in essentially any binning of the data. We will address each point in turn.

Consider first the effective area matrix, a quantity that specifies the response of the detectors to an incident neutrino of a given energy, R.A., and decl. This, of course, cannot be calculated analytically. Instead we take simulations of individual events and use these to construct ${{ \mathcal E }}_{p}$. We do this by reweighting simulated events within our energy range according to an E⁻² spectrum, as this is the spectrum that we assume our point-source population follows, see Equation (4) and the surrounding discussion. This then provides the average detector response as a function of R.A. and decl., which we can map to galactic coordinates and provide ${{ \mathcal E }}_{p}$. As claimed above, this map can vary significantly, by more than an order of magnitude between locations with highest and lowest effective area. Given IceCube's location at the geographic South Pole, this variation is exclusively in decl., which tracks whether the event is arriving above from through the atmosphere, or below through the Earth. In light of this we used a relatively large number of 50 exposure regions to convert from flux to counts.

Turning next to the PSF, recall as discussed already that the NPTF likelihood derived in the last subsection assumed perfect angular reconstruction of every event. This assumption was invoked in writing down the number of sources contributing m neutrinos in a pixel p, denoted ${x}_{p,m}$, in Equation (6). By moving directly from the expected number of sources in the pixel to the expected number of neutrino counts, implicit is the assumption that the source deposits all of its flux into that pixel. Yet detector effects will smear the flux of a real source among a number of pixels. As in the NPTF, we do not keep track of which pixels are adjacent, what we want is the distribution for how a given source deposits its flux among the pixels on the map, a quantity denoted ρ(f). Here, $f\in [0,1]$ is the fraction of the point source's flux; the case of near perfect angular reconstruction, as compared to the pixel size, corresponds to a ρ(f) peaked near f = 1, as most of the flux tends to be distributed in one pixel (the pixel where the source is located). Indeed in the limit of exact angular reconstruction, we have ρ(f) = 2δ(f − 1).⁶⁷ More generically, however, imperfect angular reconstruction leads to a distribution peaked nearer f = 0 as most often a pixel will only get a small fraction of the flux. As a concrete example, consider a 3 × 3 grid with a point source at the center. Imagine the source deposits 60% of its flux in the central pixel that it inhabits, and then 5% in each of the eight pixels surrounding it. In this case, many more pixels experience a small amount of flux, and so ρ(f) would still be peaked toward smaller values of f. Further note that ρ(f) itself is not a probability density function, instead as the point source must deposit all of its flux somewhere, the distribution is normalized so that

$$\begin{eqnarray}&&{\int }_{0}^{1}{df}\,f\,\rho (f)=1.\end{eqnarray} \tag{ 15 }$$

Imagine we have the appropriate ρ(f)—we will describe how to derive this shortly—consider how this modifies x_p,m. Previously, we used the fact that dN_p/dS provides the number of sources that contribute an expected number of counts S, to then reweight this quantity by the probability of fluctuating from S to m, given by the Poisson distribution ${S}^{m}{e}^{-S}/m!$. Now, however, a source is only expected to deposit a fraction f of that flux in the pixel under consideration, and so instead, the reweighting to obtain m neutrinos is ${({fS})}^{m}{e}^{-{fS}}/m!$. Furthermore, the probability that a given value of f is chosen is dictated by the distribution ρ(f), and integrating over all possible flux fractions, we arrive at the following modification for x_p,m:

$$\begin{eqnarray}&&{x}_{p,m}({\boldsymbol{\theta }})={\int }_{0}^{\infty }{dS}\,\displaystyle \frac{{{dN}}_{p}}{{dS}}(S;{\boldsymbol{\theta }})\times {\int }_{0}^{1}{df}\,\rho (f)\displaystyle \frac{{\left({fS}\right)}^{m}{e}^{-{fS}}}{m!}.\end{eqnarray} \tag{ 16 }$$

Another way to understand the above modification is to compare this result to Equation (6). In doing so, the modification induced by the finite PSF is seen to be equivalent to substituting in a modified source-count function,

$$\begin{eqnarray}&&\displaystyle \frac{{{dN}}_{p}}{{dS}}(S)\to \displaystyle \frac{d{\tilde{N}}_{p}}{{dS}}(S)={\int }_{0}^{1}{df}\,\displaystyle \frac{\rho (f)}{f}\displaystyle \frac{{{dN}}_{p}}{{dS}}(S/f).\end{eqnarray} \tag{ 17 }$$

Propagating the modification in Equation (16) through to the NPTF likelihood then gives a full accounting for the effect of the finite PSF, once we have ρ(f). Note that taking ρ(f) = 2δ(f − 1), i.e., perfect angular reconstruction, the above expression reduces to Equation (6), as it must, and this can also be seen clearly in Equation (17).

All that remains then in order to incorporate the IceCube PSF is an algorithm for determining ρ(f). Commonly, the instrument PSF is stated as a probability distribution for a given event to be located at some radius from the center of a source. The median reconstruction angle for ${\nu }_{\mu }+{\bar{\nu }}_{\mu }$ events can be seen, for example, in Figure 2 of Aartsen et al. (2017a), where the angular error can vary from half a degree to several degrees depending on the energy and event type. Nevertheless, this is just the median reconstruction angle; the tails of this distribution are considerably non-Gaussian and can extend out to very wide angles in the case of poorly reconstructed events. Modeling the tails correctly is critical for an NPTF analysis. As an example of why this is important, if there is a true population of sources all with identical fluxes, a mismodeled PSF can lead to the fit preferring an additional population of lower flux sources associated with mis-reconstructed neutrinos.

The event-by-event determination of the angular reconstruction can be exploited in an unbinned analysis, as done for example in Aartsen et al. (2017a). However, as the NPTF is fundamentally a binned method, we will instead consider quantities averaged within the data set of interest. There is no simple analytic expression for converting from the known PSF to ρ(f), as, for example—unlike the PSF—ρ(f) depends critically on the binning of the underlying map.⁶⁸ We can, however, determine this distribution using the following algorithmic prescription:

• 1.  
  Simulate N equal-flux point sources on a blank map with the same pixelization as the NPTF will be applied to;
• 2.  
  Determine the fraction of the total flux in each pixel, f_p, defined such that ${\sum }_{p}{f}_{p}=1$;
• 3.  
  Define a flux binning Δf, and as a function of flux f, define Δn(f) the number of pixels that have a flux between f and f + Δf; and
• 4.  
  Combine these quantities to define ρ(f) as follows:

  $$\begin{eqnarray}&&\rho (f)\approx \displaystyle \frac{{\rm{\Delta }}n(f)}{N{\rm{\Delta }}f}.\end{eqnarray} \tag{ 18 }$$

The above relation is approximate and only becomes exact in the limit $N\to \infty$ and ${\rm{\Delta }}f\to 0$.

In order to simulate this in practice, we deposit a large number of sources on the sky, and for each one, we model the neutrino distribution according to the appropriate PSF at each location. To account for the energy dependence inherent in the PSF of IceCube, following Equation (4), we assume each source has an E⁻² spectrum and draw events for the source according to this distribution. In this way, we can exactly build up ρ(f) as defined in Equation (18), and the result is shown in Figure 1 for our default full-sky analysis. In that figure, we have zoomed in on the small f values where the distribution is peaked. That the flux fraction distribution is peaked at small values is indicative of the fact that the IceCube PSF has tails that extend significantly more broadly than the size of a pixel on the map but is also quite generic of ρ(f) unless the angular reconstruction is significantly better than the pixel size. For comparison, note that the linear size of our pixels is ∼092, which is comparable to the median reconstruction angle of our events, which is ∼1°. As the PSF of IceCube varies across the sky, ρ(f) must be determined for each of our spatial distributions separately.⁶⁹

Zoom In Zoom Out Reset image size

In light of these results, the NPTF likelihood can now be applied to the IceCube data set. As an important validation of the method, we performed extensive tests based upon Monte Carlo studies where we injected and then recovered point-source populations. An example of this is shown in Figure 6, discussed in detail below. We emphasize again that, in particular, the details discussed in this section were critical to the successful recovery of injected populations. With the method validated, we can now look to calculate the expected sensitivity at IceCube, which we turn to in Section 4. Before doing so, however, in the next two subsections, we will discuss how we incorporate backgrounds that are not associated with point sources into our model and then how these various likelihoods will be combined into an inference framework we can use to test for the presence of point sources.

At the very least, due to the presence of irreducible backgrounds, we know that point sources cannot be the only contribution to the IceCube data set, and in this section, we discuss how to augment the NPTF likelihood to account for these. These additional contributions are generally expected to be described by the Poisson distribution, following an underlying spatial map. To incorporate both the Poissonian statistics and the spatial variation, we will use the language of Poissonian templates, following the language from a recent application of this topic in Lisanti et al. (2018).

To begin with, as in the NPTF case, we imagine that our model follows a spatial distribution that once pixelized can be described by a map T_p. Unlike for an NPTF model, where T_p specified the spatial distribution of point sources, in the Poisson case, we require T_p to be proportional to the expected distribution of counts, not flux. As a concrete example, in the case where our model is for isotropic extragalactic neutrino emission, the appropriate T_p would still have a spatial variation inherited from the effective area matrix described in the previous subsection. As in the non-Poissonian case, we assume for a given model that the Poissonian template, T_p, is specified ahead of time. What we fit for in this case is the overall normalization of this template, in terms of which our Poissonian model prediction is given by

$$\begin{eqnarray}&&{\mu }_{p}({\boldsymbol{\theta }})=A\,{T}_{p},\end{eqnarray} \tag{ 19 }$$

so that ${\boldsymbol{\theta }}=\{A\}$, where we emphasize that A has no pixel dependence. Given that the sum of two Poisson distributions with means μ₁ and μ₂ is again a Poisson distribution of mean μ₁+μ₂, we can readily extend this formalism to account for multiple Poisson models. For example, if we had n of these, described by template ${T}_{p}^{1},\ldots ,{T}_{p}^{n}$, then our combined model prediction in each pixel would be

$$\begin{eqnarray}&&{\mu }_{p}({\boldsymbol{\theta }})=\sum _{{\ell }=1}^{n}{A}^{{\ell }}\,{T}_{p}^{{\ell }},\end{eqnarray} \tag{ 20 }$$

where now ${\boldsymbol{\theta }}=\{{A}^{1},\ldots ,{A}^{n}\}$. To provide a concrete example, we may want to model the observed flux using a model that combines three sources: 1. terrestrial backgrounds such atmospheric neutrinos; 2. diffuse extragalactic emission; and 3. diffuse emission from the Milky Way. In such a scenario, we would have three Poissonian templates, and for each of these, ${T}_{p}^{{\ell }}$ would describe the pixel dependence of the flux, and A^ℓ the overall normalization. For the case of extragalactic emission, as the flux is expected to be isotropic, ${T}_{p}^{{\ell }}$ would then be a map of the IceCube detector's response to a uniform incident flux.

In terms of this, if all we had was a set of Poissonian models, then we could write down our likelihood according to the Poisson distribution,

$$\begin{eqnarray}&&{p}_{p,k}({\boldsymbol{\theta }})=\displaystyle \frac{{\mu }_{p}{\left({\boldsymbol{\theta }}\right)}^{k}{e}^{-{\mu }_{p}({\boldsymbol{\theta }})}}{k!}.\end{eqnarray} \tag{ 21 }$$

Nonetheless, we want to construct a likelihood incorporating both Poissonian and non-Poissonian models. For this purpose, we can use the property of generating functions we exploited earlier, specifically that the generating function that describes the sum of two independent random variables is given by the product of the generating functions for the individual variables. For the Poisson case, given in Equation (21), the associated generating function according to Equation (8), is given by

$$\begin{eqnarray}&&{P}_{p}(t;{\boldsymbol{\theta }})=\exp \left[{\mu }_{p}({\boldsymbol{\theta }})(t-1)\right].\end{eqnarray} \tag{ 22 }$$

Combined with the generating function for the non-Poissonian case given in Equation (12), we arrive at

$$\begin{eqnarray}\begin{array}{rcl}{P}_{p}(t;{\boldsymbol{\theta }}) & = & \exp \left[\displaystyle \sum _{m=1}^{\infty }{x}_{p,m}({\boldsymbol{\theta }})\left({t}^{m}-1\right)\right.\\ & & \quad \left.+\,{\mu }_{p}({\boldsymbol{\theta }})(t-1)\right],\end{array}\end{eqnarray} \tag{ 23 }$$

where now ${\boldsymbol{\theta }}=\left\{A,{F}_{b},{n}_{1},{n}_{2},{A}^{1},\ldots ,{A}^{n}\right\}$. Through the use of Equation (9), this generating function can be used to derive a combined likelihood that includes both Poissonian and non-Poissonian models, accomplishing one of the main aims of this subsection.

The main application for the Poissonian template formalism in our work will be to model the known backgrounds arising from atmospheric neutrinos and muons. For this purpose, we need to derive an appropriate T_p describing the spatial distribution of these contributions. Determining this from first principles is at present out of reach. Fortunately, however, we can estimate the distribution from the data alone. The reason for this is if we assume the data are made up of predominantly background events and a subset of point sources, then we can remove the point sources in the following way. Given the approximate azimuthal symmetry of the effective area of detector, we can take the data collected by IceCube and scramble the events by assigning them a random R.A. value. This process removes any point-source hotspots, as they will be smeared out along bands of constant decl. Furthermore, as the background is only expected to vary with decl., this process does not degrade the spatial information pertaining to the background process. Applying this process once gives a map that is still as noisy as the data. In order to extract a more appropriate map for the mean of a Poisson distribution, we repeated this scrambling process a large number of times and take the average of the resulting maps. Finally, to remove the noise in decl., we convolve this model with a von MisesFisher distribution that has a concentration corresponding to 108, chosen as this is the median angular resolution at ∼1 TeV. This last step can be justified as the real data have been scrambled on such a scale due to the PSF.

The map resulting from this procedure is shown in Figure 2, which is the Mollweide projection of the map in galactic coordinates. The most apparent feature in this map is the strong variation away from the poles toward the equator where the largest flux is observed. Taking this map as template, we can readily incorporate the largest expected background into our likelihood. In subsequent discussions, we will refer to this map as ${T}_{p}^{\mathrm{bkg}}$.

Zoom In Zoom Out Reset image size

So far in this section, we have introduced the formalism required to calculate the likelihood for a data set in the presence of a population of sources and additional Poissonian contributions. Our goal now is to put this machinery to work in the form of a test statistic that we can use to test for the presence of point sources in the data. Our test statistic will compare between two hypotheses: that point sources, distributed according to a spatial template ${T}_{p}^{\mathrm{PS}}$, are present in the data or they are not. We will refer to these as the non-Poissonian and Poissonian hypotheses, respectively.

In detail, the Poissonian hypothesis is a model consistent of two Poisson templates: one following the dominant background contribution, given by ${T}_{p}^{\mathrm{bkg}}$, and the other accounting for the possibility that for a given spatial distribution ${T}_{p}^{\mathrm{PS}}$, the data may have a diffuse rather than an unresolved point-source origin. For this second source of emission, as the flux is diffuse, it is better described by Poisson statistics and, thus, follows a template ${{ \mathcal E }}_{p}{T}_{p}^{\mathrm{PS}}$, where the extra factor of ${{ \mathcal E }}_{p}$ is required to convert to a counts map. In this case, we can write down a likelihood function as described above, and from the data, we can construct the marginal likelihood as follows:

$$\begin{eqnarray}&&{{ \mathcal L }}_{0}(d)=\int d{\boldsymbol{\theta }}\,{{ \mathcal L }}_{{\rm{P}}}(d| {\boldsymbol{\theta }})p({\boldsymbol{\theta }}),\end{eqnarray} \tag{ 24 }$$

where the subscript 0 indicates that we are using this as our null hypothesis, and the subscript P on the likelihood identifies this as the appropriate form for the Poissonian hypothesis. In detail, ${{ \mathcal L }}_{{\rm{P}}}(d| {\boldsymbol{\theta }})$ can be determined directly from Equation (21), as we are only considering Poissonian models. In addition, we have introduced $p({\boldsymbol{\theta }})$, which represents the priors on the parameters. These are given in Table 1, where A^bkg and A^P are the normalizations of the templates ${T}_{p}^{\mathrm{bkg}}$ and ${{ \mathcal E }}_{p}{T}_{p}^{\mathrm{PS}}$, respectively, and so we see this is a two-parameter model.

The non-Poissonian hypothesis is derived from the null hypothesis, except that we append one further model: a non-Poissonian template following the spatial template ${T}_{p}^{\mathrm{PS}}$. As we will take a singly broken power law to describe the source-count distribution, this hypothesis is a six-parameter model. Once more, we can form the marginal likelihood using

$$\begin{eqnarray}&&{{ \mathcal L }}_{1}(d)=\int d{\boldsymbol{\theta }}\,{{ \mathcal L }}_{\mathrm{NP}}(d| {\boldsymbol{\theta }})p({\boldsymbol{\theta }}),\end{eqnarray} \tag{ 25 }$$

where again the priors are given in Table 1. Note that all priors are uniform, except for A and F_b, which are log uniform. Also, the prior for A^P is allowed to float negative, in order to allow the fit to conserve the total amount of flux when evaluating the non-Poissonian hypothesis. To justify this choice, recall that ${T}_{p}^{\mathrm{bkg}}$ was constructed by scrambling the real data. Accordingly, if there is a detectable point-source population within the data, the flux from these sources would be picked up by the non-Poissonian template but also be present in ${T}_{p}^{\mathrm{bkg}}$ and, therefore, double counted. A negative A^P can then be loosely thought of as subtracting that flux off, but done in such a way that the combined Poissonian template ${A}^{\mathrm{bkg}}{T}_{p}^{\mathrm{bkg}}+{A}^{{\rm{P}}}{{ \mathcal E }}_{p}{T}_{p}^{\mathrm{PS}}\geqslant 0$ in every pixel. We emphasize that the model used for non-Poissonian hypothesis includes both Poissonian and non-Poissonian templates, and thus, the full generating function in Equation (23) is required.

Model selection between these two hypotheses is considered through the use of the Bayes factor

$$\begin{eqnarray}&&{{ \mathcal B }}_{\mathrm{NP}/{\rm{P}}}(d)=\displaystyle \frac{{{ \mathcal L }}_{1}(d)}{{{ \mathcal L }}_{0}(d)}=\displaystyle \frac{\int d{\boldsymbol{\theta }}\,{{ \mathcal L }}_{\mathrm{NP}}(d| {\boldsymbol{\theta }})p({\boldsymbol{\theta }})}{\int d{\boldsymbol{\theta }}\,{{ \mathcal L }}_{{\rm{P}}}(d| {\boldsymbol{\theta }})p({\boldsymbol{\theta }})}.\end{eqnarray} \tag{ 26 }$$

From this definition, it can be seen that the Bayes factor is a summary statistic: it integrates over all possible forms of the source-count function through the integral over parameters. As such, it provides a gross evaluation as to whether a point-source population is preferred by the data, rather than singling out any particular dN/dF. As a summary statistic, it can also serve a secondary role as a test statistic. Our expectation for the value of the Bayes factor can be calibrated through use of frequentist methods such as calculating the p-value of the Bayes factor.

To compare more specific model hypotheses, we introduce the pointwise likelihood ratio, defined as

$$\begin{eqnarray}&&{ \mathcal M }\left(d;{\boldsymbol{\phi }}\right)=\displaystyle \frac{{\tilde{{ \mathcal L }}}_{\mathrm{NP}}\left(d| {\boldsymbol{\phi }}\right){\pi }_{\mathrm{NP}}}{{{ \mathcal L }}_{0}(d){\pi }_{{\rm{P}}}+{{ \mathcal L }}_{1}(d){\pi }_{\mathrm{NP}}}.\end{eqnarray} \tag{ 27 }$$

where π_P and π_NP are, respectively, the model priors for the Poissonian and non-Poissonian hypotheses, which have been chosen to be equal for this presentation of the results. Here, ${\boldsymbol{\phi }}=\left\{{N}^{\mathrm{PS}},{\bar{F}}^{\mathrm{PS}}\right\}$ represents the expected number of point sources across the whole sky and the expected flux per source at this location in model space, both of which were defined in Equation (3). In this expression, ${\tilde{{ \mathcal L }}}_{\mathrm{NP}}\left(d| {\boldsymbol{\phi }}\right)$ is similar to ${{ \mathcal L }}_{\mathrm{NP}}\left(d| {\boldsymbol{\theta }}\right)$, except that n₁ and n₂ have been marginalized over. Intuitively, ${ \mathcal M }\left(d;{\boldsymbol{\phi }}\right)$ should be thought of as the probability for the NPTF model at this particular value of $\left\{{N}^{\mathrm{PS}},{\bar{F}}^{\mathrm{PS}}\right\}$, compared to the probability for an equally weighted mixture of the Poissonian and non-Poissonian models. This definition was motivated by the need for a metric that handles both the high and low signal strength regimes. When ${{ \mathcal B }}_{\mathrm{NP}/{\rm{P}}}\ll 1$, ${ \mathcal M }$ is approximately the model posterior conditioned on ${\boldsymbol{\phi }}$; while in the ${{ \mathcal B }}_{\mathrm{NP}/{\rm{P}}}\gg 1$ regime, ${ \mathcal M }$ is approximately the ratio between the posterior and prior: $p({\boldsymbol{\phi }}| d)/p({\boldsymbol{\phi }})$.

Before concluding this section, we have said a number of times that we will exploit the power of the NPTF to test several forms for the spatial dependence of the point-source population, specified by ${T}_{p}^{\mathrm{PS}}$. In addition to isotropically distributed sources, where ${T}_{p}^{\mathrm{PS}}\propto 1$, we will consider three additional templates that consider the possibility of point sources distributed within the Milky Way, all of which are shown in Figure 3. The galactic disk template is used as a generic model for sources distributed according to the disk of the Milky Way. This is the line-of-sight integrated emission from a disk with a source density that scales exponentially in radius and distance from the plane, with a scale height of 0.3 kpc and a scale radius of 5 kpc. Next, we consider sources distributed following the Fermi bubbles (Su et al. 2010), large structures observed in gamma-rays extending perpendicularly from the galactic disk. Although the emission observed from the bubbles so far appears diffuse, neutrino emission from the recently discovered small scale gas clouds within the bubbles (Di Teodoro et al. 2018), could lead to point like sources. Finally, we consider the Schlegel, Finkbeiner, and Davis (SFD) dust map (Schlegel et al. 1998), which provides a two-dimensional distribution of dust within the Milky Way, mapped out using the reddening of starlight. This distribution is interesting to consider, as dust and gas tend to have similar spatial distributions, so the map is a proxy for the distribution of hydrogen in the Milky Way. The hydrogen is a target for cosmic-ray proton to collide with and form pions. The neutral pions then decay to photons, and indeed, the SFD dust map can be seen clearly in the Fermi gamma-ray data. If these same interactions produce higher-energy pions, then the charged variants could produce neutrinos at IceCube, making this an interesting spatial distribution to consider. As the target dust has a diffuse distribution throughout the Milky Way, in order for this template to describe a point-source population, the sources of the cosmic-ray protons would need to be point like, as this would then imprint a point-source-like distribution into the neutrino data. Note that we have chosen not to divide our templates between the northern and southern sky, even though this is commonly done for extragalactic point-source searches (see, for example, Aartsen et al. 2017a). Usually, this distinction between the northern and southern sky is imposed due to the different backgrounds that dominate in each hemisphere; in the northern sky, the main background is atmospheric neutrinos, whereas in the southern sky, instead, atmospheric muons dominate. Nonetheless, for the isotropic case, we show both the sensitivity and results for the northern sky in Appendix B. There we will see that the reach for the restricted case is only slightly enhanced, justifying our choice to focus on the full sky.

Zoom In Zoom Out Reset image size