TIME-INTEGRATED SEARCHES FOR POINT-LIKE SOURCES OF NEUTRINOS WITH THE 40-STRING IceCube DETECTOR

Neutrino astronomy is tightly connected to cosmic ray (CR) and gamma-ray astronomy, since neutrinos likely share their origins with these other messengers. With a possible exception at the highest observed energies, CRs propagate diffusively losing directional information due to magnetic fields, and both CRs and gamma rays at high energies are absorbed due to interactions on photon backgrounds. Neutrinos, on the other hand, are practically unabsorbed en route and travel directly from cosmological sources to the Earth. Neutrinos are therefore fundamental to understanding CR acceleration processes up to the highest energies, and the detection of astrophysical neutrino sources could unveil the origins of hadronic CR acceleration. Whether or not gamma-ray energy spectra above about 10 TeV can be accounted for by only inverse Compton processes is still an open question. Some observations suggest contributions from hadronic acceleration processes (Morlino et al. 2009; Boettcher et al. 2009). Acceleration of CRs is thought to take place in shocks in supernova remnants (SNRs) or in jets produced in the vicinity of accretion disks by processes which are not fully understood. Black holes in active galactic nuclei (AGNs), galactic micro-quasars and magnetars, or disruptive phenomena such as collapsing stars or binary mergers leading to gamma-ray bursts (GRBs), all characterized by relativistic outflows, could also be powerful accelerators. The canonical model for acceleration of CRs is the Fermi model (Fermi 1949), called first-order Fermi acceleration when applied to non-relativistic shock fronts. This model naturally gives a CR energy spectrum similar to dΦ/dE∝E⁻² at the source. The neutrinos, originating in CR interactions near the source, are expected to follow a similar energy spectrum. More recently, models such as those in Caprioli et al. (2010) can yield significantly harder source spectra. In the framework of these models, it is possible to account for galactic CR acceleration to energies up to the knee, at about Z × 4 × 10¹⁵ eV, where Z is the atomic number of the CR. Extragalactic sources, on the other hand, are believed to be responsible for ultra-high-energy CRs observed up to about 10²⁰ eV.

The concept of a neutrino telescope as a three-dimensional matrix of photomultiplier tubes (PMTs) was originally proposed by Markov & Zheleznykh (1961). These sensors detect the Cerenkov light induced by relativistic charged particles passing through a transparent and dark medium such as deep water or the Antarctic ice sheet. The depth of these detectors helps to filter out the large number of atmospheric muons, making it possible to detect the rarer neutrino events. The direction and energy of particles are reconstructed using the arrival time and number of the Cerenkov photons. High-energy muon–neutrino interactions produce muons that can travel many kilometers. On average, the muons scatter <01 with respect to the original neutrino direction for E_ν > 10 TeV. The first cubic-kilometer neutrino telescope, IceCube, is being completed at the South Pole. IceCube has a large target mass. This gives it excellent sensitivity to astrophysical neutrinos, enabling it to test many theoretical predictions.

Reviews on neutrino sources and telescopes can be found in Anchordoqui & Montaruli (2010), Chiarusi & Spurio (2010), Becker (2008), Lipari (2006), Bednarek et al. (2005), Halzen & Hooper (2002), Learned & Mannheim (2000), and Gaisser et al. (1995). Recent results on searches for neutrino sources have been published by IceCube in the 22-string configuration (Abbasi et al. 2009a, 2009b), AMANDA-II (Abbasi et al. 2009c), Super-Kamiokande (Thrane et al. 2009), and MACRO (Ambrosio et al. 2001).

This paper is structured as follows: Section 2 describes the detector. The data sample and cut parameters are discussed in Section 3, along with the simulation. In Section 4, the detector performance is characterized for searches. Section 5 describes the unbinned maximum likelihood search method, and in Section 6 the point-source and stacking searches are discussed. After discussing the systematic errors in Section 7, the results are presented in Section 8. Section 9 discusses the impact of our results on various possible neutrino emission models, and Section 10 offers some conclusions.

The IceCube Neutrino Observatory is composed of a deep array of 86 strings holding 5160 digital optical modules (DOMs) deployed between 1.45 and 2.45 km below the surface of the South Pole ice. The strings are typically separated by about 125 m with DOMs separated vertically by about 17 m along each string. IceCube construction started with the first string installed in the 2005–2006 austral summer (Achterberg et al. 2006a) and was completed in 2010 December. Six of the strings in the final detector use high quantum efficiency DOMs and a spacing of about 70 m horizontally and 7 m vertically. Two more strings have standard IceCube DOMs and 7 m vertical spacing but an even smaller horizontal spacing of 42 m. These eight strings along with seven neighboring standard strings make up DeepCore, designed to enhance the physics performance of IceCube below 1 TeV. The observatory also includes a surface array, IceTop, for extensive air shower measurements on the composition and spectrum of CRs.

Each DOM consists of a 25 cm diameter Hamamatsu PMT (Abbasi et al. 2010a), electronics for waveform digitization (Abbasi et al. 2009d), and a spherical, pressure-resistant glass housing. A single Cerenkov photon arriving at a DOM can produce a photoelectron, which is called a hit if the analog output of the PMT exceeds a threshold equivalent to 0.25 of the average single photoelectron (SPE) charge. The waveform of the PMT total charge is digitized and sent to the surface if hits are in coincidence with at least one other hit in the nearest or next-to-nearest neighboring DOMs within ±1000 ns. Hits that satisfy this condition are called local coincidence hits. The waveforms can contain multiple hits. The total number of photoelectrons and their arrival times are extracted with an iterative Bayesian-based unfolding algorithm. This algorithm uses the template shape representing an average hit.

Forty strings of IceCube were in operation from 2008 April 5 to 2009 May 20. The layout of these strings in relation to the final 86-string IceCube configuration is shown in Figure 1. Over the entire period the detector ran with an uptime of 92%, yielding 375.5 days of total exposure. Dead time is mainly due to test runs during and after the construction season dedicated to calibrating the additional strings and upgrading data acquisition systems.

Zoom In Zoom Out Reset image size

IceCube uses a simple multiplicity trigger, requiring local coincidence hits in eight DOMs within 5 μs. Once the trigger condition is met, local coincidence hits within a readout window ±10 μs are recorded, and overlapping readout windows are merged together. IceCube triggers primarily on down-going muons at a rate of about 950 Hz in this (40-string) configuration. Variations in the trigger rate by about ±10% are due to seasonal changes affecting development of CR showers and muon production in the atmosphere, with higher rates during the austral summer (Tilav et al. 2010).

3.1. Data Sample

