1 Introduction

Based on inferences from observations of gravitational effects, it has long been believed that a significant fraction of the Universe is made up of dark matter (DM) (see [2]). However, very little is known about its properties and interactions. A weakly interacting massive particle (WIMP), whose relic abundance from a state of thermal equilibrium can make up DM has been the subject of considerable theoretical attention and experimental focus (see [3] for a comprehensive review).

Various complementary approaches have been pursued to detect the WIMPs that may constitute the DM halo of our Galaxy. Terrestrial direct detection (DD) experiments search for nuclear recoil events from the elastic scattering of WIMPs with the target nuclei of their detectors. Neutrino and gamma ray telescopes search for directional excesses over astrophysical backgrounds that may indicate the pair-annihilation of WIMPs, while collider searches look for the signatures of WIMPs being created in high-energy interactions of Standard Model particles.

Although the different search strategies have attained the sensitivity to probe the physically-motivated WIMP parameter space over the past few decades, they have failed to detect any signal. In the absence of a convincing detection, constraints have been derived on the interaction cross-sections of these hypothetical particles with Standard Model particles. Such an inference requires knowledge both of the density of DM \(\rho _{\mathrm{DM}}\) and of its velocity distribution function (VDF) \(f(\vec {v})\).

In the Standard Halo Model (SHM) [4], the DM of the halo is a collisionless gas in hydrostatic equilibrium with the stars, retaining the velocity distribution obtained during the formation of our Galaxy. An isotropic Maxwell–Boltzman velocity distribution in the Galactic rest frame is usually adopted.

Meanwhile, N-body simulations have hinted that a Maxwell–Boltzmann distribution does not accurately represent even the smooth component of the halo [5,6,7]. Recent observations point to the possibility that a dominant fraction of the DM in the Solar neighbourhood [8] may not yet have achieved dynamical equilibrium, perhaps due to the infalling tidal debris of a disrupted massive satellite galaxy of the Milky Way. New data also suggest that a substantial fraction of our stellar halo may lie in a strongly radially anisotropic population, the ‘Gaia sausage’ [9].

If so, constraints on WIMP-nucleon interactions derived assuming the SHM (from both direct and indirect searches) may be weakened. Direct detection experiments are preferentially sensitive to nuclear recoils from high velocity DM particles, while capture in the Sun is more likely for the slower fraction of the DM population. In this work we use the method of  [1], which is independent of the velocity distribution of the halo model to exploit this complementarity and derive conservative, upper limits on the spin-dependent DM-nucleon scattering cross-section by combining the results from [10, 11]. Here the DM velocity distribution is taken to be a completely general superposition of individual ’streams’ (delta functions in velocity), similarly to the halo-independent analysis of direct detection experiments  [12]. Although constraints from individual searches will now be dependent on the stream velocity, by exploiting the complementarity of the IceCube and PICO searches, constraints independent of the stream velocity can be obtained. This method also improves on previous assessments of halo model uncertainties on indirect DM detection  [13], by allowing the velocity distribution to be anisotropic. The resulting constraints are a factor of 2–4 worse than the PICO SHM constraints at low DM masses and up to an order of magnitude worse at high DM masses, depending upon the annihilation channel, but are independent of the halo model.

2 Detectors and data samples

2.1 IceCube 3 year Solar WIMP search

IceCube is a cubic-kilometer neutrino detector installed in the ice at the geographic South Pole between depths of 1450 and 2450 m. It relies on photomultiplier tubes housed in pressure vessels known as digital optical modules (DOM) for the optical detection of Cherenkov photons emitted by charged particles traversing the ice. The principal IceCube array is sensitive to neutrinos down to \(\sim \)100 GeV in energy [14,15,16]. The central region of the detector is an infill array known as DeepCore optimized in geometry and DOM density for the detection of neutrinos at lower energies, down to \(\sim \)10 GeV [17].

Over a detector uptime of 532 days corresponding to the austral winters between May 2011 and May 2014, two non-overlapping samples of upgoing track-like events, dominated by muons from charged current interactions of atmospheric \(\nu _{\mu }\) and \(\bar{\nu }_{\mu }\), were isolated [10]. During austral summers, the Sun being above the horizon, is a source of downgoing neutrinos and the signal is overwhelmed by a background of muons originating in cosmic ray interactions in the upper atmosphere.

The first sample, consisting of events that traverse the principal IceCube array, is sensitive to neutrinos in the 100 GeV–1 TeV range in energy, while the second sample is dominated by events starting in and around the DeepCore infill array, and is sensitive down to neutrinos of \(\sim \)10 GeV in energy.

