Searches for Neutrinos from Large High Altitude Air Shower Observatory Ultra-high-energy γ-Ray Sources Using the IceCube Neutrino Observatory

Galactic PeV cosmic-ray accelerators (PeVatrons) are Galactic sources theorized to accelerate cosmic rays up to PeV in energy. The accelerated cosmic rays are expected to interact hadronically with nearby ambient gas or the interstellar medium, resulting in γ-rays and neutrinos. Recently, the Large High Altitude Air Shower Observatory (LHAASO) identified 12 γ-ray sources with emissions above 100 TeV, making them candidates for PeVatrons. While at these high energies the Klein–Nishina effect exponentially suppresses leptonic emission from Galactic sources, evidence for neutrino emission would unequivocally confirm hadronic acceleration. Here, we present the results of a search for neutrinos from these γ-ray sources and stacking searches testing for excess neutrino emission from all 12 sources as well as their subcatalogs of supernova remnants and pulsar wind nebulae with 11 yr of track events from the IceCube Neutrino Observatory. No significant emissions were found. Based on the resulting limits, we place constraints on the fraction of γ-ray flux originating from the hadronic processes in the Crab Nebula and LHAASO J2226+6057.


INTRODUCTION
The origin of the high-energy cosmic rays (CRs) is a longstanding question in particle astrophysics.The cosmic-ray spectrum is typically parametrized with an E −2.7 power law up to its "knee", E CR ∼ 3 PeV, above which the spectrum softens (Abbasi et al. 2018).The CRs with energy up to the "knee" are generally believed to be confined and accelerated in the Milky Way, implying the presence of PeV accelerators known as PeVatrons that accelerate the particles to energy > 1 PeV within our Galaxy.In addition to propelling CRs to PeV energies, Galactic PeVatrons are expected to emit γ-rays and neutrinos through associate pion decays (Ahlers et al. 2016).High energy γ-rays may also originate from leptonic processes like inverse Compton scattering and bremsstrahlung; however, at E γ 50 TeV, inverse Compton scattering is inefficient due to Klein-Nishina suppression (Cristofari 2021;Hinton & Hofmann 2009;Moderski et al. 2005).As inverse Compton scattering is the primary mechanism for electrons to emit γ-rays (Blumenthal & Gould 1970), γ-ray emission from leptonic processes is dramatically suppressed.Ultra-high-energy (UHE, > 100 TeV) γ-ray sources, especially those with observed spectra following a simple power law to at least tens of TeV, are suspected to be indicators of PeVatrons.Testify that the acceleration of hadrons is taking place for those TeV sources is of importance to pinpoint the PeV CR sources.
Supernova remnants (SNRs) are widely considered candidate Galactic PeVatrons.SNRs are generally believed to be the primary candidates of hadronic sources in the Galaxy, with particles being accelerated by diffusive shock acceleration (DSA) in expanding shocks.According to the SNR hypothesis, 10% of the total explosion energy of the supernova progenitor is converted into CRs, compatible with the observed intensity of CRs and their estimated confinement time in the Galaxy (Blandford & Ostriker 1978).The DSA mechanism and strong magnetic field amplification close to the shocks in SNRs also can explain the observed power-law spectrum of CRs (Cristofari et al. 2018).Apart from that, molecular clouds located within 100 parsec of the SNRs could be effective sites for > TeV γ-ray production (Gabici & Aharonian 2007).The highest-energy particles would quickly escape SNRs and diffuse faster in the interstellar medium (ISM), arriving at nearby molecular clouds earlier than the lower energy CRs and producing neutrinos and γ-rays there.CR might also be accelerated in young massive star clusters (YMCs) and pulsar wind nebulae (PWNe).Although the γ-ray emission from PWNe is widely expected to be leptonic, several studies suggest that the presence of relativistic hadrons in the pulsar winds might also be possible (Di Palma et al. 2017).From the Hillas criterion, the magnetic fields of PWNe could effectively trap particles with energy E 3 PeV (Hillas 1984).Therefore, a large amount of energy of PWNe could be stored in relativistic protons if they are continuously produced.In hadronic scenarios, CRs might be accelerated at the termination shocks (Di Palma et al. 2017) or by the shocks of SNRs which confined the diffuse nebulae (Ohira et al. 2018).Albert et al. (2021) reported that the TeV γ-ray spectrum of HAWC J1825-134 has no indication for a cut-off up to 200 TeV, implying that the source might be a possible PeVatron.HAWC J1825-134 may be associated with a giant molecular cloud, and the YMC located inside the molecular cloud is thought to accelerate CRs to PeV energies.Evidence of a PeVatron at the center of the Galaxy has also been reported by the HESS Collaboration by observing γ-rays from collisions of a young population of CRs with close by molecular clouds, with a hard γ-ray spectrum observed up to E γ ∼ 40 TeV (Abramowski et al. 2016).Detection of astrophysical neutrinos (E ν 100 TeV) would be a smoking gun for the existence of Galactic PeVatrons.The IceCube Neutrino Observatory has previously searched for Galactic neutrino emission from individual sources, stacked similar objects (PWNe, SNRs, etc.), and the Galactic plane (Aartsen et al. 2017(Aartsen et al. , 2019(Aartsen et al. , 2020a,b),b).Joint analyses with data from the ANTARES Neutrino Telescope were also performed to search for neutrinos from the Galactic plane and the southern sky (Albert et al. 2018(Albert et al. , 2020a)), exploiting the advantageous field of view of ANTARES and large volume and high statistics of IceCube.No significant Galactic neutrino sources were found from those searches, and the contribution of Galactic sources to the all-sky isotropic astrophysical neutrino flux has been constrained to ≤ 14% (Aartsen et al. 2017).
Recently, the LHAASO Collaboration reported the detection of a dozen Galactic E γ ≥ 100 TeV γ-ray sources with > 7σ detection (Cao et al. 2021a).This is the first LHAASO catalog and the largest list of UHE γ-ray sources published to date.Searching for neutrino emission from those LHAASO UHE sources will enable differentiation between leptonic and hadronic mechanisms that contribute to UHE γ-ray flux.In this letter, we perform a stacking analysis and a catalog search to probe the hadronic nature of these potential PeVatrons.We divide the LHAASO sources into 3 groups according to their possible association (all UHE sources, all possible PWNe, and all possible SNRs) in the stacking analysis, providing flux limits to constrain the total fraction of Galactic CRs that come from the UHE γ-ray SNRs and PWNe above 100 TeV.Section 2 briefly introduces the IceCube Observatory, the LHAASO detector arrays, the neutrino dataset, and the source list of 12 UHE sources.Analysis methods are shown in Section 3. In Section 4 and Section 5, we discuss and summarize the results.

IceCube neutrino sample
The IceCube Neutrino Observatory is a cubic kilometer array of 5160 digital optical modules (DOMs) located in the geographic South Pole (Abbasi et al. 2009).These DOMs are arranged along 86 readout and support cables ("strings") that instrument the glacial ice at a depth between 1.45 to 2.45 km.Each DOM is equipped with a 10-inch photomultiplier tube in a pressure-resistant sphere along with the necessary electronics for read-out of signal (Abbasi et al. 2009).These DOMs are designed to observe the Cherenkov light induced by energetic charged particles traveling in ice.These energetic charged particles are produced by charged-current (CC) interaction or neutral-current (NC) interaction from neutrinos of all flavors (Abbasi et al. 2010).The track-like (from CC interaction from muon neutrino) and cascade-like (from other NC and CC interaction) events are recorded by DOMs and analyzed in terms of the direction and energy of the primary neutrino.IceCube also includes DeepCore, a GeV-scale neutrino sub-detector (Abbasi et al. 2012) and IceTop, an air shower array (Abbasi et al. 2013), which are not used in this study.
The dataset used in this analysis comprises 11 years of IceCube track-like data from April 6, 2008 to May 29, 2020.Track-like events have a better angular resolution (a typical ∼ TeV neutrino has angular uncertainty ∼ 1 • ) compared to cascade-like events because of the long lever arm, providing better sensitivity to astrophysical point sources.The data taken from April 6, 2008 to May 13, 2011 are from the partially built detector with 40, 59 and 79 strings (IC40, IC59, and IC79) (Aartsen et al. 2020a) and the rest of the data are taken from the the full detector configuration of 86 strings (Aartsen et al. 2014).The dominant source of background is atmospheric neutrinos and atmospheric muons originating from air showers induced by cosmic-ray interactions in the atmosphere.Appendix A shows the details of the dataset used and the relevant energy range.
The dataset also incorporates the improved in-situ calibration of the single-photoelectron charge distributions (Aartsen et al. 2020c) and other improvements in reconstruction, further boosting the sensitivity relative to the previous analyses (Aartsen et al. 2020a).Sensitivity and 5σ discovery potential of the dataset are shown in Figure 1.

LHAASO ultra-high-energy γ-ray sources
The Large High Altitude Air Shower Observatory (LHAASO) is a complex of extensive air shower detector arrays aims to study γ-rays and cosmic rays.Located in the Sichuan province of China with an elevation of 4410 m a.s.l., LHAASO consists of the three components Kilometer Square Array (KM2A), Water Cherenkov Detector Array (WCDA), and Wield Field-of-view Cherenkov Telescope Array (WFCTA).LHAASO-KM2A focuses on Galactic γ-ray sources above 30 TeV in the Northern sky and has been operating since December 2019.The KM2A comprises uniformly distributed electromagnetic particle detectors (ED) and muon detectors (MD) over 1.3 km 2 .The vast area of scintillation counters together with the MD array would effectively suppress the background from CR-induced hadronic showers for γ-rays above ∼ 50 TeV, yielding a sensitivity of the full KM2A approaching ∼ 1% for a Crab-like source with a 1-year observation at 100 TeV (Liu et al. 2016;Cui et al. 2014).The angular and energy resolution of the KM2A are 0.4 • and 28% at 30 TeV, respectively.We note that at 1 PeV, the angular resolution of the KM2A would reach 0.2 • (Cao et al. 2019).
The half-completed LHAASO-KM2A has already detected twelve UHE γ-ray sources above 100 TeV with statistical significance ≥ 7σ, from December 2019 to October 2021 (Cao et al. 2021a, Table 1).All these potential PeVatrons appear to be located in our Galaxy.Among these twelve sources, three sources (eHWCJ1825-134, eHWCJ1907+063, and eHWCJ2019+368) are also in the HAWC's high energy source list (Abeysekara et al. 2019a) with conclusive evidences of ≥ 100 TeV emission.Except for LHAASO J0534+2202, which is the Crab Nebula, all the other LHAASO UHE sources have no firm association, with several known TeV counterparts nearby.Table 2 in Cao et al. (2021a) lists all the possible TeV associations close to these twelve UHE sources, suggesting that majority of them might be associated with PWNe or SNRs.

Catalog search
A maximum likelihood ratio test is performed following the technique described in Abbasi et al. (2011) and Braun et al. (2008).The same method was also used in the previous IceCube 10 years time-integrated point source search (Aartsen et al. 2020a).The first LHAASO catalog consists of 12 UHE γ-ray sources detected above 100 TeV with > 7σ (Cao et al. 2021a).Only two of the LHAASO sources (LHAASO J0534+2202 and LHAASO J2108+5157) are considered point sources with regard to LHAASO's angular resolution.For the remaining ten sources, significant angular extensions in γ-ray are detected by LHAASO and a 0.3 degree Gaussian emission template is used to measure the γ-ray flux.Hence, a 0.3 degree Gaussian extent is assumed for the remaining ten sources to match assumptions used in the LHAASO analysis (Cao et al. 2021a).We model the spatial signal likelihood as a 2D Gaussian in tangent plane coordinates which is valid for a sufficiently small angular extension.The angular uncertainty of each event is obtained directly from the reconstruction.Source extensions and estimated angular uncertainties for events are added in quadrature when estimating spatial likelihoods for each event.We assume a power-law spectrum with spectral index Γ for the signal hypothesis.The normalization and spectral index 1 < Γ < 4 of the signal flux are optimized in a maximum likelihood analysis.A detailed description of the likelihood can be found in Appendix B.
Figure 1 shows the sensitivity and 5σ discovery potential as a function of declination of this search method for Γ = 2 and Γ = 3. Sensitivity is the neutrino flux required for 90% of the trials to yield a p-value less than 0.5.The 5 σ discovery potential is the neutrino flux required for 50% of trials with simulated signal to yield a p-value less than 2.87 × 10 −7 (5 σ significance in one-tail Gaussian test).In our analyses, we use 50 TeV as the pivot energy when calculating the neutrino flux, given that 100 TeV γ-rays from pion decay would be accompanied by neutrinos at 50 TeV (Ahlers & Murase (2014), see Appendix C for more detail).

Stacking search
A total of six stacking analyses are performed combing three catalogues -all LHAASO sources, sources potentially associated with SNRs, and sources potentially associated with PWNe -and two weighting schemes -equal neutrino flux at Earth and neutrino flux proportional to the γ-ray flux at 100 TeV reported by LHAASO (Cao et al. 2021a).The source hypothesis in the stacking analysis is the same as in the catalog search, with sources considered to have 0.3 degree Gaussian extension except for LHAASO J0534+2202 and LHAASO J2108+5157.According to the association proposed by Cao et al. (2021a), nine sources are included in the PWN stacking search, and six sources are included in the SNR stacking search.There are four gamma-ray sources in common between the SNR and PWN stacking searches.Table 1 shows the source type of the potential association of the LHAASO sources.We assume that all the sources in each stacking scheme have a universal unbroken power-law energy spectrum with the same spectral index.The spectral index and the total number of signal neutrinos are fitted as free parameters during the maximization of the likelihood.For the equal-weight cases, the relative weight of each source is 1 N where N is the number of sources in the catalog.For the flux-weighted cases, the relative weight of each source i is φi j φj where φ i is the LHAASO γ-ray flux at 100 TeV (Cao et al. 2021a).We generate 10 5 background trials by scrambling the right ascension coordinates of the track sample.The stacking search is performed on each background trial, identifying the source with the lowest pre-trial p-value in each.The distribution of these p-values is used to correct for the look-elsewhere effect.The post-trial p-value is computed by comparing the pre-trial p-value obtained from the real data and the distribution of the lowest pre-trial p-values (Aartsen et al. 2020a).

Catalog search and Stacking search results
Table 1 shows the result of the catalog search, including the best-fit number of neutrino events, best-fit power law index and pre-trial p-value.The 90% confidence interval (C.I.) flux upper limit for each source at 50 TeV assuming a spectral index Γ = 2 is shown in the last column.We find no significant correlation between the observed neutrinos and the newly identified LHAASO UHE (E γ > 100 TeV) sources in our searches.Figure 1 shows 90% flux limits for Γ = 2 and Γ = 3, comparing with the sensitivity from IceCube's 11 years dataset and ANTARES's 13 years dataset (Illuminati 2021).Appendix C explains the detail of the estimation of the neutrino flux from the observed gamma-ray flux (Ahlers & Murase 2014).Here we assume that the full γ-ray flux is produced via proton-proton interaction.TeV −1 cm −2 s −1 ) reported by LHAASO (Cao et al. 2021a).
The neutrino 90% C.I. flux upper-limit (φ 90% ) is parameterized as: The relevant energy range of the flux upper limit can be found in Appendix C. The two smallest pre-trial p-values of 0.046 and 0.045 observed for two sources correspond to a post-trial p-value of 0.42 with the assumption that those 12 sources are independent given their large separation.
The two lowest p-value sources are LHAASO J1908+0621 and LHAASO J2018+3651.LHAASO J1908+0621 is thought to be the UHE counterpart of MGRO J1908+06, which is also the lowest p-value Galactic source in the source list of the IceCube 10-year point-source search (Aartsen et al. 2020a).LHAASO J2018+3651 might be associated with MGRO J2019+37, and an HII region, Sh 2 − 104, is located to the west.The region Sh 2 − 104 comprises several YMCs with diffuse X-ray emission probably from the colliding winds of young stars (Gotthelf et al. 2016).
Ten of the LHAASO UHE sources with δ < 42 • are in the source list for the ANTARES 13-year point-source search (Illuminati 2021).LHAASO J2018+3651 is the only LHAASO source with reported pre-trial p-value of 0.14.LHAASO J1908+0621, LHAASO J1929+1745, and LHAASO J2032+4102 have p-values 0.06 in the searches of Huang & Li (2021) with 10 years of IceCube track events (Abbasi et al. 2021).The pre-trial p-values of LHAASO J1929+1745 and LHAASO J2032+4102 are ∼ 0.03 and ∼ 0.04 with an extension of 0.6 degree, while that of LHAASO J1908+0621 is ∼ 0.06 for the point-source hypothesis.LHAASO J2032+4132 is the third significant source in our catalog search with a pre-trial p-value of 0.075.The differences between the significance level for LHAASO J2018+3651 and LHAASO J1929+1745 in our analysis and that in Huang & Li (2021) might result from the differences in datasets' livetime and in track reconstruction approaches or from the limitations of the public dataset.
Table 2 shows the result of the stacking searches, including the best-fit total number of signal events, best-fit power law index, and pre-trial p-value for each stacking scheme.The 90% C.I. flux upper limit at 50 TeV for each tested weighting scheme assuming a spectral index of 2 is shown in the last column.No significant emission is found from all 6 stacking schemes, with the lowest p-value results found from all possible LHAASO PWNe and all LHAASO UHE sources weighted equally.The 90% flux upper limits for each weighting schemes tested as a function of spectral index is shown in Figure 2.
The smallest pre-trial p-value of 0.022 observed in the two stacking hypotheses corresponds to a post-trial p-value of 0.06.

