Indirect Search for Dark Matter from the Galactic Center and Halo with the Super-Kamiokande Detector

We present a search for an excess of neutrino interactions due to dark matter in the form of Weakly Interacting Massive Particles (WIMPs) annihilating in the galactic center or halo based on the data set of Super-Kamiokande-I, -II, -III and -IV taken from 1996 to 2016. We model the neutrino flux, energy, and flavor distributions assuming WIMP self-annihilation is dominant to $\nu \overline{\nu}$, $\mu^+\mu^-$, $b\overline{b}$, or $W^+W^-$. The excess is in comparison to atmospheric neutrino interactions which are modeled in detail and fit to data. Limits on the self-annihilation cross section $\langle \sigma_{A} V \rangle$ are derived for WIMP masses in the range 1 GeV to 10 TeV, reaching as low as $9.6 \times10^{-23}$ cm$^3$ s$^{-1}$ for 5 GeV WIMPs in $b\bar b$ mode and $1.2 \times10^{-24}$ cm$^3$ s$^{-1}$ for 1 GeV WIMPs in $\nu \bar \nu$ mode. The obtained sensitivity of the Super-Kamiokande detector to WIMP masses below several tens of GeV is the best among similar indirect searches to date.

We present a search for an excess of neutrino interactions due to dark matter in the form of Weakly Interacting Massive Particles (WIMPs) annihilating in the galactic center or halo based on the data set of Super-Kamiokande-I, -II, -III and -IV taken from 1996 to 2016. We model the neutrino flux, energy, and flavor distributions assuming WIMP self-annihilation is dominant to νν, µ + µ − , bb, or W + W − . The excess is in comparison to atmospheric neutrino interactions which are modeled in detail and fit to data. Limits on the self-annihilation cross section σAV are derived for WIMP masses in the range 1 GeV to 10 TeV, reaching as low as 9.6 × 10 −23 cm 3 s −1 for 5 GeV WIMPs in bb mode and 1.2×10 −24 cm 3 s −1 for 1 GeV WIMPs in νν mode. The obtained sensitivity of the Super-Kamiokande detector to WIMP masses below several tens of GeV is the best among similar indirect searches to date.

I. INTRODUCTION
There is compelling evidence that ordinary baryonic matter composes only 5% of the total energy density of the Universe, which is dominated by dark energy (68%) and dark matter (27%) whose nature is unknown [1]. Some well-motivated candidates for particle dark matter (DM) arise within supersymmetric extensions of the Standard Model [2,3]. These particles belong to a collective group referred to as Weakly Interacting Massive Particles (WIMPs). The lightest supersymmetric particle, the neutralino (χ), has been considered one of the most promising WIMP candidates. However, the analysis presented here would be also valid for any other dark matter candidate annihilating into standard model particles.
WIMPs present in the galactic halo may be observed directly via elastic scattering off nuclei in detectors [4,5] or indirectly through detection of the products of their annihilations (or decays) into standard model particles, including neutrinos [6][7][8][9]. It is expected that dark matter particles will accumulate in massive celestial objects like stars or planets and be bound by their gravitational potentials [10]. Previous searches for WIMP-induced neutrinos from the Sun and Earth cores, based on data collected with the Super-Kamiokande (SK) detector, have shown no excess of dark-matter-induced neutrinos over the atmospheric neutrino background [11][12][13].
In this study, we assume that dark matter particles concentrate in the central regions of galaxies, as predicted by many halo models [14][15][16][17][18][19]. We search for neutrinos from WIMP annihilation in the center of the Milky Way and from its halo. This search with neutrinos is complementary to other indirect searches relying on annihilation products like γ, e + /e − orp [20]. The experimental advantage of neutrinos is that they travel unimpeded and undeflected from their origin. However, their low interaction cross section puts them at a disadvantage relative to these other annihilation products.
We constrain thermally averaged self-annihilation cross section σ A V for WIMP pair annihilation for their masses from 1 GeV to 10 TeV. This is the first search for WIMPs from the Milky Way based on data acquired by the Super-Kamiokande detector and extending to dark matter masses below 10 GeV.
The Super-Kamiokande detector is nearly 100% efficient for detecting neutrino interactions above ∼ 100 MeV. The dominant background is atmospheric neutrino interactions, which have been well studied with our detector. We have made precise measurements of the flux, angular distribution and energy spectra of atmospheric neutrinos [21]. We also have a detailed model of systematic uncertainties [22].
Below we present results of two different approaches to analyzing the same data sample, a full fit employing both data and Monte Carlo templates (combined fit) and a method comparing observed event rates from a region around the galactic center to those in a region located in the opposite part of the sky (on-source off-source analysis). The first method fits for a simulated WIMP-induced neutrino contribution to the Super-Kamiokande data together with the expected atmospheric neutrino background component. Various WIMP annihilation channels and masses are tested. All event categories of detected neutrino interactions used in recent atmospheric analyses [22], binned in angle and momentum, are used in the fit. The second method provides a check that is entirely data-driven. We compare the number of events in a certain angular region around the galactic center (onsource) and in a region of the same size but shifted by 180 • in right ascension (off-source). As the atmospheric neutrino background equally affects both regions, any excess of events in the on-source data would indicate an additional source of neutrinos from the area around the galactic center.