Traditional astrophysical neutrino point-source searches have used the Earth to block all upward traveling (up-going) particles except muons induced by neutrinos, as in Abbasi et al. (2009b). There remains a background of up-going muons from neutrinos, which are created in CR air showers and can penetrate the entire Earth. These atmospheric neutrinos have a softer energy spectrum than many expectations for astrophysical neutrinos. The measurement of the atmospheric neutrino spectrum for the 40-string detector is discussed in Abbasi et al. (2010b). A large number of muons produced in CR showers in the atmosphere and moving downward through the detector (down-going) are initially misreconstructed as up-going. These mask the neutrino events until quality selections are made, leaving only a small residual of misreconstructed events.

The down-going region is dominated by atmospheric muons that also have a softer spectrum compared to muons induced by astrophysical neutrinos. At present, this large background reduces the IceCube sensitivity to neutrino sources in the southern sky in the sub-PeV energy region. While veto techniques are in development which will enable larger detector configurations to isolate neutrino-induced events starting within the detector, point-source searches can meanwhile be extended to the down-going region if the softer-spectrum atmospheric muon background is reduced by an energy selection. This was done for the first time using the previous 22-string configuration of IceCube (Abbasi et al. 2009a), extending IceCube's field of view to −50° declination. In this paper, we extend the field of view to −85° declination (the exclusion between −85° and −90° is due to the use of scrambled data for background estimation in the analysis, described in Section 5). Downgoing muons can also be created in showers caused by gamma rays, which point back to their source like neutrinos. The possibility for IceCube to detect PeV gamma-ray sources in the southern sky is discussed in Halzen et al. (2009), which concludes that a realistic source could be detected using muons in the ice only after 10 years of observing. Gamma-ray sources will not be considered further here.

Two processing levels are used to reduce the approximately 3.3 × 10¹⁰ triggered events down to a suitable sample for analysis (see Table 1). Random noise at the level of about 500 Hz per DOM is mainly due to radioactive decays in the materials in the DOMs. The contribution to triggered events by this random noise is highly suppressed by the local coincident hit requirement. To further reduce the contribution from noise, only hits within a 6 μs time window are used for the reconstructions. This time window is defined as the window that contains the most hits during the event. About 5% of down-going muons which trigger the detector are initially misreconstructed as up-going by the first stages of event processing. A persistent background that grows with the size of the detector is CR muons (or bundles of muons) from different showers which arrive in coincidence. At trigger level, they make up about 13% of the events. These coincident muon bundles can mimic the hit pattern of good up-going events, confusing a single-muon fit.

Table 1. Number of Events at Each Processing Level for the 375.5 days of Livetime

Triggered events	3.3 × 1010
L1 filtered events	8.0 × 108
Events in final sample	36,900

Download table as:  ASCIITypeset image

A likelihood-based muon track reconstruction is first performed at the South Pole (L1 filter). The likelihood function (Ahrens et al. 2004) parameterizes the probability of observing the geometry and timing of the hits in terms of a muon track's position, zenith angle, and azimuth angle. This likelihood is maximized, yielding the best-fit direction and position for the muon track. Initial fits are performed using an SPE likelihood that uses the time of the leading edge of the first photon arriving in each DOM. These reconstructions yield robust results used for the first level of background rejection. All events that are reconstructed as up-going are kept, while events in the down-going region must pass an energy cut that tightens with decreasing zenith angle. Events pass this L1 filter at an average rate of about 22 Hz and are buffered before transmission via a communications satellite using the South Pole Archival and Data Exchange (SPADE) system.

The processing done in the North includes a broader base of reconstructions compared to what is done at the South Pole. Rather than just the simple SPE fit, the multiple photoelectron (MPE) fit uses the number of observed photons to describe the expected arrival time of the first photon. This first photon is scattered less than an average photon when many arrive at the same DOM. The MPE likelihood description uses more available information than SPE and improves the tracking resolution as energy increases, and this reconstruction is used for the final analysis. The offline processing also provides parameters useful for background rejection, reconstructs the muon energy, and estimates the angular resolution on an event-by-event basis. Reducing the filtered events to the final sample of this analysis requires cutting on the following parameters:

• 1.  
  Reduced log-likelihood. The log-likelihood from the muon track fit divided by the number of degrees of freedom, given by the number of DOMs with hits minus five, the number of free parameters used to describe the muon. This parameter performed poorly on low-energy signal events. It was found that low-energy efficiency could be increased by instead dividing the likelihood by number of DOMs with photon hits minus 2.5. Both the standard and modified parameters were used, requiring events to pass one selection or the other. This kept the efficiency higher for a broader energy range.
• 2.  
  Angular uncertainty, σ. An estimate of the uncertainty in the muon track direction. The directional likelihood space around the best track solution is sampled and fit to a paraboloid. The contour of the paraboloid traces an error ellipse indicating how well the muon direction is localized (Neunhoffer 2006). The rms of the major and minor axes of the error ellipse is used to define a circular error. This parameter is effective both for removing misreconstructed events and as an event-by-event angular uncertainty estimator.
• 3.  
  Muon energy proxy. The average photon density along the muon track, used as a proxy for the muon energy. It is calculated accounting for the distance to DOMs, their angular acceptance, and average scattering and absorption properties of photons in the ice. The energy loss of a muon moving through the detector scales with the muon energy above about 1 TeV when stochastic energy losses due to bremsstrahlung, pair production, and photonuclear interactions dominate over ionization losses. The energy resolution obtained is of the order of 0.3 in the log₁₀ of the muon energy (at closest approach to the average hit location) for energies between about 10 TeV and 100 PeV. Since the interaction vertex is often an unknown distance from the detector, the muon in the detection volume has already lost an unknown amount of energy. Figure 2 shows the distribution of this energy parameter versus the true neutrino energy for a simulated spectrum dΦ/dE∝E⁻². Despite the uncertainty on the neutrino energy, for a statistical sample of events this energy estimator is a powerful analysis tool because of the wide range over which energies are measured.
• 4.  
  Zenith-weighted likelihood ratio. The likelihood ratio between an unbiased muon fit and a fit with an event weight according to the known down-going muon zenith distribution as a Bayesian prior. Applied to up-going tracks, a high likelihood ratio establishes strong evidence that the event is actually up-going and not a misreconstructed down-going event.
• 5.  
  NDir. The number of DOMs with direct photons, defined as arriving within −15 ns to +75 ns of the expectation from an unscattered photon emitted from the reconstructed muon track at the Cerenkov angle. Scattering of photons in the ice causes a loss of directional information and will delay them with respect to the unscattered expectation.
• 6.  
  LDir. The maximum length between direct photons, projected along the best muon track solution.
• 7.  
  Zenith directions of split events. The zenith angles resulting from splitting of an event into two parts and reconstructing each part separately. This is done in two ways: temporally, by using the mean photon arrival time as the split criterion, and geometrically, by using the plane both perpendicular to the track and containing the average hit location as the split criterion. This technique is effective against coincident muon bundles misreconstructed as single up-going tracks if both sub-events are required to be up-going.

Zoom In Zoom Out Reset image size

