No role for muons in the DAMA annual modulation results

This paper gathers arguments and reasons why muons surviving the Gran Sasso mountain cannot mimic the Dark Matter annual modulation signature exploited by the DAMA/NaI and DAMA/LIBRA experiments. A number of these items have already been presented in individual papers. Further arguments have been addressed here in order to present a comprehensive collection and to enable a wider community to correctly approach this point.


Introduction
The DAMA/NaI and DAMA/LIBRA experiments at the Gran Sasso underground laboratory (LNGS) of the I.N.F.N. have been and are, respectively, investigating the presence of the Dark Matter (DM) particles in the galactic halo by exploiting the model independent DM annual modulation signature, originally suggested in the middle of '80s in Refs. [1,2]. In fact, as a consequence of the Earth annual revolution around the Sun, which is moving in the Galaxy traveling with respect to the Local Standard of Rest towards the star Vega near the constellation of Hercules, the Earth should be exposed to a higher flux of Dark Matter particles around ∼2 June (when the Earth orbital velocity is added to the one of the solar system with respect to the Galaxy) and by a smaller one around ∼2 December (when the two velocities are subtracted).
This DM annual modulation signature is very distinctive since the effect induced by DM particles must simultaneously satisfy all the following requirements: a e-mail: rita.bernabei@roma2.infn.it (I) the event rate must contain a component modulated according to a cosine function; (II) with period equal to one year; (III) with a phase roughly around June 2nd in case of usually adopted halo models (slight variations may occur in case of presence of non thermalized DM components in the halo); (IV) this modulation must be present only at low energy, where DM particles can induce signals; (V) it must be present only in those events where just a single detector, in a multi-detector set-up, actually "fires" (single-hit events), since the probability that DM particles experience multiple interactions is negligible; (VI) the modulation amplitude in the region of maximal sensitivity has to be 7 % in case of usually adopted halo distributions, but it may be significantly larger in some particular scenarios.
Thus, this signature has a different origin and peculiarities than effects correlated with seasons on the Earth.
To mimic such a signature spurious effects or side reactions should be able not only to account for the observed modulation amplitude but also to simultaneously satisfy all the requirements of the signature; thus, no other effect investigated so far in the field of rare processes offers a so stringent and unambiguous signature.
Let us now briefly describe the DAMA/LIBRA experiment [3], recalling its model independent annual modulation results [4,5]. The present DAMA/LIBRA set-up, installed at the Gran Sasso underground laboratory, is made of 25 highly radiopure NaI(Tl) crystal scintillators in a 5-rows 5-columns matrix. Each NaI(Tl) detector has 9.70 kg mass and a size of (10.2 × 10.2 × 25.4) cm 3 . The scintillation light (decay time 240 ns) of each crystal is collected (through two 10 cm long highly radiopure quartz light guides, which also act as optical windows being directly coupled to the bare crystal) by two low-background photomultipliers working in coincidence at single photoelectron threshold. The software energy threshold in the present data taking is 2 keV electron equivalent (hereafter keV) and the measured light response is 5.5-7.5 photoelectrons/keV depending on the detector. In order to reject afterglows, Cherenkov pulses in the light guides and Bi-Po events, a 500 µs veto occurs after each event [3]. The detectors are housed in a low radioactivity sealed copper box installed in the center of a low-radioactive Cu/Pb/Cd-foils/polyethylene/paraffin shield; moreover, about 1 m concrete (made from the Gran Sasso rock material) almost fully surrounds (mostly outside the barrack) this passive shield, acting as a further neutron moderator. In particular, the neutron shield reduces by a factor larger than one order of magnitude the thermal neutrons flux [4]. The copper box is maintained in HP Nitrogen atmosphere in slight overpressure with respect to the external environment; it is part of the 3-levels sealing system which prevents environmental air reaching the detectors. The experiment takes data up to the MeV scale despite the optimization is made for the lowest energy region. The linearity and the energy resolution of the detectors at low and high energy have been investigated using several sources as discussed in Ref. [3]; routine calibrations are carried out in the same conditions as the production runs, by using the glovebox installed in the upper part of the apparatus [3].
The DAMA/LIBRA data released so far correspond to six annual cycles for an exposure of 0.87 ton×yr [4,5]. Considering these data together with those previously collected by the former DAMA/NaI over 7 annual cycles (0.29 ton×yr) [6,7], the total exposure collected over 13 annual cycles is 1.17 ton×yr. Several analyses on the model-independent DM annual modulation signature have been performed (see Refs. [4,5] and references therein). A clear modulation is present in the (2-6) keV single-hit events and fulfills all the requirements of the DM annual modulation signature. In particular, no modulation is observed either above 6 keV or in the (2-6) keV multiple-hits events.
The results provide a model independent evidence of the presence of DM particles in the galactic halo at 8.9σ C.L. on the basis of the investigated DM signature. In particular, with the cumulative exposure the modulation amplitude of the single-hit events in the (2-6) keV energy interval, measured in NaI(Tl) target, is (0.0116 ± 0.0013) cpd/kg/keV; the measured phase is (146 ± 7) days (corresponding to May 26 ± 7 days) and the measured period is (0.999 ± 0.002) yr, values well in agreement with those expected for the DM particles.