Hadronic constrains
In this section, we try to constrain the fraction of γ-ray emission from hadronic processes in the LHAASO sources.To obtain meaningful hadronic constraints based on the non-observation of neutrinos, we look up the observed γ-ray indexes for LHAASO sources from the literature in order to estimate a more realistic γ-ray spectrum for the 90% C.I. upper limits.Considering the observed TeV γ-ray spectra and comparing the theoretically predicted neutrino flux (the estimation detail is described in Appendix C) with 90% upper-limit flux calculated using the TeV spectral parameters, we are able to constrain the hadronic contributions to the γ-ray flux for three of the LHAASO UHE sources: LHAASO J0534+2202, LHAASO J1849-0003, and LHAASO J2226+6057.The hadronic constraints, neutrino upper limits, and expected neutrino fluxes for the three LHAASO sources are shown in Table 3.We note that the constraints illustrated here are the best possible bounds we could obtain as we assume that the γ-ray emissions are completely hadronic.The limits would be weaker if part of the UHE γ-rays from the sources originated from leptonic processes.
The strongest hadronic constraint is put on LHAASO J2226+6057 with a limit of 47%, with the assumption that the neutrino spectrum follows a simple power law with an index of 3.01 reported by LHAASO.The UHE source may be associated with SNR G106.3+02.7,suggested as a potential PeVatron by Albert et al. (2020b) and Amenomori et al. (2021).The observations of HAWC and Tibet ASγ for the SNR indicate an unbroken hard TeV spectrum extends up to at least 100 TeV.Amenomori et al. (2021) suggested that the VHE protons might be accelerated and trapped in the SNR given the hard spectral index.On the other hand, the UHE γ-ray flux of J2226+6057 detected by LHAASO from 20 to 500 TeV is better described with log-parabola rather than a single power-law (Cao et al. 2021a).However, the hadronic component of J2226+6057's γ-ray flux cannot be constrained when applying the log-parabola model to calculate the neutrino flux limits.
For LHAASO J0534+2202, which is the UHE counterpart of the Crab Nebula, we constrain the hadronic UHE γ-ray emission to 59% of the total, assuming a power law of Γ = 3.12 as reported by LHAASO.Applying the log-parabola Table 3. Table of TeV spectral parameters and the corresponding hadronic constraints, neutrino upper limits, and expected neutrino fluxes at 50 TeV.The parameter φ 90% represents the neutrino 90% C.I. flux limits parameterized as : The parameter φ Predicted is the predicted neutrino flux with the assumption of γ-ray s are entirely hadronic.Hadronic constraints correspond to the ratio φ 90% /φ Predicted at 50 TeV calculated with fixed spectral parameters.The TeV spectral and morphology information (columns 3 − 6) is taken from model by HAWC, the hadronic constraint on the Crab is 84%, compared to 70% reported previously by Huang & Li (2021).Although a 1.12 PeV γ-induced shower was detected by LHAASO-KM2A from the direction of the Crab, the hadronic origin of TeV γ-rays from it is considered to be less likely (Cao et al. 2021b).If the UHE emission from the Crab is hadronic, CRs should be accelerated in the pulsar's magnetosphere or at the wind termination shock.An alternative scenario proposed by Cheng et al. (1990) suggests that the CRs could be accelerated in the outer gap of the pulsar.
We could only marginally constrain the hadronic fraction of the γ-ray flux for LHAASO J1849−0003 to 94%, with a power-law index of 1.99 (Huang & Li 2021).