In the up-going region, all parameters are used. The zenith-weighted likelihood ratio and event splitting are specifically designed to remove down-going atmospheric muon backgrounds that have been misreconstructed as up-going while the other parameters focus on overall track quality.

In the down-going region, without a veto or Earth filter, muons from CR showers overwhelm neutrino-induced muons, except possibly at high energies if the neutrino source spectra are harder than the CR spectrum. The aim of the analysis in this region is therefore to select high-energy, well-reconstructed events. We use the first three parameters in the list above as cut variables, requiring a higher track quality than in the up-going range. Energy cuts were introduced in the down-going region to reduce the number of events to a suitable size, cutting to achieve a constant number of events per solid angle (which also simplifies the background estimation in the analysis). This technique keeps the high-energy events which are most important for discovery. Cuts were optimized for the best sensitivity using a simulated signal of muon neutrinos with spectrum dΦ/dE∝E⁻². We checked that the same cuts resulted in a nearly optimal sensitivity for a softer spectrum dΦ/dE∝E⁻³ in the up-going region where low-energy sensitivity is possible and for a harder spectrum dΦ/dE∝E^−1.5 in the down-going region.

Of the 36,900 events passing all selection criteria, 14,121 are up-going events from the northern sky, mostly muons induced by atmospheric neutrinos. Simulations of CR air showers with CORSIKA (Heck et al. 1998) show a 2.4% ± 0.8% contamination due to misreconstructed down-going atmospheric muons. The other 22,779 are down-going events from the southern sky, mostly high-energy atmospheric muons. An equatorial sky map of these events is given in Figure 3.

Zoom In Zoom Out Reset image size

3.2. Data and Simulation Comparison

Simulation of neutrinos is used for determining event selection and calculating upper limits. The simulation of neutrinos is based on ANIS (Gazizov & Kowalski 2005). Deep inelastic neutrino–nucleon cross sections use CTEQ5 parton distribution functions (Lai et al. 2000). Neutrino simulation can be weighted for different fluxes, accounting for the probability of each event to occur. In this way, the same simulation sample can be used to represent atmospheric neutrino models such as Bartol (Barr et al. 2004) and Honda (Honda et al. 2007) neutrino fluxes from pion and kaon decays (conventional flux). Neutrinos from charmed meson decays (prompt flux) have been simulated according to a variety of models (Martin et al. 2003; Enberg et al. 2008; Bugaev et al. 1989). Seasonal variations in atmospheric neutrino rates are expected to be a maximum of ±4% for neutrinos originating near the polar regions. Near the equator, atmospheric variations are much smaller and the variation in the number of events is expected to be less than ±0.5% (Ackermann & Bernardini 2005).

Atmospheric muon background is simulated mostly to guide and verify the event selection. Muons from CR air showers were simulated with CORSIKA (Heck et al. 1998) with the SIBYLL hadronic interaction model (Ahn et al. 2009). An October polar atmosphere, an average case over the year, is used for the CORSIKA simulation, ignoring the seasonal variations of ±10% in event rates (Tilav et al. 2010). Muon propagation through the Earth and ice is done using Muon Monte Carlo (Chirkin & Rhode 2004). Using measurements of the scattering and absorption lengths in ice (Ackermann et al. 2006), a detailed simulation propagates the photon signal to each DOM (Lundberg et al. 2007). The simulation of the DOMs includes their angular acceptance and electronics. Experimental and simulated data are processed and filtered in the same way.

In Figure 4, we show the cosine of zenith and in Figure 5 the muon energy proxy, reduced log-likelihood, and angular uncertainty estimator distributions of all events at trigger level, L1 filter level, and after final analysis cuts for data and Monte Carlo (MC). In these figures, the simulation uses a slightly modified version of the poly-gonato model of the galactic CR flux and composition (Hoerandel 2003). Above the galactic model cutoff at Z × 4 × 10¹⁵ eV, a flux of pure iron is used with an E⁻³ spectrum. This is done because currently CORSIKA cannot propagate elements in the poly-gonato model that are heavier than iron. Moreover, the poly-gonato model only accounts for galactic CRs and does not fully account for the average measured flux above 10¹⁷ eV, even when all nuclei are considered (see Figure 11 in Hoerandel 2003). These corrections then reproduce the measured CR spectrum at these energies. There is a 23% difference in normalization of data and CR muon events at trigger level. This normalization offset largely disappears after quality cuts are made. Generally good agreement is achieved at later cut levels.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 4 shows some disagreement between data and simulation for zenith angles around 80°. Muons created in CR showers in the atmosphere near this zenith angle must travel about 15 km to reach the bottom of IceCube. Only very high energy muons can travel such distances. For the simulation to produce the correct zenith distribution for these nearly horizontal events, CR composition can be important since protons can produce higher energy muons than iron nuclei with the same energy.

In addition to the slightly modified version of the poly-gonato model, discussed above, a simpler pure proton and iron two-component model with a much higher contribution of protons is used for comparison (Glasstetter & Hoerandel 1999). The final zenith distribution with each of these models is shown in Figure 6. The atmospheric muon simulation is not only affected by the primary composition uncertainties at high energy; it is also affected by the hadronic model, affecting the production rate of muons at the level of 15% in the region of interest for IceCube, greater than about 1 TeV, as discussed in the SIBYLL model paper (Ahn et al. 2009) and in the comparison between different hadronic models used in CORSIKA presented in Berghaus et al. (2008).

Zoom In Zoom Out Reset image size

For the up-going region, several models of atmospheric neutrino fluxes, both conventional fluxes from pion and kaon decay and prompt fluxes from charmed meson decay, are shown in Figure 6. To represent the low and high predictions, conventional and prompt models are used in pairs: Honda (Honda et al. 2007) for the conventional flux paired with Sarcevic (Enberg et al. 2008) for the prompt flux represent the low prediction, and Bartol (Barr et al. 2004) for the conventional flux paired with Naumov (Bugaev et al. 1989) for the prompt flux represent the high prediction. Additional uncertainty in the predicted atmospheric neutrino rate is estimated to be about 40% at 1 TeV (Barr et al. 2006). We conclude that our data agree with background simulation at the final level, within the range of uncertainties allowable by existing CR composition and atmospheric neutrino models.

The performance of the detector and the analysis is characterized using the simulation described in Section 3.2. For a spectrum of neutrinos dΦ/dE∝E⁻², the median angular difference between the neutrino and the reconstructed direction of the muon in the northern (southern) sky is 08 (06). Along with more severe quality selection in the southern sky, the different energy distributions in each hemisphere, shown in Figure 7, cause the difference in these two values. This is because the reconstruction performs better at higher energy due to the larger amount of light and longer muon tracks. The cumulative point-spread function (PSF) is shown in Figure 8 for two energy ranges and compared with simulation of the complete IceCube detector using the same quality selection, as well as the median PSF versus energy for the two hemispheres.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The neutrino effective area A^eff_ν is a useful parameter to determine event rates and the performance of a detector for different analyses and fluxes. The expected event rate for a given differential flux dΦ/dE is

