SEARCH FOR PROMPT NEUTRINO EMISSION FROM GAMMA-RAY BURSTS WITH ICECUBE

While cosmic rays have been observed with energies up to 10²⁰ eV, their sources remain unknown. Gamma-ray bursts (GRBs) have been proposed (Vietri 1995) as promising candidate sources of ultra-high-energy cosmic rays (UHECRs) because of their extremely large energy release over timescales of only $\sim {{10}^{-3}}-{{10}^{3}}$ s. In the popular fireball model (e.g., Shemi & Piran 1990; Piran 2004; Mészáros 2006), gamma-rays are produced by the dissipation of kinetic energy in an ultra-relativistic fireball flowing outward from a cataclysmic stellar collapse or merger. If GRBs accelerate protons with comparable efficiency to electrons, then they could account for most or all of the UHECR flux (Waxman 1995). In this case, protons and gamma-rays in the fireball interact through channels such as the Δ-resonance process $p+\gamma \to {{{\Delta }}^{+}}\to n+{{\pi }^{+}}$. The charged pions decay leptonically via ${{\pi }^{+}}\to {{\mu }^{+}}+{{\nu }_{\mu }}$ followed by ${{\mu }^{+}}\to {{e}^{+}}+{{\nu }_{e}}+{{\bar{\nu }}_{\mu }}$. Waxman & Bahcall (1997) noted that this neutrino flux could be measured on Earth by a sufficiently large detector. Neutrinos correlated with GRBs would be a "smoking-gun" signal for UHECR acceleration in GRBs. To date, however, neither IceCube (Abbasi et al. 2011b, 2012) nor ANTARES (Adrián-Martínez et al. 2013) have observed such a signal.

IceCube is a km³ scale neutrino detector deployed deep in the south polar ice cap. The completed detector consists of 5160 digital optical modules (DOMs), with 60 DOMs mounted on each of 86 strings. Construction was performed during Austral summers, with the final strings deployed in 2010 December. Photomultiplier tubes (PMTs) in the DOMs detect Cherenkov light emitted by energetic charged particles produced in neutrino–nucleon interactions in the ice. When a DOM collects sufficient charge, digitized PMT waveforms are transmitted to the data acquisition system at the surface of the ice. When eight DOMs initiate such launches within 5 μs, a trigger is formed which results in initial processing, filtering, and further transmission of data via satellite to servers in the north. In previous publications, the PMTs (Abbasi et al. 2010), data acquisition methods (Abbasi et al. 2009), and overall detector operations (Achterberg et al. 2006) have been discussed in detail. Data sets were collected during construction using the partially completed detector configurations, each of which was active for approximately one year. The results presented here are derived from the first year of data from the completed 86 string detector in addition to data from the 40, 59, and 79 string configurations.

While IceCube is sensitive to neutral and charged-current interactions of all neutrino flavors coming from any direction, in this analysis, we restrict our focus to up-going charged-current ${{\nu }_{\mu }}$ interactions at energies above 1 TeV. Product muons from such a signal can travel several kilometers through the ice, providing high detection efficiency and good angular resolution that both improve with increasing neutrino energy. By selecting up-going muons with declination greater than −5°, we use the Earth (and, near the horizon, the ice cap itself) as a shield to attenuate the large flux of muons produced by cosmic-ray interactions in the atmosphere. The search will be extended to all interaction channels and the entire sky in separate papers.

The originating direction of muons passing through IceCube is reconstructed using a maximum likelihood method (Ahrens et al. 2004) to fit the spatial and temporal Cherenkov light pattern observed by the DOMs. IceCube is sensitive to muons with sufficiently high energy that the interaction frame is highly boosted with respect to the detector frame so that the muon trajectory is nearly collinear with the neutrino. Neutrino angular resolution is affected by both the deviation angle of the product muon, which decreases with increasing neutrino energy, and the accuracy of the reconstruction of the muon track, which is limited by light timing uncertainties due to photon scattering in the ice. Including both of these effects, the median neutrino angular error for simulated neutrinos surviving the quality cuts used in this analysis is 1° at ∼TeV energies; at ∼PeV energies, this value improves to 0.5° and the muon deviation angle is negligible. For each neutrino individually, the angular uncertainty (${{\sigma }_{\nu }}$) is estimated using the width of the optimum in the fit likelihood space (Neunhöffer 2006).

