Seasonal Variation of Multiple-Muon Cosmic Ray Air Showers Observed in the NOvA Detector on the Surface

We report the rate of cosmic ray air showers with multiplicities exceeding 15 muon tracks recorded in the NOvA Far Detector between May 2016 and May 2018. The detector is located on the surface under an overburden of 3.6 meters water equivalent. We observe a seasonal dependence in the rate of multiple-muon showers, which varies in magnitude with multiplicity and zenith angle. During this period, the effective atmospheric temperature and surface pressure ranged between 210 K to 230 K and 940mbar to 990mbar, respectively; the shower rates are anti-correlated with the variation in the effective temperature. The variations are about 30% larger for the highest multiplicities than the lowest multiplicities and 20% larger for showers near the horizon than vertical showers.

We report the rate of cosmic ray air showers with multiplicities exceeding 15 muon tracks recorded in the NOvA Far Detector between May 2016 and May 2018. The detector is located on the surface under an overburden of 3.6 meters water equivalent. We observe a seasonal dependence in the rate of multiple-muon showers, which varies in magnitude with multiplicity and zenith angle. During this period, the effective atmospheric temperature and surface pressure ranged between 210 K to 230 K and 940 mbar to 990 mbar, respectively; the shower rates are anti-correlated with the variation in the effective temperature. The variations are about 30% larger for the highest multiplicities than the lowest multiplicities and 20% larger for showers near the horizon than vertical showers.

I. INTRODUCTION
Several experiments have observed seasonality in the rate of muons from cosmic ray air showers which depends on the density profile of the earth's atmosphere. As the density changes, the relative numbers of pions and kaons which interact or decay in the air showers change, altering the observed rate of muons. During the summer months, when temperature is high and the density is lowest, meson decays are more probable, which leads to an expected peak in the muon rates. This summer peak has been confirmed by many experiments in single-muon air showers [1][2][3][4] and is explained by existing models [5].
However, existing models fail to fully explain the observations of multiple-muon air showers. The DECOR [6] and GRAPES [7] experiments have reported peak rates of multiple-muon showers in the winter in detectors close to the surface, opposite to expectations and observations for single-muons. DECOR attributed its observation to geometric effects arising from altitude differences in meson production. The MINOS experiment [8] observed a winter peak in two underground detectors, with minimum muon energies of 60 GeV and 700 GeV, and showed that at a depth of at least 225 meters water equivalent the effect from altitude differences suggested by DECOR * Now at Oak Ridge National Laboratory was too small to fully explain the observation. NOvA also previously reported a peak rate in winter of multiplemuon cosmic ray air showers using its Near Detector [9] located at the same depth as one of the MINOS detectors with a threshold energy of 60 GeV.
The NOvA Far Detector is located near the surface where no seasonal variation is expected for low energy, single-muon air showers [10]. However, the Far Detector has a top surface area which is 15 times larger than the Near Detector making it sensitive to much higher multiplicity showers.
In this paper, we report the observation of a winter maximum of multiple-muon air showers using NOvA's Far Detector. Since no quantitative models for multiplemuon air showers reproduce the effects we observe, the seasonal effect will be quantified using two different methods. First, the rate of multiple-muon air showers is compared to the temperature and surface pressure of the atmosphere above the detector site and, second, by fitting the observed muon rate to a cosine function. We also show how the strength of this observation varies with observed muon multiplicity and arrival direction in a surface detector for the first time.