$$\begin{equation} N_{{\rm events}}(\delta _\nu) = \int {\it dE}_{\nu } A^{\rm eff}_{\nu }(E_{\nu }, \delta _\nu) \frac{d\Phi _\nu (E_{\nu }, \delta _\nu)}{{\it dE}_{\nu }} \,, \end{equation} \tag{ 1 }$$

and is calculable using simulation. The A^eff_ν represents the size of an equivalent detector if it were 100% efficient. Figure 9 shows the A^eff_ν for fluxes of $\nu _{\mu }+\bar{\nu }_{\mu }$ and $\nu _{\tau }+\bar{\nu }_{\tau }$, for events at final selection level. Neutrinos arriving from the highest declinations must travel through the largest column depth and can be absorbed: this accounts for the turnover at high energies for nearly vertical up-going muon neutrinos. Tau neutrinos have the advantage that the tau secondary can decay back into a tau neutrino before losing much energy.

Zoom In Zoom Out Reset image size

Although tau (and electron) neutrino secondaries usually produce nearly spherical showers rather than tracks, tau leptons will decay to muons with a 17.7% branching ratio (Amsler et al. 2008). At very high energy (above about 1 PeV), a tau will travel far enough before decaying that the direction can be reconstructed well, contributing to any extraterrestrial signal in the muon channel. For the upper limits quoted in Section 8, we must make an assumption on the flavor ratios at Earth, after oscillations. It is common to assume $\Phi _{\nu _e}\hbox{{:}}\Phi _{\nu _{\mu }}\hbox{:}\Phi _{\nu _{\tau }} = 1\hbox{:}1\hbox{:}1$. This is physically motivated by neutrino production from pion decay and the subsequent muon decay, yielding $\Phi _{\nu _e}\hbox{:}\Phi _{\nu _{\mu }}\hbox{:}\Phi _{\nu _{\tau }} = 1\hbox{:}2\hbox{:}0$. After standard oscillations over astrophysical baselines, this gives an equal flux of each flavor at Earth (Athar et al. 2000). Under certain astrophysical scenarios, the contribution from muon decay may be suppressed, leading to an observed flux ratio of $\Phi _{\nu _e}\hbox{:}\Phi _{\nu _{\mu }}\hbox{:}\Phi _{\nu _{\tau }} = 1\hbox{:}1.8\hbox{:}1.8$ (Kashti & Waxman 2005), or the contribution of tau neutrinos could be enhanced by the decay of charmed mesons at very high energy (Enberg et al. 2009). For a spectrum dΦ/dE∝E⁻² and equal muon and tau neutrino fluxes, the fraction of tau neutrino-induced events is about 17% for vertically down-going, 10% for horizontal, and 13% for vertically up-going. Because the contribution from tau neutrinos is relatively small, assuming only a flux of muon neutrinos can be used for convenience and to compare to other published limits. We have tabulated limits on both $\Phi _{\nu _{\mu }}$ and the sum $\Phi _{\nu _{\mu }}+\Phi _{\nu _{\tau }}$, assuming an equal flux of each, while in the figures we have specified that we only consider a flux of muon neutrinos. Limits are always reported for the flux at the surface of the Earth.

An unbinned maximum likelihood ratio method is used to look for a localized, statistically significant excess of events above the background. We also use energy information to help separate possible signal from the known backgrounds.

The method follows that of Braun et al. (2008). The data are modeled as a two-component mixture of signal and background. A maximum likelihood fit to the data is used to determine the relative contribution of each component. Given N events in the data set, the probability density of the ith event is

$$\begin{equation} \frac{n_s}{N}\mathcal {S}_i + \left(1 - \frac{n_s}{N}\right)\mathcal {B}_i, \end{equation} \tag{ 2 }$$

where $\mathcal {S}_i$ and $\mathcal {B}_i$ are the signal and background probability density functions (PDFs), respectively. The parameter n_s is the unknown contribution of signal events.

For an event with reconstructed direction $\vec{x}_i=(\alpha _i,\delta _i)$, where α_i is the right ascension (R.A.) and δ_i is the declination, we model the probability of originating from the source at $\vec{x}_s$ as a circular two-dimensional Gaussian,

$$\begin{equation} \mathcal {N}(\vec{x}_i|\vec{x}_s,\sigma _i) = \frac{1}{2\pi \sigma _i^2}{\rm exp}\left({-{\frac{|\vec{x}_i - \vec{x}_s|}{2\sigma _i^2}}}\right), \end{equation} \tag{ 3 }$$

where σ_i is the angular uncertainty reconstructed for each event individually (Neunhoffer 2006) and $|\vec{x}_i - \vec{x}_s|$ is the space angle difference between the source and reconstructed event. While the average angular uncertainty decreases with increasing energy, the individual σ_i values are estimated from the reconstruction likelihood shape itself, and therefore the PSF incorporates this dependence without explicitly being a function of energy. The PSFs for different ranges of σ_i are in Figure 10, showing the correlation between the estimated angular uncertainty and actual track reconstruction error.

Zoom In Zoom Out Reset image size

The energy PDF $\mathcal {E}(E_i|\gamma,\delta _i)$ describes the probability of obtaining a reconstructed muon energy E_i for an event produced by a source of a given neutrino energy spectrum E^−γ at declination δ_i. We describe the energy distribution using 22 declination bands. Twenty bands, spaced evenly by solid angle, cover the down-going range where the energy distributions are changing the most due to the energy cuts in the event selection, while two are needed to sufficiently describe the up-going events, with the separation at δ = 15°. We fit the source spectrum with a power law E^−γ; γ is a free parameter. The probability of obtaining a reconstructed muon energy E_i for an event produced by a source with spectral index γ, for spectral indices 1.0 < γ < 4.0, is determined using simulation. Two examples of these energy PDFs are shown in Figure 11.

Zoom In Zoom Out Reset image size

The full signal PDF is given by the product of the spatial and energy PDFs:

$$\begin{equation} \mathcal {S}_i = \mathcal {N}(\vec{x}_i|\vec{x}_s,\sigma _i) \cdot \mathcal {E}(E_i|\gamma,\delta _i). \end{equation} \tag{ 4 }$$

The background PDF $\mathcal {B}_i$ contains the same terms, describing the angular and energy distributions of background events:

$$\begin{equation} \mathcal {B}_i = \mathcal {N}_{\rm Atm}(\vec{x}_i) \cdot \mathcal {E}(E_i|\rm Atm,\delta _i), \end{equation} \tag{ 5 }$$

