Searching for Dark Matter Annihilation with IceCube and P-ONE

We present a new search for weakly interacting massive particles utilizing ten years of public IceCube data, setting more stringent bounds than previous IceCube analysis on massive dark matter to neutrino annihilation. We also predict the future potential of the new neutrino observatory, P-ONE, showing that it may even exceed the sensitivities of gamma-ray searches by about 1-2 orders of magnitude in 1-10 TeV regions. This analysis considers the diffuse dark matter self-annihilation to neutrinos via direct and indirect channels, from the galactic dark matter halo and extra-galactic sources. We also predict that P-ONE will be capable of pushing these bounds further than IceCube, even reaching the thermal relic abundance utilizing a galactic center search for extended run-time.

Here we consider Majorana weakly interacting massive particles (WIMPs) as dark matter candidates.In the case of dark matter self-annihilation to neutrino pairs, the spectrum would possess a distinct shape: a peak at the DM rest mass.This feature differs fundamentally from the measured diffuse astrophysical neutrino power-law spectrum [8,19], as well as the atmospheric background.Due to this, we perform energy-binned likelihood analyses, searching for these exotic energy distributions.
In standard WIMP freeze-out scenarios, DM particles are assumed to be in thermal equilibrium, with the expanding universe, the dark sector decouples due to the expansion rate increasing greater than the interaction rate, ceasing DM production and self-interaction leading to the relic density [20][21][22].In this scenario, to account for current DM observations, the thermally averaged annihilation cross section would then be ⟨σν⟩ = 3×10 −26 cm 3 s −1 [20].Thus, should DM not be found, the searches aim to push the constraints on the cross-section below this value, excluding these WIMP scenarios.Note that for DM with mass, m χ > 10 GeV the thermal relic abundance constraint on the cross-section shows minor mass dependence [23].
We follow [23] when calculating the direct neutrino annihilation channel.In this case, the WIMP spin plays an insignificant role and hence can be neglected [19].For indirect channels, we rely on the simulations performed in [24] for their neutrino production spectrum.The data released in [24] include various intermediate particle states, such as W s, bs, and τ s, and their consequent decay products, e.g.neutrinos, electrons, etc.In section 2 and section 3, we give a thorough discussion of the signal and background modeling for P-ONE [1] and IceCube [25].The final limits on the thermally averaged cross-section are shown in the section 4.

