Solar Atmospheric Neutrinos and the Sensitivity Floor for Solar Dark Matter Annihilation Searches

TL;DR

This paper calculates the flux of solar atmospheric neutrinos, identifying a fundamental background that limits the sensitivity of indirect solar dark matter searches, and suggests this flux could be observed with current and future neutrino telescopes.

Contribution

It provides a detailed new calculation of solar atmospheric neutrino flux and establishes the sensitivity floor for indirect WIMP annihilation searches.

Findings

Sensitivity floor is less than an order of magnitude beyond current limits.

High-energy solar atmospheric neutrinos may be detectable.

The flux calculation includes systematic uncertainties in shower evolution.

Abstract

Cosmic rays interacting in the solar atmosphere produce showers that result in a flux of high-energy neutrinos from the Sun. These form an irreducible background to indirect solar WIMP co-annihilation searches, which look for heavy dark matter particles annihilating into final states containing neutrinos in the Solar core. This background will eventually create a sensitivity floor for indirect WIMP co-annihilation searches analogous to that imposed by low-energy solar neutrino interactions for direct dark matter detection experiments. We present a new calculation of the flux of solar atmospheric neutrinos with a detailed treatment of systematic uncertainties inherent in solar atmospheric shower evolution, and we use this to derive the sensitivity floor for indirect solar WIMP annihilation analyses. We find that the floor lies less than one order of magnitude beyond the present experimental limits on spin-dependent WIMP-proton cross sections for some mass points, and that the high-energy solar atmospheric neutrino flux may be observable with running and future neutrino telescopes.

Tables

Table 1. Table 1 : The expected event rates in each of the detectors considered in this work. The rows marked with a * are for detector configurations where only the averaged effective area is provided, and so the neutrino and antineutrino effective areas set equal. Because the better estimate is provided where the separate effective areas are available, the IceCube+PINGU rate is constructed using the IC79 and PINGU effective areas.

Experiment	Expected ${\nu }_{\mu }$ rate $R$ (evts / yr)	Expected ${\overline{\nu }}_{\mu }$ rate $R$ (evts / yr)
IceCube (IC79)
IceCube (IC86)*
ANTARES*
IceCube+PINGU
KM3NeT*

Equations

$$\frac{d{N}_{\nu }}{d{E}_{\nu }}=\frac{{\Gamma }_{\odot }^{ann}}{4\pi d^2}\sum \frac{d{N}_f}{d{E}_{\nu }},$$

$${\Gamma }_{\odot }^{ann}=\frac{{\Gamma }_{\odot }^{capt}}{2}{\tanh }^2(t_0/{\tau }_A)$$

$${\Gamma }^{capt}=1.3\times 10^{25}{\mathrm{s}}^{-1}\left(\frac{{\rho }_{\chi }}{0.3\mathrm{GeV}\mathrm{c}{\mathrm{m}}^{-3}}\right)\left(\frac{270km/s}{v_{sun}}\right)\left(\frac{\mathrm{GeV}}{m_{\chi }}\right)\left(\frac{{\sigma }_{\chi p}}{10^{40}{\mathrm{cm}}^2}\right)S(\frac{m_{\chi }}{m_H}),$$

$$\tau (E_{\nu };b)=\exp \left[-{\int }_{l(b)}{n}_N(x;b){\sigma }_{\nu N}(E_{\nu })dx\right],$$

$$\Phi (E_{\nu }{)}_{\alpha }=2\int dbb{\Phi }_{b,\alpha }(E_{\nu }){\Omega }_{\odot },$$

$${\bar{\Phi }}_{\alpha }(E_{\nu })=\frac{1}{T}\sum _{\beta }{\int }_0^{T}dt{P}_{osc}^{\beta ,\alpha }(r(t),E_{\nu }){\Phi }_{b,\beta }(E_{\nu }),$$

Full text

1 Introduction: Direct and indirect dark matter searches

The existence of gravitationally interacting dark matter has been well established through a wide variety of cosmological and astrophysical observations [1]. However, the experimental observation of particle dark matter remains elusive [2]. Among techniques used to search for particle dark matter are direct detection [3], accelerator production [4] and indirect astrophysical searches [5, 6]. These techniques are highly complementary, allowing access to different regions of dark matter model space.

Direct detection is based on the principle that, if dark matter particles have interactions beyond gravitational ones with standard model particles, low-background particle detectors may observe the recoils of nuclei scattering from dark matter particles in the galactic halo. The standard paradigm is to assume that dark matter particles are WIMPs (weakly interacting massive particles), and following assumptions about the velocity and density distribution of the Galactic halo of dark matter [7], non-observation of nuclear recoils in such detectors have allowed limits to be placed on the WIMP-nucleus cross section. These limits presently extend to $\sigma=10^{-45}\mathrm{cm^{2}}$ at the most sensitive mass of $40\,\mathrm{GeV}$ [8]. The limits are considerably weaker if the dark matter-nucleus interaction is spin dependent, extending only to $5\times 10^{-40}\mathrm{cm^{2}}$ in a similar mass range [9, 10].