where $\mathcal {N}_{\rm Atm}(\vec{x}_i)$ is the spatial PDF of atmospheric background and $\mathcal {E}(E_i|{\rm Atm},\delta _i)$ is the probability of obtaining E_i from atmospheric backgrounds (neutrinos and muons) at the declination of the event. These PDFs are constructed using data and, for the energy term, in the same 22 declinations bands as the signal PDF. All non-uniformities in atmospheric background event rates caused by the detector acceptance or seasonal variation average out in the time-integrated analysis. Therefore, $\mathcal {N}_{\rm Atm}(\vec{x}_i)$ has a flat expectation in R.A. and is only dependent on declination. Because the data are used in this way for background estimation, the analysis is restricted from −85° to 85° declination, so that any point-source signal will still be a small contribution to the total number of events in the same declination region.

The likelihood of the data is the product of all event probability densities:

$$\begin{equation} \mathcal {L}(n_s, \gamma) = \prod _{i=1}^{N} \Big [\frac{n_s}{N}\mathcal {S}_i + \left(1 - \frac{n_s}{N}\right)\mathcal {B}_i\Big ]. \end{equation} \tag{ 6 }$$

The likelihood is then maximized with respect to n_s and γ, giving the best-fit values $\hat{n}_s$ and $\hat{\gamma }$. The null hypothesis is given by n_s = 0 (γ has no meaning when no signal is present). The fit has been restricted to the physical signal region n_s ⩾ 0. The likelihood ratio test statistic is

$$\begin{equation} {\rm TS} = \left\lbrace \begin{array}{@{}lr}-2 \log \dsty\frac{\mathcal {L}(n_s=0)}{\mathcal {L}(\hat{n}_s,\hat{\gamma }_s)} & n_s \ge 0, \\ 0 & n_s < 0. \end{array} \right. \end{equation} \tag{ 7 }$$

The significance of the analysis is determined by comparing the TS from the real data with the distribution of TS from the null hypothesis (events scrambled in R.A.). We define the p-value as the fraction of randomized data sets with equal or higher test statistic values than the real data. Since we do not allow negative values of n_s, all underfluctuations result in TS = 0, the lowest possible value. This yields a p-value of 100%, which happens in approximately half of the searches. We evaluate the median sensitivity and upper limits at a 90% confidence level (CL) using the method of Feldman & Cousins (1998) and calculate the discovery potential as the flux required for 50% of trials with simulated signal to yield a p-value less than 2.87 × 10⁻⁷ (i.e., 5σ significance if expressed as the one-sided tail of a Gaussian distribution). The distributions of TS and the corresponding p-value for 10 million trials are shown in Figure 12 for a fixed point source at δ = 25°. Distributions with simulated signal events injected following a spectrum dΦ/dE∝E⁻² are included, as well as a χ² distribution with 2 degrees of freedom, which is used to estimate the 5σ significance threshold for calculating the discovery potential since simulating enough scrambled data sets requires a large amount of processing time.

Zoom In Zoom Out Reset image size

Although sensitivities and limits for sources with dΦ/dE∝E⁻² have become a useful benchmark for comparing performance, a wide range of other spectral indices are possible along with cutoffs over a wide range of energy. To understand the ability of the method to detect sources with cutoff spectra, typically observed in gamma rays to be in the range 1–10 TeV for galactic sources, Figure 13 shows the discovery potentials for a wide range of exponential cutoffs, demonstrating the ability of the method to detect sources with cutoff spectra. Typically, cutoffs observed in gamma rays are in the range 1–10 TeV for galactic sources. The likelihood fit is still performed using a pure power law.

Zoom In Zoom Out Reset image size

The likelihood analysis is not only more sensitive than binned methods, but it can also help extract astrophysical information. Figure 14 shows our ability to reconstruct the spectral index for power-law neutrino sources at a declination of 6°. The effective area is high for a broad range of energies here, and the spectral resolution is best. For each spectrum shown, the statistical uncertainty (1σ CL) in the spectral index will be about ±0.3 when enough events are present to claim a discovery. Spectral resolution worsens to ±0.4 at both δ = −45° and δ = 45° when enough events are present for a discovery in each case.

Zoom In Zoom Out Reset image size

Stacking multiple sources in neutrino astronomy has been an effective way to enhance discovery potential and further constrain astrophysical models (Achterberg et al. 2006b; Abbasi et al. 2009c). We can consider the accumulated signal from a collection of sources using a method similar to Abbasi et al. (2006). Only a modification to the signal likelihood is necessary in order to stack sources, breaking the signal hypothesis into the sum over M sources:

$$\begin{equation} \mathcal {S}_i \Rightarrow \mathcal {S}_i^{{\rm tot}} = \frac{\sum _{j=1}^{M} W^j R^j(\gamma) \mathcal {S}_i^j}{\sum _{j=1}^{M} W^j R^j(\gamma) }, \end{equation} \tag{ 8 }$$

where W^j is the relative theoretical weight, R^j(γ) is the relative detector acceptance for a source with spectral index γ (assumed to be the same for all stacked sources), and $\mathcal {S}_i^j$ is the signal probability density for the ith event, all for the jth source. As before, the total signal events n_s and collective spectral index γ are fit parameters. The W^j coefficients depend on our prior theoretical assumptions about the expected neutrino luminosity. They are higher for sources that are, on theoretical grounds, expected to be brighter. Tables for R^j(γ), given as the mean number of events from a source with dΦ/dE∝E^−γ, are calculated using simulation. The flexibility built into the method by the relative detector acceptance and theoretical weights allows us to use source catalogs covering the whole sky and with large variations in source strengths, as well as to directly test model predictions.

We would also like to consider sources that are spatially extended (with respect to the PSF). For an example of how important this can be, the significance observed by the Milagro experiment in the location of the Fermi source J0634.0+1745 (associated with the Geminga pulsar) rises from 3.5σ to 6.3σ by fitting for an extended source (Abdo et al. 2009). The only modification to the method required is to convolve the source distribution with the PSF. Since we model our PSF as a circular two-dimensional Gaussian distribution, it is easy to also model a source as a circular two-dimensional Gaussian of width σ_s. The convolution results in a broader two-dimensional Gaussian of width $\sqrt{\vphantom{A^A}\smash{\hbox{${\sigma _i^2 + \sigma _{\mathrm{s}}^2}$}}}$ and the likelihood uses this distribution for the signal spatial term. The discovery potential flux for a range of source extensions is shown in Figure 15 and compared to the (incorrect) hypothesis of a point source. For a source with true extent σ_s = 2°, the point-source hypothesis requires nearly a factor of two times more flux for discovery compared to the correct extended-source hypothesis.

Zoom In Zoom Out Reset image size

We have performed five searches:

• 1.  
  a scan for the most significant point source in the entire sky;
• 2.  
  a search over an a priori defined list of 39 interesting astrophysical objects;
• 3.  
  a stacking search for 16 Milagro TeV gamma-ray sources, some seen only in coincidence with the Fermi-LAT, and one unconfirmed hot spot (17 total sources);
• 4.  
  a stacking search for 127 local starburst galaxies (Becker et al. 2009);
