Searches for neutrinos from cosmic-ray interactions in the Sun using seven years of IceCube data

The first theoretical calculations for SA$\nu$s date back to 1991 [1, 2]. The authors modeled gamma-ray, neutrino, antiproton, neutron, and antineutron fluxes that are initiated by the interactions of cosmic rays with the solar atmosphere. The flux originates from the solar disk as cosmic rays that are mirrored in the solar atmosphere are expected to contribute significantly to the flux [1]. While these early predictions were based on semi-analytical calculations, full numerical simulations of the interactions based on the Monte Carlo method have been performed in Ref. [4]. In more recent publications [7, 8], uncertainties in the predicted neutrino energy spectra from the choice of the primary cosmic-ray flux models, particle interaction models, solar density models, and neutrino oscillation parameters are discussed. Neutrinos generated in the atmosphere of a non-magnetic Sun are propagating through the Sun and their attenuation due to absorption is included. The impact of mirroring of charged particles in solar magnetic flux tubes anchored at the bottom of the photosphere is not considered as it is expected to be sub-dominant at energies above 100 GeV [1]. We will discuss the energy dependent spatial distribution of the solar disk flux in Sec. 5.3.1.

Solar atmospheric neutrino fluxes have been implemented in the simulation framework, WIMPSim [28], which we used for our signal prediction. We obtained neutrino fluxes for all neutrino flavor channels. However, for our analysis we only consider the $\nu_{\mu}+\bar{\nu}_{\mu}$ channel to benefit from IceCube’s excellent angular resolution O(1°). The impact of additional flavor contributions are discussed as part of our systematic uncertainty studies (see Sec. 4 and 5.3).

In Fig. 1, the $\nu_{\mu}+\bar{\nu}_{\mu}$ neutrino flux predictions as well as their uncertainties are shown as the shaded regions. The range of the shaded gray area spans the energy spectra of the results published in [8, 30]. The red region represents the simulation results obtained by running the built-in codes in WIMPSim [7, 28]. Neutrino oscillations are fully taken into account when propagating the neutrinos from the Sun to the Earth and appear as wiggles in the theoretical flux predictions. The oscillations of these high-energy neutrinos, however, cannot be resolved due to the limited energy resolution of IceCube. In Fig. 1 we compare the energy spectra from a set of representative models [7, 8]. For the visualization, we sample the flux predictions with coarse energy bins, that were chosen to reflect IceCube’s energy response function. Only the parametrized energy flux from Ref. [4] (IT1996) did not include neutrino oscillations. Ref. [6] has shown that if the primary flavor ratio of SA$\nu$s $(\nu_{e}:\nu_{\mu}:\nu_{\tau})$ is $(1:2:0)$, it would be roughly close to $(1:1:1)$ at Earth. The flavor ratio of the IT1996 fluxes integrated in the range $(10^{2.0},10^{7.0})$ GeV is $(0.92,2.08,0)$. For simplicity, IT1996 fluxes for $\nu_{\mu}+\bar{\nu}_{\mu}+\nu_{e}+\bar{\nu}_{e}$ are divided by a factor of 3 to apply neutrino oscillation effect, shown in Fig. 1 as the black line. Newer reference fluxes [7, 8] already include the effect of the oscillations.

We measure the flux normalization of the SA$\nu$s in this analysis. A comparison of signal predictions [7, 8, 4] shows that the SA$\nu$ spectral shapes are similar enough that we are not expected to be sensitive to individual models. We therefore choose one representative baseline energy spectrum (shown as the blue line in Fig. 1). The baseline energy spectrum is chosen from [7] and uses the Hillas-Gaisser 3-generation model [31] for the primary cosmic-ray spectrum, a combination of the Serenelli [32] and the Stein et al. [33] models for the solar density profile, and the normal mass ordering.

Finally, we note that the current leading models neglect solar magnetic field effects. These effects influence cosmic-ray propagation and the cascade development, which in turn influence the neutrino signal. The effect of magnetic fields on cosmic-ray propagation can be indirectly measured through the absorption of cosmic rays in the Sun, which in turn makes a corresponding deficit of cosmic rays in the direction of the Sun. The so-called cosmic-ray Sun shadow has been observed by the Tibet air shower array, including a variation of the intensity correlated with the solar cycle [34]. IceCube also observed the Sun shadow and found a correlation with the sunspot number with a likelihood of 96% [35]. The Sun shadow is sensitive to magnetic field models [36, 37, 38, 39] and recent works with numerically computed trajectories of charged cosmic rays confirm the observationally established correlation between the magnitude of the shadowing effect and both the mean sunspot number and the polarity of the magnetic field during a solar cycle [17]. In general, however, high-energy cosmic rays are expected to be energetic enough not to be influenced by magnetic fields. Therefore, only for neutrino production below 200 GeV [1] or 1 TeV [40] is it expected to become significant. Theoretical works using HAWC’s Sun shadow observation predict a factor of about two difference in SA$\nu$ flux between solar minimum and maximum at 200 GeV [40].