The precision of many direct dark matter detection experiments is now nearing sensitivity to coherent scatters of low-energy solar neutrinos. These neutrinos represent a source of background that is effectively irreducible in traditional detectors. The WIMP scattering cross section at which the dark matter signal would be buried under coherent neutrino scattering events is commonly referred to as a “sensitivity floor” for direct dark matter detectors. Although a WIMP signal could still be observed over well-predicted background models, working in the large-background regime implies a significant impediment to the rate of experimental progress [11], which can make sensitivity improvements in traditional dark matter searches prohibitively slow. On the other hand, it may be still possible to exploit additional observables, including directional detection [12] and annual modulation [13], to move beyond the limitations imposed by the coherent neutrino background, if suitable detectors can be devised.

Indirect detection is an independent and complementary method of searching for dark matter. The technique relies on searching for annihilation products of dark matter particles in astrophysical environments [14, 15], where the local dark matter density is expected to be high enough to facilitate copious WIMP self-annihilation. Potentially detectable annihilation products include neutrinos, gamma rays, positrons, antiprotons and anti-nuclei [1]. A channel that has drawn particular attention is the production of high-energy neutrinos from WIMP self-annihilation in the solar core. The Sun would act as a concentrator of dark matter particles if and only if WIMPs interact sufficiently often with regular matter to effectively transfer their kinetic energy to solar material as they traverse the Sun, causing them to pool near the center. Thus the strength of this signature dependends on both the self-annihilation cross section (calculable in WIMP models) and the WIMP-proton cross section.

Large neutrino detectors, including SuperKamiokande [16, 17], ANTARES [18], and IceCube [19], have searched for fluxes of high-energy neutrinos produced through solar WIMP self-annihilation and set limits on both spin-dependent and spin-independent WIMP-proton couplings. The limits assume annihilation to an intermediate particle that decays into neutrinos, with the strength of the limit depending on the final state and its kinematics. Three channels of particular interest are $\chi\chi\rightarrow b\bar{b}$, $\chi\chi\rightarrow W^{+}W^{-}$, and $\chi\chi\rightarrow\tau^{+}\tau^{-}$. In all cases, the limits on the spin-dependent cross section exceed the strength of those from direct searches. At the time of writing, IceCube holds the strongest published spin-dependent cross section limit, using 3 years of data from the full 86-string detector to set limits that peak at $3\times 10^{-39}\,\mathrm{cm^{2}}$ for $\chi\chi\rightarrow b\bar{b}$, $5\times 10^{-41}\,\mathrm{cm^{2}}$ for $\chi\chi\rightarrow W^{+}W^{-}$, and $2\times 10^{-41}\,\mathrm{cm^{2}}$ for $\chi\chi\rightarrow\tau^{+}\tau^{-}$ modes [20]. IceCube and ANTARES will continue collecting data and improving their limits, and future neutrino telescopes, including PINGU [21] and KM3NeT [22], will also be sensitive to dark matter annihilation signatures.

The experimental technique involves treating the Sun as a point source and searching for an unexpected flux of high-energy neutrinos from that direction–we discuss this signature in more detail in Section 2. The backgrounds in the leading searches today are dominated by the flux of atmospheric neutrinos produced by cosmic ray interactions in the Earths atmosphere. This background is in principle reducible by improving the angular resolution of large neutrino detectors, either by improved analysis methods in existing experiments or by design improvements for future detectors.

A second source of background, presently sub-dominant but effectively irreducible, is the flux of high-energy neutrinos produced by cosmic ray induced showers in the solar atmosphere. These neutrinos have been discussed before, for example in [23, 24, 25, 26, 27]. In this paper we present a new, systematically detailed calculation of this flux111We will refer to our calculation as the Fedynitch-Jones-Argüelles-Wasseige solar flux, or FJAWs, and we discuss its implication as a background to indirect WIMP annihilation searches. We predict the event rates in existing and proposed neutrino telescopes, and calculate the effective “sensitivity floor” imposed by solar atmospheric neutrinos on indirect solar WIMP annihilation searches.

This paper is organized as follows. Section 2 discusses the experimental signature of WIMP self-annihilation and the method of searching for it with neutrino telescopes. Section 3 details our calculation of the solar atmospheric neutrino flux, from primary modeling, through shower simulation, oscillation and propagation effects. In Section 4 we predict the event rates in various running and proposed neutrino telescopes, and in Section 5 we consider the effects of this flux as a background to WIMP searches, deriving the sensitivity floor in terms of spin-dependent cross section for the $\chi\chi\rightarrow b\bar{b}$, $\chi\chi\rightarrow W^{+}W^{-}$, and $\chi\chi\rightarrow\tau^{+}\tau^{-}$ self-annihilation channels. Finally, in Section 6 we present our conclusions.

2 The indirect dark matter signature

2.1 Dark matter annihilation channels

In order to compare the high energy solar atmospheric neutrino flux to the proposed signals from dark matter self-annihilation, we must specify the dark matter self-annihilation rate and neutrino yield. With these elements, the neutrino spectrum is given by

[TABLE]

where $d$ is the mean Sun-Earth distance and the neutrino spectrum, $\frac{dN_{f}}{dE_{\nu}}$, for dark matter annihilation after solar propagation is taken from [28]. The sum is over the dark matter annihilation channels $f$. The annihilation rate is related to the dark matter capture rate by [29]