II. THE NOVA FAR DETECTOR
The NOvA Far Detector is a 14 kt sampling calorimeter, 15.5 m × 15.5 m × 59.8 m in size, segmented into 4 cm × 6 cm × 15.5 cm channels. The channels are arranged into alternating horizontal and vertical planes. The detector was designed to detect neutrino interactions in the NuMI beam from Fermi National Accelerator Laboratory [11]. It is located on the surface in Minnesota near the U.S.-Canada Border at (48.4 • N, 92.8 • W). The detector has been operating with more than 97% up-time efficiency since 2014. This analysis samples 15% of the total cosmic ray data set collected between May 2016 and May 2018. The detector design and detection mechanism are described in [12].
The detector sits just below the surface level beneath an overburden to shield the detector from cosmic ray photons and electrons. It consists of 1.2 m of concrete and 15 cm of barite rock giving a total of 3.6 meters water equivalent. Three sides of the building are surrounded by a sloped berm of granite rock at 30 • to the surface. This shielding is not present north of the detector where the detector assembly hall is located. Above the horizon, this overburden adds an additional muon energy threshold of E thr cos θ ≈ 1.5 GeV, where θ is the zenith angle, to reach the detector. On average, 10 billion cosmic ray muons traverse the detector each day.
Data from multiple-muon showers are recorded for analysis [13] whenever the detector records total visible energy in excess of approximately 20 GeV of energy deposited in a 50 µs readout distributed among at least 120 of the detector's total of 343,968 channels. A typical muon with θ = 30°traversing the center of the detector will deposit around 2.5 GeV of visible energy.

III. ATMOSPHERIC AND MUON DATA
The signal of muon air shower events in the detector is a large number of coincident, parallel tracks. Fig. 1 shows the signal topology of a multiple-muon shower recorded in the detector. Reconstruction of these showers begins by isolating the time range containing the activity of interest from other detector activity and suppressing isolated detector hits which do not contribute to tracks. A Hough transform determines the overall shower angle in each view of the detector. These angles seed the construction of individual muon tracks. These algorithms produce a zenith angle, azimuthal angle, and multiplicity assignment for each air shower. The multiplicity reported here is the observed multiplicity within the detector. Because air showers can be much larger than the surface area of the NOvA detector, no attempt was made to estimate the true multiplicity.
The reconstruction was optimized and validated using air-shower simulations based on COsmic Ray SImulations for KAscade (CORSIKA) [14] with a range of primary cosmic ray energies in order to explore performance on showers with a variety of multiplicities. The reconstructed multiplicity is within ±1 of the true multiplicity in the detector for 70% of showers and within ±4 for 95% of showers. The most common reconstruction failure is due to muon tracks that are nearly overlapping in space in one view and are treated as a single track.
We apply a number of data selection requirements to ensure uniform detector acceptance for air showers. While it is possible for the NOvA detectors to operate with only a subset of components active, we only analyze data-taking periods when the detector was completely active to ensure continuity of the muon tracks. Air showers travelling nearly parallel to the detector planes (up-down, east-west) are removed as the reconstruction cannot produce complete tracks for these orientations. Candidate events with very large energy deposits per hit likely contain large hadron showers from the overburden. These events are either not associated with air showers or hinder our ability to reconstruct the air shower direction or multiplicity and are removed from the sample. Showers with a reconstructed multiplicity less than 15 are removed to avoid trigger inefficiencies at low multiplicities and to ensure a uniform efficiency over the analysis sample. From a CORSIKA simulation, the typical primary energy to make 15 muons in the detector is 30 TeV to 100 TeV. The reconstructed multiplicity and zenith angle of all selected showers can be seen in Fig. 2.
The livetime used to compute the shower rates is recorded by data acquisition processes, which monitor the trigger data streams. Fig. 3 shows the rate of cosmic ray air showers for the full analysis period. The rates reach their maximum values during the winter months. In total, 7.64 × 10 6 multiple-muon showers with an average multiplicity of 28 are analyzed. The showers have an average rate of R µ = 1.09 Hz.
The measured rate of multiple-muon showers will be compared to atmospheric conditions. We use atmospheric data provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) [15]. These global data are provided four times per day with a resolution of 0.75°in latitude and longitude. In this analysis, we average the temperature at the four points nearest to the detector location that are closest in time to the air shower.
The effective temperature of the atmosphere above the detector is a weighted average of the temperature recorded at pressure levels ranging from 1 hPa to 1000 hPa over the depth of the atmosphere with higher weight given to altitudes with muon production [5]. This model is only valid for muons from leading pions and kaons produced in the primary interaction and is applicable for single-muon events. However, this computation approximates what the effective temperature would be for multiple-muon events. The average effective temperature for this analysis is T eff = 223 K. The surface pressure data is also reported with an average value of P = 968 mbar. Fig. 3 shows the variations in these quantities over the analysis period. Fig. 3 shows a clear seasonal variation in the rate of multiple-muon air showers with peak rates in both winters. We employ two methods to quantify the significance of these peaks.

