Modulations of the Cosmic Muon Signal in Ten Years of Borexino Data

We have measured the flux of cosmic muons in the Laboratori Nazionali del Gran Sasso at 3800\,m\,w.e. to be $(3.432 \pm 0.003)\cdot 10^{-4}\,\mathrm{{m^{-2}s^{-1}}}$ based on ten years of Borexino data acquired between May 2007 and May 2017. A seasonal modulation with a period of $(366.3 \pm 0.6)\,\mathrm{d}$ and a relative amplitude of $(1.36 \pm0.04)\%$ is observed. The phase is measured to be $(181.7 \pm 0.4)\,\mathrm{d}$, corresponding to a maximum at the 1$^\mathrm{st}$ of July. Using data inferred from global atmospheric models, we show the muon flux to be positively correlated with the atmospheric temperature and measure the effective temperature coefficient $\alpha_\mathrm{T} = 0.90 \pm 0.02$. The origin of cosmic muons from pion and kaon decays in the atmosphere allows to interpret the effective temperature coefficient as an indirect measurement of the atmospheric kaon-to-pion production ratio $r_{\mathrm{K}/\pi} = 0.11^{+0.11}_{-0.07}$ for primary energies above $18\,\mathrm{TeV}$. We find evidence for a long-term modulation of the muon flux with a period of $\sim 3000\,\mathrm{d}$ and a maximum in June 2012 that is not present in the atmospheric temperature data. A possible correlation between this modulation and the solar activity is investigated. The cosmogenic neutron production rate is found to show a seasonal modulation in phase with the cosmic muon flux but with an increased amplitude of $(2.6 \pm 0.4)\%$.


Introduction
Cosmic muons are produced mainly in the decays of kaons and pions that originate from the interaction of primary cosmic rays with nuclei in the upper atmosphere [1]. For detectors situated deep underground, the flux of cosmic muons is strongly reduced. Only muons surpassing a certain threshold energy E thr contribute, while lower-energetic muons are absorbed in the rock overburden. At great depths, the residual high energy muons are bound to be produced by parent mesons that decay in flight without any inelastic interactions and without elastic interactions of large momentum transfer before the decay. As a consequence, the density and temperature variations of the upper atmosphere that alter the mean free path of the decaying mesons introduce, in first approximation, a seasonal modulation of the underground muon flux, which has been investigated for many decades [2]. Several underground experiments located at the Laboratori Nazionali del Gran Sasso (LNGS) in Italy such as MACRO [3], LVD [4,5], Borexino [9], and GERDA [6] and at other underground sites, e.g. IceCube [7] or MINOS [8], have studied this phenomenon. Compared to the previously published investigation based on four years of Borexino data acquired between 2007 and 2011 [9], the present analysis of ten years of data from 2007 to 2017 achieves a significantly better precision on the muon flux, on the modulation parameters, and on the effective temperature coefficient. In addition, we expand the former analysis by measuring the atmospheric kaon-to-pion production ratio and observe a long-term modulation of the cosmic muon flux, hints for a correlation between this modulation and the solar activity, and the seasonal modulation of the cosmogenic neutron production rate.
Borexino is an organic liquid scintillator detector situated at the LNGS, covered by a limestone overburden of 3, 800 m w.e. [10]. It is designed for the spectroscopy of low energy solar neutrinos that are detected via elastic scattering off electrons. Based on the data acquired after the start of data taking in May 2007, Borexino accomplished measurements of the solar 7 Be [11][12][13][14], 8 B [15, 16], pep [13,17], and pp neutrino fluxes [13,18], a limit on the flux of solar neutrinos produced in the CNO cycle [13,17], and a spectroscopic measurement of antineutrinos produced in radioactive decays in the Earth, so-called geo-neutrinos [19][20][21]. Investigating the background to the neutrino analyses, Borexino further performed detailed studies of high energetic cosmic muons as well as of cosmogenic neutrons and radioactive isotopes from muon spallation on the detector materials [22]. The Borexino detector geometry allows to identify muons passing through a spherical volume with a cross section of 146 m 2 . The detection efficiency is virtually independent of the muon's incident angle, resulting in minimum systematics when measuring the muon flux and its variations. Detailed air temperature data is provided by weather forecasting centers [23] for the location of the laboratory and can be used to investigate the correlation between the flux of high energetic cosmic muons and the atmospheric temperature to determine the atmospheric temperature coefficient. In this article, we present an analysis of the cosmic muon flux as measured by Borexino based on ten years of data. In section 2, we briefly introduce the Borexino detector. In section 3, we report on the measured flux of cosmic muons and its seasonal modulation. In section 4, we introduce a model describing the expected relation between the flux of cosmic muons and the atmospheric temperature. In section 5, we present the modulation of the atmospheric temperature. In section 6, we analyze the correlation between the flux of cosmic muons and the atmospheric temperature. In section 7, we use the inferred effective temperature coefficient to measure the kaon-to-pion production ratio in the upper atmosphere. In section 8, we further analyze both the cosmic muon flux and the effective atmospheric temperature using a Lomb-Scargle periodogram. In section 9, we report the evidence for a long-term modulation with a period in agreement with the solar cycle. In section 10, we report on the seasonal modulation of the cosmogenic neutron production rate in Borexino. In section 11, we summarize and conclude our results.