[TABLE]

where $t_{0}=4.5~{}{\rm Gyr}$ is the Sun age and $\tau_{A}$ is the capture and annihilation time-scale. In our calculation, which is analogous to the procedure described in [30], we will assume thermal equilibrium so $\Gamma^{ann}_{\odot}\approx\frac{\Gamma^{capt}_{\odot}}{2}$. The dark matter capture rate, $\Gamma^{capt}_{\odot}\propto\sigma_{\nu p}\rho_{\rm DM}$, as a function of the dark matter mass and the dark matter proton scattering cross section, is calculated using formulae from [31, 32] assuming a spin dependent interaction. We assume a local dark matter density of 0.3 ${\rm GeV}/{\rm cm^{3}}$ with Maxwell-Boltzmann distribution with a Sun local velocity of 270 km/s [33]. Thus, following [32], the capture rate is

[TABLE]

where $m_{H}$ is the hydrogen atom mass and $S$ is kinematic supression factor given in [34].

In this work we consider three different independent annihilation scenarios: ${\rm DM}-{\rm DM}\to b\bar{b}$, ${\rm DM}-{\rm DM}\to W^{+}W^{-}$, and ${\rm DM}-{\rm DM}\to\tau^{+}\tau^{-}$. In the $\tau^{+}\tau^{-}$ and $W^{+}W^{-}$ scenarios, neutrinos are produced promptly by decay, thus producing a hard spectrum. On the other hand, in the $b\bar{b}$ channel the fermion pair hadronizes into $B$ mesons, which at the high solar densities interact before decaying, thus producing a softer spectrum [29, 28].

2.2 The WIMP annihilation signature in neutrino detectors

Neutrinos are ideal candidates to search for WIMP annihilations in celestial bodies such as the Sun or the Earth as their small cross section allows them to escape high density medium without significantly interacting. Since the interaction length becomes smaller than the solar radius for neutrino energies above 1 TeV, the solar WIMP searches are typically performed in an energy range from a few GeV up to a few TeV.

The event selection criteria are typically optimized to extract charged-current (CC) muon (anti-) neutrino interactions. This is because alternative channels, including neutral current and CC electron/tau neutrino interactions suffer from lower angular resolution. Muons from CC muon neutrino interactions produce track-like events, which in this energy range are effectively collinear with the initial neutrino direction. It is therefore possible to pinpoint the direction of the Sun by reconstructing the muon track direction, thus allowing the rejection of the terrestrial atmospheric neutrino background.

A night-search, i.e. when the Sun is below the horizon, comes with an advantage of further reducing this background as the Earth would play the role of a shield for atmospheric muons. Day-searches can also be performed by the use of harder cuts in the event selection to limit the atmospheric muon contamination, leading inevitably to a smaller neutrino acceptance compared to night-searches. Terrestrial atmospheric neutrinos reach the detector in both day and night searches, being rejected using the angular reconstruction described above.

Statistical methods (e.g. likelihood or counting experiment) can then be performed to highlight an excess of neutrino events in the direction of the Sun compared to the rest of the observable sky. A significant excess of events in the angular distribution, sometimes coupled with an energy distribution, would therefore highlight a potential dark matter signal.

To date, no significant excess of events with respect to the expected background has been found. Limits on the solar neutrino flux can therefore be set. Several assumptions are needed in order to express the results in a limit on the WIMP-proton scattering cross section, comparable with direct search results. Among these are the local dark matter density, the velocity distribution of the halo, and the assumption of equilibrium between capture and annihilation. These determine the density of WIMPs acreted in the solar core.

3 The solar atmospheric neutrino flux

Just as cosmic rays impinging on the Earth’s atmosphere produce air showers, leading to hadrons which decay to atmospheric neutrinos [35], cosmic rays interacting in the solar atmosphere also produce a high-energy neutrino flux. This flux was studied in references [25, 26, 27]. The main conclusion was that this flux of neutrinos is small, and therefore unlikely to be useful for detailed study of (for example) neutrino oscillations. However, as the precision of neutrino telescopes improves and increasingly high statistics samples are collected, this flux will assume a new importance as the limiting background to indirect solar WIMP annihilation searches. In this section we present a new calculation of the solar atmospheric flux, incorporating detailed treatments of systematic uncertainties inherent in neutrino production in the solar atmosphere.

The production of neutrinos in solar showers is different from their production in terrestrial air showers in a few key ways. First, the region of solar atmosphere where the majority of production is localized is significantly less dense and further extended than its terrestrial counterpart. This allows for longer decay lengths of high-energy hadrons before they are absorbed through inelastic interactions, reducing the suppression of the high-energy neutrino flux observed in the Earths atmosphere. On the other hand, the solar core is very large and dense relative to the Earth, so more high-energy neutrinos are lost through interactions when propagating across it. Finally, the path lengths in the solar atmosphere are long enough (thousands of kilometers) that high-energy muons decay and produce a sizeable contribution to the solar atmospheric neutrino flux, whereas in terrestrial air showers these would be stopped abruptly in the Earths crust.