CONCLUSIONS AND OUTLOOK
In this study, we have performed catalog and stacking searches on twelve LHAASO UHE γ-ray sources.We have found no significant neutrino emission, with the smallest post-trial p-value of 0.06 in PWN equal-weight stacking search.The resulting 90% neutrino flux upper limits for an E −2 spectrum are given.For LHAASO J1849-0003, LHAASO J0534+2202, and LHAASO J2226+6057, we are able to constrain the hadronic γ-ray emission tailored to the available spectral information of these sources.Joint fits between γ-rays and neutrinos (Fan et al. 2021) and physics-motivated lepto-hadronic modeling of the spectrum of each individual source may provide a better and more accurate limit.Future detectors should also improve our understanding of these potential PeVatrons.The planned IceCube-Gen2 detector (Aartsen et al. 2021), with an increased effective area relative to IceCube and improved angular resolution, is expected to become up to five times more sensitivity for Galactic neutrino sources compared to IceCube.
The IceCube collaboration acknowledges the significant contributions to this manuscript from Yu-Ling Chang, Kwok Lung Fan and Michael Larson.We also acknowledge support from: USA -U.S.   4 shows the event counts, livetime, start and end date of each detector configuration of the track-like neutrino dataset used in the analysis.Figure 3 shows the distribution of events as a function of reconstructed declination and estimated energy and also the 90% energy range for data and MC signal with a spectral index of 2 and 3.In pp interactions, protons interact with nearby material and produce charged and neutral pions with a ratio of about 2 to 1 (Ahlers & Murase 2014).Charged pion decay to neutrinos and other products, and the resulting neutrinos carry roughly 1/4 the energy of the pion.Neutral pion decay into two γ-rays, each with 1/2 of the energy (Ahlers & Murase 2014).On Earth, the flux observed is where J γ (E γ ) and J να (E ν ) are the gamma-ray differential particle flux and neutrino differential particle flux for neutrino favor α .Assuming the absorption of gamma-rays is negligible and equal flux from each neutrino flavor due to oscillations, relationship between γ-rayand neutrino fluxes is then )