The muon flux at LNGS
The muons surviving the coverage of the Gran Sasso laboratory either can have direct interactions in the experimental set-up or can produce in the surroundings and/or inside the set-up secondary particles, such as fast neutrons, γ 's, electrons, spallation nuclei, hypothetical exotics, etc., possibly depositing energy in the detectors. Such direct or indirect events are a potential background for low count rate experiments, as DAMA is. In this paper, the muon induced background in DAMA/LIBRA will be investigated and any possible role in the DAMA results will be quantitatively ruled out.
The surviving muon flux (Φ μ ) has been measured in the deep underground Gran Sasso Laboratory (3600 m w.e. depth) by various experiments with very large exposures [14][15][16][17][18][19]; its value is Φ μ 20 muons m −2 d −1 [14], that is about a factor 10 6 lower than the value measured at sea level. The measured average single muon energy at the Gran Sasso laboratory is 270 ± 3 (stat) ± 18 (syst) GeV [20]; this value agrees with the predicted values using different parametrizations [21]. A 2 % yearly variation of the muon flux was firstly measured years ago by MACRO; when fitting the data of the period January 1991-December 1994 all together, a phase around middle of July was obtained [14]. It is worth noting that the flux variation of the muons is attributed to the variation of the temperature in the outer atmosphere, and its phase changes each year depending on the weather condition. Recently, other measurements have been reported by LVD, quoting a lower amplitude (about 1.5 %) and a phase, when considering the data of the period January 2001-December 2008 all together, equal to (5 July ± 15 days) [15]. Finally, Borexino, has quoted a phase of (7 July ± 6 days), still considering the data taken in the period May 2007-May 2010 all together [16,17]. More recently, the Borexino collaboration presented a modified phase evaluation (29 June ± 6 days), 1 with a still lower modulation amplitude: about 1.3 % [18,19], by adding the data collected in a further year; the appreciable difference in the fitted values further demonstrates the large variability of the muon flux feature year by year.

Why muons cannot play any role
The measured muon variation at LNGS has no impact on the DAMA annual modulation results, recalled in Sect. 1. In the following sections we summarize the key items. It is worth noting the arguments reported in Ref. [22], where no evidence for a correlation between cosmic rays and DAMA result has been found and it is shown that the two phenomena differ in their power spectrum, phase, and amplitude.

Intrinsic inability of muons to mimic the DM annual modulation signature
Let us firstly recall [3][4][5][6][7][8][9][10][11][12][13] that a muon flux variation cannot mimic the Dark Matter annual modulation signature in DAMA/LIBRA (and even less in the smaller DAMA/NaI) set-up, not only because it may give rise just to quantitatively negligible effects (see later for details), but also because it is unable to mimic the DM signature. In fact, it would fail some requirements of the signature; namely e.g.: (i) it would induce variation in the whole energy spectrum.
(ii) it would induce variation in the multiple-hits events (events in which more than one detector "fires"), (iii) it would induce variation with a phase and amplitude distinctively different from the DAMA measured one (see later).
3.2 Inconsistency between the phase of muons and of the muon-induced effect and the DAMA phase The phase of muons surviving the Gran Sasso coverage, measured deep underground at LNGS, and the phase of the (2-6) keV single-hit events measured by DAMA are distinctively different (see Fig. 1). In particular, the values quoted by MACRO, LVD and Borexino experiments for the muon phase have to be regarded as mean values of the muon phases among the analyzed years and the associated errors are not simply due to statistical fluctuation of the data, but rather to the variations of the muon phase in the different years. The phase of the DAMA observed effect has instead a stable value in the different years [4,5] and is 5.7(5.9, 4.7)σ from the LVD (Borexino, first and recently modified evaluations, respectively) "mean" phase of muons (7.1σ from the MACRO one). This simple approach does not consider that the experimental errors in the muon flux are not completely Gaus- Fig. 1 The phase of the DAMA annual modulation signal [5] and the muon phases quoted by Borexino in two analyses (May 2007-May 2010 [16,17], and May 2007-May 2011 [18,19]), by LVD (January 2001-December 2008 [15]), and by MACRO (January 1991-December 1994 [14]). The muon phases quoted by those three experiments have to be regarded as mean values among the muon phases in all the considered years since the muon phase depends on the temperature of the outer atmosphere and, thus, it changes each year. The phase of the DAMA observed effect has instead a stable value in the different years [4,5]. The horizontal dashed line corresponds to 2nd June (date around which the phase of the DM annual modulation is expected). The middle of June is also marked as an example; in fact, the maximum temperature of the T eff at the LNGS location (see text) cannot be as early as the middle of June (and for several years), date which is still 3σ far away from the phase of the DAMA observed effect sian; however, it gives the right order of magnitude of the confidence level for the incompatibility between the DAMA phase and the phase of muons and of the muon-induced effects. Analyses carried out by different authors confirm these outcomes; for example, a disagreement in the correlation analysis between the LVD data on muon flux and the DAMA residuals with a confidence level greater than 99.9 % is reported in Ref. [22].
It is also worth noting that the expected phase for DM is significantly different than the expected phase of muon flux at Gran Sasso: in fact, while the first one is always about 152.5 day of the year, the second one is related to the variations of the atmospheric temperature above the site location, T eff . In particular, the atmosphere is generally considered as an isothermal body with an effective temperature T eff ; the behaviour of T eff at the LNGS location as function of time has been determined e.g. in Refs. [18,19]. As first order approximation T eff was fitted with a cosinusoidal behaviour and the phase turned out to be (24 June ± 0.4 days) [18,19]; this is later than e.g. the middle of June, date which is still 3σ far away from the DAMA measured phase (see Fig. 1). In addition, fitting the temperature values at L'Aquila in the Fig. 2 Time difference, (t side − t μ ), as function of τ . The t side parameter is the expected phase in case of any contribution to the modulation amplitude arising from the decay or the de-excitation (or whatever else) of hypothetical cosmogenic product (even exotic) produced by muons and having mean-life time τ . The muon phase of each considered year is indicated as t μ . Any exotic contribution related to muons would have a phase larger than the muon phase and, thus, even more distant from the expected DM phase and from the DAMA measured phase. T is the 1 year period; see text years 1990-2011 2 with a cosinusoidal function, a period of (365.1 ± 0.1) days and a phase of (25 July ± 0.6 days) are obtained.
Thus, in conclusion, the phase of the DAMA annual modulation signal [5] is significantly different than the phases of the surface temperature and of the T eff , on which the muon flux is dependent, and than the phases of the muon flux measured by MACRO, LVD and Borexino experiments.
The above argument also holds for every kind of cosmogenic product (even hypothetical exotics) due to muons.
In particular, when the decays or the de-excitations of any hypothetical cosmogenic product have mean-life time τ , the expected phase, t side , would be (much) larger than the muon phase (of each considered year) t μ , as shown in Fig. 2, and even more different from the one measured by the DAMA experiments and expected from the DM annual modulation signature ( June 2nd). In fact, the number of the cosmogenic products, N(t), satisfies the following equation: given by the sum of two contributions, the former due to the decay of the species and the latter due to their production, showing the typical pattern of muon flux with b/a 0.015, and period T = 2π/ω = 1 year; a is the mean production rate. Solving this differential equation, one has: 2 The values of the temperature at L'Aquila are taken from [23].
where A is an integration constant, and t side = t μ + arctg(ωτ ) ω (see Fig. 2). In condition of secular equilibrium (obtained for time scale greater than τ ), the first term vanishes and the third term shows an annual modulation pattern with phase t side . The relative modulation amplitude of the effect is: Two extreme cases can be considered: if τ T /2π , one gets t side t μ + τ ; else if τ T /2π , one gets t side t μ + T /4 ( t μ + 90 days) and the relative modulation amplitude of the effect is 1.5 %.
In conclusion, the phase of muons and of whatever (even hypothetical) muon-induced effect is inconsistent with the phase of the DAMA annual modulation effect.