To calculate the solar atmospheric neutrino flux, a model of the solar atmospheric density profile is required. This is discussed in Section 3.1. The calculation of the solar atmospheric neutrino flux within this density profile is best performed by solving a system of cascade equations [36]. For this purpose we use the general purpose MCEq hadronic cascade solver [37, 38], which has previously been used for detailed studies of neutrino production in terrestrial air showers. We derive not only a prediction of the solar atmospheric flux production in the solar atmosphere but also a suite of systematic uncertainties on this flux in Section 3.2. The decay of muons which punch through the end of the MCEq geometry is then simulated using a numerical method, yielding a neutrino flux to be propagated across the higher density solar matter.

Having produced the shower flux, oscillation and propagation effects for neutrinos travelling across the Sun must be treated. The oscillation calculation accounts fully for the details of the solar matter potential and its coherent and incoherent effects on neutrino oscillations, using the nuSQuIDS software framework [39], discussed in Section 3.3. After propagation to the Earth, our predicted flux at Earth is presented in Section 3.4.

3.1 Modelling of the solar atmosphere

One of our motivations to improve on the solar neutrino flux calculations presented in past works (e.g. Ingelman and Thunman [25]), was the significant improvements made on solar atmosphere modelling in the past years. Ingelman and Thunman used an exponential density profile fit on the data of [40] for the atmosphere and of [44] for the deeper layers of the Sun. An update of [40] has been produced by Fontenla et al. [41]. As its predecessor, this model is a 1D semi-empirical model of the quiet-Sun. A 3D model (e.g. [45]) including the movements of the convection cells is not needed in the work presented here as the speed of these cells ($\sim$ 1km/s) can be neglected compared to the speed of the particles we consider.

An additional requirement for our solar model is a continuous density profile from the top of the atmosphere to the interior of the Sun. While the interaction region will be located in the low chromosphere and the photosphere, the density and composition of the core will be of importance for the neutrino propagation across the Sun.

A hybrid model based on the following combination allows to fulfill the two requirements described above:

• •

An up-to-date 1D atmosphere model [41],

• •

A brief use of the Harvard-Smithsonian atmosphere model [43], extensively used in high-energy physics, for the deeper layers of the atmosphere,

• •

A convection zone modeled by [46]. This maintains continuity between the deep layers of the atmosphere and the core.

• •

A core based on the BS05_agsop model [42].

The complete hybrid model density profile is shown in Figure 1, along with the delineation of the zones of our calculation: first MCEq air-shower cascade evolution in high solar atmosphere, followed by numerical muon decay simulation in the convection zone, and finally nuSQuIDS neutrino density matrix evolution across the core to the opposite side of the Sun. Figure 1, bottom shows how our hybrid model is assembled in terms of the parts described above, and also how it compares to the previous approximation for the solar density profile used in [25].

In this calculation we follow previous works such as [25] by considering neutrino production in the atmosphere of a non-magnetic Sun. It has been suggested [23] that at the lowest cosmic ray energies, corrections to the neutrino flux due to various atmospheric and heliospheric effects may be non-negligible. Among these effects, the most relevant in the current framework are the propagation of low energy cosmic-ray through large-scale, such as the interplanetary magnetic field or the coronal field, and small-scale magnetic fields taking place in the low atmosphere of the Sun [23]. The large scale fields will reduce the fraction of the cosmic ray flux that are absorbed by the Sun, leading to a smaller neutrino yield from the low energy cosmic rays. The latter changes the cascade evolution through charged particle reflection from solar magnetic flux tubes anchored in the bottom of the photosphere. The cosmic ray particle, or one of its subsequent charged products, can be mirrored and trapped by the magnetic structure, modifying the density profile seen by the cascade. Because of the complexity of accurately modelling the effects of solar magnetic fields, such corrections are not included in the present work. This may introduce some additional uncertainty at the lowest energies (below 100 GeV) though over most of the energy range of our calculation the effects are expected to be sub-dominant to other sources of uncertainty.

3.2 Numerical cascade equation solution in the solar atmosphere

The software package MCEq [37, 38] is a linear, one-dimensional cascade equation solver. We use it to calculate the production neutrino fluxes from hadrons and muons produced in cosmic ray induced showers in the solar atmosphere.

We consider profiles at different cosmic ray impact parameters across the solar sphere, and these trajectories define chords through the solar atmosphere. Density profiles along the chords are calculated at discrete impact parameter values of $0<b<1$, with $b=0$ corresponding to a directly core-crossing trajectory and $b=1$ a glancing trajectory at the solar radius. Our calculation is divided into distinct regions as a function of radial depth. The cascade equation for neutrino production is solved for each chord to predict the total flux crossing a surface approximately 6500 km inside the conventionally defined solar radius (see Figure 1). The section of solar density profile used for the cascade calculation thus starts at a solar density of $10^{-14}\,\mathrm{g\,cm}^{-2}$ and ends at $3\times 10^{-3}\,\mathrm{g\,cm}^{-2}$. Shorter geometries were also tested with negligible difference observed in the final $\pi$- and $K$-induced neutrino flux, validating that neutrino production by hadrons is complete by this radius and in our energy range. The most important region for neutrino production in the solar atmosphere lies at densities of $1\times 10^{-4}\mathrm{g\,cm^{-3}}$. This reflects that distance scales of the solar atmosphere are sufficiently long that by this point there has been ample column depth for the initial flux of primaries to have fully interacted, and for the resultant boosted hadrons to decay.