• 5.  
  a stacking search for five nearby clusters of galaxies (CGs), testing four different models for the CR spatial distribution (Murase et al. 2008).

The analyses and event selection procedure were determined before unblinding the R.A. of the data. We require a 5σ significance for discovery. Final p-values are calculated for each search individually.

6.1. All-sky Scan

6.2. A Priori Source List

6.3. Milagro TeV Source Stacking

6.4. Starburst Galaxy Stacking

6.5. Galaxy Cluster Stacking

CGs are another potential source of high-energy protons and, through interactions with intracluster material (ICM), neutrinos. CGs are the largest gravitationally bound objects in the universe and continue to grow through merging and accretion of dark matter and baryonic gas, generating shock fronts on megaparsec scales. The possibility for CGs to be sources of ultra-high-energy CRs above 3 × 10¹⁸ eV is described in, e.g., Norman et al. (1995) and Kang et al. (1997). Murase et al. (2008) discuss the possibility of CGs being a significant contribution to the CR spectrum between the second knee at about 3 × 10¹⁷ eV and the ankle at about 3 × 10¹⁸ eV. They give predictions for neutrinos from five nearby (z < 0.03) CGs: Virgo, Centaurus, Perseus, Coma, and Ophiuchus. Information on location, distance, and size of CGs (virial radii) was taken from Reiprich & Boehringer (2002). These nearby CGs appear to us as spatially extended objects with virial radii subtending 13–69, so an extended spatial distribution of neutrinos is possible. Whereas the distribution of the ICM is well known from X-ray observations (Pfrommer & Ensslin 2004), the distribution of CRs is highly uncertain. The distribution of neutrinos is given by the product of the CR and ICM distributions. Four CR models have been considered for neutrino production, discussed in Murase et al. (2008) and references therein (e.g., Colafrancesco & Blasi 1998; Berezinsky et al. 1997):

• 1.  
  Model A. CRs are uniformly distributed within the cluster shock radius, taken to be 0.56 of the virial radius for the dynamical parameters considered.
• 2.  
  Model B. CRs are uniformly distributed within the virial radius, yielding the most conservative neutrino flux distributed over the largest area.
• 3.  
  Isobaric. CRs follow the distribution of thermal gas.
• 4.  
  Central AGN. In a two-step acceleration scenario, CRs are accelerated in the central AGN up to a maximum energy before diffusing throughout the cluster and possibly undergoing further acceleration. For the purposes of IceCube searches, this model can be treated as a point source. This model is discussed in detail by Kotera et al. (2009).

Signal neutrinos were simulated according to each of the four models. We modeled the source extensions in the likelihood as two-dimensional Gaussian distributions with the width for each source and each model determined by optimizing for the best discovery potential. Although the modeling of the source extension as a Gaussian in the likelihood is not ideal, it is straightforward and computationally fast. The exact shape of the sources is not important for small signals; we may be able to analyze the shape with more detail depending on the intensity of any signal.

We performed a stacking search for five nearby CGs mentioned above following the model predictions of Murase et al. (2008) as weights in the likelihood. The size of the clusters in the likelihood fit was allowed to vary discretely between the optimal widths for each CR distribution model. The optimal width and ν_μ differential flux for each source and each model are given in Table 2. The differential fluxes are parameterized as broken power laws:

$$\begin{equation} \frac{d\Phi }{{\it dE}} (\mathrm{TeV}^{-1}\,\mathrm{cm}^{-2}\,\mathrm{s}^{-1}) = \left\lbrace \begin{array}{@{}lr}A \cdot (E / \mathrm{TeV})^{-\gamma _1} & E \le E_{\mathrm{break}}, \\ B \cdot (E / \mathrm{TeV})^{-\gamma _2} & E > E_{\mathrm{break}}. \end{array} \right. \end{equation} \tag{ 9 }$$

The parameter B = A · Eγ2 − γ1break after enforcing continuity at the break energy.

Table 2. Galaxy Cluster Parameters

Source	R.A. (°)	Decl. (°)	Model	σs (°)	A (TeV−1 cm−2 s−1)	γ1	γ2	Ebreak (TeV)
Virgo	186.63	12.72	Model A	2.0	1.42 × 10−12	−2.14	−4.03	2.16 × 106
			Model B	4.0	1.18 × 10−12	−2.14	−4.03	2.16 × 106
			Isobaric	3.0	7.57 × 10−13	−2.14	−4.03	2.16 × 106
			Central AGN	0.0	6.47 × 10−12	−2.42	−4.24	2.13 × 106
Centaurus	192.20	−41.31	Model A	0.25	2.78 × 10−13	−2.14	−4.03	2.15 × 106
			Model B	0.5	2.20 × 10−13	−2.14	−4.03	2.15 × 106
			Isobaric	0.25	1.09 × 10−13	−2.15	−4.07	2.33 × 106
			Central AGN	0.0	5.10 × 10−13	−2.45	−4.28	2.39 × 106
Perseus	49.95	41.52	Model A	0.0	5.83 × 10−14	−2.15	−4.07	2.32 × 106
			Model B	0.5	4.60 × 10−14	−2.15	−4.07	2.32 × 106
			Isobaric	0.0	6.17 × 10−13	−2.15	−4.07	2.32 × 106
			Central AGN	0.0	5.97 × 10−13	−2.40	−4.20	1.88 × 106
Coma	194.95	27.94	Model A	0.25	2.14 × 10−14	−2.14	−4.03	2.12 × 106
			Model B	0.25	1.34 × 10−14	−2.14	−4.03	2.12 × 106
			Isobaric	0.25	1.83 × 10−13	−2.15	−4.07	2.30 × 106
			Central AGN	0.0	2.13 × 10−13	−2.41	−4.20	1.89 × 106
Ophiuchus	258.11	−23.36	Model A	0.0	4.87 × 10−14	−2.15	−4.07	2.29 × 106
			Model B	0.5	1.50 × 10−14	−2.15	−4.07	2.29 × 106
			Isobaric	0.0	5.50 × 10−13	−2.15	−4.11	2.49 × 106
			Central AGN	0.0	2.55 × 10−13	−2.43	−4.24	2.12 × 106

Notes. σ_s is the optimized sigma of a two-dimensional Gaussian distribution used in the likelihood. Numerical calculations of the differential fluxes (K. Murase 2009, private communication) for each model described in Murase et al. (2008) are fit well to broken power laws, parameterized in Equation (9).

Download table as:  ASCIITypeset image

The analyses described in Section 6 give reliable statistical results (p-values) due to the ability to generate background-only data sets by scrambling the data in R.A. By using the data to estimate background, the systematic errors come only from signal and detector simulation used to calculate flux upper limits, sensitivities, and discovery potentials.