No role for the muons interacting in the detectors directly
In addition to the previous arguments, the direct interaction of muons crossing the DAMA set-ups cannot give rise to any appreciable variation of the measured rate. In fact, the exposed NaI(Tl) surface of DAMA/LIBRA is about 0.13 m 2 (and smaller in the former DAMA/NaI); thus the muon flux in the 250 kg DAMA/LIBRA set-up is about 2.5 muons/day. In addition, the impinging muons give mainly multiple-hits events and over the whole energy spectrum.
The order of magnitude of such a contribution has been estimated by a Monte Carlo calculation [24, 25] which takes into account the measured muon intensity distribution as function of the incident direction and energy, the Gran Sasso overburden map [26], and the geometry of the set-up [3]. The direct interaction of muons is simulated according to the impinging direction of muons on the detector and to the energy loss per specific path in the detector. Three topologies of the detectors' locations in the 5 × 5 DAMA/LIBRA detectors' matrix have been considered. The results are shown in Fig. 3. One can easily infer that the contributions from muons interacting in the detectors directlyfor single-hit events in the (2-6) keV energy region-to the DAMA total counting rate (that is around 1 cpd/kg/keV) and to the observed annual modulation amplitude (that is around 10 −2 cpd/kg/keV) are negligible (many orders of magnitude lower). In addition, as mentioned above, this contribution would also fail some requirements of the DM annual modulation signature such as III, IV and V.