An exception to the above principle is the large flux of highly boosted muons, many of which remain at the end of the MCEq simulation region, and have to be decayed subsequently in the next region. After the MCEq calculation is complete, the flux of muons penetrating out of the cascade equation geometry is extracted and fed to a second post-decay simulation step, which propagates them through a further radial 28000 km of solar atmosphere, to a region with density around $10^{-2}\,\mathrm{g\,cm^{-3}}$. This region has sufficient integrated density that the total energy loss of the most energetic muons generated in the MCEq stage moves them below the lowest energy bin. In each bin of the calculated muon flux, energy losses are applied according to tables from [47] in steps, and at each step a fraction of the muon flux is converted into a neutrino population via three-body decay according to the polynomial parametrization from  [48]. The propagation is evaluated in distance steps corresponding to whichever is the smallest between 1% of the decay length of the boosted muon and 1% of the energy loss length, and repeated until the muon energy falls below the threshold of 10 GeV relevant to this study. The effect of the additional post-decayed neutrino flux is most relevant for the electron flavor neutrinos, which are dominated over most of the energy range of our study by production from muon decay. At energies below 100 GeV, this post-decayed flux makes a small ($<3$%) contribution, since at low energies most muons have decayed in the MCEq volume. As the energy increases and the muons become more boosted, the additional flux becomes a more substantial component, contributing an additional 40% to the electron neutrino flux and 9% to the muon neutrino flux at 1 TeV. At still higher energies, the relative flux of the re-decayed muons falls again as muons become subdominant to hadrons like $\pi$, $K$, and eventually charmed hadrons.

The predicted flux of $\nu_{e}$, $\nu_{\mu}$, and $\nu_{\tau}$ in the solar atmosphere after both MCEq and post-decay steps is shown in Figure 2. For illustration, the results at impact parameter $b=0$ are shown, although the impact parameter dependence of the flux is weak at this stage. To obtain the total flux from the whole Sun, the flux at each $b$ must be propagated across the solar body, which introduces much larger b dependencies that must be accounted for before integrating. This step will be discussed in Section 3.3.

In Figure 2, the total flux, shown as a black line, is broken into neutrino and antineutrino components, and separated by parent particle. The electron neutrino / antineutrino component is dominated by muon decay until the highest energies where this becomes subleading to kaon decays. The muon neutrinos / antineutrinos are dominated by muon decays at the lowest energies, then briefly by pions at around 1 TeV, and finally by kaon decay at the highest energies. Neutrinos from charmed hadrons, often called “prompt” neutrinos, are sub-dominant everywhere for the muon and electron flavors. However, since neither pions, kaons or muons have sufficient mass to produce a $\tau$ lepton, this charmed contribution is the only source of tau neutrinos, as can be seen in Figure 2, right.

The inputs to the calculation described above are:

1. I

The solar density profile, as discussed in Section 3.1. 2. II

The composition / hadronic interaction model, which are supplied together as MCEq yield and decay tables, used in solving the cascade equation. For the calculation shown in Figure 2 we used the SIBYLL-2.3 [49, 50] model assuming a purely protonic atmosphere and MRS prompt model [51]. 3. III

The primary cosmic ray spectrum, which is injected as the initial condition. For the calculation shown in Figure 2 we used the cosmic ray flux model H4a from [52], also called the Hillas-Gaisser H4a model.

To assess the scale of systematic uncertainty on our prediction we explore a variety of primary models and composition / hadronic models. In addition to the H4a parametrization of [52], the primary models studied include the “Gaisser-Honda” model, a widely used five-mass-group model tuned to balloon data [35]; the “combined Gaisser-Honda with H4a” model (cHGP) [53]; the “Thunman model”, a broken power law parametrization constructed in [54]; the “Polygonato model” [55], constructed to explain the knee of the cosmic ray spectrum in terms of sequential cut-offs of different mass components; and the “Zatsepin-Sokolskaya” model [56], constructed in terms of populations of cosmic rays accelerated by discrete types of nova and supernova events. The effects of varying the primary model on each of the flavor fluxes in the lower solar atmosphere, relative to the prediction in Figure 2 is shown in the top row of Figure 3.

The set of hadronic models we consider are SIBYLL-2.3 [49, 50]; the previous SIBYLL version, SIBYLL-2.1 [57]; the EPOS-LHC [58] model; and QGSJET-II-04 [59, 60]. In each case we couple the model for conventional hadron production to the MRS prompt model [51]. With the exception of SIBYLL-2.3, these models do not have MCEq tables prepared for purely protonic environments. Hence to estimate the effect of changing the hadronic model we take the ratio of the flux predictions calculated using an air atmosphere of the same nucleon density profile as the Sun, $\phi_{air}^{model}(E)$. We then correct this flux using the SIBYLL-2.3 flux predicted in a protonic atmosphere according to the approximation $\phi_{p}^{model}(E)\approx\phi_{air}^{model}(E)/\phi_{air}^{SIBYLL2.3}(E)\times\phi_{p}^{SIBYLL2.3}(E)$. The effects of varying the hadronic model is shown in the central panel of Figure 3.