Figure 1 .
Figure1.The blue solid and blue dashed lines are the per-source IceCube 5σ discovery potential and sensitivity for 0.3 degree extension and Γ = 2 spectrum (a), and the orange solid and dashed lines are those for 0.3 degree extension and Γ = 3 spectrum (b) of this analysis.The green solid and green dashed lines on the left plot are the 5σ discovery potential and sensitivity of a E −2 spectrum of ANTARES point source analysis(Illuminati 2021).The red stars represent the neutrino flux predicted from the LHAASO measurement if photon flux is assumed to be 100% of hadronic origin(Ahlers & Murase 2014).The blue and orange triangles represent the 90% flux limits obtained from this analysis for E −2 and E −3 spectrum, respectively.The green crosses are the ANTARES 90% flux limits(Illuminati 2021).

Figure 2 .
Figure 2. Our 90% flux upper limits on neutrino emission at 50 TeV from (a) all 12 LHAASO sources, (b) sources that could be associated with SNRs, and (c) sources that could be associated with PWNe.The orange line represents the equal-weight case and the blue line represents the flux-weighted case.The hadronic flux prediction assumes that the flux measured by LHAASO is of purely hadronic origin.

Figure 3 .
Figure 3.The 2D distribution of events in data taken in 2019 with full 86 strings as a function of reconstructed declination and estimated energy.The 90% energy range for the data is shown in the black solid line.Simulated astrophysical signal Monte Carlo (MC) for an E −2 and an E −3 spectrum is shown in purple and orange respectively as a guide for the relevant energy range of IceCube

Table 1 .
Table of best-fit parameters with corresponding test statistic (TS) and p-value of catalog search.The γ-ray flux is the γ-ray flux of the source at 100 TeV in Crab units (CU) (6.1 × 10 −17

Table 2 .
Table of best-fit parameters for the stacking analyses.The 90% flux upper-limit (φ 90% ) is parameterized as : National Science Foundation-Office of Polar Programs, U.S. National Science Foundation-Physics Division, U.S. National Science Foundation-EPSCoR, Wisconsin Alumni Research Foundation, Center for High Throughput Computing (CHTC) at the University of Wisconsin-Madison, Open Science Grid (OSG), Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS), Frontera computing project at the Texas Advanced Computing Center, 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; European Union -EGI Advanced Computing for research; 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, Carlsberg Foundation, and European Commission; 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.