Most events in IceCube are downward-going atmospheric muons from cosmic-ray air showers in the Earth atmosphere. These muons can be efficiently rejected by selecting events reconstructed upward, i.e. with declination $\delta>-5\degree$. The well-established event selection for the upward-going neutrino events has achieved a purity of 99.7% [19]. In the remaining sample, the main background arises from terrestrial atmospheric neutrinos produced by decays of mesons within cosmic-rays air showers. Another irreducible, but sub-dominant background is due to isotropic astrophysical neutrinos. They can be described by an unbroken power-law with a spectral index of $2.19\pm 0.1$ and a flux normalization, $\Phi_{\mathrm{100~TeV}}=1.01^{+0.26}_{-0.23}\cdot 10^{-18}\mathrm{GeV^{-1}cm^{-2}s^{-1}sr^{-1}}$ at 100 TeV, obtained from fits to the data [19]. The astrophysical neutrinos are included as a background in this analysis, but the uncertainties of the best-fit parameter values are negligible due to it’s small contribution to the background rate.

Neutrinos from dark matter annihilations in the Sun could result in a competing signal that has been extensively searched for at neutrino telescopes [41, 42, 43, 44, 45, 46]. The expected neutrino spectra from solar dark matter strongly depend on the dark matter mass and annihilation channels. As dark matter annihilations are expected to occur in the center of the Sun, neutrino absorption becomes important for energies above 100 GeV and fluxes are significantly attenuated above that energy. As a result, spectra are expected to be significantly different from that of SA$\nu$s [47]. Purely based on event rate expectations at neutrino detectors, one can compute a sensitivity floor for indirect dark matter searches from the Sun [47, 8, 7, 40]. Past dark matter searches were not sensitive enough to have significant backgrounds from solar atmospheric neutrinos. However, in the near future they are expected to reach the neutrino floor from SA$\nu$s.

Another competing signal may arise from the interactions of cosmic rays with thermal solar photons. These can interact to form $\Delta^{+}$ baryons which quickly decay, producing muons and neutrinos from subsequent pion decays [48]. The expected flux from $\Delta^{+}$ is small and few events are expected in IceCube, so we assume no contributions from the process in this analysis. Larger active volumes, like those proposed for IceCube-Gen2 [49, 50], may be needed to observe events from these interactions.

A good angular resolution is necessary to search for SA$\nu$s because the angular size of the Sun is $\theta_{\odot}\sim 0.27\degree$. Muons traversing the entire detector are reconstructed with good angular resolution as kilometer-long tracks, so-called “through-going muons.” We restrict ourselves to IceCube’s neutrino sample of predominantly through-going muons [21] providing $1.0\degree$ and $0.6\degree$ median angular resolutions at 1 TeV and 10 TeV neutrino energy, respectively. As the events are not fully contained in the detector volume, the energy resolutions are limited to $\Delta\log_{10}(E/\mathrm{1GeV})\sim 0.3$ and $\sim 0.5$ at 1 TeV and 10 TeV, respectively.

The data samples consist of three sub-samples covering a total of seven years. There are three time periods: IC79-2010, IC86-2011, and IC86-(2012-2016). An optimized event selection has been used for each configuration. The ranges of the reconstructed energies are ($10^{2.2}$, $10^{7.2}$) GeV and ($10^{2.0}$, $10^{7.0}$) GeV for IC79-2010 and IC86-(2011,2012-2016). Events below the horizon (declination, $\delta>-5\degree$) are selected to exclude atmospheric muon events. Unlike Ref. [21], we only consider events where the Sun is below the horizon, resulting in a total analysis livetime of 1406.62 days.

Angular separation ($\theta_{\odot}$) is defined as an angular distance between the reconstructed position and the center of the Sun at the event trigger time. We use an IceCube internal software implementation that relies on the Positional Astronomy Library (PAL) [51] to obtain the position of the Sun based on the Modified Julian Date and the IceCube detector location. The tracking of the Sun was cross checked with data from NOAA to verify that it agrees within 0.01 degrees. We define a Region of Interest (RoI) as a $5\degree$ window around the center of the Sun. The size of the RoI was optimized taking into account the signal sensitivity and the computational time. The sensitivity only marginally improves for larger windows as the RoI contains 96% of the reconstructed signal events.

Simulations are used to obtain probability density functions (PDFs) of the signal and background in the muon neutrino and muon anti-neutrino channels. Simulations samples are taken from the well established and tested IceCube point source analysis [21]. We reuse these simulations but apply selection criteria on the angular separations within the RoI and reweight them with the effective analysis livetime. The background expectations are constructed using simulation, weighted to best-fit parameters for atmospheric and astrophysical neutrino backgrounds found from previous fits to data [19]. In Fig. 2, the comparisons between the simulations for terrestrial atmospheric neutrinos and the total data samples are shown for the reconstructed zenith angle and energy distributions. The simulation samples and the data samples are well-matched within 4% differences.