Although we use the MRS model in all cases above, the SIBYLL models also have a built-in charm production parametrization available. We compare the effects of switching to this charm model in Figure 4. This comparison is made using SIBYLL-2.3 coupled to either its intrinsic charm model or to the MRS model, in both cases using a protonic atmosphere. The charm model uncertainties have a drastic effect on the prompt flux prediction, though this is a small effect for the $\nu_{\mu}$ and $\nu_{e}$ fluxes where the prompt contribution is strongly sub-dominant. On the other hand, $\nu_{\tau}$ are produced exclusively through decays of charmed hadrons in this energy range, so the charm model effects are the leading source of uncertainty at production for $\nu_{\tau}$ in the solar atmosphere. The scale of the variation induced by changing the charm model on the prompt fluxes is shown in Figure 4.

Finally we investigate the effects of changing the atmospheric density and composition. We find that extreme outliers are required to produce flux distortions comparable with the hadronic and primary model uncertainties. Figure 3, bottom shows a comparison between an assumed protonic ($\langle A\rangle=1$) or air ($\langle A\rangle=14.5$) composition for the solar atmosphere. With this highly unrealistic composition model, an order 10% distortion is observed. Similarly, switching to the much less accurate Ingelman Thunman atmospheric density model allows us to test the effect of a gross distortion in the solar density profile [25]. Again, with such a highly exaggerated distortion, the flux variations are of comparable scale to those obtained from switching between primary and hadronic models. Because the true uncertainties on the composition and density profile are much smaller than those investigated here, which generate effects of comparable scale to the primary and hadronic model variations, we conclude that these uncertainties can be neglected in our final uncertainty budget.

We compile the effects of all variations described above into a total uncertainty band for each neutrino and antineutrino flavor, shown in Figure 5. The band is drawn between a flux evaluated at the extremal interaction model, charm model, and primary model combination at each energy. We conclude that the $\nu_{\mu}$ and $\nu_{e}$ fluxes in both neutrinos and antineutrinos are predictable within a factor of approximately 2 below 1 TeV, with slightly larger uncertainties at higher energies. The $\nu_{\tau}$ flux is dominated by charm production uncertainties and is much less predicable than $\nu_{\mu}$ or $\nu_{e}$. While the $\nu_{\tau}$ flux at production is difficult to estimate, as will be shown in the next section, oscillations during solar propagation strongly mix the $\nu_{\mu}$ and $\nu_{\tau}$ flavors, resulting in the charm model uncertainty not being significant in our final results.

3.3 Oscillation and absorption effects

The neutrino flux calculated in Section 3.2 at the outer most layers of the Sun needs to be propagated through the dense solar medium. For each impact parameter, $b$, we assume a one dimensional propagation through the Sun. As neutrinos travel through it, two important phenomena will modify the flavor composition and flux strength: neutrino oscillations and solar absorption effects, respectively. To first order, the opacity induced by traversing the Sun is given by

[TABLE]

where $l(b)$ is the chord traversed, $n_{N}(x;b)$ is the nucleon number density along the trajectory, and $\sigma_{\nu N}(E)$ is the total neutrino nucleon cross section. In the absorption the neutrino electron induced opacity is ignored since, for $10~{}{\rm GeV}<E_{\nu}<1~{}{\rm PeV}$, $\sigma_{\nu e}(E)\ll\sigma_{\nu N}(E)$ [61, 62]. The nucleon number density is obtained from [42]. We use the neutrino nucleon cross section as calculated in [63]. Equation (3.1) has some deficiencies that prevent us from using it directly. First, it neglects neutrino neutral current re-population and also tau regeneration, which make the Sun more transparent to high-energy neutrinos [64, 65]. Furthermore, it omits neutrino oscillations. A consistent formalism to treat both neutrino interactions together with oscillation, and including tau regeneration, is given in [66, 67, 68]; which relies on the density matrix formalism to represent the neutrino flux and describe the neutrino evolution by a quantum Boltzmann equation. This formalism is implemented in the publicly available nuSQuIDS package [39, 69], which is used in this work. We set the neutrino propagation medium according to the solar model given in [42] and set the neutrino oscillation parameters to the best fit values of [70]. This implementation has been tested by reproducing similar propagation calculations shown in [29, 71, 30].

The result of the neutrino propagation through the Sun is shown in Figure 6 for three different impact parameters. As expected due to oscillations the $\nu_{\mu}$ and $\nu_{\tau}$ fluxes are at similar intensity, where as oscillations with the $\nu_{e}$ flavor are suppressed by matter effects in the inner part of the Sun [72, 73, 74]. At 1 TeV the matter suppression is significant until approximately $0.9R_{\odot}$, whereas at a 100 GeV its only dominant until $\sim 0.5R_{\odot}$. When matter effects are no longer dominant the oscillation scale is comparable to the solar radius at 100 GeV and $\sim$10% at 1 TeV for $\nu_{e}\to\nu_{x}$ transitions. Furthermore, the solar opacity increases with energy and is larger for neutrinos than antineutrinos.

3.4 Predicted solar atmospheric neutrino flux at Earth