The main systematic uncertainties on the flux limits come from photon propagation in ice, absolute DOM sensitivity, and uncertainties in the Earth density profile as well as muon energy loss. All numbers are for a spectrum dΦ/dE∝E⁻² of muon neutrinos. We evaluate the systematic uncertainty due to photon propagation by performing dedicated signal simulations with scattering and absorption coefficients varied within their uncertainties of ±10% (Ackermann et al. 2006). The maximum difference was between the case where both scattering and absorption were increased by 10% and the case where both were decreased by 10%. The deviation in the observed number of events between these two cases was 11%. The range of uncertainty in the DOM sensitivity is taken as ±8%, based on the measured uncertainty in the PMT sensitivity (Abbasi et al. 2010a). Similarly, another dedicated simulation where we varied the DOM sensitivity inside the uncertainty leads to a maximum systematic uncertainty on the number of events of 9%. These uncertainties on the flux varied by only about 2% between the northern and southern sky, so only averages over the whole sky are reported. Finally, uncertainties in muon energy losses, the neutrino–nucleon cross section, and the rock density near the detector introduce an 8% systematic uncertainty for vertically up-going events (Achterberg et al. 2007). These events are the most affected, and this uncertainty is applied to all zeniths to be conservative. A sum in quadrature of the systematic uncertainties on the flux gives a total of 16% systematic uncertainty in the signal simulation. These systematic uncertainties are incorporated into the upper limit and sensitivity calculations using the method of Conrad et al. (2003) with a modification by Hill (2003).

The results of the all-sky scan are shown in the map of the pre-trial p-values in Figure 16. The most significant deviation from background is located at 11375 R.A., 1515 decl. The best-fit parameters are n_s = 11.0 signal events above background, with spectral index γ = 2.1. Since the best-fit spectral index is substantially less than the expectation from background, much of the significance comes from the higher energies of the associated events. The pre-trial estimated p-value of the maximum log-likelihood ratio at this location is 5.2 × 10⁻⁶. In trials using data sets scrambled in R.A., 1817 out of 10,000 have an equal or higher significance somewhere in the sky, resulting in the post-trial p-value of 18%. Upper limits for a flux dΦ/dE∝E⁻² of $\nu _{\mu } + \bar{\nu }_{\mu }$ are presented in Figure 17. In all cases, an equal flux of neutrinos and antineutrinos is assumed.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The results of the point-source search in the direction of 39 source candidates selected a priori are given in Table 3 and also shown in Figure 18 with the IceCube median sensitivity. Since the fit was restricted to physical signal values n_s ⩾ 0, approximately half of the results have n_s = 0 exactly, corresponding to p-values equal to 100% and upper limits equal to the median upper limit (i.e., the sensitivity). The most significant source on the list was PKS 1622−297 with a pre-trial estimated p-value of 5%. The post-trial p-value of 62% was again determined as the fraction of scrambled data sets with at least one source with an equal or higher significance. The mean number of events at the final cut level required for the discovery of a point source is also shown in Figure 19, along with the average background in a circular bin with 1° radius. Included in Figure 18 is a preliminary comparison to the ANTARES experiment (Coyle 2010). ANTARES is primarily sensitive to GeV–TeV energy neutrinos in the southern sky, so the coverage in energy is quite complementary to this IceCube analysis.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3. Results for the A Priori Source Candidate List

Object	R.A. (°)	Decl. (°)	$\Phi _{\nu _{\mu }}^{90}$	$\Phi _{\nu _{\mu }+\nu _{\tau }}^{90}$	p-value	ns	γ	$N_{1^{\circ }}$	$B_{1^{\circ }}$
Cyg OB2	308.08	41.51	6.04	10.54	1.00	0.0	...	2	1.8
MGRO J2019+37	305.22	36.83	7.50	13.3	0.44	1.0	2.8	2	1.9
MGRO J1908+06	286.98	6.27	3.73	6.82	0.43	1.5	3.9	4	3.1
Cas A	350.85	58.81	9.04	15.92	1.00	0.0	...	1	1.8
IC443	94.18	22.53	3.80	6.62	1.00	0.0	...	1	2.0
Geminga	98.48	17.77	3.91	6.66	0.48	0.7	2.1	1	2.3
Crab Nebula	83.63	22.01	3.70	6.58	1.00	0.0	...	1	2.0
1ES 1959+650	300.00	65.15	10.74	19.18	1.00	0.0	...	0	2.0
1ES 2344+514	356.77	51.70	7.24	12.96	1.00	0.0	...	0	1.8
3C66A	35.67	43.04	10.89	19.70	0.24	3.4	3.9	3	1.9
H 1426+428	217.14	42.67	6.14	10.94	1.00	0.0	...	3	1.8
BL Lac	330.68	42.28	10.80	18.70	0.25	2.6	3.9	3	1.8
Mrk 501	253.47	39.76	8.11	14.14	0.41	1.3	3.9	3	2.0
Mrk 421	166.11	38.21	11.71	20.14	0.15	2.6	1.9	2	2.0
W Comae	185.38	28.23	4.46	8.06	1.00	0.0	...	0	1.9
1ES 0229+200	38.20	20.29	6.89	12.06	0.19	4.0	2.8	4	2.1
M87	187.71	12.39	3.42	5.98	1.00	0.0	...	2	2.5
S5 0716+71	110.47	71.34	13.28	23.56	1.00	0.0	...	0	1.6
M82	148.97	69.68	19.14	32.84	0.40	2.0	3.9	4	1.8
3C 123.0	69.27	29.67	5.59	10.66	0.44	1.3	2.7	1	1.9
3C 454.3	343.49	16.15	3.42	5.92	1.00	0.0	...	1	2.3
4C 38.41	248.81	38.13	6.77	11.86	0.48	0.9	3.9	2	2.0
PKS 0235+164	39.66	16.62	6.77	11.62	0.15	5.3	3.0	5	2.3
PKS 0528+134	82.73	13.53	3.63	6.72	1.00	0.0	...	2	2.4
PKS 1502+106	226.10	10.49	3.26	5.78	1.00	0.0	...	0	2.5
3C 273	187.28	2.05	3.61	6.54	1.00	0.0	...	3	3.4
NGC 1275	49.95	41.51	6.04	10.54	1.00	0.0	...	2	1.8
Cyg A	299.87	40.73	7.84	13.44	0.46	1.0	3.5	3	1.9
IC-22 maximum	153.38	11.38	3.26	5.86	1.00	0.0	...	1	2.5
Sgr A*	266.42	−29.01	80.56	139.26	0.41	1.1	2.7	4	3.3
PKS 0537−441	84.71	−44.09	113.90	201.82	1.00	0.0	...	3	3.5
Cen A	201.37	−43.02	109.51	191.56	1.00	0.0	...	4	3.5
PKS 1454−354	224.36	−35.65	92.56	156.74	1.00	0.0	...	4	3.5
PKS 2155−304	329.72	−30.23	105.41	182.90	0.28	1.7	3.9	3	3.4
PKS 1622−297	246.53	−29.86	152.28	263.86	0.048	3.0	2.6	4	3.3
QSO 1730-130	263.26	−13.08	24.83	43.30	1.00	0.0	...	4	3.5
PKS 1406−076	212.24	−7.87	16.04	28.72	0.42	1.3	2.3	4	3.3
QSO 2022-077	306.42	−7.64	12.18	21.78	1.00	0.0	...	2	3.3
3C279	194.05	−5.79	11.94	21.36	0.33	3.6	3.0	7	3.5