An unbinned maximum likelihood ratio analysis of the directions and energies of the events that make up the two samples was unable to identify a statistically significant excess of neutrinos from the direction of the Sun. This enabled 90% CL upper limits on the DM annihilation induced neutrino flux to be computed according to the prescription of  [18] as presented in [10].

This can be interpreted as both a constraint on the annihilation rate of DM particles in the Sun, as well as on the scattering cross-section of DM with nucleons, although this has been usually done under the SHM assumption. In particle physics models where the DM couples to the spin of the nucleus and annihilates preferentially into SM particles that decay to produce a large number of high energy neutrinos (such as \(\tau ^+\tau ^-\)), the resultant constraints are the most stringent for DM mass above \(\sim 80~\hbox {GeV}\) [19].

2.2 PICO

The PICO collaboration searches for WIMPs using superheated bubble chambers operated at temperature and pressure conditions which lead to being virtually insensitive to gamma and beta radiation [20]. Events in PICO consist of the transition from liquid to gas phase, signalled by the nucleation of a bubble in the target material. This phase change is imaged by the cameras surrounding the active area, which trigger upon detecting the formation of a pocket of gas. Additional background suppression is achieved through the measurement of the acoustic signal generated by the event, allowing alpha particles to be discriminated from nuclear recoils. Details of the apparatus are available in [21]. The data used in this study were obtained from the PICO-60 detector, consisting of a \(52.2\pm 0.5~\hbox {kg C}_{3}~\hbox {F}_{8}\) target, operated roughly 2 km underground at SNOLAB in Sudbury, Ontario, Canada. The results used here come from an efficiency-corrected exposure of 1167 kg-days taken between November 2016 and January 2017 [11].

The response of the detector to WIMPs is dependent on the thermodynamic conditions, and is calibrated using in situ nuclear and electronic recoil sources. Additionally, the Tandem Van de Graaff facility at the University of Montreal was used to determine the detector response, using well-defined resonances of the \(^{51}\hbox {V}\)(p,n)\(^{51}\hbox {Cr}\) reaction to produce mono energetic neutrons at 61 and 97 keV. The combination of these measurements is simulated using differential cross-sections for elastic scattering on fluorine to produce the detector response.

3 DM velocity distributions and impact on constraints: the method

Following the method of [1], the velocity distribution of the DM (WIMP) population in the Solar system, \(f(\vec {v})\) can be expressed as the superposition of streams with fixed velocity \(\vec {v}_0\) with respect to the Solar frame.

$$\begin{aligned} f(\vec {v}) = \int _{|\vec {v}_0| \le v_{\mathrm{max}}} d^3 v_0 \delta ^{(3)}(\vec {v} - \vec {v}_0)f(\vec {v}_0) \end{aligned}$$
(1)