In order to predict the flux at the Earth from the results of Section 3.3 we must aggregate the neutrino fluxes calculated at several impact parameters, and also account for oscillation effects between Sun and Earth. The aggregated atmospheric solar neutrino flux is given by

[TABLE]

where $\phi_{b,\alpha}(E_{\nu})$ is the flux of flavor $\alpha$ at impact parameter $b$ after solar propagation and $\Omega_{\odot}$ is the total solid angle of the Sun. The aggregated flux is dominated by the larger impact parameters due to a geometrical enhancement and because the fluxes are less suppressed as $b$ increases (see Fig. 6). We evaluate the above integral over $b$ using a trapezium rule implementation with support points at intervals of $\delta b$=0.1. Some uncertainty is introduced by this discretization, particularly in the outer layers of the Sun where the $b$ dependence is largest and where most of the solid angle is contained. This was assessed by repeating the calculation with ten extra support points between $b=0.90$ and $b=0.99$ for one flux model. The integration uncertainty was found to be $<5$% below 10 TeV, growing to  10% at 100 TeV where the absorption effects are largest, and thus sub-dominant to other uncertainties treated in our calculation.

To account for neutrino oscillations from the Sun to the Earth we average the flux over one orbit, namely,

[TABLE]

where $T$ is one year, $t$ is the time through the year, $r(t)$ is the distance from the Sun to the Earth, and $P_{osc}$ is the vacuum oscillation probability. The distance difference between the aphelion and perihelion is $4\times 10^{6}\,{\rm km}$. This distance is $\sim 100\times$ larger than the neutrino production region (see Fig. 1), which implies we can ignore effects of oscillation averaged due to the production region uncertainties. The neutrino oscillation length in vacuum at $E_{\nu}=100~{}{\rm GeV}$ is approximately 1% of the Earth-Sun distance yearly variation, while they are comparable at 10 TeV. Thus neutrinos with energies less than $\sim 1$ TeV will exhibit oscillations which are averaged by the changing Earth-Sun distance, while oscillations at the highest neutrino energies remain observable. The result of this process is shown in Fig. 7.

4 Solar atmospheric neutrinos as a signal

Having predicted the solar atmospheric neutrino flux at Earth, it remains for us to consider its observability by realistic experiments. To calculate the number of observable events and their spectrum, we require both a flux and a detector effective area. The effective areas we consider are from a range of existing and proposed neutrino telescopes, shown in Figure 8. In all cases, we consider only the $\nu_{\mu}$ / $\bar{\nu}_{\mu}$ channel, since only this sample is expected to have sufficient pointing resolution to perform a high quality point source search in the direction of the Sun.

The IceCube effective area is taken from a dedicated WIMP search analysis [19] which has three discrete and independent sub-samples, and was provided in a data release available at [75]. We add the three effective areas to form a whole-experiment effective area for the purposes of this study. In the data release, effective areas are only provided up to an energy of 1775 GeV, and so we calculate the neutrino floor for IceCube-like effective areas up to this upper energy limit. PINGU extends IceCube’s reach to lower energies [21], and we assume that when PINGU runs, IceCube will also be running and maintaining its existing effective area at high-energy. Thus we construct the combined effective area as the published PINGU effective area below 40 GeV joined onto the IceCube effective area at higher energies. This is likely to be a slight underestimate of the PINGU+IceCube combined effective area since PINGU will contribute slightly to improving IceCube’s sensitivity above 40 GeV. The more recent IceCube result [20] provides only effective areas averaged over neutrinos and antineutrinos, and so can only be included approximately in our calculations. Because the sensitivity floor is only weakly sensitive to the details of the effective area, we use the more complete effective area model of [75] to calculate it, but will later estimate event rates for the new effective area from [20] by assuming equal effective areas for neutrinos and antineutrinos. Similarly, the ANTARES effective area reported in [76] and the KM3Net effective area from [77] are both reported averaged over neutrinos and antineutrinos, and we make similar approximations in these cases. The shape and location of the sensitivity floor, where solar atmospheric neutrinos overwhelm the expected dark matter signal, is independent of the normalization of the effective area and depends only on its shape. We will later show that in all cases described, the location of the sensitivity floor is largely independent of the experiment under consideration.

The expected rate of detectable solar atmospheric neutrinos, on the other hand, does depend on the effective area magnitude as well as shape. Considering the neutrino telescopes discussed above, we calculate the expected number of events per year detectable given each effective area. Our results are shown tabulated in Table 1. If sufficient pointing resolution can be achieved, according to our predictions this flux should already be detectable above the terrestrial atmospheric neutrino background in IceCube’s multi-year muon neutrino sample.

5 Solar atmospheric neutrinos as a background

In this section we compare the detectable WIMP self-annihilation signals described in Section 2 with the background from the flux of high-energy, solar atmospheric neutrinos. Figure 9 shows the comparison between the detectable flux of dark matter from self-annihilation modes $\chi\chi\rightarrow b\bar{b}$, $\chi\chi\rightarrow W^{+}W^{-}$ and $\chi\chi\rightarrow\tau^{+}\tau^{-}$, for various values of the cross section near the sensitivity floor, and WIMP masses 100 GeV, 250 GeV and 1000 GeV. This plot uses the effective area of the IceCube detector from [75]. Our next goal is to establish quantitatively which regions of WIMP parameter space fall below the sensitivity floor imposed by the solar atmospheric neutrino background.