IV. SEASONAL EFFECTS
The first method correlates the atmospheric temperature and the muon rates. We compute a temperature correlation coefficient, α T [5]: where ∆R µ is the difference from the mean rate of multiple-muon showers and similarly for the effective temperature, T eff . This effective temperature model has been shown to be closely correlated with the rate of single-muon showers. The value of α T is dependent on the threshold energy of detected muons and thus on the depth of the detector. Surface detectors are expected to exhibit no temperature dependence with α T ≈ 0 for single-muon showers. However, the model does not accurately explain the development of multiple-muon showers [10] where many competing effects contribute to the observed rate. Despite this limitation, the magnitude of α T still demonstrates a correlation between temperature and observed multiple-muon rate. Applying a linear fit to the multiple-muon rates as a function of the effective temperatures from Fig. 3, results in a temperature correlation coefficient of α T = −1.14 ± 0.02. The atmospheric pressure at the surface can affect the survival probability of low-energy muons as they approach the detector and alter the observed rate [16]. The barometric coefficient, β, is measured by: where ∆P is difference in pressure from the mean. The barometric coefficient has been shown to accurately relate rates of single-muon showers to the pressure [2]. A fit between the multiple-muon rates and the surface pressures from Fig. 3, yields a barometric correlation coefficient of β = (−0.08 ± 0.01) %/mbar.
The value of α T is statistically different from zero and negative, as expected for an anti-correlation between temperature and multiple-muon shower rates. β is also negative, which is as expected for the single-muon case for a surface detector [2]. However, the relationship between the pressure and shower rate is found to be nonlinear, suggesting additional corrections are needed for multiple-muon showers. The measurement of α T was repeated after correcting the observed shower rates for the changes in pressure, but the difference in the new value was negligible.
The second method used to quantify the significance of the seasonal effect is the amplitude of a cosine fit to the data. The temperature of the atmosphere is not expected to strictly track a cosine and, thus, neither does the observed rate of events. However, the modulations are expected to be periodic, with a period close to, but not exactly, one year, with higher temperatures in the summer and colder in the winter. The amplitude of the fit demonstrates the strength of the effect independently from the used model. Fig. 4 shows the average percent change in the multiple-muon rate as a function of month of the year, which is found to be more sinusoidal than the more finely binned version. The function used for the fit in Fig. 4 is where V 0 is the function average, V is the relative change, T is the period of modulation fixed at one year, and the phase φ is the time of maximum. The best fit values are V 0 = 0.00 ± 0.01%, V = 5.86 ± 0.05%, and φ = 0.88 ± 0.02 months. The phase implies a peak rate around January 27. The value of V gives a qualitative measurement for how much the multiple-muon showers rate varies throughout any given year. The models used to measure the above quantities provide imperfect descriptions of the relationship between atmospheric conditions and the rate of multiple-muon air showers. However, we can still consider how systematically changing the detector observables affects the measured quantities to determine if a systematic effect could give rise to the observed variations. Here, we discuss the largest such effects.
The temperature and pressure measurements made by the ECMWF have an associated systematic uncertainty of ±0.31 K and ±1 mbar, respectively [3]. These uncertainties are the largest known systematic uncertainties in the measurement and are already included in the measurement of α T and β in the fits above.
The detector data acquisition system writes data files in periods of up to 2.5 minutes, depending on the data triggers operating at the time. The detector events are written to the files in the order the trigger issues verdicts, so the data may appear out of order, and the files are closed when they reach a fixed file size [13]. Occasionally, data at the end of one file and the start of the next will be misordered. The livetime in such cases may be overestimated by as much as 1 s. Systematically increasing all livetimes by 1 s has less than 0.5% effect on the values of V and α T and 1.5% on the value of β.
The detector is a rectangular prism with length about four times its width and height. Showers directed at the two smaller faces of the detector will have a smaller visible cross section of detector and will, thus, have lower multiplicities. To account for this geometric effect, such showers have both their multiplicities and rate systematically increased by a factor of four. The values of α T and V decrease by less than 0.2% and the value of beta increases by only 1%. Since the effect is small, this is not used to correct angular effects in the data.
The detector electronics are sensitive to the temperature and humidity of the atmosphere. As a result, the electronics are noisier when the operating temperature is warmer, and up to 5% more noise hits are observed in the summer months. However, noise hits have much lower ADC counts than those made by the signal muon tracks and were not found to have any impact on the reconstruction of multiple-muon events.
None of the considered effects are large enough to have artificially created the observed winter maximum in the multiple-muon rate. Additionally in the following section, these effects cancel when considering the relative measurements between the bins of multiplicity and the bins of zenith angle.