No role for fast neutrons produced by muon interaction
The surviving muons and the muon-induced cascades or showers can be sources of neutrons in the underground laboratory. Such neutrons produced by cosmic rays are substantially harder (extending up to several hundreds MeV energies [21]) than those from environmental radioactivity; their typical flux is about 10 −9 neutrons/cm 2 /s [27], that is three orders of magnitude smaller than the neutron flux produced by radioactivity.
In particular, the fast neutron rate produced by muons interaction is given by: where M eff is the effective mass where muon interactions can give rise to events detected in the DAMA set-up and Y is the integral neutron yield, which is normally quoted in neutrons per muon per g/cm 2 of the crossed target material.
The integral neutron yield critically depends on the chemical composition and on the density of the medium through which the muons interact. The dependence on atomic weight is well described by a power law [21]: Y = 4.54 × 10 −5 A 0.81 neutrons per muon per g/cm 2 ; alternatively, it can also be expressed as [21]: Y = 1.27 × 10 −4 (Z 2 /A) 0.92 neutrons per muon per g/cm 2 . Thus, the integral yield of neutrons produced by muons deep underground at LNGS is Y (1-7) × 10 −4 neutrons per muon per g/cm 2 [21,28] for relatively light nuclei and Y 4.5 × 10 −3 neutrons per muon per g/cm 2 [21] for lead.
Consequently, the modulation amplitude of the singlehit events in the lowest energy region induced in DAMA/ LIBRA by a muon flux modulation can be estimated according to: where g is a geometrical factor, is the detection efficiency for neutrons, f E is the acceptance of the considered energy window (E ≥ 2 keV), f single is the single-hit efficiency and 1.5 % is the muon modulation amplitude. Since: assuming the very cautious values: and taking for M eff the total mass of the heavy shield, 15 ton, one obtains: 0.5 % of the observed single-hit events modulation amplitude. Even assuming in the calculation the yield Y for lead, the upper limit of S (μ) m (<1.5 × 10 −4 cpd/kg/keV) remains still lower than 1 % of the observed single-hit events modulation amplitude.
We stress that-in addition-the latter value has been overestimated by orders of magnitude both because of the extremely cautious values assumed in the calculation and because of the omission of the effect of the neutron shield of the set-up.
In conclusion, any appreciable contribution from fast neutrons produced by the muon interactions can be quantitatively excluded. In addition, it also would fail some of the requirements of the DM annual modulation signature such as III, IV and V.
For completeness, in the next two sub-sections we will address the case of environmental neutrons of whatever origin (and, thus, also including those induced by muons). In the first sub-section the outcomes in Refs. [4][5][6][7]29] will be recalled without entering in details, while in the second one the case of the neutron capture on Iodine will be analysed in depth.
3.5 . . . and no role for environmental neutrons Environmental neutrons cannot give any significant contribution to the annual modulation measured by the DAMA experiments [4][5][6][7]29]. In fact, the thermal neutron flux surviving the multicomponent DAMA/LIBRA shield has been determined by studying the possible presence of 24 Na from neutron activation of 23 Na in NaI(Tl). In particular, 24 Na presence has been investigated by looking for triple coincidences induced by a β in one detector and by the two γ 's in two adjacent ones. An upper limit on the thermal neutron flux surviving the multicomponent DAMA/LIBRA shield has been derived as: <1.2 × 10 −7 n cm −2 s −1 (90 % C.L.) [4], that is at least one order of magnitude lower than the value of the environmental neutrons measured at LNGS. The corresponding capture rate is: <0.022 captures/day/kg. Fig. 4 The decay schema of 128 I [31] Even assuming cautiously a 10 % modulation (of whatever origin 3 ) of the thermal neutrons flux, and with the same phase and period as for the DM case, the corresponding modulation amplitude in the lowest energy region would be [4,6]: <0.01 % of the DAMA observed modulation amplitude. Similar outcomes have also been achieved for the case of fast neutrons; the fast neutrons have been measured in the DAMA/LIBRA set-up by means of the inelastic reaction 23 Na(n, n ) 23 Na * (2076 keV) which produces two γ 's in coincidence (1636 keV and 440 keV). An upper limitlimited by the sensitivity of the method-has been found: <2.2 × 10 −7 n cm −2 s −1 (90 % C.L.) [4], well compatible with the value measured at LNGS; a reduction at least an order of magnitude is expected due to the neutron shield of the set-up. Even when cautiously assuming a 10 % modulation (of whatever origin) of the fast neutrons flux, and with the same phase and period as for the DM case, the corresponding modulation amplitude is <0.5 % of the DAMA observed modulation amplitude [4,6].
Moreover, in no case the neutrons can mimic the DM annual modulation signature since some of the peculiar requirements of the signature would fail, such as III, IV and V.