The Borexino Detector
A schematic drawing of the Borexino detector [10] is shown in figure 1. In the present analysis, we consider muons passing through the Inner Detector (ID). It consists of a central organic scintillator target of 278 t composed of the solvent PC (1,2,4-trimethylbenzene) doped with the wavelength shifter PPO (2,5-diphenyloxazole) at a concentration of 1.5 g/l. The scintillator mixture is contained in a spherical and transparent nylon Inner Vessel (IV) with a diameter of 8.5 m and a thickness of 125 µm. To shield this central target from external γ-ray backgrounds and to absorb emanating radon, the IV is encompassed by two layers of buffer liquid in which the light quencher DMP (dimethylphthalate) is added to the scintillator solvent. A Stainless Steel Sphere (SSS) of 13.7 m diameter holding 2212 inward-facing 8" E.T.L. 9351 photomultiplier tubes (PMTs) that detect the scintillation light caused by particle interactions in the central region completes the ID. The ID is embedded in a steel dome of 18 m diameter and 16.9 m height that is filled with 2.1 kt of ultra-pure water. With the instrumentation of the outer surface of the SSS and the floor of the water tank with 208 PMTs, this Outer Detector (OD) provides an extremely efficient detection and tracking of cosmic muons via the Cherenkov light that is emitted during their passage through the water [25].

Seasonal Modulation of the Cosmic Muon Flux
The upper atmosphere is affected by temperature variations on the scale of seasons that alter the mean free path of the muon-producing mesons at the relevant production heights. These fluctuations are expected to be mirrored in a seasonal modulation of the underground muon flux since the high energies necessary for muons to pass through the rock overburden require that the parent mesons decay in flight without any former virtual interaction. The present analysis is based on ten years of Borexino data acquired between the 16 th of May 2007 and the 15 th of May 2017. Besides cosmic muons, the CERN Neutrino to Gran Sasso (CNGS) beam [26] that was operational between 2008 and 2012 introduces muon events in the Borexino detector [27]. These events have been carefully removed from the data sample via a comparison of the event time at Borexino and the beam extraction times as in [25]. To prevent statistical instabilities in the data sample, only data acquired on 3, 218 days for which a minimum detector livetime of eight hours was provided are considered. Besides a phase in 2010 and 2011 during which the liquid scintillator target underwent further purification, no prolonged downtime of the detector is present in the data set. Borexino features three different methods for muon identification of which two rely on the detection of the Cherenkov light generated in the OD. The Muon Trigger Flag (MTF) is set if a trigger is issued in the OD when the detected Cherenkov light surpasses a threshold value. The Muon Clustering Flag algorithm (MCF) searches for clusters in the OD PMT hit pattern. Further, muons can be identified via their pulse shape in the ID (IDF). The mean detection efficiencies have been measured to be 0.9925(2), 0.9928 (2), and 0.9890(1), respectively, and were found to remain stable. For details on the muon identification methods and the calculation of the efficiencies, we refer to [25]. In the present analysis, we define muons as events that are identified by the MCF. To account for small fluctuations of the muon identification efficiency, we estimate this efficiency for each bin and correct the measured muon rate. We discard events that do not trigger the ID to select tracks penetrating both the ID and OD volumes. The relevant detector cross section is, thus, 146 m 2 as given by the radius of the SSS, independent of the incident angle of the muon. The resulting effective exposure of the data set is ∼ 4.2 · 10 5 m 2 · d, in which ∼ 1.2 · 10 7 muons were detected. Most of the muons arriving at the Borexino detector are produced in decays of kaons and pions in the upper atmosphere. In the stratosphere, temperature modulations mainly occur on the scale of seasons, while short-term weather phenomena usually only affect the temperature of the troposphere, with the exception of stratospheric warmings that may lead to extreme temperature increases in the polar stratosphere during winter [28]. Since the higher temperature in summer lowers the average density of the atmosphere, the probability that the muon-producing mesons decay in flight before their first virtual interaction is increased due to their longer mean free paths. Only muons produced in these decays obtain enough energy to penetrate the rock coverage and reach the Borexino detector. As a consequence, the cosmic muon flux as measured by Borexino is expected to follow the modulation of the atmospheric temperature. In first order, the muon flux I µ (t) may be described by a simple sinusoidal behavior as with I 0 µ the mean muon flux, δI µ the modulation amplitude, T the period, and t 0 the phase. Short-or long-term effects are expected to perturb the ideal seasonal modulation. Moreover, temperature and flux maxima and minima will occur at different dates in successive years. The cosmic muon flux measured with Borexino is shown in figure 2 together with a fit according to eq. 3.1. For better visibility, the measured average muon flux per day is shown in weekwide bins while the presented results are inferred applying a fit to the muon flux in a day-wide binning. The lower panel shows the residuals (Data − Fit)/σ. We measure an average muon flux of phase of the strictly seasonal modulation is found to be t 0 = (181.7 ± 0.4) d, corresponding to a maximum on the 1 st of July. We consider this as our final estimate of the phase of the seasonal modulation. Especially in winter and spring, clear deviations from the sinusoidal assumption of the fit may be observed that can be attributed to a more turbulent environment of the upper atmosphere due to, e.g., stratospheric warmings [28]. Thus, the reduced χ 2 of the fit is χ 2 /NDF = 13702/362. To check the result, we selected a sample of muons as identified by the MTF and performed the same analysis steps. Similar results were obtained and we conclude that no systematic effects based on the muon definition are introduced. The flux of cosmic muons and the seasonal modulation have formerly been investigated by several experiments located at the LNGS, namely by MACRO [3], LVD [4], GERDA [6], and Borexino [9]. The results are summarized and compared to the present analysis in table 1. The LNGS consist of three experimental halls labelled A, B, and C. Note that Borexino reports a slightly higher rate value than all other experiments but GERDA. While Borexino is located in Hall C, the other measurements of the cosmic muon flux at the LNGS were performed in halls A and B. Since the halls are separated by approximately 80 m, the different locations of the experiments are expected to result in a slightly altered effective coverage due to the shape of the rock overburden, which might affect the muon rate at the respective experimental sites. Unlike Borexino, the acceptance of the other experiments carried out at the LNGS contains a dependence on the muon incident angle that must be carefully modelled.   flux may be affected by variations of the mean temperature or a long-term modulation of the cosmic muon flux. The seasonal modulation is found by all experiments and the phases agree well with that determined in the present work. Only GERDA reports a later maximum of the cosmic muon flux but their analysis is based on three years of data only.