Muon energy is reconstructed by measuring the charge collected by the DOMs as the muon traverses the detector. Very good neutrino energy resolution is possible for analyses requiring the interaction vertex to be contained within the instrumented volume (Aartsen et al. 2014a). In this search, most of the sensitivity comes from neutrinos interacting outside of the instrumented volume. Since the location of the interaction vertex is generally not known, muons can lose significant energy before reaching the instrumented volume. Therefore, the reconstructed muon energy must be interpreted as an approximate lower bound on the neutrino energy.

Down-going cosmic-ray-induced muons trigger the completed detector at a rate of over 2 kHz. A large fraction of these events are correctly reconstructed as down-going and are easily excluded from this analysis. The dominant remaining backgrounds are muons passing near the boundary of the instrumented volume and emitting light upwards and multiple independent muons traversing the detector at the same time. These backgrounds, which often yield incorrect up-going reconstructions, are rejected using parameters described in previous work (Abbasi et al. 2011c) including (1) fit quality parameters from a progression of reconstructions that apply increasingly detailed ice and DOM response modeling; (2) comparison of the fit quality for unbiased and down-going-biased reconstructions; (3) reconstruction results for time- and geometry-based split subsets of the event data; and (4) topology variables related to the distribution of DOM pulses about the reconstructed muon path. Event selection criteria were optimized separately for each detector configuration. For the 40 and 59 string configurations, previously published event selection criteria were re-used. For the 40 string configuration, a simple set of cuts selected events which performed well in several quality criteria (Abbasi et al. 2011b), while for subsequent configurations, Boosted Decision Tree forests (Freund & Schapire 1997) were used to synthesize a single quality parameter from all available event information. The final sample has a data rate of ∼3.8 mHz in the completed detector and consists primarily of atmospheric muon neutrinos from the northern hemisphere with ∼15% contamination from misreconstructed cosmic-ray-induced muons. Atmospheric neutrinos constitute an irreducible background which can only be separated statistically from astrophysical neutrinos based on reconstructed energy and temporal and directional correlation with a GRB.

Between 2008 April 5 and 2012 May 15, 592 GRBs were observed at declinations greater than −5° and reported via the GRB Coordinates Network⁵¹ and the Fermi GBM catalogs (Gruber et al. 2014; von Kienlin et al. 2014). Bursts during commissioning and calibration phases are excluded. This analysis includes 506 bursts which occurred during stable IceCube data collection. The search window is determined by the time of gamma emission and the location in the sky for each burst. When multiple satellites observed a given burst, the gamma emission time (${{T}_{100}}$) is defined by the most inclusive start and end times (T₁ and T₂) reported by any satellite. The angular window is determined by the direction and angular uncertainty (${{\sigma }_{{\rm GRB}}}$) given by the satellite reporting the smallest angular uncertainty. Fermi GBM, which observes the most bursts, typically has a total statistical plus systematic uncertainty of a few degrees or more, but for bursts observed by other satellites, the uncertainty is generally $\ll 1{}^\circ$ (Winkler et al. 2003; Gehrels et al. 2004; Feroci et al. 2007; Hurley et al. 2010). When an asymmetric error ellipse is reported, the larger axis is used. The small GRB time and space windows, along with the low atmospheric neutrino rate, make this a nearly background-free search, with a sensitivity that improves nearly linearly with the number of bursts observed. For modeling neutrino fluence predictions, gamma-ray fluence parameters are taken from satellite measurements, and unmeasured model inputs are assumed as in our previous work (Abbasi et al. 2010). We catalog burst information in a publically accessible online database.⁵²

We use an unbinned maximum likelihood analysis based on Braun et al. (2008) to test for a correlation between GRBs and neutrino events. The likelihood ${\mathcal{S}}$ that a given event is a signal event and ${\mathcal{B}}$ that it is a background event are the products of separately normalized time, direction, and energy probability distribution functions (PDFs):

$$\begin{eqnarray}&&{\mathcal{S}}/{\mathcal{B}}={{(S/B)}_{{\rm time}}}\ {{(S/B)}_{{\rm dir}}}\ {{(S/B)}_{{\rm energy}}}.\end{eqnarray} \tag{ 1 }$$

For a given burst, the signal time PDF is constant during gamma emission. Before and after gamma emission, the signal time PDF falls smoothly to zero with Gaussian tails that have a width parameter given by