No role for 128 I decay
Reference [30] has claimed that environmental neutrons (mainly thermal and/or epithermal), occasionally produced by high energy muon interactions, once captured by Iodine might contribute to the modulation observed by DAMA 3 For the sake of correctness, it is worth noting that a variation of the neutron flux in the underground Gran Sasso laboratory has never been suitably proved. In particular, besides few speculations, there is just an unpublished 2003 short internal report of the ICARUS collaboration, TM03-01, that seemingly reports a 5 % environmental neutron variation in hall C by exploiting the pulse shape discrimination (PSD) in commercial BC501A liquid scintillator. However, the stability of the data taking and of the applied PSD procedures over the whole data taking period and also the nature of the discriminated events are not fully demonstrated. Anyhow, even assuming the existence of a similar neutron variation, it cannot quantitatively contribute to the DAMA observed modulation amplitude [4,6,7] as well as satisfy all the peculiarities of the DM annual modulation signature. through the decay of activated 128 I (that produces-among others-low energy X-rays/Auger electrons). Such an hypothesis is already excluded by several arguments given above (as e.g. those in Sects. 3.1 and 3.5), moreover it has already been rejected in Refs. [11,12]; anyhow, in the following we will focus just on its main argument avoiding to comment on several other wrong statements present in Ref. [30].
The 128 I decay schema is reported in Fig. 4. When 128 I decays via the EC channel (6.9 %), it produces low energy X-rays and Auger electrons, totally contained inside the NaI(Tl) detectors; thus, the detectors would measure the total energy release of all the X-rays and Auger electrons, that is the atomic binding energy either of the K-shell (32 keV) or of the L-shells (4.3 to 5 keV) of the 128 Te. The probability that so low-energy gamma's and electrons would escape a detector is very small. In Ref. [30] it is claimed that such low-energy gamma's and electrons from the L shells may contribute to the DAMA observed annual modulation signal; but: 1. considering the branching ratios of the EC processes in the 128 I decay, the K-shell contribution (around 30 keV) must be about 8 times larger than that of L-shell; while no modulation has been observed by DAMA above 6 keV (see [4,5] and references therein) and, in particular, around 30 keV; 2. the 128 I also decays by β − with much larger branching ratio (93.1 %) than EC (6.9 %) and with a β − end-point energy at 2 MeV. Again, no modulation has instead been observed in DAMA experiments at energies above 6 keV [4,5]; 3. the L-shell contribution would be a gaussian centered around 4.5 keV; this shape is excluded by the behaviour of the measured modulation amplitude, S m , as a function of energy (see Fig. 6 (bottom)). The efficiencies to detect an event per one 128 I decay are: 2 × 10 −3 , 6 × 10 −3 , and 2 × 10 −3 in (2-4) keV, (4-6) keV and (6-8) keV respectively, as calculated by the Monte Carlo code. Thus, the contribution of 128 I in the (2-4) keV would be similar to the one in the (6-8) keV, while the data exclude that.
Moreover, the data collected by DAMA/LIBRA allow the determination of the possible presence of 128 I in the detectors. In fact, neutrons would generate 128 I homogeneously distributed in the NaI(Tl) detectors; therefore studying the characteristic radiation of the 128 I decay and comparing it with the experimental data, one can obtain the possible 128 I concentration. The most sensitive way to perform such a measurement is to study the possible presence of the 32 keV peak (K-shell contribution) in the region around 30 keV. This was already done and published by DAMA-for other purposes-in Ref. [25] before Ref. [30] appeared. As it can be observed in Fig. 5, there is no evidence of such a peak in the DAMA/LIBRA data; hence an upper limit on the area of a peak around 32 keV can be derived to be: 0.074 cpd/kg (90 % C.L.) [25]. Considering the branching ratio for Kshell EC, the efficiency to detect events in the energy interval around 30 keV for one 128 I decay has been evaluated by the Monte Carlo code to be 5.8 %. Hence, one can obtain a limit on possible activity of 128 I (a 128 ): This upper limit allows us to derive the maximum counting rate which may be expected from 128 I in the keV region; it is reported in Fig. 6 (top) together with the cumulative low-energy distribution of the single-hit scintillation events measured by DAMA/LIBRA [3]. It can be noted that any hypothetical contribution from 128 I would be negligible. Moreover, even assuming the hypothetical case of a 10 % environmental neutron flux modulation, and with the same phase and period as the DM signal, the contribution to the DAMA   6 Top-Data points: cumulative low-energy distribution of the single-hit scintillation events measured by DAMA/LIBRA [3] above the 2 keV software energy threshold of the experiment. Histogram: maximum expected counting rate from 128 I decays corresponding to the measured upper limit on 128 I activity in the NaI(Tl) detectors: <15 µBq/kg (90 % C.L.); see the data in Ref. [25] and the text. Bottom-Data points: the DAMA measured modulation amplitude as a function of the energy. Histogram: maximum expected modulation amplitude multiplied by a factor 30 as a function of the energy from 128 I decays corresponding to the measured upper limit on 128 I activity given above when assuming an hypothetical 10 % neutron flux modulation, and with the same phase and period as a DM signal (Color figure online) measured (2-6) keV single-hit modulation amplitude would be <3×10 −4 cpd/kg/keV at low energy, as reported in Fig. 6 (bottom), that is <2 % of the DAMA observed modulation amplitudes. In conclusion, any single argument given in this section excludes a role played by 128 I.

No role for hypothetical phosphorescence induced by muons
In Ref. [32] it is argued that delayed phosphorescent pulses induced by the muon interaction in the NaI(Tl) crystals might contribute to the (2-6) keV single-hit events. Many wrong statements are put forward in that reference. We have already critically addressed Ref. [32] in our Ref. [13]. We will just focus on the main argument.
Thus, the number of muons is too low to allow a similar effect to contribute to the DAMA observed amplitude; in fact, to give rise to the DAMA measured modulation amplitude each μ should give rise to about 270 (∼10 counts/day/(0.0375 μ/day), see above) single-hit correlated events in the (2-6) keV energy range in a relatively short period. But: (i) such a hypothesis would imply dramatic consequences for every NaI(Tl) detector at sea level (where the μ flux is 10 6 times larger than deep underground at LNGS), precluding its use in nuclear and particle physics; (ii) phosphorescence pulses (as afterglows) are single and spare photoelectrons with very short time decay (∼10 ns); they appear as "isolated" uncorrelated spikes. On the other side, scintillation events are the sum of correlated photoelectrons following the typical time distribution with mean time equal to the scintillation decay time (∼240 ns). Pulses with short time structure are already identified and rejected in the noise rejection procedure described in detail in [3] (the information on the pulse profile is recorded). Thus, in addition, phosphorescence pulses are not present in the DAMA annual modulation data; (iii) because of the poissonian fluctuation on the number of muons, the standard deviation of the (2-6) keV singlehit modulation amplitude due to a similar effect would be 13 times (see Appendix) larger than that measured by DAMA, and therefore no statistically-significant effect, produced by any correlated events, could be singled out. Even just this argument (that will be further illustrated in the following) is enough to discard the hypothesis of Ref. [32] (similar considerations are also reported in Ref. [33]); (iv) the muon phase is inconsistent with the phase measured by DAMA (see Sect. 3.2).
Thus, the argument regarding a possible contribution from delayed phosphorescent pulses induced by muons can be safely rejected.