This task is non-trivial, given that there is no unique definition of the sensitivity floor. We will define that the sensitivity floor has been reached when, in a representative toy analysis, the number of expected background events becomes equal in size to the signal being sought. As in the case for the sensitivity floor in direct dark matter searches, the caveats that small signals can still be detected over large backgrounds with sufficiently large statistics, and that additional degrees of freedom can open up new parameter space, will also apply here.

We thus make the assumption that the detector is ideal in the sense that events can be identified as being of solar origin unambiguously through pointing information, and the only additional information is the event energy, which is known precisely. In reality, both direction and energy will suffer from imperfect reconstructions. The uncertainty on the neutrino energy, in particular, will depend on energy reconstruction of the detected muon as well as the distribution of the initial energy among interaction products. It will result in a smearing of the events around the true initial neutrino energy with a wider spread at low energies. The neutrino floor derived in this work is therefore likely to be too deep, especially for low mass WIMPs.

As a representative dark matter search analysis we consider each WIMP model and define an energy region of interest (ROI) in which to count events. This ROI is defined such that 90% of the signal from the WIMP model, after convolution with the relevant effective area, is contained between the upper and lower thresholds of the ROI, with 5% escaping above and 5% escaping below. We then predict the expected signal (S) and background (B) counts in this region. We define the sensitivity floor as being surpassed when S=B. This analysis method is likely to be slightly weaker than a template likelihood based approach as described in [19], though it is stronger than a simple full-sample counting analysis described in [78].

Using the procedure described above, we derive the position of the neutrino floor for each detector and each channel. These are shown in Figure 10. Since the effective areas differ in their magnitude but not significantly in their shape, the neutrino floor always falls in a similar range in regions where the experiments have sensitivity. Figure 11 shows the comparison between these sensitivity floors, defined as the minimum accessible spin-dependent cross section given the effective areas considered, and the present leading experimental limits from IceCube from [20]. We see that the solar atmospheric neutrino background will effectively obscure solar dark matter signals less than one order of magnitude below existing spin-dependent cross section limits for some mass points. It is notable that PINGU does not reduce the sensitivity floor beyond IceCube in all channels, because its primary role is to extend the search region to lower energies where the solar atmospheric background is larger. On the other hand, augmenting IceCube with PINGU will slightly increase the detectable rates of both the solar atmospheric neutrinos and of the WIMP signal for some mass points, thus allowing this floor to be reached somewhat sooner. In all cases, our study shows that the sensitivity floor lies approximately one order of magnitude beyond the existing limit in the strongest exclusion regions.

6 Conclusions

We have derived the sensitivity floor for indirect solar WIMP annihilation searches, which arises due to interactions of cosmic rays in the atmosphere of the Sun. Using a new calculation of this flux, called the FJAWs model, we predict detectable event rates of around one event per year in existing and proposed neutrino experiments, which introduces a sensitivity floor for indirect dark matter searches using $\chi\rightarrow W^{+}W^{-}$, $\chi\rightarrow\tau^{+}\tau^{-}$, and $\chi\rightarrow b\bar{b}$ self-annihilation modes.

Our analysis assumes perfect energy and directional resolution and performs a counting experiment in the energy ROI. The resolution assumptions are generous, making our statement of the sensitivity floor slightly deep. On the other hand, the use of a full likelihood analysis may lead to a small improvement in sensitivity, though this is unlikely to change our main conclusions or to compensate for the generous assumptions about detector resolution.

We thus conclude that present limits on spin-dependent dark matter interactions from solar WIMP self-annihilation searches are less than one order of magnitude in sensitivity above the floor imposed by the irreducible, high-energy solar neutrino background at some mass points. Next-generation neutrino telescopes running for multiple years should be able to reach this sensitivity level within their operating lifetimes, thus saturating this particular method of searching for dark matter.

Note added

Shortly after the appearance of V1 of this paper on the arXiv pre-print server, two pre-prints that also propose the neutrino floor were submitted, [79] and [80], based on similar work developed independently. Disagreements in V1 were scrutinized and, through collaboration with the authors of [79], a minor error was identified and corrected in our work, bringing the proposed neutrino floors into closer agreement.

Also during final preparation of this manuscript, an updated IceCube solar WIMP result was published [20] which updates the originally cited work [19] using three years of data from the IC86 detector. We include here rate estimates and limits based on this new effective area. However, since the new paper does not provide separated $\nu$ and $\bar{\nu}$ effective areas, we retain the separate effective areas from [19] for calculating the neutrino floor. Since the details of the effective area model have only a minor impact on the location of the floor, this detail does not substantially affect our conclusions.

Acknowledgments

The authors thank David Seckel and Ali Kheirandish for their comments on the text, and Jean DeMerit and Niyousha Davichi for careful proof-reading. We thank the authors of [79] and [80] for rewarding scientific discussions. CAD is supported by NSF grants No. PHY-1505858 and PHY-1505855; GdW is supported by the Belgian Inter-University Attraction Pole (IUAP) network ‘fundamental interactions’; and BJPJ is supported by the University of Texas at Arlington.