Notes. $\Phi _{\nu _{\mu }}^{90}$ and $\Phi _{\nu _{\mu }+\nu _{\tau }}^{90}$ are the upper limits of the Feldman–Cousins 90% confidence intervals for a dΦ/dE∝E⁻² flux normalization of ν_μ and ν_μ + ν_τ, respectively, in units of 10⁻¹² TeV⁻¹ cm⁻² s⁻¹ (i.e., dΦ/dE ⩽ Φ⁹⁰ · (E/TeV)⁻²). n_s is the best-fit number of signal events; the (pre-trial) p-value is also calculated and the spectral index γ is given when n_s > 0. N$_{1^{\circ }}$ is the actual number of events observed in a bin of radius 1°. The background event density at the source declination is indicated by the mean number of background events $B_{1^{\circ }}$ expected in a bin of radius 1°.

Download table as:  ASCIITypeset image

The Milagro TeV source stacking search resulted in a final p-value of 32% with best-fit n_s = 7.6 (total number of signal events above background) and spectral index γ = 2.6. The starburst galaxy stacking search resulted in an underfluctuation with best-fit n_s = 0 and a final p-value of 100%, since we do not allow negative values of n_s. Finally, the CGs stacking search yielded a final p-value of 78% with n_s = 1.8. These results, along with sensitivities and upper limits, are summarized in Table 4.

The IceCube Neutrino Observatory aims to further our understanding of astrophysical phenomena, constraining models even in the absence of a detection. Figure 20 shows our sensitivity to some specific predictions. The model of Morlino et al. (2009) is for SNR RX J1713.7-3946. This analysis is largely insensitive to spectra which cut off below 100 TeV in the southern sky. Applying this emission model at the location of the Crab Nebula (δ = 2201), we obtain an upper limit that rules out a flux 3.2 times higher than the prediction.

Zoom In Zoom Out Reset image size

The Milagro hot spot MGRO J1852+01 is the brightest source of six SNR associations considered in Halzen et al. (2008). The stacking results were already shown in Table 4. Our upper limit for just this one brightest source is a factor of 7.9 away from excluding this model at 90% CL. The best fit for MGRO J1852+01 is to 7.0 events with γ = 2.9, which increases the upper limit compared to the average background-only case.

The nearest AGN, Centaurus A (Cen A), has been discussed as a potential source of ultra-high-energy CRs in the context of results from the Pierre Auger Observatory (PAO). The point-source likelihood fit at the location of Cen A resulted in zero signal events in this analysis, setting an upper limit that is 5.3 times higher than the ν_μ prediction by Koers & Tinyakov (2008) for the most optimistic case where the protons have a spectral index α_p = 3.

Starburst galaxies were already presented as sources of interest in Section 6.4. Recent detections (Abdo et al. 2010; Acciari et al. 2009; Aharonian et al. 2005) of very high-energy photons from the nearest luminous starburst galaxies M82 and NGC 253, each characterized by star-forming regions with high supernova rates in the core, support the belief that the observed enhanced gamma-ray emission is due to CR interactions. Under the assumption that the GeV–TeV photons originate from the decay of neutral pions produced when protons that are shock-accelerated by SNRs in the starburst core inelastically scatter against target hydrogen atoms with densities of the order of 10² cm⁻³ (de Cea del Pozo et al. 2009; Persic et al. 2008), an order-of-magnitude calculation of the resulting flux of muon neutrinos based on Kelner et al. (2006) can be made. The muon neutrino upper limit from M82 is about 400 times higher than the rough prediction. For NGC 253 in the southern sky, the muon neutrino upper limit is about 6000 times higher than the prediction.

A search for sources of high-energy neutrinos has been performed using data taken during 2008–2009 with the 40-string configuration of the IceCube Neutrino Observatory. Five searches were performed: (1) a scan of the entire sky for point sources, (2) a predefined list of 39 interesting source candidates, (3) stacking 16 sources of TeV gamma rays observed by Milagro and Fermi, along with an unconfirmed hot spot (17 total sources), (4) stacking 127 starburst galaxies, and (5) stacking five nearby CGs, testing four different models for the CR distribution. The most significant result of the five searches came from the all-sky scan with a p-value of 18%. The cumulative binomial probability of obtaining at least one result of this significance or higher in five searches is 63%. This result is consistent with the null hypothesis of background only.

The sensitivity of this search using 375.5 days of 40-string data already improves upon previous point-source searches in the TeV–PeV energy range by at least a factor of two, depending on declination. The searches were performed using a data set of up-going atmospheric neutrinos (northern sky) and higher energy down-going muons (southern sky) in a unified manner. IceCube construction is now complete, with 86 strings in operation. The full IceCube detector should improve existing limits by at least another factor of two with one year of operation. Additional improvement is foreseeable in the down-going region by developing sophisticated veto techniques and at lower energies by using the new dense sub-array, DeepCore, to its fullest potential.

We thank K. Murase for helpful discussions regarding clusters of galaxies.

We acknowledge the support from the following agencies: U.S. National Science Foundation–Office of Polar Programs, U.S. National Science Foundation–Physics Division, University of Wisconsin Alumni Research Foundation, the Grid Laboratory Of Wisconsin (GLOW) grid infrastructure at the University of Wisconsin-Madison, the Open Science Grid (OSG) grid infrastructure; U.S. Department of Energy, and National Energy Research Scientific Computing Center, the Louisiana Optical Network Initiative (LONI) grid computing resources; National Science and Engineering Research Council of Canada; Swedish Research Council, Swedish Polar Research Secretariat, Swedish National Infrastructure for Computing (SNIC), and Knut and Alice Wallenberg Foundation, Sweden; German Ministry for Education and Research (BMBF), Deutsche Forschungsgemeinschaft (DFG), Research Department of Plasmas with Complex Interactions (Bochum), Germany; Fund for Scientific Research (FNRS-FWO), FWO Odysseus Programme, Flanders Institute to encourage scientific and technological research in industry (IWT), Belgian Federal Science Policy Office (Belspo); University of Oxford, UK; Marsden Fund, New Zealand; Japan Society for Promotion of Science (JSPS); the Swiss National Science Foundation (SNSF), Switzerland; A. Groß acknowledges support by the EU Marie Curie OIF Program; J. P. Rodrigues acknowledges support by the Capes Foundation, Ministry of Education of Brazil.