Absence of long term modulation
In Ref. [33] it is argued that high-energy muons measured by LVD might show a long term modulation with a period of about 6 years, suggesting that a similar long term modulation might also be present in the DAMA data. We avoid to comment here on the other arguments reported in Ref. [33] that are actually already addressed elsewhere in the present paper, and already discard such a suggestion. However, for completeness such a long term modulation has also been looked for in the DAMA data.
For each annual cycle of DAMA/NaI and DAMA/ LIBRA, we calculated the annual baseline counting ratesthat is the averages on all the detectors (j index) of flat j Fig. 7 Power spectrum of the annual baseline counting rates for (2-4) keV (dashed blue curve) and (2-6) keV (dotted red curve) single-hit events. Also shown for comparison is the power spectrum obtained by considering the 13 annual cycles of DAMA/NaI and DAMA/LIBRA for the single-hit residuals in (2-6) keV (solid black line) [5]. The calculation has been performed according to Ref. [34,35], including also the treatment of the experimental errors and of the time binning. As can be seen, a principal mode is present at a frequency of 2.735 × 10 −3 d −1 , that corresponds to a period of 1 year. The 99.7 % confidence lines for excluding the white noise hypothesis are also shown (see text). No statistically-significant peak is present at lower frequencies and, in particular, at frequency corresponding to a 6-year period. This implies that no evidence for a long term modulation in the counting rate is present (Color figure online) (i.e. the single-hit rate of the j -th detector averaged over the annual cycle, see e.g. [4])-for the (2-4) keV and (2-6) keV energy intervals, respectively. Their power spectra (dashed (blue online) and dotted (red online) curves, respectively) in the frequency range 0.0002-0.0018 d −1 (corresponding to a period range 13.7-1.5 year) are reported in Fig. 7; the power spectrum (solid black line) above 0.0022 d −1 , obtained when considering the (2-6) keV single-hit residuals of Fig. 1 of Ref. [5], is reported for comparison. To evaluate the statistical significance of these power spectra we have performed a Monte Carlo simulation imposing constant null expectations for residuals; from the simulated power spectra the probability that an apparent periodic modulation may appear as a result of pure white noise has been evaluated. The 99.7 % confidence lines for excluding the white noise hypothesis are shown in Fig. 7. A principal mode is present in the power spectrum of the experimental data for a frequency equal to 2.735 × 10 −3 d −1 (black solid curve), corresponding to a period of 1 year, while no statistically-significant peak is present at lower frequencies and, in particular, at frequency corresponding to a 6-year period.
A further investigation of any hypothetical 6-year period has been performed, taking into account that the LVD muon data have a 1-year period modulation amplitude equal to 1.5 % [15], and-according to the claim of Ref.
[33]-a 6-year period modulation amplitude equal to 1 % (actually Fig. 8 A detail of Fig. 7, with superimposed power spectrum (solid black curve), expected in the hypothesis of a contribution from muons. The shaded region is the 1σ (68 % C.L.) band. The peak at 6 year period would be in this hypothesis well evident and above the threshold of detectability at 3σ C.L. On the contrary, the power spectrum of the experimental data (dashed blue and dotted red curves) is completely outside the 1σ band. For simplicity, the calculations are shown just for the cumulative (2-6) keV energy interval. This further shows that no evidence for a long term modulation in the counting rate is present, as it should already be expected on the basis of the many other arguments discussed in this paper. See text (Color figure online) 1 %, as reported in Ref. [33]). Thus, in the case that muons might contribute to the DAMA effect, a 6-year period modulation would be present in the DAMA data with amplitude 0.008 cpd/kg/keV, that is a 1 %/1.5 % fraction of the 1-year period modulation amplitude measured by DAMA. A simulation of 10 6 experiments has been performed and the power spectrum of each simulation has been derived. The average of all the simulated power spectra is reported in Fig. 8 as solid black curve; the shaded region is the 1σ (68 % C.L.) band. The hypothetical peak at 6-year period would be under these assumptions well evident and above the threshold of detectability at 3σ C.L. On the contrary, the power spectrum of the experimental data is completely outside the 1σ band. For simplicity, these calculations have been reported just for the cumulative (2-6) keV energy interval.
This further shows that no evidence for a long term modulation in the counting rate is present, as-on the other hand-it should already be expected on the basis of the many other arguments (and just one suffices) discussed in this paper, further demonstrating that there is no role for muons.