Atmospheric Model and Effective Atmospheric Temperature
Since the mesons and consequently also the muons from their decays are produced at various heights in the atmosphere, it is extremely difficult to determine the point in the temperature distribution of the atmosphere where an individual muon was produced. To facilitate the investigation of the correlation between fluctuations of the atmospheric temperature and the cosmic muons flux observed underground, the atmosphere is modelled as an isothermal mesonproducing entity with the effective temperature T eff [2]. T eff is defined as the temperature of an isothermal atmosphere that produces the same meson intensities as the actual atmosphere.
Properly chosen weighting factors must be assigned to the respective depth levels accounting for the physics that determine the meson and muon production. A common parametrization is given by [8,30] where the approximation considers that the temperature is measured at discrete levels X n . The temperature coefficients α π (X) and α K (X) relate the atmospheric temperature to the muon flux considering pion and kaon contributions, respectively. These coefficients are translated into the weights W π n and W K n via numerical integration over the atmospheric levels ∆X n to allow the approximation. The weights are defined in appendix A. Figure 4 shows the ten year average temperature for different pressure levels at the closest point to the LNGS provided by the European Center for Medium-range Weather Forecasts (ECMWF) [23] and the assigned weights to the respective altitude levels. Higher layers of the atmosphere are assigned higher weights since muons possessing sufficient energy to penetrate the rock coverage of the LNGS are mainly produced at these altitudes. On the contrary, muons produced at lower altitudes are usually less energetic with the majority not having the threshold energy E thr to reach the detector. A so-called effective temperature coefficient may be defined as where W n (X) = W π n (X) + W K n (X). Thus, fluctuations of the cosmic muon flux may be related to fluctuations of the effective temperature via and α T quantifies the correlation between these two observables as discussed in section 6.

Seasonal Modulation of the Effective Atmospheric Temperature
To verify the correlation between the observed modulation of the cosmic muon flux and fluctuations of the atmospheric temperature, we analyze atmospheric temperature data provided by the ECMWF [23] for the time period corresponding to the muon flux measurement. This  Figure 4. The ten year average temperature [23] at the location of the LNGS is shown by the red line and the normalized weighting factors W π n + W K n by the black line, as functions of the pressure levels. The right vertical axis shows the altitude corresponding to the pressure level on the left vertical axis. data is generated by interpolating several atmospheric observables based on different types of observations (surface measurements, satellite data, or upper air sounding) and a global atmospheric model. For this analysis, we use the temperature for the location at 42.75 • N and 13.5 • E, which is the closest grid point to the LNGS available. The model provides atmospheric temperature data at 37 discrete pressure levels in the range from [0-1000] hPa four times per day at 00.00 h, 06.00 h, 12.00 h, and 18.00 h. Based on these data, we calculated T eff for each of the temperature sets based on eq. 4.1. The effective atmospheric temperature T eff of the respective day was computed as the average of the four values calculated during the day, their variance was used to estimate the uncertainty. Figure 5 shows the mean effective atmospheric temperature in a week-wide binning. Analogously to the cosmic muon flux, the modulation parameters were inferred by a fit to the data in a day-wide binning. A fit similar to eq. 3.1 returns an average effective atmospheric temperature of T 0 eff = (220.89 ± 0.01) K, a modulation amplitude of δT eff = (3.41 ± 0.01) K = (1.54 ± 0.06)%, a period of τ = (365.42 ± 0.04) d, and a phase of t 0 = (182.8 ± 0.2) d. While period and phase of the temperature modulation clearly show the leading seasonal behavior of the effective atmospheric temperature and agree well with the results of the muon flux discussed in section 3, the slightly higher modulation amplitude indicates that not all mesons relevant for the production of muons penetrating the LNGS rock coverage are affected by the density variations of the atmosphere. With a χ 2 /NDF = 106392/3214, a sinusoidal is only a very poor reproduction of the fine-grained temperature data. Similar to the flux of cosmic muons (see section 3), short term variations of the effective atmospheric temperature and, especially, additional secondary maxima in winter and spring are observed. These maxima may be ascribed to stratospheric warmings [28]. Sudden Stratospheric Warmings (SSW) sometimes even feature amplitudes comparable to the leading seasonal modulation [31], as visible e.g. in winter 2016/2017.