V. MULTIPLE-MUON OBSERVABLES
The two methods in the previous section demonstrate a clear seasonal variation with a peak during the winter.  As in the Near Detector analysis [9], we observe that the strength of this effect changes under different shower observables that act as a proxy for the primary cosmic ray energy. Here we examine the changes in the strength of the seasonal effect with the multiplicity and zenith angle of the shower. The multiple-muon showers are divided into five sets depending on their multiplicities and zenith angles, respectively. These divisions, shown in Fig. 2, contain nearly equal numbers of showers. Fig. 5 shows the multiple-muon data and sinusoidal fit for each multiplicity bin as a function of month of year. The sinusoidal fit amplitude, temperature coefficient, and barometric coefficient are measured within each multiplicity bin and reported in Table I. Both the cosine fit amplitude and temperature coefficient demonstrate a stronger seasonal dependence at higher multiplicities.
We perform a similar analysis for each zenith angle bin. The results are reported in Table II. The bins of zenith angle nearest the horizon exhibit the strongest seasonal effect. The table also shows the average multiplicity of showers within each bin. The showers coming from nearest the horizon also have the lowest average multiplicities, which would be expected to exhibit the weakest seasonal variation from Table I. We observe that the most ver-

VI. SUMMARY
We observed that the rate of multiple-muon cosmic ray air showers in a detector near the surface presents a seasonal variation with a peak rate in the winter. Additionally, we showed that this effect is dependent on the primary cosmic ray energy by looking at two detector observables, the multiplicity and zenith angle. A stronger seasonal effect is seen for air showers with higher multiplicities or zenith angles near the horizon. The amplitude of the seasonal modulation grows by 30% from the lowest multiplicities to the highest multiplicities and by 20% from the most vertical to the most horizontal showers considered.
For surface detectors where the threshold energy for detection is low, the production altitude of muons in cosmic ray air showers is of the same magnitude as the muon decay length. For example, a typical muon reaching the surface begins with 5 GeV of energy at production and will traverse on average 31.1 km before decaying. Muon production occurs around 15 km to 20 km, with higher altitudes in the summer months when the atmosphere is expanded. Thus, the longer muon path length in summer months would give the produced muons a higher chance to decay before reaching the surface and reduce the number of observed muons. The effects of particle decay on the seasonal rate of muons could be confirmed using Monte Carlo simulation such as CORSIKA [14]. However, this effect is negligible in underground detectors where the muon energies are at least ten times larger and cannot explain all observations of seasonal variations for multiple-muon showers.
An interesting continuation of this study will be the inclusion of low multiplicity air showers from another detector trigger. The two datasets could be combined to see if there is a threshold where the seasonal behavior for multiple-muons inverts as in underground detectors or flattens as expected for single-muons. Additionally, comparisons to Monte Carlo simulation could be used to trace detector observables back to the primary cosmic ray energy.
This document was prepared by the NOvA collaboration using the resources of the Fermi National Accelera-