No role for muons from statistical considerations
In addition to the previous arguments, let us finally demonstrate that any hypothetical effect-even with high-multiplicity production-due to muons crossing the NaI(Tl) de- Fig. 9 Expected standard deviation (solid line), σ (A), of the (2-6) keV single-hit annual modulation amplitude, A, as a function of the effective area, A eff , in the hypothetical case that muons might produce in A eff "something" (gamma's, beta's, neutrons, phosphorescence pulses, exotics, . . .) able to contribute to the DAMA (2-6) keV single-hit modulation amplitude (see text and the Appendix). The value experimentally observed by DAMA is shown as (red) dashed line; the σ (A) curve approaches this value just for A eff 50 m 2 while even in the most cautious case the A eff value is much smaller. The four interesting cases described in the text are shown; as it can be seen, in all four cases the fluctuation is much larger than that observed by DAMA (Color figure online) tectors and/or the surroundings of the set-up cannot give any appreciable contribution to the observed (2-6) keV singlehit event rate, just owing to statistical considerations. In fact, because of the poissonian fluctuation on the number of muons, the standard deviation, σ (A), of any hypothetically induced (2-6) keV single-hit modulation amplitude, A, would be much larger than measured by DAMA, thus, giving rise to no statistically-significant effect. To explain this argument, Fig. 9 (see Appendix) reports the expected σ (A) (solid line) as a function of an effective area, A eff , around the set-up. This effective area is defined as the area crossed by the muons that would give rise to the bulk contribution (say 90 % of the total) of the DAMA measured annual modulation amplitude through their "products" (gamma's, beta's, neutrons, phosphorescence pulses, whatever exotics, . . .). Therefore, muons outside this area can safely be neglected for the purpose of this description. As a matter of fact, the size of the effective area depends on the considered "product" of muons, on its nature and on its interactions with the materials around the DAMA/LIBRA set-up. The smaller is A eff , the smaller would be the number of muons that might hypothetically contribute (either directly or indirectly) to the DAMA annual modulation amplitude, and the larger would be the fluctuation σ (A).
The value experimentally observed by DAMA for σ (A), 0.0013 cpd/kg/keV [5], is shown as (red online) dashed line in Fig. 9; the σ (A) curve approaches this value just for A eff 50 m 2 , while-as we will demonstrate-even in the most cautious case the A eff value is much smaller. In fact, to have an idea on the possible sizes of A eff , four interesting cases are shown in Fig. 9: (i) the muons interacting directly in the NaI(Tl) DAMA/LIBRA detectors (hypothetically producing either very short range particles or phosphorescence pulses as discussed in Sect. 3.7, . . .), corresponding to A eff equal to the DAMA/LIBRA exposed surface: 0.13 m 2 ; (ii) the effective area equal to the one calculated (by a Monte Carlo simulation) just considering the 1/r 2 dependence for the flux of the "products", without including any shield effect, corresponding to A eff 1.65 m 2 ; (iii) the effective area equal to the area of the heavy passive shield, A eff 1.75 m 2 ; (iv) the effective area equal to the area of the heavy passive shield plus the neutron moderator, A eff 3 m 2 . In all the four cases the fluctuation, driven by the small number of muons, is much larger than that observed by DAMA. Since it is extremely safe to consider that any hypothetical mechanism would have a corresponding A eff within the previous considered cases (that is A eff 50 m 2 ), we can conclude that all (standard and exotic) mechanisms, because of the low number of the involved muons, provide too high fluctuations of the data, not observed in DAMA. Even just this argument is enough to discard any kind of hypothesis about muons.

Conclusion
In this paper we have compiled some of the main scientific and quantitative arguments which demonstrate that there is no room for any hypothetical contribution from muons to the (2-6) keV single-hit annual modulation amplitude measured by DAMA experiments. Some comments about incorrect arguments reported in recent papers [30,32,33,36,37] have also been addressed. In conclusion, the hypothesis of a role for muons in the DAMA observed (2-6) keV single-hit annual modulation can be safely rejected for many scientific reasons.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.
T is the period (1 year); (iv) t 0 is the phase of the annual modulation (t 0 = 152.5 day of the year). Moreover, it turns out useful to define the following quantities: the total exposure W = Σ i w i , the mean value of the cosine function β = Σ i w i c i W , and the variance of the cosine function For a data taking involving the whole year, one can expect β 0 and (α − β 2 ) 0.5.
The variable X allows us to obtain information about the parameters of the annual modulation. In fact, from the experimental data one can obtain the value of X and its error σ X = √ Var(X); the variance can be calculated directly from the experimental data hypothesizing that N i is distributed as poissonian variable (that is Var(N i ) = N i ): These two quantities are directly connected with the parameters of the annual modulation; see later.
In the following two cases will be considered: (i) the case of an annual modulation induced by DM particles; (ii) the case of an annual modulation due to "generic products" of the muons crossing the set-up and/or the surroundings.