Correlation Between Muon Flux and Temperature
As expected, the modulation parameters inferred for the cosmic muon flux in section 3 and the effective atmospheric temperature in section 5 point towards a correlation of the two observables. Figure 6 shows the measured muon flux in Borexino and the effective atmospheric temperature scaled to percent deviations from their means I 0 µ and T 0 eff for ten years in a daywide binning. I 0 µ and T 0 eff were determined via sinusoidal fits to the respective data sets. Besides the consistency of the leading seasonal modulations of both observables, we find short-term variations of the temperature to be promptly mirrored in the undergound muon flux. Exemplarily, the short-term and non-seasonal temperature increase around January 2016 generates a secondary maximum of the muon flux. To quantify the correlation of the two observables, we plot ∆I µ /I 0 µ versus ∆T eff /T 0 eff for each day as shown in figure 7. Indeed, we find a positive correlation coefficient (R-value) of 0.55. Based on eq. 4.3, we determine the effective temperature coefficient by performing a linear regression using a numerical minimization method and accounting for error bars on both axes. We obtain α T = 0.90 ± 0.02 stat. in agreement with the former Borexino result of α T = 0.93 ± 0.04 stat. [9] but with the statistical uncertainties reduced by a factor ∼ 2.
To check for systematic uncertainties, we computed I 0 µ and T 0 eff by fitting a constant function to the data sets and performed the analysis for a two-year moving subset. We also selected muons based on an identification by the Borexino MTF and repeated the analysis. Since we find the result to be stable within the statistical expectations in all cases, we assume the systematic uncertainty to be small compared to the statistical uncertainty obtained from the fit. In table 2, the result of this analysis is compared to several further measurements performed at the LNGS. The results agree well within their uncertainties. The GERDA experiment [6] reported two values of α T using two different sets of temperature data. The theoretical expectation of α T at the location of the LNGS considering muon production from both kaons and pions is 0.92 ± 0.02 [9], also in agreement with our measurement.