where \(v_{\mathrm{max}}\) is the maximum velocity at which WIMPs can be found, typically the escape velocity of the Galaxy. For every stream with velocity \(\vec {v}_0\) with respect to the Sun, upper limits can be derived from the null results of IceCube by requiring that the capture rate for the stream \(C_{\vec {v}_0}\) be less than or equal to \( C_{\mathrm{max}}\), the upper limit on the capture rate from the results of the experiment. For a direct detection experiment, which sees the same stream with velocity \(\vec {v}_0 - \vec {v}_{\mathrm{E}}(t)\) with respect to the Earth, similar constraints can be derived for each stream velocity by requiring that the event rate for the stream \(R_{\vec {v}_0}\) be less than or equal to \(R_{\mathrm{max}}\), the upper limit on the event rate from the results of the experiment. \(C_{\vec {v}_0}\) and \(R_{\vec {v}_0}\) are computed by evaluating the integrals of equations 2 and 3 of [1]. Since the PICO exposure period was too short for the Earth’s velocity \(\vec {v}_{\mathrm{E}}(t)\) to average out to zero, velocities are conservatively shifted by 30.29 km s\(^{-1}\) (the velocity of the Earth around the Sun at perihelion [22]) when computing \(R_{\vec {v}_0}\). For the capture rates in the Sun, the integrals were evaluated using the density profile and nuclear abundances in the Sun for protons and nitrogen nuclei (the second most abundant species with nuclear spin) in the standard Solar model [23] as implemented in sunpy [24]. Nuclear form factors as implemented in dmdd [25] for spin-dependent scattering, corresponding to the \(\Sigma '_{1M}\) (Axial transverse electric response) and \(\Sigma ''_{1M}\) (Axial longitudinal response), Table 1 of  [26] were employed for the event rate calculations in PICO.

Figure 1 demonstrates the evolution of the constraints on the spin-dependent DM-proton scattering cross-section from both IceCube and PICO as \(|v_0|\) is varied. The individual constraints on the cross section are computed from the constraints on the capture rate in the Sun already derived in [10] as well as the constraint on the event rate within PICO presented in [11]. For a WIMP of mass M scattering off a nucleus of mass m, the maximum stream velocity at which capture is allowed is given by [1]:

$$\begin{aligned} v_{\mathrm{max}} = 2 v_{\mathrm{esc}} \sqrt{\frac{Mm}{|M-m|}} \end{aligned}$$
(2)

where \(v_{\mathrm{esc}}\) is the escape velocity. Consequently, above certain threshold values of the stream velocity, capture by scattering off protons is kinematically impossible and only nitrogen nuclei contribute to the capture rate.

Fig. 1
figure 1

Constraints at \(\gtrsim 90\%\) CL on the spin-dependent DM-proton scattering cross-section from both IceCube and PICO for different values of \(|v_0|\), for 40 (700) GeV WIMPs annihilating to \(b\bar{b}\) are shown on the left (right). For 40 GeV WIMPs, as the efficiency of PICO falls off below stream velocities of c (the speed of light) \(\times 10^{-3}\), Solar capture by scattering off hydrogen nuclei provides a complementary bound, while for 700 GeV WIMPs, a bound is provided only by the much less abundant nitrogen nuclei in the Sun

Full size image
Fig. 2
figure 2

DM velocity distribution independent constraints on the SD DM-nucleon interaction cross-section \(\gtrsim 90\%\) CL. Systematic uncertainties are presented as shaded regions. The traditional SHM upper limits at \(90\%\) CL from IceCube and PICO are shown as dashed and dash dotted lines. The kinks in the constraints at \(\sim 250~\hbox {GeV}\) (for \(W^{+}W^{-}\)and \(\tau ^{+}\tau ^{-}\)) and \(\sim 700~\hbox {GeV}\) (for \(b\bar{b}\)) are explained in Section 4

Full size image

Subsequently, the largest value of the scattering cross-section allowed by both IceCube and PICO, \(\sigma ^{\mathrm{HI}}\), can be determined at the velocity of least constraint, \(v_{\mathrm{LC}}\), where \(\sigma ^{\mathrm{PICO}}_{\mathrm{max}}(v_{\mathrm{LC}}) = \sigma ^{\mathrm{IceCube}}_{\mathrm{max}}(v_{\mathrm{LC}})\). This procedure is illustrated in Figure 1 for two specific models, 40 GeV and 700 GeV WIMPs annihilating to \(b\bar{b}\).

4 Results and conclusions

The resultant DM velocity independent constraints are illustrated in Fig. 2 and presented in Table 1. For the “hard” channels ( \(W^{+}W^{-}\)and \(\tau ^{+}\tau ^{-}\)), which produce a relatively large number of neutrinos at energies just below the DM mass, the DM-velocity-independent constraints are in general worse only by a factor of 2–4 compared to the PICO SHM constraints. However, at a DM Mass of \(\sim \)250 GeV (\(\sim \)700 GeV for \(b\bar{b}\)), the constraints are significantly worse because the DM particle velocities just below the PICO threshold are still too high to be captured by scattering off protons in the Sun (see Fig. 1). At immediately higher masses, the constraints improve because the IceCube sensitivity improves with the DM mass in this range. The constraints are in agreement with the findings by [27]. The IceCube constraints were recomputed with Monte-Carlo data sets under varying assumptions of all systematic uncertainties as described in [10]. The dominant uncertainties were found to originate in the photodetection efficiency of the photomultiplier tubes that make up the DOMs, as well as the optical properties of the ice. Since these constraints correspond to the same annihilation rates of DM particles in the Sun reported in [10], capture-annihilation equilibrium continues to be a valid assumption. The dominant uncertainties in the detector acceptance of PICO originate in the uncertainties of the neutron beam used in the calibration process. These are propagated to the final level and shown as shaded regions. Conservatively, the pessimistic efficiencies of PICO have been used to derive the constraints. While these constraints are robust with respect to any uncertainties in the velocity distribution of DM particles, they are still susceptible to uncertainties and/or fluctuations in the local density of DM, and are presented for the benchmark local density of \(\rho _{\mathrm{DM}} = 0.3~\hbox {GeV cm}^{-3}\), and scale inversely with this quantity.

Table 1 Constraints on the SD DM-nucleon cross-section. SHM constraints from PICO and IceCube, as well as the DM velocity distribution independent constraint are presented at \(\gtrsim 90\%\) CL. The velocity of DM particles at which the cross-section is least constrained (\(V_{\mathrm{LC}}\)) is also presented for each point. The constraints are conservative with respect to systematic uncertainties
Full size table

IceCube Collaboration:

USA – U.S. National Science Foundation-Office of Polar Programs, U.S. National Science Foundation-Physics Division, Wisconsin Alumni Research Foundation, Center for High Throughput Computing (CHTC) at the University of Wisconsin-Madison, Open Science Grid (OSG), Extreme Science and Engineering Discovery Environment (XSEDE), U.S. Department of Energy-National Energy Research Scientific Computing Center, Particle astrophysics research computing center at the University of Maryland, Institute for Cyber-Enabled Research at Michigan State University, and Astroparticle physics computational facility at Marquette University; Belgium – Funds for Scientific Research (FRS-FNRS and FWO), FWO Odysseus and Big Science programmes, and Belgian Federal Science Policy Office (Belspo); Germany – Bundesministerium für Bildung und Forschung (BMBF), Deutsche Forschungsgemeinschaft (DFG), Helmholtz Alliance for Astroparticle Physics (HAP), Initiative and Networking Fund of the Helmholtz Association, Deutsches Elektronen Synchrotron (DESY), and High Performance Computing cluster of the RWTH Aachen; Sweden – Swedish Research Council, Swedish Polar Research Secretariat, Swedish National Infrastructure for Computing (SNIC), and Knut and Alice Wallenberg Foundation; Australia – Australian Research Council; Canada – Natural Sciences and Engineering Research Council of Canada, Calcul Québec, Compute Ontario, Canada Foundation for Innovation, WestGrid, and Compute Canada; Denmark – Villum Fonden, Danish National Research Foundation (DNRF), Carlsberg Foundation; New Zealand – Marsden Fund; Japan – Japan Society for Promotion of Science (JSPS) and Institute for Global Prominent Research (IGPR) of Chiba University; Korea – National Research Foundation of Korea (NRF); Switzerland – Swiss National Science Foundation (SNSF). United Kingdom – Department of Physics, University of Oxford.

PICO Collaboration:

The PICO collaboration wishes to thank SNOLAB and its staff for support through underground space, logistical and technical services. SNOLAB operations are supported by the Canada Foundation for Innovation and the Province of Ontario Ministry of Research and Innovation, with underground access provided by Vale at the Creighton mine site. We wish to acknowledge the support of the Natural Sciences and Engineering Research Council of Canada (NSERC) and the Canada Foundation for Innovation (CFI) for funding. We acknowledge the support from the National Science Foundation (NSF) (Grant 0919526, 1506337, 1242637 and 1205987). We acknowledge that this work is supported by the U.S. Department of Energy (DOE) Office of Science, Office of High Energy Physics (under award DE-SC-0012161), by the DOE Office of Science Graduate Student Research (SCGSR) award, by DGAPA-UNAM (PAPIIT No. IA100118) and Consejo Nacional de Ciencia y Tecnoloǵıa (CONACyT, México, Grant No. 252167 and A1-S-8960), by the Department of Atomic Energy (DAE), Government of India, under the Centre for AstroParticle Physics II project (CAPP-II) at the Saha Institute of Nuclear Physics (SINP), European Regional Development Fund- Project “Engineering applications of microworld physics” (No. CZ.02.1.01/0.0/0.0/16 019/0000766), and the Spanish Ministerio de Ciencia, Innovación y Universidades (Red Consolider MultiDark, FPA2017–90566–REDC). This work is partially supported by the Kavli Institute for Cosmological Physics at the University of Chicago through NSF grant 1125897 and 1806722, and an endowment from the Kavli Foundation and its founder Fred Kavli. We also wish to acknowledge the support from Fermi National Accelerator Laboratory under Contract No.DE-AC02-07CH11359, and from Pacific Northwest National Laboratory, which is operated by Battelle for the U.S. Department of Energy under Contract No. DE- AC05-76RL01830. We also thank Compute Canada (www.computecanada.ca) and the Centre for Advanced Computing, ACENET, Calcul Quebeć, Compute Ontario and WestGrid for computational support.