$$\begin{eqnarray}{{\sigma }_{{\rm time}}}=\left\{ \begin{array}{ccccccccccccccc} 2\;{\rm s} & {{T}_{100}}\lt 2\;{\rm s}, \\ {{T}_{100}} & 2\;{\rm s}\leqslant {{T}_{100}}\lt 30\;{\rm s}, \\ 30\;{\rm s} & 30\;{\rm s}\leqslant {{T}_{100}}. \\ \end{array} \right.\end{eqnarray} \tag{ 2 }$$

The burst time window is truncated at $4{{\sigma }_{{\rm time}}}$ before and after the gamma emission, and the background time PDF is constant throughout this time window. The signal direction PDF is a two-dimensional circular Gaussian:

$$\begin{eqnarray}&&{{S}_{{\rm dir}}}(\nu ,{\rm GRB})=\frac{1}{2\pi \sigma _{{\rm dir}}^{2}}\;{\rm exp} \left( -\frac{{\Delta }{{{\Psi }}^{2}}}{2\sigma _{{\rm dir}}^{2}} \right),\end{eqnarray} \tag{ 3 }$$

where $\sigma _{{\rm dir}}^{2}=\sigma _{{\rm GRB}}^{2}+\sigma _{\nu }^{2}$ and ${\rm \ \Delta\ \ \Psi\ }$ is the angular separation between the burst and the reconstructed muon direction. The background direction PDF is constructed from off-time data, accounting for the declination-dependent atmospheric neutrino event rate. The energy PDFs are computed from the reconstructed muon energy. While this reconstruction only provides a lower bound on the neutrino energy, it is nevertheless useful for probabilistically distinguishing a possible astrophysical flux from the atmospheric background, which has a softer spectrum. The background energy PDF is taken from off-time data in the energy range where we have good statistics; at higher energies, this PDF is extended using simulated atmospheric neutrinos. The signal energy PDF is computed using simulated signal events with an ${{E}^{-2}}$ spectrum, which provides good sensitivity to a wide range of GRB model spectra.

In this search, the observed number of events N in the on-time window is not known a priori. For supposed signal and background event rates n_s and n_b, respectively, the probability of observing N events is given by the Poisson distribution:

$$\begin{eqnarray}&&P({{n}_{s}},{{n}_{b}})=\frac{{{({{n}_{s}}+{{n}_{b}})}^{N}}}{N!}{\rm exp} [-({{n}_{s}}+{{n}_{b}})].\end{eqnarray} \tag{ 4 }$$

Without knowledge of the signal and background PDFs, the probabilities of an observed event representing signal or background are ${{n}_{s}}/({{n}_{s}}+{{n}_{b}})$ and ${{n}_{b}}/({{n}_{s}}+{{n}_{b}})$, respectively. These probabilities are combined with the per-event signal and background likelihoods ${{{\mathcal{S}}}_{i}}$ and ${{{\mathcal{B}}}_{i}}$ to obtain a single likelihood for each event i:

$$\begin{eqnarray}&&{{{\mathcal{L}}}_{i}}({{n}_{s}},{{n}_{b}})=\frac{{{n}_{s}}{{{\mathcal{S}}}_{i}}+{{n}_{b}}{{{\mathcal{B}}}_{i}}}{{{n}_{s}}+{{n}_{b}}}.\end{eqnarray} \tag{ 5 }$$

The product of the Poisson probability and the per-event likelihoods give an ensemble likelihood. We replace the background rate hypothesis n_b with the measured rate $\langle {{n}_{b}}\rangle$, which is well-measured in off-time data. Because the background rate varies with detector configuration due to the increasing size of the instrumented volume after each construction season, an ensemble likelihood is calculated for each configuration c. The overall likelihood is a function of the per-configuration signal rates $\left\{ {{({{n}_{s}})}_{c}} \right\}$ and is given by the product of the per-configuration likelihoods:

$$\begin{eqnarray}&&{\mathcal{L}}(\left\{ {{({{n}_{s}})}_{c}} \right\})=\mathop{\prod }\limits_{c}P({{({{n}_{s}})}_{c}}){{\mathop{\prod }\limits_{i=1}^{{{N}_{c}}}}^{{}}}{{{\mathcal{L}}}_{i}}({{({{n}_{s}})}_{c}}).\end{eqnarray} \tag{ 6 }$$

Our test statistic is the log-likelihood-ratio $T={\rm ln} \left[ {\mathcal{L}}(\left\{ {{({{{\hat{n}}}_{s}})}_{c}} \right\})/{\mathcal{L}}(\{0\}) \right]$, where the values $\left\{ {{({{{\hat{n}}}_{s}})}_{c}} \right\}$ maximize the likelihood and ${\mathcal{L}}(\{0\})$ is the likelihood for background-only. The test statistic can be written as

$$\begin{eqnarray}&&T=\mathop{\sum }\limits_{c}\left\{ -{{({{{\hat{n}}}_{s}})}_{c}}+\mathop{\sum }\limits_{i=1}^{{{N}_{c}}}{\rm ln} \left[ \frac{{{({{{\hat{n}}}_{s}})}_{c}}\;{{{\mathcal{S}}}_{i}}}{{{\langle {{n}_{b}}\rangle }_{c}}\;{{{\mathcal{B}}}_{i}}}+1 \right] \right\}.\end{eqnarray} \tag{ 7 }$$

We use a frequentist method to derive statistical significance and fluence upper limits from actual observations. The significance of an observed test statistic ${{T}_{{\rm obs}}}$ is the probability p of finding $T\geqslant {{T}_{{\rm obs}}}$ given background alone. To find this probability, pseudo-experiments are performed in which background-like data samples are generated by drawing from the reconstructed energy, direction and angular error distributions observed in off-time data. The resulting T distribution sets the significance of any single observation. We calculate fluence upper limits using a Feldman–Cousins approach (Feldman & Cousins 1998). Simulated events weighted to a given spectrum and normalization are added to pseudo-experiments; the exclusion confidence level (CL) is the fraction of pseudo-experiments which yield $T\geqslant {{T}_{{\rm obs}}}$.

When expressing constraints in terms of a quasi-diffuse flux, we assume that the 506 northern hemisphere bursts included in our four-year analysis are representative of ${{n}_{{\rm GRB}}}$ bursts per year that are potentially observable by existing satellites. Potentially observable bursts can go unseen because they are hidden by the Sun or Moon; they occur outside the field of view of any satellite or during satellite downtime; or, in this analysis, because they are in the southern sky. The extrapolation from actually observed bursts to potentially observable bursts is uncertain due to the differing fields of view and sensitivities of existing satellites, but here we assume ${{n}_{{\rm GRB}}}=667$—the same approach used in our previous publications (Abbasi et al. 2011a, 2012). Our results can be reinterpreted for a different supposed burst rate $n_{{\rm GRB}}^{\prime }$ by multiplying our reported flux values by $n_{{\rm GRB}}^{\prime }/667$. A potentially large population of nearby, low-luminosity GRBs (Liang et al. 2007) may contribute to an observable diffuse neutrino flux (Murase et al. 2006), but because they rarely trigger gamma-ray detectors, these bursts are not directly constrained by our analysis.

Our results are subject to systematic uncertainties in our neutrino signal simulation. Detector response and ice property uncertainties are accounted for by repeating the simulation with varied values for these inputs. Uncertainties due to muon propagation, Earth model parameters, and neutrino interaction cross sections have been studied in detail in previous work; these effects give a maximum uncertainty of ∼8% (Achterberg et al. 2007). The cumulative amplitude of these effects, which are included in all results presented in this paper, is spectrum-dependent, but generally the fluence corresponding to a given exclusion CL is increased by ∼10%.

In four years of data, we find a single neutrino candidate event correlated with a GRB, yielding a significance of p = 0.46. The burst and neutrino properties are listed in Table 1. Because this observation is not significant, we are able to improve upon our previously published upper limits (Abbasi et al. 2012). First, we consider a simple class of models for which each burst produces the same flux with a doubly broken power law spectrum in the Earth's frame, such that the total quasi-diffuse flux takes the form:

$$\begin{eqnarray}{{{\Phi }}_{\nu }}(E)={{{\Phi }}_{0}}\cdot \left\{ \begin{array}{ccccccccccccccc} {{E}^{-1}}\varepsilon _{b}^{-1} & E\lt {{\varepsilon }_{b}}, \\ {{E}^{-2}} & {{\varepsilon }_{b}}\leqslant E\lt 10{{\varepsilon }_{b}}, \\ {{E}^{-4}}{{(10{{\varepsilon }_{b}})}^{2}} & 10{{\varepsilon }_{b}}\leqslant E. \\ \end{array} \right.\end{eqnarray} \tag{ 8 }$$

We show exclusion contours for such models in Figure 1. Our treatment here is similar to that in Abbasi et al. (2012), but with the following modifications: (1) the inclusion of the second spectral break at $10{{\varepsilon }_{b}}$, and (2) the use of an updated Waxman-Bahcall prediction which accounts for more recent measurements of the UHECR flux (Katz et al. 2009) and typical gamma break energy (Goldstein et al. 2012) in accordance with the original prescription from Waxman & Bahcall (1997). The model by Ahlers et al. (2011) assumes that only neutrons escape from the GRB fireball to contribute to the UHECR flux; this scenario is strongly excluded by our limit. The Waxman-Bahcall model allows protons to escape the fireball as UHECRs directly without producing neutrinos, so it is not yet strongly excluded by our observations.

Zoom In Zoom Out Reset image size

Table 1.  GRB and Neutrino Properties for the Single Coincidence Observed in Four Years of Data

	GRB100718A	IceCube ν
Time	T100 = 39 s	T1 + 15 s
Angular separation		16°
Angular uncertainty	10.2°	1.3°
GRB fluence	$2.5\times {{10}^{-6}}\;{\rm erg}\;{\rm c}{{{\rm m}}^{-2}}$
ν energy		$\gtrsim 10\;{\rm TeV}$

Notes. The quoted GRB angular uncertainty is the Fermi GBM statistical error for this burst. In our analysis, the statistical error for GBM bursts is added in quadrature with a two-component estimated systematic error: 2.6° with 72% weight plus 10.4° with 28% weight (Paciesas et al. 2012). No GCN circular was produced for this burst; however, its observation was reported in the second Fermi GBM catalog (von Kienlin et al. 2014). The reconstructed energy of the product muon is 10 TeV. As discussed above, the neutrino energy may be larger

Download table as:  ASCIITypeset image

In models that predict per-burst neutrino spectra based on the details of the measured gamma-ray spectra, the fluence normalization scales linearly with the baryonic loading ${{f}_{p}}=1/{{f}_{e}}$, where f_e is the ratio of the kinetic energy in electrons to the total energy in protons within the fireball. In response to our previously published model-dependent limits (Abbasi et al. 2012), Baerwald et al. (2014) and others have observed that the relevant parameter space for f_p in the context of UHECR production depends on the energy range over which the baryonic loading is defined. We adopt the convention that f_p is defined over all proton energies—not just energies relevant to cosmic-ray production. Additional modeling corrections have also been studied. More detailed treatment of the $p+\gamma \to {{{\Delta }}^{+}}$ process leads to a fluence reduction while the use of numerical simulation to include other standard model $p\gamma$ interaction channels gives a fluence enhancement (Hümmer et al. 2012).

Using a wrapper for SOPHIA (Mücke et al. 2000) to calculate per-burst spectra, we evaluate exclusion contours in three scenarios. One is the standard fireball picture (Hümmer et al. 2012). Another is a photospheric model which moves the neutrino production to the photosphere, where the fireball transitions from optically thick to optically thin for $\gamma \gamma$ interactions (Rees & Mészáros 2005; Murase 2008; Zhang & Kumar 2013). Finally, we consider a Poynting-dominated flux model—Internal Collision-induced MAgnetic Reconnection and Turbulence, or ICMART (Zhang & Yan 2011)—in which internal shocks and particle acceleration take place at a much higher radius, typically ${{10}^{16}}\;{\rm cm}$ (Zhang & Kumar 2013).

For each model, we scan the parameter space for the bulk Lorentz factor of the fireball Γ and the baryonic loading ${{f}_{p}}=1/{{f}_{e}}$. In each case, we consider $1\lt {{f}_{p}}\lt 200$. For the standard and photospheric models, we test $100\lt {\Gamma }\lt 950$ while for ICMART, which varies more strongly with Γ, we test $50\lt {\Gamma }\lt 400$. The predicted spectra, summed over all analyzed bursts, are shown in Figure 2; the resulting exclusion contours are shown in Figure 3. Our results rule out some of the parameter space for f_p and Γ in regions that allow GRBs to be dominant UHECR sources. For very large values of Γ, IceCube would require a very long exposure to constrain the models. However, this region can be probed in other ways, such as by improved energy calibration of cosmic-ray measurements (Baerwald et al. 2014). We note that the constraints calculated here do not account for a possible enhancement to the high energy neutrino flux due to acceleration of secondary particles (Winter et al. 2014) or a distribution of differing Γ (He et al. 2012); nor do we attempt to account for a possible reduction of the neutrino flux if the brightest GRBs (in gamma-rays) have a smaller baryonic loading (Asano & Mészáros 2014).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

IceCube has recently established (Aartsen et al. 2014b, 2015) the existence of an astrophysical neutrino flux whose sources, like those of the UHECRs, are not yet known. This flux is established by neutrino events above expected backgrounds in the 10 TeV to few PeV range. The observed signal is consistent with an isotropic flux and can be parameterized as ${{{\Phi }}_{\nu }}(E)={{{\Phi }}_{0}}{{(E/{{E}_{0}})}^{-\gamma }}$. If E₀ is taken to be 100 TeV, then the best fit gives a per-flavor $\nu +\bar{\nu }$ normalization $E_{0}^{2}{{{\Phi }}_{0}}=2.06_{-0.3}^{+0.4}\times {{10}^{-8}}\;{\rm GeV}\;{\rm c}{{{\rm m}}^{-2}}\;{{{\rm s}}^{-1}}\;{\rm s}{{{\rm r}}^{-1}}$ and spectral index $\gamma =2.46\pm 0.12$ (Aartsen et al. 2015). To constrain the contribution to this flux from GRBs, we follow the prescription applied above for doubly broken power law spectra, except this time the simulation is weighted to unbroken spectra with $2\lt \gamma \lt 2.6$. Only simulated events above 10 TeV are considered; at very high energies, where the flux is already much smaller, no explicit cutoff is made. We find that the allowed GRB per-flavor $\nu +\bar{\nu }$ normalization, at 90% CL, is $E_{0}^{2}{{{\Phi }}_{0}}\sim 2\times {{10}^{-10}}\;{\rm GeV}\;{\rm c}{{{\rm m}}^{-2}}\;{{{\rm s}}^{-1}}\;{\rm s}{{{\rm r}}^{-1}}$. This constraint weakens only slightly with increasing γ. Thus potentially observable GRBs, as defined in this paper, contribute no more than ∼1% of the observed diffuse flux.

In this work, we have only considered a handful of possible neutrino spectra. In recognition of the large space of possible models to test, we now provide an online tool for calculating limits on alternative spectra. The subset of analyzed bursts to include as well as the per-burst spectra must be provided by the user. These choices are applied to our full analysis chain, and the results are sent back to the user via e-mail. Calculating limits in this way accounts for the details of our unbinned likelihood analysis, most importantly including the energy PDF; it also accounts for the one low-significance event which has been observed so far. See http://icecube.wisc.edu/science/tools for more details.

Using four years of IceCube data, we set the most stringent limits yet on GRB neutrino production, with a sensitivity improvement of ∼2× relative to our previous results. We constrain parts of the parameter space relevant to the production of UHECRs in the latest models. In addition to the work presented here, complementary analyses are underway. We are improving our acceptance with a search in the cascade channel, which is sensitive to the whole sky and to all neutrino interactions other than muon charged-current, as well a search for GRB-correlated high energy starting events, which has an extremely low background rate and therefore is sensitive to very early precursor or late afterglow neutrinos. Results from these searches will soon be published separately. In the absence of an emerging signal in the coming years, IceCube limits will increasingly constrain GRBs as dominant sources of UHECRs.

We acknowledge the support from the following agencies: U.S. National Science Foundation-office of Polar Programs, U.S. National Science Foundation-physics Division, University of Wisconsin Alumni Research Foundation, the Grid Laboratory Of Wisconsin (GLOW) grid infrastructure at the University of Wisconsin—Madison, the Open Science Grid (OSG) grid infrastructure; U.S. Department of Energy, and National Energy Research Scientific Computing Center, the Louisiana Optical Network Initiative (LONI) grid computing resources; Natural Sciences and Engineering Research Council of Canada, WestGrid and Compute/Calcul Canada; Swedish Research Council, Swedish Polar Research Secretariat, Swedish National Infrastructure for Computing (SNIC), and Knut and Alice Wallenberg Foundation, Sweden; German Ministry for Education and Research (BMBF), Deutsche Forschungsgemeinschaft (DFG), Helmholtz Alliance for Astroparticle Physics (HAP), Research Department of Plasmas with Complex Interactions (Bochum), Germany; Fund for Scientific Research (FNRS-FWO), FWO Odysseus programme, Flanders Institute to encourage scientific and technological research in industry (IWT), Belgian Federal Science Policy Office (Belspo); University of Oxford, United Kingdom; Marsden Fund, New Zealand; Australian Research Council; Japan Society for Promotion of Science (JSPS); the Swiss National Science Foundation (SNSF), Switzerland; National Research Foundation of Korea (NRF); Danish National Research Foundation, Denmark (DNRF).