Atmospheric Kaon-to-Pion Production Ratio
Since kaons and pions are affected differently by atmospheric temperature variations due to their distinct properties like mass, lifetime, or attenuation length, the strength of the correlation between the cosmic muon flux underground and the atmospheric temperature
depends on the production ratio of kaons and pions in the atmosphere. In the following, we infer an indirect measurement of the atmospheric kaon-to-pion production ratio based on the measurement of the effective temperature coefficient reported in section 6. For a properly weighted temperature distribution, the effective temperature coefficient α T is theoretically predicted to [2] with T being the temperature. The differential muon spectrum at the surface may be parametrized as [1] dI µ dE µ A × E −(γ+1) µ 1 1 + 1.1E µ cos θ/ π + 0.38 · r K/π 1 + 1.1E µ cos θ/ K (7.2) with r K/π the atmospheric kaon-to-pion ratio, θ the zenith angle, γ = 1.78 ± 0.05 [38] the muon spectral index, and π = (114±3) GeV and K = (851±14) GeV [1] the critical pion and kaon energies, respectively. The critical meson energy separates the decay from the interaction regime and mesons with an energy below this energy are more likely to decay while mesons with a higher energy most probably interact in the atmosphere before decaying. As shown in [2], equation 7.1 may be transformed into with the threshold energy E thr . The muon intensity underground may be approximated for the muon surface spectrum described by eq. 7.2 as [2, 3] With this approximation, the predicted α T may be calculated as and A K = 0.38×r K/π describing the kaon contribution to the cosmic muon flux [8]. E thr cos(θ) is the product of the threshold energy for a muon arriving from a zenith angle θ at the detector and the cosine of this angle. The mean value of this product allows to properly parametrize and compare the depths of various underground sites taking into account that the threshold energy is direction dependent due to the shape of the respective rock overburden. Figure 8 shows the weighted mean of α T for measurements performed at the LNGS together with measurements at other underground laboratories from Barrett [2], IceCube [7], MINOS [8], Double Chooz [32], and AMANDA [33]. The experimental results are plotted as a function of E thr cos θ , which constitutes an effective value of the measurement depths. The insert shows the LNGS based measurements from MACRO [3], GERDA [6], LVD [5], and the two Borexino measurements from 2012 [9] and from this work. For the LNGS, a value of E thr cos θ = (1.34 ± 0.18) TeV has been calculated based on a Monte Carlo simulation (see below). The red line shows the expected α T as a function of E thr cos θ considering muon production using the literature value of the atmospheric kaon-to-pion ratio of r K/π = 0.149 ± 0.06 [24], the dashed and dotted lines illustrate the extreme cases of pure pion or pure kaon production, respectively. The green line indicates the result of a fit to the measurements according to eq. 7.5 with r K/π as a free parameter. We obtain r K/π = 0.08 ± 0.02 stat. at a χ 2 /NDF = 2/6. Note, however, that systematic uncertainties like the exact value of E thr cos θ for the respective experimental sites are not fully determined and the result only constitutes an indication. Also, the measured values of α T depend on the assumed kaon-to-pion ratio since this quantity is included in the computation of T eff . We do not take into account this inter-value dependency here. The value of r K/π can also be inferred indirectly from the combination of a theoretical calculation of α T with the measurement from Borexino. In this case, no further experimental data has to be included. We performed a Monte Carlo simulation to calculate the expected value of α T at the location of the LNGS depending on r K/π . For muons with an energy E µ π , which is true for the muons arriving at the LNGS, the zenith angle distribution can be approximated by sec θ for θ < 70 • instead of the usual cos 2 θ [34]. We generated a toy Monte Carlo set of muons by randomly drawing a zenith angle from this distribution and an energy from the distribution given by eq. 7.2 for the respective θ. Moreover, a random azimuth angle φ was selected and the rock coverage D(θ, φ) for muons arriving from this direction was calculated based on an altitude profile of the Gran Sasso mountains obtained from the Google Maps Elevation API [35] and the density of the Gran Sasso rock of ρ = (2.3 ± 0.1) g/cm 3 [36]. We converted this into a direction dependent threshold energy E thr (θ, φ) for surface muons to arrive at the underground LNGS. For muons with E µ > E thr (θ, φ), we compute α T for r K/π values from 0 to 0.3 in steps of 0.01. For each value of r K/π , we calculated the mean value of the effective temperature coefficient α T for a sample of 10,000 muons surpassing E thr (θ, φ). Additionally, we performed the same calculations using the MUSUN [37] simulation code for the LNGS and obtained similar results. We compared the zenith angle distribution predicted by the simulation with the measured distribution and found them in good agreement. To estimate the systematic uncertainty of α T , we varied the input parameters of the simulation. We considered contributions from a 5% uncertainty of the altitude profile, the uncertainty of the measured rock density of the Gran Sasso rock of 0.05 g/cm 3 , the uncertainty of the measurement of the muon spectral index of 0.05 [38], the uncertainties of the critical meson energies ∆ π = 3 GeV and ∆ K = 14 GeV, and a 10% uncertainty of the drawn zenith angle. For the combined systematic uncertainty of α T , we found ∼ 0.015. However, the strength of several contributions of the above factors depends on r K/π . Especially the larger uncertainty of K compared to π leads to an increasing uncertainty of α T with rising r K/π . This simulation was used as well to calculate E thr cos θ = (1.34 ± 0.18) TeV for the location of the LNGS. Figure 9 shows the experimental and theoretical values of α T as functions of r K/π . Also the experimental value of α T has a weak dependence on r K/π since it enters into the calculation of the effective temperature T eff . To investigate this dependence, we calculated the daily T eff for the same range of r K/π values as above and redetermined α T for each set of T eff values via the correlation to the measured muon flux as in section 6. The resulting dependence is very weak and strongly overpowered by the statistical uncertainties of the measurements. Finally, to determine the kaon-to-pion production ratio, we estimate the intersection of the two allowed α T bands to obtain a value of r K/π = 0.11 +0.11 −0.07 . The allowed region in r K/π and α T has been determined by adding the χ 2 profiles of the Borexino measurement and the theoretical prediction. Former indirect measurements of the kaon-to-pion ratio were presented by the MINOS [8] and IceCube [7] experiments using a similar approach. Direct measurements have been carried out at accelerators, e.g. by STAR for Au+Au collisions at RHIC [39], by NA49 for Pb+Pb collisions at SPS [40], and by E735 forp + p collisions at Tevatron [41]. Results of  (   Figure 10. Comparison of several measurements of the kaon-to-pion production ratio. The STAR measurement was performed using Au+Au collisions at RHIC [39], the NA49 using Pb+Pb collisions at SPS [40], and the E735 usingp + p collisions at Tevatron [41]. The MINOS [8], IceCube [7], and Borexino measurements were performed indirectly via a measurement of the effective temperature coefficient. many older measurements using various reactions are summarized and referred to in [42]. The theoretical uncertainty of the kaon-to-pion ratio in current cosmic ray models is of the order of 40% [43]. Even though the indirect measurements do not directly compare with the accelerator experiments since the latter are performed with fixed beam energies, the central values are consistent as shown in figure 10. We place the Borexino data point in figure 10 at a center-of-mass energy √ s = (190 ± 28) GeV. Since cosmic muons obtain on average one tenth of the energy of the primary cosmic ray particle, we calculate the center-of-mass energy for Borexino assuming an average collision of a primary proton of 18 TeV (ten times the threshold energy for a muon arriving at the LNGS from straight above) on a fixed nucleon target. Due to the broad energy range of contributing muons, uncertainties on the center-of-mass energy need to be considered for the indirect measurements. Our result agrees with former indirect and direct measurements. Note that while the indirect measurements feature larger uncertainties than the accelerator experiments, they benefit from measuring the atmospheric kaon-to-pion ratio in situ. Due to the smaller muon statistics at greater depths, our measurement uncertainty is larger than for the MINOS and IceCube results. However, Borexino contributes the data point at the highest center-of-mass energy for indirect as well as fixed target measurements.

Lomb-Scargle Analysis of Muon Flux and Temperature
Besides the seasonal modulation of the cosmic muon flux underground, further physical processes might affect the cosmic muon flux and cause modulations of different periods. To investigate the presence of such non-seasonal modulations in the cosmic muon flux with Borexino, we perform a Lomb-Scargle analysis of the muon flux data. Lomb-Scargle (LS) periodograms [44,45] constitute a common method to identify sinusoidal modulations in a binned data set described by where N (t) is the expected event rate at time t given the data set is modulated with a period T , a relative amplitude A, and a phase φ. The LS power P for a given period T in a data set containing n data points may be calculated via is the difference between the data value in the j th bin and the weighted mean of the data set N 0 and σ 2 is the weighted variance. The weight w j = σ −2 j / σ −2 j of the j th bin is computed as the inverse square of the statistical uncertainty of the bin divided by the average inverse square of the uncertainties of the data set. The phase τ satisfies [46] tan 4π .

(8.3)
Since the quadratic sums of sine and cosine are used to determine the LS power, the latter is unaffected by the phase of a modulation as long as its period is short compared to the overall measurement time. Figure 11 shows a LS periodogram for the ten year cosmic muon data acquired with Borexino. To estimate the significance at which a peak in LS power exceeds statistical fluctuations, we use the known detector livetime distribution and mean muon rate to produce 10 4 white noise spectra distributed equally to the data. We define a modulation of period T to be significant if it surpasses a LS power P thr that is higher than 99.5% of the values found for white noise spectra at that period. This threshold is indicated by the red line in figure 11. Besides the leading peak of the seasonal modulation at 365 d, several secondary peaks are visible in figure 11. The second most significant peak is found at a long-term period of ∼ 3, 000 d. Further, additional peaks on the verge of significance are present at periods slightly larger than seasonal. These peaks might be caused by a peculiarity of the LS method that often characterizes harmonics of the leading modulation as significant. Finally, a significant peak is present at ∼ 180 d caused by the minor maxima in winter and spring as described in section 3. Due to these maxima, a bent sinusoidal is a better description of the muon data than the pure seasonal modulation and a peak in the LS periodogram arises.
Since only the peak of highest power in a LS periodogram allows a rigid statistical conclusion [47], we filter out the seasonal modulation from the data set by subtracting a sinusoidal function as described by eq. 3.1 with the fit parameters obtained in section 3. In this mode,  Figure 12. The left side shows the LS periodogram for the ten year effective atmospheric temperature data at the location of the LNGS [23]. On the right side, the LS periodogram of the effective atmospheric temperature data after the seasonal modulation was subtracted statistically is shown. The red lines indicate the sigificance level of 99.5%.
we repeat the LS analysis to probe specifically the significance of the striking long-term modulation. As shown on the right side of figure 11, the modulation with a period of ∼ 3, 000 d remains significant. This is also true for the 180 d peak indicating its physical origin. However, the peaks corresponding to harmonics of the seasonal modulation are now beneath the significance level. With the period of the long-term modulation being close to our overall measurement time, the phase of the modulation is expected to affect the LS power. To investigate this, we artificially generated data samples including a seasonal and a long-term modulation of 3, 000 d period equally binned as the muon flux data. For each sample, the phase of the long-term modulation was altered and we computed a LS periodogram. We found the location of the peak to vary between ∼ 2, 550 d and ∼ 3, 750 d, which indicates the absolute uncertainty of the period. However, the long-term modulation appears as a significant peak in the LS periodogram independent of the inserted phase. On the left panel of figure 12, we show the LS periodogram of the effective atmospheric temperature. Here, only the seasonal modulation is found to be significantly present in the data. No further period surpasses the threshold power.
To ensure that no long-term modulation might be inserted by the statistical subtraction of the seasonal modulation from the data, we repeated this procedure for the effective temperature data. As illustrated in the right panel of figure 12, no further significant peaks are introduced by this approach. Hence, we conclude that the significant long-term modulation in the cosmic muon flux at ∼ 3, 000 d is not present in and, hence, not related to the effective atmospheric temperature. We determine the phase and amplitude of the observed long-term modulation by fitting a function accounting for both the seasonal and the long-term modulation of the form to the day-widely binned data. We obtain a long-term modulation with a period T long = (3, 010 ± 299) d = (8.25 ± 0.82) yr, a phase t long 0 = (1, 993 ± 271) d, and an amplitude δI long µ = (14.7 ± 1.8) d −1 = (0.34 ± 0.04)%. T long is in good agreement with the period inferred from the LS periodogram and the phase of the long-term modulation indicates a maximum of the modulation in June 2012 for the investigated time frame. The χ 2 /NDF reduces from 3921/3214 when only a single modulation according to eq. 3.1 is fitted to the data to 3855/3211. Figure 13 shows our recorded muon data in year-wide bins after correcting for residual effects of the seasonal modulation caused by the uneven detector livetimes in different years. The red

Long-Term Modulation of the Cosmic Muon Flux and the Solar Activity
A long-term modulation of the cosmic muon flux has been observed before e.g. in the LVD data [48] and in [47], where data from MACRO, LVD, and Borexino were combined and a long-term modulation not noticed in the atmospheric temperature but in agreement with the solar activity was found.
To investigate the possibility of such a correlation, we perform a LS analysis of the daily sunspot data provided by the World Data Center SILSO, the Royal Observatory of Belgium in Brussels [49] for the time frame corresponding to the cosmic muon data acquired by Borexino as shown in figure 14. Since individual solar cycles are known to vary significantly in their period, it is not sensible to use a data set including earlier sunspot data. In figure 15, the most significant peak in the LS periodogram occurs at a period of ∼ 3, 000 d, in coincidence with the long-term modulation of the cosmic muon flux. The significance level of 99.5% was calculated following the procedure outlined in section 8. Further modulation periods are found to be significant in the sunspot data. This is expected since several authors observed minor modulations in the solar activity besides the solar cycle (see [50] and refs. therein). Also the increase in LS powers towards very high periods is expected since solar activity modulations larger than the solar cycle have been observed. An individual solar cycle can be described by a cubic power law and a Gaussian decline as [51]  The large uncertainties of the parameters of the long-term modulation of the cosmic muon flux for both fits as well as the large reduced χ 2 of the fit to the sunspot data indicate that the results need to be treated with care. A correlation between the solar activity and the flux of high energetic cosmic muons can neither be ruled out nor clearly be proven. However, we find indications encouraging further investigation of this phenomenon, especially considering the striking agreement between the modulation periods observed in the LS analysis. To eventually prove or negate a correlation, longer measurement times for the muon flux underground across several solar cycles may be necessary. A further indication of a positive correlation between solar activity and the flux of high energetic cosmic rays was found by the Tibet AS array [52] in observations of the size of the shade the Sun casts on ∼ 10 TeV cosmic rays. This shade was observed to be modulated depending on the solar activity with the shade shrinking by ∼ 50% at the maximum activity. Also Monte Carlo simulations performed by this collaboration for several solar surface models predicting the coronal magnetic field gave consistent results. However, the amplitude δI long µ = (0.34 ± 0.04)% we measure for the long-term modulation of the cosmic muon flux seems too high to leave a modulation of the solar shade as the sole explanation of a possible correlation.

Modulation of the Cosmogenic Neutron Production Rate
Cosmic muons may produce cosmogenic neutrons through various spallation processes on carbon in the Borexino organic scintillator target [53]. Neutrons are detected via the emission of a 2.2 MeV γ-ray following the capture on hydrogen or a total γ-ray energy of 4.9 MeV after the capture on 12 C. The capture time is τ = (259.7 ± 1.3 stat ± 2.0 syst ) µs [22]. As a secondary product of cosmic muons, the number of cosmogenic neutrons is expected to also undergo a seasonal modulation. Cosmogenic neutrons have been discussed as a possible background for the expected modulation in direct DM searches [54]. We investigate here the amplitude and phase of the cosmogenic neutron production rate. We select cosmogenic neutrons in a special 1.6 ms acquisition gate that is opened after each ID muon [25]. Events with a visible energy corresponding to at least 800 keV are selected. The efficiency of the neutron selection has been measured to be det = (91.7 ± 1.7 stat. ± 0.9 syst. )%  Figure 16. Rate of neutron-producing muons per day (left) and the cosmogenic neutron production rate (right) selecting only events with neutron multiplicity n ≤ 10. Ten years of data are projected to a single year with monthly binning. The red line depicts a sinusoidal fit to the data with the period fixed to twelve months.
after the stabilization of the electronics baselines ∼ 30 µs after the passage of a muon. Due to the stable muon detection efficiency and no significant changes of the detector, we expect this efficiency to be stable. We fail to observe the seasonal modulation of the cosmogenic neutron production rate, which we attribute to the occurrence of showering muons producing extremely high neutron multiplicities of up to ∼ 1, 000 [22] and following a non-Poissonian probability distribution. This hypothesis is sustained by the fact that an annual modulation can indeed be seen in the LS periodogram for the rate of neutron-producing muons. However, in order for the modulation to be significant in the neutron production rate, we need to remove those neutrons that are produced in high multiplicity showers from the sample. Figure 16 shows on the left the monthly binned data of neutron-producing muons projected to one year with a sinusoidal fit similar to eq. 3.1 and the period fixed to twelve months. Without efficiency correction, we obtain an average rate of neutron-producing muons R 0 µn = (36.8 ± 0.1) d −1 , an amplitude δR µn = (0.9 ± 0.2) d −1 = (2.3 ± 0.5)%, and a phase t 0 = (6.3 ± 0.4) months. The reduced χ 2 of the fit is χ 2 /NDF = 11/9. For the cosmogenic neutron production rate including neutrons produced in showers featuring up to 10 neutrons shown in figure 16 on the right, we observe a rate of R 0 n = (61.9 ± 1.5) d −1 , an amplitude δR n = (1.6±0.2) d −1 = (2.6±0.4)%, and a phase t 0 = (6.0±0.3) months at a χ 2 /NDF = 42/9. Beyond a neutron multiplicity of 10, we are unable to observe the seasonal modulation in the applying LS periodogram. The observed phase is in good agreement with the seasonal modulation of the entire cosmic muon flux. However, we find the amplitude of the modulation to be higher with a difference of about 2σ compared to the relative amplitude of (1.36 ± 0.04)% measured for the entire cosmic muon flux in section 3. Table 3 lists the phase and amplitude of the seasonal modulation measured for the number of neutron-producing muons as well as for the neutron production rate applying increasingly high neutron multiplicity cuts. Consistent results are found for all samples with neutron multiplicities n ≤ 10 beyond which the modulation is no longer significant in the LS periodogram. We find the phase of the cosmogenic neutron production rate to agree with the muon flux. Correspondingly, the maximum occurs approximately one month later than expected for dark matter particles. However, the amplitude of the modulation is higher compared to the muon flux independent of the actually applied multiplicity cut. To further probe this effect, we  Table 3. Parameters of the seasonal modulation observed for the number of neutron-producing muons and the neutron production rate applying increasingly high neutron multiplicity cuts. The second and third column show the phase and the amplitude of the seasonal modulation observed in the fit to the projected data, respectively. The last column shows the relative amplitude of the modulation following the residual approach.
computed the average cosmogenic neutron production and neutron-producing muon rate in three months in summer and three months in winter for each data set. We inferred the amplitude by assuming a sinusoidal modulation around the mean of the two values. The amplitudes observed following this residual approach are listed in the last column of table 3. We find consistent values to the ones obtained from the fit in each data sample confirming the increased modulation amplitude. The modulation of the cosmogenic neutron production rate has formerly been measured by the LVD experiment reporting an even higher modulation amplitude of δR n = (7.7 ± 0.8)% [55].
Since the neutron production depends on the muon energy, the larger amplitude of the modulation of the cosmogenic neutron production rate was interpreted as an indirect measurement of a seasonal modulation of the mean energy E µ of cosmic muons observed in the LNGS. The values measured by LVD implied a modulation amplitude of the mean muon energy of ∼ 10% or ∼ 28 GeV [55]. Following this interpretation and the calculation outlined in [56], our measurements of the amplitudes of the modulations of the cosmic muon flux and of the cosmogenic neutron production rate would indicate a modulation of the mean energy of cosmic muons of ∼ 4.5 GeV. However, we find a modulation of the mean energy for muons of several GeV to be unlikely. Based on the MUSIC/MUSUN [37] simulation codes, we altered the muon spectral index by 0.01 for the muon energy spectrum at the surface. For this input spectrum, we observed an increase of the mean muon energy of only 1 GeV but the muon flux increased by ∼ 14%, in strong contrast to our measurement. Additionally, we performed simulations of the muon production using the MCEq code [57] for various atmospheric models and inferred the muon surface spectra at Gran Sasso in winter and summer. We used these spectra as inputs for MUSIC/MUSUN to predict the corresponding underground spectra. Based on these simulations, we obtained a difference of the cosmic muon flux between the summer and winter spectra of about 1.4% in accordance to our measurement but observed only slight deviations of the mean muon energy of less than ∼ 0.1 GeV.
With the exclusion of a modulation of the mean muon energy, a more complex energy dependence of the cross section for neutron production than the conventionally assumed E α µ law [56] must be hypothesized to explain our observations. However, relatively large uncertainties of the atmosphere models and of the neutron production processes in the scintillator make it difficult to further investigate this percent-level effect.

Conclusions
We have presented a new precision measurement of the flux of cosmic muons in the LNGS under a rock coverage of 3, 800 m w.e. using ten years of Borexino data acquired between May 2007 and May 2017. We have measured a cosmic muon flux of (3.432 ± 0.001) · 10 −4 m −2 s −1 with minimum systematics due to the spherical geometry of the detector. The seasonal modulation of the flux of high energetic cosmic muons is confirmed and we have observed an amplitude of (1.36 ± 0.04)% and a phase of (181.7 ± 0.4) d corresponding to a maximum on the 1 st of July. We have used data from global atmospheric models to investigate the correlation between variations of the muon flux and variations of the atmospheric temperature and showed that the seasonal modulation is also present in the effective atmospheric temperature. The correlation coefficient between the two data sets is 0.55 indicating a positive correlation and we have measured the effective temperature coefficient α T = 0.90 ± 0.02 reducing the statistical uncertainties of our former measurement by a factor ∼ 2. The measurement is in good agreement with theoretical estimations and previous measurements carried out at the LNGS. We have performed a Monte Carlo simulation to calculate the theoretical expectation of α T at the location of the LNGS as a function of the atmospheric kaon-to-pion ratio r K/π . By calculating the intersection region of the expected value and our measurement of α T in dependence on r K/π , we have indirectly measured r K/π = 0.11 +0.11 −0.07 . This measurement is compatible with former indirect and accelerator measurements and constitutes a determination of r K/π in a new energy region for fixed target experiments. Based on a Lomb-Scargle periodogram, we have found evidence for a long-term modulation of the flux of high energetic cosmic muons with a period of ∼ 3, 000 d that is not present in the effective atmospheric temperature data. The amplitude of the long-term modulation is measured to be (0.34 ± 0.04)% and the maximum occurs around June 2012. We have found indications of an agreement between this modulation and the solar activity. However, the physical reason for a correlation between the high energetic part of the cosmic muon flux and the solar activity remains unclear. We have analyzed the production rate of cosmogenic neutrons in the Borexino detector as well as of the number of neutron-producing muons and found a seasonal modulation in phase with the cosmic muon flux but increased amplitudes of ∼ (2.6 ± 0.4)% and ∼ (2.3 ± 0.5)%, respectively. We have shown that a strong modulation of the mean muon energy underground as an explanation of this phenomenon is disfavored by performing simulations of the muon surface spectrum in summer and winter using the MCEq software code [57] and simulating the corresponding underground spectra using the MUSIC/MUSUN codes [37].