II. DARK MATTER ANNIHILATION IN THE MILKY WAY
The energy spectrum of neutrinos from the annihilation of dark matter WIMPs and the branching ratio (BR) for their production is unknown and must be modeled. In the following analysis we consider direct annihilation into pairs of neutrinos, χχ → νν, as well as annihilation into pairs of µ + µ − , bb and W + W − . Each annihilation mode is considered separately assuming 100% BR irrespective of the actual nature of DM annihilation. In the bb and W + W − channels, neutrinos are mainly created in decays of mesons produced during the hadronization of the primary annihilation products. In the µ + µ − channel neutrinos are produced directly in the decays of the muons. In the νν channel, a monoenergetic spectrum and equal fluxes of DM-induced neutrinos of every flavor are assumed, and the energy of the neutrinos equals the mass of the annihilating dark matter particles.
Neutrinos that travel over galactic distances experience multiple oscillations and arrive at the detector with a flavor composition that is predictable based on the values of the neutrino mixing parameters [23,24]. In applying oscillations to DM-induced neutrinos we followed the approach presented in [25], adopting parameter values of sin 2 2θ 12 = 0.86, sin 2 2θ 23 = 1.0 (maximal mixing) and θ 13 0 to check how long-distance oscillations affect the neutrino production spectra and their flavor ratios at Earth. The resulting neutrino fluxes at the detector can be described with the following effective formulas: and where φ 0 νe,µ,τ are the initial fluxes, and s 2 is defined as [25]. Note that the above formulas lead to an equal flux of neutrino flavors at the Earth, that is flux ratios of 1:1:1 for ν e :ν µ :ν τ , for the typical scenario of neutrino production via the decays of light mesons and muons in which the initial flavor ratio is 1:2:0 at the production point. The DarkSUSY [26] simulation package was used to obtain the differential neutrino energy spectra for the DM annihilation channels considered as shown in Fig. 1 after taking into account neutrino oscillations.
The expected flux of DM annihilation products depends on the density distribution of DM particles in the Milky Way. There are various models describing the structure of the DM halo obtained on the basis of N-body simulations [27] and gravitational lensing observations [28]. The predictions of the Navarro-Frenk-White (NFW) [14], the Kravtsov et al. [15] and the Moore et al. [16] models of the expected DM density ρ in the halo are shown in Fig. 2, assuming a flat DM density profile for the inner most regions of the galaxy. The Moore model anticipates a high density cusp in the center of the Milky Way, while the Kravtsov model yields a flatter density profile. The NFW profile is between these two extreme predictions and is similar to other commonly considered models such as Einasto [17,18]. In the analysis presented in this paper the NFW model is used as a benchmark.
The DM density distribution in the NFW model can be written as: where r is the distance from the center of the galaxy and r s = 20 kpc is the scale radius. The normalization is set so as to match the DM density at the radius of the solar system (R sc = 8. Moore [16], NFW [14] and Kravtsov [15]. curves and corresponds to ρ(R sc ) = 0.3 GeV cm −3 for the NFW profile [8,9]. The annihilation intensity J a (numerical flux per solid angle) at an angle ψ with respect to the galactic center (GC) direction is proportional to the line of sight (l) integral of the DM density squared (as two DM particles are required for annihilation) [8,9] J a (ψ) = 1 R sc ρ 2 where ρ 2 sc = ρ 2 (R sc ) and R sc are used as scaling parameters to make J a dimensionless. The upper limit of integration, depends on the assumed size of the Milky Way halo, R M W . However, contributions beyond the size of the visible stellar halo (20-30 kpc) are negligible. The dependence of J a (ψ) on the angular distance from the GC is shown in Fig. 3. The average value of the intensity of DM annihilation products received at the Earth from a cone with halfangle ψ around the GC, that spans a field of view of ∆Ω = 2π(1 − cos ψ), can be cast as Then, the neutrino flux at the Earth can be related to the average intensity of DM annihilation products as [8,9] where M χ is the assumed mass of the DM particle, the factor 1/2 is needed for self-conjugate WIMPs, and 1/4π is for isotropic emission. Here dN/dE is the differential neutrino multiplicity per one annihilation (cf. Fig. 1).

III. SUPER-KAMIOKANDE DETECTOR AND DATA SAMPLES
Super-Kamiokande is a 50 kton water Cherenkov detector located at the Kamioka Observatory operated by the Institute for Cosmic Ray Research of the University of Tokyo [29]. It is used to search for proton decay as well as to investigate atmospheric, man-made and extraterrestrial neutrinos, including solar neutrinos and those produced in supernovae bursts. The detector consists of an inner cylindrical volume (inner detector, ID) and an outer part (outer detector, OD), which plays the role of a veto region for penetrating particles. The detection of neutrinos relies on observation of charged particles, primarily leptons, produced in ν interactions inside or outside the ID's 22.5 kton fiducial volume. Cherenkov radiation from charged particle in the water is emitted in a characteristic conical pattern, is projected onto the detector walls, and is recorded by photomultipliers. The detected light pattern allows for the reconstruction of the particle energy and direction and provides differentiation between electromagnetic showers (e-like) and muon-type non-showering particles (µ-like).
The sample of neutrino interactions in which we search for WIMP-induced neutrinos consists mainly of interactions of atmospheric neutrinos, which are the main background of this search. Based on the topology and energy of the detected events, atmospheric neutrinos can be assigned to three main categories: fully-contained (FC), partially-contained (PC), and upward-going muons (UPµ). The FC events have a reconstructed neutrino interaction vertex inside the fiducial volume and particles produced in the parent neutrino interaction stop inside the ID. The true neutrino energy of FC events is in the range of hundreds of MeV up to several GeV. The true neutrino energy of PC events, which have particles exiting into the OD, ranges from 1 GeV to 100 GeV with median energy 10 GeV. Neutrinos of higher energies, 10 GeV to a few TeV, are detected as upward-going muons due to muon neutrino interactions in the surrounding rock. Downward-going muons entering from outside the detector are removed from the analysis samples as they cannot be separated from downward-going cosmic ray muons.
The FC, PC and UP-µ event classes can be further divided into more specific subcategories as listed in Table I. The main criteria for dividing FC events are the type of the primary ring (e-like vs. µ-like), the number of reconstructed rings (single-ring or multi-ring), the number of electrons from muon decays, the presence of π 0 , and the likelihood to beν e . If the reconstructed momentum of the primary lepton is below 1.33 GeV, the event is classified as sub-GeV, and multi-GeV otherwise. The PC events either stop in or escape the veto region and are thus divided into PC-stopping and PC-throughgoing. Similarly, UP-µ events are categorized as "stopping" in the detector or "through-going". The most energetic muons induce showers while passing through the detector. Therefore, "through-going" UP-µ events are further categorized as "showering" and "non-showering". A detailed description of the event classification can be found elsewhere [22].
Since this search covers a wide range of neutrino energies, all of the SK atmospheric neutrino data and their corresponding subsamples are taken into account in the analysis. The search is based on data from the SK-I The Monte Carlo (MC) simulation of the detector response is based on a GEANT3 model [30]. Neutrino interactions are modeled with the NEUT generator [31] and the background flux is taken from Honda et al. [32]. In this analysis, MC corresponding to 500 years of livetime separately for SK-I, -II, -III and -IV sets (2000 years of livetime in total) is used.

IV. COMBINED FIT
In this search it is assumed that the neutrino data collected with the SK detector can be described by two components, WIMP-induced neutrinos (the signal) and atmospheric neutrinos (the background). In the com- bined fit analysis we search for a signal excess in the atmospheric neutrino data sample that is compatible with a neutrino source from the GC by introducing DM signal templates and fitting them together with the atmospheric background assuming a variety of WIMP masses. The spectra and directional distributions of signal events are correlated across all of the event subsamples. This is different from previous SK searches for DM-induced neutrinos that were based only on the angular information of UP-µ events [12,13]. Monte Carlo event samples are used to simulate both signal and background. The standard simulation of atmospheric neutrino interactions used in SK's oscillation studies is used to estimate the background and includes simulated tau neutrino interactions [33]. In order to simulate the signal, a separate and independent sample of events from the atmospheric MC is reweighted to the angular and energy distributions expected for DM-induced neutrinos for a given WIMP mass and annihilation channel.
The Monte Carlo samples allow DM-induced signal simulation for ν µνµ with energies 100 MeV≤ E ν ≤ 10 TeV, for ν eνe in 100 MeV≤ E ν ≤ 80 GeV and for ν τντ in 100 MeV≤ E ν ≤ 150 GeV. These upper limits are determined by event statistics in the atmospheric MC. For WIMP masses (M χ ) greater than several hundreds of GeV, the detected event rate is dominated by interactions in the rock rather than the detector volume, therefore only muon neutrinos are used to simulate the high energy part of the DM-induced neutrino energy spectrum. Although the analysis in this energy range lacks contributions from other neutrino flavors, it has the advantage of the best pointing resolution.
The signal and background contributions are compared against the data for various WIMP masses and annihilation channels using the "pull" χ 2 [34] where n indexes the analysis bins, N AT M n is the atmospheric neutrino background expectation including oscillations, N DM n is the simulated DM-induced neutrino contribution and N data n is the number of observed events in the n th bin. The parameter β represents the normalization of the simulated signal. Determining the "best-fit" β is the main goal of this analysis. Systematic errors are incorporated into the fit via parameters i , where i is the systematic error index. Here f i n ( g i n ) represents the fractional change in the background (signal) MC in bin n for a σ i change in the i th systematic error. Note that systematic uncertainties related to the detector, data selection and neutrino interactions are fully correlated between the signal and background. That is, f i n = g i n for those errors. Atmospheric neutrino flux errors, on the other hand, have no impact on the signal and g i n = 0 for those errors, accordingly. This analysis only considers DM-induced neutrinos from the GC and halo and does not include other sources. During the fit, Eq. 8 is minimized with respect to the i at each point in parameter space according to ∂χ 2 ∂ i = 0 [34]. The "best-fit point" is defined as the global minimum χ 2 on the grid of all tested points. The index i covers systematic uncertainties from ref. [22] considered in the atmospheric neutrino analysis of SK data, but adjusted for the GC angular binning scheme used in this analysis.
The atmospheric neutrino background, N AT M n , depends on the values of neutrino oscillation parameters, θ 23 , θ 12 , θ 13 , ∆m 2 23 , ∆m 2 12 , and δ CP . The values of these parameters are not determined in this analysis as a given GC coordinate does not correspond to a fixed travel path length in the Earth and therefore provides little sensitivity to neutrino oscillations. Accordingly, the oscillation parameters are set to δ CP = 0, sin 2 2θ 23 = 1.0, ∆m 2 23 = 2.44 × 10 −3 eV 2 [35], sin 2 2θ 13 = 0.092 [36], sin 2 θ 12 = 0.32, and ∆m 2 12 = 7.46 × 10 −5 eV 2 [37]. However, the uncertainty on each oscillation parameter is included in the fit via systematic error parameters which effectively change the N AT M n prediction. There are 19 data samples used in this analysis, including both e-like and µ-like event categories. All the samples listed in Table I are used in the search. In total, there are 520 analysis bins (n = 520) as each sample is binned in momentum and the cosine of the angle between the event direction and the direction of the galactic center (cos θ GC ). The momentum binning is the same as in other SK atmospheric ν analyses [22]. The definition of the event direction depends on the number of rings in the event. In the case of the single ring events, the angle to the GC is calculated using the observed Cherenkov ring direction. For events with multiple rings, the direction of the event is obtained as the momentum-weighted average direction taken over all ring directions. The angular resolution of sub-GeV events strongly depends on the parent neutrino energy and is roughly tens of degrees on average (Fig. 32 in [38]). At higher energies, the true neutrino direction is more accurately determined due to the high Lorentz boost of the interaction products.
An illustration of the samples used in the fit is shown in Fig. 4. The samples in the first column use only one angular bin and several momentum bins due to the poor angular resolution at low neutrino energy and lepton momentum. All other samples use multiple angular and momentum bins, though the latter have been merged together in the figure. The data are simultaneously fit in all neutrino flavors. Tau neutrinos end up mainly categorized as multi-GeV single-ring and multi-ring e-like events due to the complex final state of secondary particles produced in ν τ interactions and τ decay. Depending on M χ and the annihilation mode, the signal appears in different samples of Fig. 4. Even samples containing no signal contribution constrain the background rate and associated systematic errors, which impacts all bins in the analysis.

V. RESULTS OF COMBINED FIT ANALYSIS
For each tested WIMP mass and annihilation channel, we try to determine the value of β that governs the allowed contribution to the SK data from the simulated signal. The fitted number of DM-induced neutrino events accounting for the "best-fit" signal normalization β is shown in Fig. 5 along with the expected limits as calculated by MC under the assumption of no WIMP signal (denoted as 90% and 99% C.L. sensitivities).
It should be noted that the points in Fig. 5 are not independent, as the same set of data is used for each M χ hypothesis. Accordingly, a correlated rise is seen in all annihilation channels at low WIMP masses as parts of the annihilation spectra from neighboring M χ overlap in some analysis bins. At higher WIMP masses the results fall into the unphysical region, where the fitted number of WIMPs is negative. A strong excess across all annihilation channels that exceeds the background-only expectation around a given WIMP mass is expected for a real signal. However, there is no such signal contribution allowed by the data.
The final fit result is translated into an upper 90% C.L. limit on the fitted number of DM-induced neutrinos using a Bayesian approach [39] to adjust results which fall into the unphysical region. Figure 6 shows the result together with the sensitivity assuming no WIMP contri-  8. Upper limits at 90% C.L. on the DM self-annihilation cross section σAV (the region above the lines is excluded) obtained in the combined fit. SK limits for νν, bb, W + W − , µ + µ − obtained for the NFW halo profile are compared to results from IceCube [41] and ANTARES [42].
bution. The 1 − 2σ uncertainty was derived from pseudo MC data sets constructed without any WIMP-induced neutrinos.
Prior to fitting the data, MC studies assuming a WIMP signal were used to determine a significance threshold for the analysis. An excess above the atmospheric background with a local p-value greater than 3σ was chosen as the criterion for a possible WIMP signal. Figure 7 shows the p-value distributions for the obtained fit results. There is a ∼ 2σ excess observed in the M χ = 5 GeV bb channel (2.08σ) , for the M χ = 1 GeV µ + µ − (1.74σ) channel, and for the M χ = 2 GeV νν (1.82σ) channel. Though all p-values are consistent with no WIMP contribution, additional checks were performed for the most significant result.
A fit in which the position of the GC was treated as a free parameter moving in right ascension (RA), but fixed in declination (DEC), was performed for the bb annihilation channel and M χ = 5 GeV. Fixing the declination allows us to use the same systematic uncertainties, as they depend on the zenith angle, which itself is proportional to the declination. This additional analysis found a similar signal excess of approximately 2σ for a WIMP signal, but for a GC position of RA in range from 210 • to 260 • , though the true GC position is at 266 • . Accordingly, we find no indication of a signal consistent with the expectation for DM from the GC halo in this channel.
Based on the limit on the number of DM-induced neutrinos, the corresponding limit on the diffuse flux is de-rived as a function of M χ using Eq. 7 and translated into a limit on the DM self-annihilation cross section σ A V . The latter is shown in Fig. 8 and compared to limits obtained by other neutrino experiments. The line at the σ A V = 3 · 10 −26 cm 3 s −1 is the expectation for WIMPs produced thermally during the evolution of the Universe [40]. Despite the smaller effective area of the SK detector when compared to the IceCube detector [41], the limits obtained in this analysis are stronger due to the fact that the SK detector can probe the GC with upward-going events. At the location of the SK detector, the GC is below the horizon for ∼ 64% of the year, while for IceCube it is always above the horizon and can only be directly probed with downward-going events which typically suffer from more backgrounds from cosmic ray muons. ANTARES is operating in the same hemisphere as SK and its limits [42] are stronger than the ones obtained here for M χ > 50 − 500 GeV (depending on the annihilation mode) due to the larger effective area of its detector. Weaker constraints from ANTARES observed for M χ < 500 GeV for bb and for the M χ < 100−150 GeV for W + W − /µ + µ − annihilation channels are due to the different detection thresholds between ANTARES and SK. In this M χ range, a substantial part of the DMinduced neutrino signal is expected below several tens of GeV (cf. Fig 1), a region that is well covered by the SK detector.

VI. ON-SOURCE OFF-SOURCE ANALYSIS
The on-source off-source analysis provides a datadriven cross check of the analysis shown in the previous section, albeit with weaker sensitivity. It has the advantage of being able to estimate the background directly from data. Equally sized on-and off-source regions are defined in right ascension and declination as shown in Fig. 9. Most of the signal is expected to come from the on-source region centered around the GC (266 • RA,−29 • DEC) while an independent background estimation is obtained from the off-source region, which is offset 180 • in right ascension but at the same declination. Note that the atmospheric neutrino background is expected to be the same in both regions as they correspond to identical zenith angles in SK's local coordinate system. Therefore, the expected number of events in the on-source region can be interpreted as N ON = N bkg on + N sig on , while for the off-source it is N OF F = N bkg of f + N sig of f . Here N sig (N bkg ) stands for the number of signal (background) neutrinos in the on-source (off-source) regions. In this analysis we subtract the number of off-source events from the onsource observation: N ON − N OF F = N sig on − N sig of f . This number effectively equals N sig on as N sig of f is expected to be significantly smaller than N sig on assuming a true source from the GC halo. The result of this subtraction is directly proportional to σ A V . Systematic uncertainties related to the background should equally effect the onand off-source regions and therefore cancel in the subtraction.
The angular size of the on-and off-source regions is determined to maximize S/ √ B, where S stands for the number of expected signal events and B for the number of background events. This optimization is performed with the same signal and background Monte Carlo used in the combined fit. As the angular resolution of an event depends on the neutrino energy, the optimal size of the onand off-source regions differs between the FC sub-GeV, FC multi-GeV, PC and UP-µ event samples as shown in Table II for the three considered halo profiles. The number N ON − N OF F is obtained for each subcategory of events and similar numbers of events are observed in both the on-and off-source regions for all classes. Figure 10 shows the asymmetry, A = (N ON − N OF F )/(N ON + N OF F ), for the benchmark NFW model and all event categories.
As no asymmetry in the event rate is observed, we constrain σ A V by introducing information on the halo model, the mass of the DM particles and the annihilation channel. Limits are calculated for M χ in the range from 1 GeV to 10 TeV for the bb, µ + µ − , W + W − and νν annihilation channels as shown in Fig. 11. Differences between the NFW model and the Moore and Kravtsov models can reach an order of magnitude, due to large discrepancies in the predicted densities of the DM particles in the inner parts of the galaxy. Figure 12 compares the results of this analysis with the constraints from the combined fit presented in Section V. The combined fit analysis yields limits that are roughly one order of magnitude stronger than the on-source offsource approach for bb, µ + µ − and W + W − . However,  Upper limits at 90% C.L. on the DM selfannihilation cross section σAV as a function of Mχ for the combined fit (solid) and on-source off-source (dashed) analyses for the bb (blue), W + W − (maroon), µ + µ − (purple) and νν (orange) annihilation channels assuming NFW profile.
for low WIMP masses, the obtained limits are of similar strength, as the signal is mostly present in the sub-GeV samples. The pointing resolution of these samples is poorer than their higher energy counterparts which reduces their contribution to the sensitivity of the combined fit.
For the νν mode the limits in both analyses are similar due to the fact that in the on-source off-source method the number N ON − N OF F used in the limit calculation is based only on the sample with the expected DM contribution for the assumed mass. This effectively resembles the situation in the combined fit where the νν signal also appears in only a limited number of samples. In contrast, the limit calculation for bb, µ + µ − and W + W − is based on all samples together since the signal is broadly distributed for these channels.

VII. SUMMARY
The analyses presented above show no excess of DMinduced neutrinos from the galactic center or its halo; the data are consistent with the expectation from atmospheric neutrino backgrounds. The strongest exclusion is obtained for WIMP masses below several tens of GeV in the bb, νν, µ + µ − channels in the combined fit. Limits on σ A V reach as low as 9.6 × 10 −23 cm 3 s −1 for 5 GeV WIMPs in the bb channel and 1.2 × 10 −24 cm 3 s −1 for 1 GeV WIMPs in the νν channel. These are the strongest indirect limits among similar WIMP-induced neutrino searches to date. Strong constraints on GeV-scale M χ reflect the long exposure and the discrimination power of multiple subsamples reconstructed using the Super-Kamiokande detector. A complementary and completely data-driven result which compared event rates from the angular region near the GC and one of the same size but located in the opposite direction also found no evidence of an excess beyond the background expectation.