Signal simulations are obtained by re-weighting the simulated events with the given SA$\nu$ energy spectrum for muon neutrinos (see Sec. 3.1). The angular separations between the center of the Sun and the events are calculated. The azimuthal directions of the signal events are uniformly scrambled. Events are also randomized in zenith using the probability distribution as a function of angular distance from the center of the Sun for the given source hypothesis. We account for the movement of the Sun in zenith by weighting events using the fraction of livetime spent by the Sun in 30 zenith bins from $85\degree$ to $113.4\degree$. In Fig. 3, two-dimensional angular distributions are shown for the baseline signal and background assumptions.

An unbinned likelihood method [52] is applied to find evidence of SA$\nu$s in seven years of the data sample. The likelihood function, $L_{j}$, for each sub-sample, j, is defined by

where $j$ is the index of the sub-sample, $i$ is the event index, $n_{\textrm{tot},\,j}$ is the total number of events and $n_{\textrm{s},\,j}$ is a number of signal events. For each event, $\theta_{i}$ is the angular separation to the Sun and $E_{i}$ is the reconstructed muon energy [53]. The function of $p_{\textrm{sig},\,j}$ and $p_{\textrm{bkg},\,j}$ are the signal and background PDFs evaluated at the location of each event, respectively. In Fig. 4, PDFs of the IC86-(2012-2016) sub-sample are shown. The PDFs are obtained from the simulations and the corresponding likelihood functions are used to study a particular energy spectrum $M_{\textrm{sig}}$. We combine different sub-samples with a uniform signal emission and use the maximum likelihood estimator to estimate the signal strength. The total likelihood function, $L$, is a multiplication of the likelihood functions, $L_{j}$, for the three sub-samples mentioned in Sec. 4.1. The fractions ($f_{j}$) of the total expected signal events for each sub-sample are calculated: $f_{j}=\bar{n}_{\textrm{s},\,j}/\Sigma_{k}\bar{n}_{\textrm{s},\,k}$ where $\bar{n}_{s,\,j}$ is an expected number of signal events from the simulations. The total likelihood function is redefined as a function of the total signal strength $\mu$ with converting $n_{\textrm{s},\,j}$ to $\mu f_{j}$:

The theoretical distribution of flux across the solar disk is expected to depend on neutrino energy via the energy dependence of IceCube reconstructions. To include this correlation, two dimensional PDFs, shown in Fig. 4, are used to model the signal and background distributions in the likelihood functions.

We define the test statistic (TS)

as the likelihood ratio between the best-fit value and the null hypothesis. The range of $\hat{\mu}$ is not restricted to positive values. Therefore, we can track the sign of $\hat{\mu}$ to separately determine sensitivities for a positive or negative signal strength. The negative signs of $\hat{\mu}$ can appear when the alternate hypothesis represents an under-fluctuation relative to the background prediction, especially that the under-fluctuation can be enhanced by the Sun shadow, see Sec. 5.3.3.

Pseudo-experiments are conducted to obtain the TS distribution for a given hypothesis. Each pseudo-experiment consists of mock samples generated by random sampling based on each PDF of a certain hypothesis. The number of signal and background events are random variables that are Poisson distributed. The mean of the Poisson distribution for the number of background events is given by the expected number of events from the simulations, $\bar{n}_{\textrm{bkg}}=1147.4$ in the RoI. Depending on the hypotheses, the mean for the signal $\bar{\mu}$ is scaled, e.g. $\bar{\mu}$=0 for the null hypothesis and $\bar{\mu}=C_{s}\cdot\bar{n}_{\textrm{sig}}$, where $C_{s}$ is a scale factor to test $C_{s}$ times larger signal hypotheses. $\bar{n}_{\textrm{sig}}$ is the expected number of signal events determined by combining a given signal model with the simulated detector response. The expected number of background events, $\bar{n}_{\textrm{bkg}}$, and signal events, $\bar{n}_{\textrm{sig}}$, also change with the PDFs according to the hypotheses chosen.

In Fig. 5, the blue histogram is the TS distribution for the background-only hypothesis. Negative TS values appear when the likelihood function is maximized with a negative signal strength, due to under-fluctuations in the background rate. The median of the histogram, indicated by the vertical dashed blue line in Fig. 5, is close to zero. The 90% confidence interval (C.I.) is obtained with the Feldman-Cousins method [54] for each alternate hypothesis. $\mu_{90}$ is defined by $\bar{\mu}$ of the Poisson mean when the minimum of the 90% C.I. is larger than the median of the TS distribution for the null hypothesis. The 90% confidence level (C.L.) upper sensitivities are set with $\mu_{90}=C_{s,90}\cdot\bar{n}_{\textrm{sig}}$ for each SA$\nu$ flux model given by Refs. [7, 8, 4]. PDFs are used in the likelihood functions and the random sampling for the mock samples of the pseudo-experiments. The sensitivities to each flux model are calculated with the corresponding PDFs. The sensitivity to the baseline signal spectrum is shown as part of our final result, see the red solid line in Fig. 8. It is 12.8 times larger than the theoretical expectation from the baseline model flux [7].

We investigate how the sensitivity of our analysis depends on different choices for flux distributions on the solar disk, oscillation parameters, the effect of the Sun shadow on the backgrounds and detector uncertainties. The differences between the sensitivities are quantified relative to the baseline model as systematic uncertainties.