Galactic Contribution
In the case of direct annihilation, the galactic dark matter halo contributes trivially due to the negligible redshift.This means the flux spectrum arriving at Earth is almost identical to the spectrum at the production sites, given by Equation 2.1, assuming an equal neutrino flavor decomposition 1:1:1 after long-distance propagation ⟨σν⟩ is the thermally averaged self-annihilation cross-section.We have used κ = 2 for Majorana DM with mass m χ for the DM mass.The factor of 1/3 is due to equal distribution among neutrino flavors.dN/dE is the number spectrum of the neutrinos.Equation 2.2 depicts the number spectrum of neutrinos produced via DM direct annihilation to neutrinos that have a typical delta peak 2) The spectrum shape varies depending on the annihilation channels.DM annihilation to W-boson pairs, bquark pairs, and τ -lepton pairs are conventional choices.These heavy Standard Model (SM) particles can again decay into stable SM particles such as electrons, gamma-rays, or neutrinos.Modeling of the branching ratios is required to include different channels.In [24], thorough modeling has been performed using PYTHIA [29] and HERWIG [30] for the DM mass range from 100 GeV to 100 TeV.
For the direct annihilation to neutrinos, we neglect the Electroweak (EW) corrections.These would generate tails for the annihilation spectrum, increasing the yield of low-energy neutrinos and broadening the peak.As shown in [13,26,[31][32][33][34], the broadening of the peak is less than 10%, which is well below the typical energy resolution for track-like events in IceCube [35].These corrections would be significant in the case of cascade-like events, where the energy reconstruction is better.Of note, is that these corrections are relevant when considering other final states, such as γ, e + , and p, since these corrections induce a significant spectrum [33].Contrary to the direct production channel, for the indirect channels, the EW corrections have been included in the spectrum from PPPC4DM [24].In Appendix A, we show the spectra for different annihilation channels.The J in Equation 2.1, is the 'J-factor', which stands for the three-dimensional integration of the host galaxy's DM density ρ χ over a solid angle Ω, within the line of sight l.o.s.
In the case of P-ONE, we use three distinct regions.The J-factor values are listed in Table 1.In Table 1 we give three J-factors: s-wave, p-wave, and d-wave.The s-wave corresponds to the thermally averaged cross section independent of velocity, p-wave to [⟨σν⟩ ∝ (v/c) 2 ] and d-wave [⟨σν⟩ ∝ (v/c) 4 ] each denoted as J s , J p , and J d respectively.The p-wave and d-wave contributions are suppressed due to the typical Virial velocity of 100 km/s, hence, they can be neglected.The sensitivity region of the sky for different neutrino detectors affects the value of the integrated J-factor, as presented in [13,26].
We modify the J-factor in IceCube's case, splitting it into northern (up-going) and southern (down-going) sky.We then focus on the up-going events, suppressing most of the atmospheric muon background.While this effectively removes most of this background, it also removes the galactic center from the field of view.Depending on the core or cuspy nature of the dark matter distribution, this can make a significant difference.Here, we use a cusped Navarro-Frenk-White (NFW) profile [36], which leads to a J-factor ∼ 5 times smaller than the All-sky one shown in Table 1 and [13,26].This results in a sensitivity drop of ∼ 2. Note that using a different profile, such as an Einasto [37] or Burkert [38], leads to an approximately three times more or less stringent constraint, respectively.This has been widely studied in various previous works [19,39,40].

Experiment
J s /10 0.13 0.12 0.18 Table 1.J-factors calculated using the NFW density profile [36], for IceCube (All-sky) and P-ONE [26].The units of these J-factors are GeV 2 cm −5 sr.θ is the zenith angle.The columns are J-factors for the s-wave, p-wave, and d-wave contributions.
For our analysis, we neglect the p-wave and d-wave contributions, due to their strong suppression.
For the calculation of these J-factors, we assume the sun is located at a distance of R 0 = 8.127 kpc from the galactic center (GC), as determined by [11].The dark matter density used to calculate the J-factor in Equation 2.3 is parameterized as (2.4) Here, r is the distance from GC.We use the best-fit values from [41], with a local density of ρ 0 = 0.4 GeVcm −3 , a slope parameter γ = 1.2, and a scale density ρ s at scale radius r s = 20 kpc.Inverting Equation 2.4 and setting ρ χ (R 0 ) ≡ ρ 0 we obtain ρ s .In Figure 1 we show the differential neutrino flux at the earth, produced by DM annihilation with various masses from the galactic halo.There, we adopt ⟨σν⟩ = 1 × 10 −23 cm 3 s −1 and consider the neutrino production via DM annihilation to τ -pairs.The production spectra are taken from PPPC4DM [24].  1, are used.The differential fluxes are calculated with the spectrum data from PPPC4DM [24] for DM annihilation to neutrino via the indirect τ -channel.

Extra-Galactic Contribution
There are several approximation models for Equation B.3.Here, we use the one described in [27].We give a brief discussion of the different approximations in Appendix B and Appendix C along with a few intermediate results of our calculations.In Figure 2, we show the flux results for two different DM indirect annihilation channels to neutrinos, τ -lepton(upper) and b-quark(bottom).Here we have used ⟨σν⟩ = 3 × 10 −26 cm 3 s −1 for comparing the flux results with Fig. 5 in reference [42].There is a factor of 10 difference between our results and the [42] flux, which can be explained due to differences in the approximation methods.Additionally, [42] mentioned they were using redshifts up to z = 49 in their simulations, which differ from our calculations.Even using the higher fluxes from [42], does not change the final CL limits significantly.For the direct annihilation channel to neutrinos, [26] showed that the extra-galactic flux would, at maximum contribute 10% to the final CL limit.We consider the 10% as an upper bound on contributions by extra-galactic sources.Although the possible extra-galactic contribution could vary between 0.0001% to 10% of the galactic one, it would at best minimally improve the final constraint value on the cross-section.Therefore, we neglect the extra-galactic component in the following analysis.
More analyses were done on extra-galactic contribution from the halo boost factor and halo substructures in [43] and [44].According to [43], the total sub-halo contributions are dominant compared to those of smooth halo structures in the high mass regime, approximately above 10 4 GeV masses.In such scenarios, the contributions from extra-galactic sources boosted by substructures could be comparable with the galactic contribution.Future directional searches targeting nearby galaxies would be ideal for probing such substructure contributions further.

Effective Areas and Background
This section is dedicated to the effective areas of the IceCube and P-ONE observatories and the corresponding background fluxes.We describe the flux-to-count conversion with the help of effective areas since the later statistical analysis requires an event rate prediction.The published IceCube effective areas [25] are binned with a zenith angle grid between 0 • and 180 • .In comparison, the simulated effective areas for P-ONE are differentiated into trimesters in the sky between 0 • to 180 • ( −1 < cos θ < 1).In Figure 3, we show the effective areas for both P-ONE and IceCube.The solid line style corresponds to −1 < cos θ < −0.5, and the dashed-dotted to −0.5 < cos θ < 0.5.Note that the effective area line for 0.5 < cos θ < 1 would overlap with the −0.5 < cos θ < 0.5 line.In addition to the effective area, in IceCube's case, we use the provided mixing matrices, to move from true neutrino energies to smeared energies.For P-ONE, we use the IceCube public energy smearing to approximate the energy reconstruction.The model is described at the end of this section.Observed and projected background counts are used for the likelihood analysis of the signal.To model the atmospheric background for P-ONE, we use MCEq [45], with the primary model H4a [46] and interaction model SYBILL2.3c[47].Since the published effective areas are designed for track-like events, we will exclusively use ν µ charged-current events.
Along with the atmospheric background, we have to include the astrophysical background for both detectors.We assume a power-law spectrum can model the astrophysical flux where we set the parameters to the best-fit values measured by IceCube, with a spectral index of 2.53±0.07 and ϕ 0 = 1.66 +0.25 −0.27 × 10 −18 GeV cm 2 s sr −1 [48].The differential fluxes for the background along with the signal fluxes for the galactic (Equation 2.1) signals are then convolved with the effective areas of the neutrino detectors to produce the signal and background counts represented as N events with t run (run time of the detector) in Equation 3.2 Using this energy-smearing procedure, we obtain the expected signal counts, depending on the dark matter mass and thermally averaged cross-section.Details on the energy smearing are given in Appendix D. In Figure 4 we compare ten years of observed IceCube events [25] with our predicted dark matter signals.Similarly, we show the expected background and dark matter to neutrino signal events after ten years for P-ONE in Figure 5.There we give the results for an IceCube-like energy smearing (top) and a perfect detector (bottom).The energy smearing procedure is done with a log-normal distribution of events.This procedure is briefly explained in Appendix D. The parameters used for the smearing are similar to the IceCube parameterization depicted in [35].In Figure 5, the smeared flux curves (upper graph) show changes in the energy distribution shape for different mass regimes.This is due to the different parameters according to their true energy region (resonant peaks with the DM mass) as discussed in Appendix D.

Analysis and Results
In this section, we perform a log-likelihood ratio test on IceCube data and predicted events for P-ONE.Both of these analyses were done with the NFW profile and combining the signal fluxes as discussed in section 2 with galactic contributions.We define the binned likelihood function for the signal hypothesis, H 1 , to be The binning used for the background and signal distributions are the same.Here, P is a Poisson distribution, i runs over the energy bins, d i is the measured event count taken from data and µ s,i is the expected event count given the DM parameters θ. µ atmos,i and µ astro,i are the expected atmospheric and astrophysical events depending on nuisance parameters η 1 and η 2 .For the astrophysical flux, η 2 is set to reflect the uncertainties on the parameters given in Equation 3.1, while for the atmospheric background, we assume a flat 20% uncertainty.A comparison between the CR, astrophysical flux uncertainties, and η i , is shown in Appendix F. We define the null hypothesis, H 0 , to have zero signal events (µ i,s (θ) = 0).We then define the test statistic q µ as the log-likelihood ratio Since this is a diffuse analysis, we can disregard the right ascension of the events and create a background model by scrambling data in the right ascension.Note, that this reduces the sensitivity of this analysis since we remove The solid red curve represents IceCube's high energy sensitivity estimate from [26] (red).The solid dark green and yellow curves are the current bounds set by KM3NeT [50] and ANTARES [51].We also show light green curve for Bergmez et.al.
any possible benefits from potential DM overabundance in the galactic center.To create the background model, we scrambled 10 years of IceCube data 10 6 times and constructed the mean binned in energy.Then we fit our atmospheric and astrophysical neutrino fluxes to the mean by scaling our calculations, see section 3. The scaling values are 1.1 for the atmospheric flux and 0.98 for the astrophysical flux.We then add the previously mentioned nuisance parameters η 1 and η 2 as floating normalizations to these best fits.We then follow [49] and marginalize over these parameters, to minimize the resulting limits.This is to account for the signal possibly contributing to the background.
In Figure 6, we show the resulting confidence level (C.L.) limits for direct annihilation of DM pairs to ν µ νµ using ten years of public IceCube data.The red solid curve represents the 95% C.L. sensitivity predicted in [26] for the high m χ region with a Background Agnostic method.The green line shows the current bound set by KM3NeT [50] using a fraction of the planned detector.In yellow, we show the bounds set by ANTARES [51].In Appendix G, we compare these results to the expected ones, as well as a comparison of best-fit and injected background and signal events.
The strong bound set by ANTARES is due to its location and the analysis method.While far smaller than IceCube, the galactic center lies in ANTARES' most sensitive region.This allows ANTARES to perform a directional study of the galactic center, boosting its sensitivity.[52] shows that the location leads to approximately a factor of 20 sensitivity penalty for IceCube when performing similar analyses.
In [39] the potential of a full KM3NeT-like [53] detector is shown calculated with an Angular Power Spectrum (APS) method assuming a run time of ten years.This would push the KM3NeT line below the ANTARES line.
We show the resulting predicted sensitivities for P-ONE in Figure 7.The stark change in shape between 1 and 10 TeV mass is due to the energy-smearing model applied here.This change is especially relevant in the direct annihilation channel to neutrinos due to the expected peak at the mass resonance.If one removes our energy reconstruction assumptions, the P-ONE sensitivities will improve greatly, following a similar shape as the red line.The P-ONE sensitivity could be improved by performing a spatial analysis similar to [39].Purely from a scaling perspective, when comparing the effective areas, we would expect the P-ONE sensitivity to improve further by a factor between five and ten when utilizing spatial information.
In Appendix E, we show example test statistic distributions and the resulting confidence limits for P-ONE and IceCube for the scanned model parameters, m χ and ⟨σν⟩.
We now analyze DM pair annihilation to W-bosons and consequent decay to muon neutrinos, depicted in Figure 8.This analysis is done for the P-ONE detector with the numerical data from [24] assuming equilateral distribution amongst the flavors.The low energy cutoff was lowered to 500 GeV for this analysis.Our new limits (in the case of IceCube) are in energies above 1 TeV more stringent than previous analyses, e.g., IC86 by IceCube [54].At the same time, we predict P-ONE's sensitivity using ten years of data to be even greater.
Similarly, we have analyzed the τ -lepton channel for IceCube and P-ONE, with an NFW profile for ν µ νµ pair production.Figure 9 shows the improvement between our results and the previous IC86 Galactic halo with all sky cascade signal [54] study.In both annihilation scenarios, we expected P-ONE to outperform IceCube by an order of magnitude.The "jump" in the 10 3 to 10 3.5 energy region which one can see in Figure 7, is still present in Figure 8 and Figure 9.However, the jump is not as intense as it is for the direct channel since the neutrino spectrum is not a mass resonant peak rather, it is a distribution over the energy up to the dark matter mass discussed in Figure 9 shows that the projected P-ONE sensitivity to dark matter annihilation will exceed the constraints set by Fermi+MAGIC γ-ray limits for 15 dwarf satellite galaxies [55] in the high-mass region as well as exceeds that of predicted by ANTARES collaboration for 4532 days [51].[26], the green curve represents the results from KM3NeT [50], and the yellow those from ANTARES [51].We also show light green curve for Bergmez et.al.
[39](KM3NeT-like).The constraint linear shows a jump in the 10 3 GeV to 10 3.7 GeV energy region.This is due to the log-normal distribution used in our P-ONE energy-smearing model (see Appendix D).
Figure 8. Thermally-averaged WIMP annihilation cross section as a function of the WIMP mass for W + W − annihilation channels.The limits for IceCube (pink band) and estimated sensitivity P-ONE (purple band) are shown.Two W + W − results from ANTARES(4532 days) [51] with yellow curve and KM3NeT/ARCA6 [50] are also presented as a comparison with dark green.

Conclusion
We performed dark matter annihilation searches on ten years of public IceCube data, setting the most stringent constraints on DM self-annihilation to neutrinos in the high-mass regime.Compared to previous analyses, the constraints set here also show almost one order of magnitude improvement to previous neutrino studies for the galactic halo, galactic center, and extra-galactic diffused sources in both the direct and indirect annihilation channels.We also modeled sensitivities for a new proposed neutrino telescope, P-ONE.These show even greater potential than the constraints derived in this work for IceCube and can compete with constraints set by Fermi-LAT gammaray experiments.This indicates that P-ONE will play an important role in future DM searches, especially in the 10-100 TeV range.We expect its detection potential could be pushed further toward the thermal relic abundance when performing an analysis similar to [39], which would include spatial information, unlike the diffuse analysis presented here.Directional information may become especially relevant when analyzing individual galaxies and studying sub-halo contributions [56] to the neutrino flux.

B Extra-Galactic Contributions
The extra-galactic flux is particularly interesting and rarely studied since the neutrinos produced in the extragalactic sources would suffer non-negligible redshift effects.The neutrino flux from extra-galactic sources is calculated by Eq.Equation B.1 for DM direct annihilation to the neutrino pair channel Where is the time-dependent Hubble parameter, ρ c is the critical density of the Universe, and Ω m , Ω r , and Ω Λ are, respectively, the fractions of ρ c made up of matter, radiation, and dark energy [58].The number spectrum of neutrino pair production is given in Equation B.2 which is similar to Equation 2.2.However, the redshift needs to be considered, which results in an energy transformation.The energy parameterization at Earth is E = E ′ (1 + z), with the energy at the source E ′ and the redshift, z, of the extra-galactic source.Similarly, the numerical spectra from [24] also require the transformation factor of 1/(1 + z) for each energy bin.G(z) is the halo boost parameter at redshift z describing the clustering effect of matter in a galaxy's halo, given by With the halo boost, the DM annihilation rate has been parameterized with respect to redshift z and halo mass M .dn/dM describes the number distribution of halo masses and is strongly related to the halo mass function (HMF) [28].A detailed discussion about its calculation is provided in Appendix C. The selection of minimum halo mass affects the uncertainties of the final integral result.The smaller halos are more concentrated and contribute more to the neutrino flux, so choosing the lower integration bound is important.In this calculation, we set M min = 10 −3 M ⊙ as a conservative lower limit [60], [61].

C Halo Annihilation Boost-factor
The halo boost factor G(z) depends on the halo mass function, which itself has a dependence on the variance of linear field density σ[59] Equation C.1 defines the σ.The z dependence of σ is from the growth factor D(z) [59].The W (kR) is the top-hat filter function and P(k) is the power spectrum.These two parts have only M dependence but not z.Therefore, we The results from [72] and [27] have been compared here.The deviation is hard in the low halo mass region, whereas in the higher mass region, both approximations converge.
can treat M,z-dependence of σ separately.For σ(M, z = 0) we use the parameterization described in the appendix of [27] and as in Equation C.2: ln σ −1 = 0.2506 M 0.07536 − 2.6M 0.001745 .(C.2) The Equation C.2 does not include the z dependence of σ.Its z dependence is introduced with the growth factor D(z).In Figure C.2 we compared the approximation from [27] and from [72] at a wide mass range.At the higher mass, the two approximations converge.However, at the low mass region, there seems a significant divergence.The mass function is non-trivially dependent on σ via function f(σ) as shown in Equation C.3 [59] f (σ) function is called the differential mass function, which has various fitting methods developed over the years, as in [28], [67], and [66], etc.In our analysis, we used the approximation described in [27] because it is comparable with the results from [28], as shown in  The integral with ρ χ (DM halo mass density) can be approximated by its concentration parameter c = r ∆ r as in Equation C.4: The integral has been described as a clumsiness factor in [93] or as an enhancement factor described in [82] as well as in [27].The enhancement factor is dependent on the concentration parameter, c ∆ , and the halo mass, M .The g has been approximated as Equation C.5 [27] g (C.5) The concentration parameter can be bounded by the upper limit of 100, hence, one has a limit on the enhancement factor.There are several extensive models to approximate the concentration parameter described in [72], [79], [83], [84] and [91].The concentration parameter comparison at z=3 throughout the halo mass for various models from [72], [79] and [92].The approximations vary in the lower mass region by a greater margin.However, in the high-energy region, the variation is small.
In our analysis, we have opted for the approximation of the concentration parameter c from [72].The validity ranges over halo mass and redshift for various approximations change drastically, increasing the deviations even more.This results in larger differences in the estimations of halo boost factors and consequently affects the flux estimations.In

D Smearing to approximate the Energy Reconstruction
The energy of an incoming particle (true energy) is converted into the optical module-registered photon.The reconstruction procedure converts the registered photon into an energy distribution.The process is specific to each detector.Since P-ONE is still in its early developing stage, we have used a log-normal distribution for the P-ONE energy reconstruction with the parameters close to the IceCube smearing matrix parameters published in [35].
[TeV] E true < 1 1 ≤ E true ≤ 10 10 < E true µ 0.7 linear spline E true σ 0.45 linear spline 0.35  With the help of Equation D.1 we "reconstruct" the energy distribution of the detected neutrinos.
f  2. We are using three different "Energy Reconstruction" parameters from low to high energy regions, i.e. below 1 TeV between 1-10 TeV and above 10 TeV.This is modeled in such a way to emulate "IceCube-like" energy smearing behavior as parameterized in [35].In the

E Confidence Level Example
In     F Uncertainties due to atmospheric and astrophysical neutrino fluxes In Figure F.11 we compare the uncertainties caused by CR and astrophysical neutrino uncertainties to the normalization uncertainties we introduce with η 1 and η 2 .In the analysis range we consider here, CR and astrophysical uncertainties are smaller than those from η 1 and η 2 .The uncertainties from the CR flux are calculated by injecting different primary cosmic ray models, from [96].We compare these uncertainties to the effects of η 1 and η 2 (gray) used in the analysis here.

G Analysis Tests
In this section, we test the analysis method employed here.In Figure G.12 we compare the injected counts with the obtained best-fit values.In Figure G.13, we compare the expected limit to the one set using IceCube data.

Figure 1 .
Figure 1.Differential neutrino fluxes for the galactic halo for various DM masses between 100 GeV and 100 TeV.Here ⟨σν⟩ = 1 × 10 −23 cm 3 s −1 , and the J-factors for P-ONE, see Table1, are used.The differential fluxes are calculated with the spectrum data from PPPC4DM[24] for DM annihilation to neutrino via the indirect τ -channel.

Figure 2 .
Figure 2. Differential neutrino fluxes of extra-galactic sources via in-direct annihilation of 10TeV-DM pairs into a τ + τ − pair with ⟨σν⟩ = 3 × 10 −26 cm 3 s −1 .The difference in the flux results (dashed) compared to [42] (dotted) is due to uncertainties associated with the approximation methods.The blue line is MCEq simulated atmospheric neutrino background.Here we have parameterized the red-shift grid with the energy of neutrino at production E prod and the mass of DM m χ .

Figure 5 .
Figure 5. Expected counts for the P-ONE detector with ten years run-time for the direct DM to neutrino annihilation channel.The upper graph depicts expected counts with energy smearing.The bottom graph represents an ideal scenario with perfect energy reconstruction at the detector.The smearing parameters applied are similar to the IceCube public data energy smearing matrix from[35].The 10 4 GeV peak seen in the fluxes is due to the parametrization we chose here.

Figure 6 .
Figure 6.IceCube 95% C.L. limits using ten years of public data on the ν µ νµ annihilation channel (black).The dark pink and light pink bands represent uncertainties due to the dark matter density and further model-related uncertainties.The solid red curve represents IceCube's high energy sensitivity estimate from[26] (red).The solid dark green and yellow curves are the current bounds set by KM3NeT[50] and ANTARES[51].We also show light green curve for Bergmez et.al.[39](KM3NeT-like),thermal relic abundance[23] (black dots), and the unitarity bound (black dashed).

Figure 7 .
Figure 7. Sensitivity estimation for P-ONE.The solid black line shows the 95% C.L. limit.The 90% (dark blue) and 68% (light blue) bands represent uncertainties due to the dark matter density.The sensitivity is estimated with ten years of run-time, including energy smearing.The solid red curve represents the previous results from[26], the green curve represents the results from KM3NeT[50], and the yellow those from ANTARES[51].We also show light green curve for Bergmez et.al.[39](KM3NeT-like).The constraint linear shows a jump in the 10 3 GeV to 10 3.7 GeV energy region.This is due to the log-normal distribution used in our P-ONE energy-smearing model (see Appendix D).

Figure C. 2 .
Figure C.2. Dependence of the linear density field σ, on the Mass, M. The results from[72] and[27] have been compared here.The deviation is hard in the low halo mass region, whereas in the higher mass region, both approximations converge.

Figure C. 3 .
The ∆ symbol means the density described in terms of critical density ρ c multiplied with a constant ∆.

Figure C. 3 .
Figure C.3.The fitted functions for the differential mass function f compared to σ at z=0 for both ∆ parameters.The results from Watson et al.[28] and Lopez et al.[27] are within a difference of one order of magnitude.

Figure C. 4 .
Figure C.4.The concentration parameter comparison at z=3 throughout the halo mass for various models from[72],[79] and[92].The approximations vary in the lower mass region by a greater margin.However, in the high-energy region, the variation is small.
Figure C.5 we compared the results of the halo boost factor with the approximation from [82] and simulation from [72].The halo boost factor calculated from Figure C.5 will then be used in Equation B.1 to calculate the differential flux.

Figure C. 5 .
Figure C.5.The halo boost factor G with respect to redshift z range (0, 20).A comparison between Prada et al.[72] simulated result and A. Moline et al.[82] simulated result is shown here.Both curves deviate to an even larger throughout the larger redshift range.

Figure D. 6 .
Figure D.6."Energy reconstruction" PDF used for P-ONE.It is a combination of log-normal distributions, with parameters shown in Table2.The color gradient of gray pixels presents the possibility of 0 (light) to 1 (dark) for a true energy to be reconstructed as the respective smeared energy.The red line indicates the 68% confidence interval.

2 , (D. 1 )
The PDF describes the smeared energy(E recon ) distribution probability of each True energy (E true ) bin throughout the energy grid.The Figure D.6 shows the smeared energy distribution (E Recon ) with respect to the true energy (E true ) with the parameters described in Table

Figure D. 6 ,
we see the log-normal distribution's behavior changing drastically from 10 3 GeV.This change results in the peak observed in the P-ONE limits between the energy band of 1 to 10 TeV.The decline from the peak also occurs due to a change in log-normal distribution at 10 TeV.

Figure D. 7 .
Figure D.7.Example test statistic distributions for the background and signal hypothesis for IceCube.The red lines denote the 90% (dashed) and 10% (solid) quantiles for the signal and background hypothesis, respectively.
Figure E.8 and Figure D.7, we show an example distribution of the test statistic for P-ONE and IceCube.

Figure E. 8 .
Figure E.8.Example test statistic distributions for the background and signal hypothesis for P-ONE.The red lines denote the 90% (dashed) and 10% (solid) quantiles for the signal and background hypothesis, respectively.

Figure E. 9 .
Figure E.9.The confidence level plot for P-ONE.We have drawn the 95% confidence level contour for reference.

Figure E. 10 .
Figure E.10.The confidence level plot for IceCube.We have drawn the 95% confidence level contour for reference.

Figure E. 9
Figure E.9 and Figure E.10 show the 95 % confidence levels (solid line) for P-ONE and IceCube respectively.This illustrates the distribution of confidence levels.

Figure F. 11 .
Figure F.11.Here we show the uncertainties due to the cosmic ray flux (green) and the astrophysical neutrino flux (orange).We compare these uncertainties to the effects of η 1 and η 2 (gray) used in the analysis here.

Figure G. 12 .
Figure G.12. Comparing the injected signal and background (red) to the best-fit values for the background (black dashed) and signal (green dashed).

Figure G. 13 .
Figure G.13.A comparison of the expected limit (dashed) to the one set using data (solid).

Table 2 .
Reconstruction parameters for P-ONE for different true energies given in TeV.The linear spline describes the linear extrapolation between the 1 TeV to 10 TeV boundaries.