Case of annual modulation induced by DM particles
In such a case the expectation of the stochastic variable N i is given by: Here b is the background rate, S 0 is the unmodulated component of the DM signal and S m is the modulation amplitude. The expected value of the X variable is connected with the modulation amplitude by that is, since the first term is null because of the definition of β: Therefore, an estimate of X through the experimental data gives the S m evaluation. The sensitivity can be determined by considering that Var(N i ) (b + S 0 )w i , since the number of events have a poissonian distribution and the S m c i terms can be safely neglected. Hence, the variance of X can be written from eq. (A.2) as: The sensitivity reachable on the modulation amplitude is given by: . Let us now consider the case of the DAMA/NaI (target mass 100 kg) and DAMA/LIBRA (target mass 250 kg) results: M · T = 425428 kg day, E = 4 keV, η 0.7, W (α − β 2 ) = 5.96 × 10 5 kg day keV. The counting rate of the set-ups (b + S 0 ) is around 1 cpd/kg/keV. Thus, σ (S m ) 0.0013 cpd/kg/keV, well comparable with the one obtained by DAMA. Therefore, such a sensitivity has allowed the measurement of a modulation amplitude S m of the order of 10 −2 cpd/kg/keV. We remind that the singlehit modulation amplitude in the (2-6) keV energy region measured by DAMA/NaI and DAMA/LIBRA is (0.0114 ± 0.0013) cpd/kg/keV (here period and phase fixed in the fit).
In Fig. 9 the experimental value of σ (S m ) and the expected σ (S m ) are reported as a (red online) dashed line; actually, they are overlapped.
In conclusion, the experimental value of σ (S m ) is well compatible with the expectation, giving further support to the evidence for an annual modulation effect induced by DM particles.
Case of annual modulation directly or indirectly due to muons crossing the set-up and/or the surroundings We consider here the case of hypothetical effects correlated with muons crossing the NaI(Tl) detectors and/or the surroundings of the set-up.
To easily describe the model, an effective area, A eff , can be defined as the area around the DAMA/LIBRA setup where whatever (even hypothetical) product of muons (gamma's, beta's, neutrons, phosphorescence pulses, any exotics, . . .) would give the bulk contribution (we would say 90 % of the total) to the measured (2-6) keV single-hit modulation; muons outside this area are for simplicity not considered in the following, because of their slight contribution.
The rate of muons inside the effective area is r μ = Φ μ × A eff muons/day and its relative annual modulation is Δ 0.015 [15][16][17]. Let us assume for a while as true the scenario: (i) where each muon, crossing the NaI(Tl) detectors and the A eff area, might produce during the incoming period (minutes, hours, days) ε single-hit events in the considered low energy bin in all the DAMA detectors; (ii) the period and the phase of the muon flux is equal to those of the DM signal. By the way, ε is distributed as a poissonian variable with expectation and variance equal toε.
The expected value of N i would be in this case (for simplicity η ∼ 1 is assumed): here b is still the background and R μ = r μ M E . Now considering the results obtained by DAMA and the assumed scenario, an estimation of theε parameter can be obtained from the measured modulation amplitude, A 10 −2 cpd/kg/keV: . For example, if A eff = 0.13 m 2 , as for the case of direct interaction of muons producing e.g. hypothetical delayed phosphorescence pulses [32] (hypothesis already discarded above), each muon should have to give rise toε 255 single-hit (2-6) keV events in all the DAMA/LIBRA set-up, value in agreement with that given in Sect. 3.7. Moreover, R μ ·ε A Δ 0.68 cpd/kg/keV; considering that the total counting rate of the DAMA detectors is around ∼1 cpd/kg/keV, one can derive b 0.32 cpd/kg/keV.
Following the same procedure as above, one gets: Finally, the sensitivity reachable on the modulation amplitude, A, is given by: where σ 2 B = b/[W (α − β 2 )] takes into account the fluctuation of the number of the background events, σ 2 μ = (R με 2 )/[W (α −β 2 )] the fluctuation of the number of muons crossing the A eff area and, finally, σ 2 ε = (R με )/[W (α − β 2 )] takes into account the fluctuation of the number of the hypothetically produced low energy events after a muon.
Considering for the parameters the values given above, one has: Therefore, the major contribution to the error for A eff < 100 m 2 is from the fluctuation of the number of muons crossing the NaI(Tl) detectors and the A eff area; in addition, this contribution becomes larger for smaller A eff , since the smaller would be the number of muons hypothetically able to directly or indirectly contribute to the DAMA (2-6) keV single-hit modulation amplitude, the larger would be the fluctuation σ (A). The standard deviation σ (A), calculated as in Eq. (A.13), is reported in Fig. 9 as function of the effective area A eff .
For example, if A eff is the exposed NaI(Tl) surface of DAMA/LIBRA, the standard deviation σ (A) is in this case equal to 0.017 cpd/kg/keV, that is 13 times larger than the one measured by DAMA.
Other enlightening cases, reported in the text and in Fig. 9, show effective areas at level of few m 2 ; in all these cases the fluctuations are much larger than the observed value by DAMA. Since it is extremely safe to assume that any hypothetical mechanism would have a corresponding A eff within the previous considered cases, we can conclude that all (standard and exotic) mechanisms, because of the low number of the involved muons, provide too high fluctuations of the data, not observed in DAMA.
In conclusion any hypothetical and quantitatively appreciable effect correlated with muons crossing either the NaI(Tl) detectors or the surroundings of the set-up can be further excluded, even already just owing to statistical considerations about the poissonian fluctuation on the number of muons.
In other words, because of the low number of the involved muons, all (standard and exotic) muon-induced mechanisms provide fluctuations of the data larger than those observed by DAMA, and therefore-in addition to all the other arguments given in this paper-cannot give rise to any evidence of modulation in DAMA.