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 ref. [1]. 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 crossed by a larger flux of Dark Matter particles around ∼ 2 June (when the Earth orbital velocity is summed 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: 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 2 nd 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 in a well-defined low energy range, 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 [2], recalling its model independent annual modulation results [3,4]. 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 [2]. 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 [3]. The copper box is maintained in HP Nitrogen atmosphere in slightly overpressure with respect to the external environment; it is part of the 3-levels sealing system which excludes the detectors from environmental air. 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. [2]; routine calibrations are carried out in the same conditions as the production runs, by using the glove-box installed in the upper part of the apparatus [2].
The DAMA/LIBRA data released so far correspond to six annual cycles for an exposure of 0.87 ton×yr [3,4]. Considering these data together with those previously collected by the former DAMA/NaI over 7 annual cycles (0.29 ton×yr) [5,6], 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. [3,4] 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.
Careful investigations on absence of any significant systematics or side reaction able to account for the measured modulation amplitude and to simultaneously satisfy all the requirements of the signature have been quantitatively carried out (see e.g. ref. [3,4,2,5,6,7,8,9,10,11], refs therein); none has been found or suggested by anyone over more than a decade. In particular, the case of muons has been deeply investigated. This paper will further demonstrate that neither muons nor muon-induced events can significantly contribute to the DAMA observed annual modulation signal. Some already-published arguments [3,4,2,5,6,7,8,9,10,11] are also recalled.

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, as fast neutrons, γ's, electrons, spallation nuclei, hypothetical exotics, etc., possibly releasing energy in the detectors. Such direct or indirect events are a potential background for lowcounting 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 [12,13,14,15]; its value is Φ µ ≃ 20 muons m −2 d −1 [12], 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 [16]; this value agrees with the predicted values using different parametrizations [17]. A ≃ 2% yearly variation of the muon flux was firstly measured years ago by MACRO; when fitting the data of several years all together, a phase around middle of July was obtained [12]. 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 several years all together, equal to (5 July ± 15 days) [13]. Finally, Borexino, almost for the same period, has quoted a phase of (7 July ± 6 days), still considering the data taken in various years all together [14]; more recently, the Borexino collaboration presented a modified phase evaluation on the same set of data, (29 June ± 6 days), with a still lower modulation amplitude, about 1.3% [15].

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. [18], 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 possibly in amplitude.

Intrinsic inability of muons to mimic the DM annual modulation signature
Let us firstly remind [3,4,5,6,7,8,9,10,2,11] 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 well different from the DAMA measured one (see later).

Inconsistency between the phase of muons and of the muoninduced 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 well 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 [3,4] and is 5.7 (5.9, 4.7) σ far from the LVD (Borexino, first and recently modified evaluations, respectively) "mean" phases of muons (7.1 σ far from the MACRO one). In addition, considering the seasonal weather, the maximum temperature of the outer atmosphere (on which the µ flux variation is dependent) at the LNGS location cannot be as low as e.g. the middle of June, date which is still 3 σ far away from the DAMA phase (see Fig. 1). 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 muons 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: MACRO [12] surviving µ the mountain Figure 1: The phase of the DAMA annual modulation signal [4] and the muon phases quoted by Borexino in two analyses [14,15], by LVD [13], and by MACRO [12]. 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 [3,4]. 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 outer atmosphere at the LNGS location cannot be as low as the middle of June (and for several years), date which is still 3 σ far away from the phase of the DAMA observed effect.
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: 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 , 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 much distant from the expected DM phase and from the DAMA measured phase. T is the 1 year period; see text.
term vanishes and the third term shows an annual modulation pattern with phase t side . The relative modulation amplitude of the effect is: b/a √ 1+(ωτ ) 2 . 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) muoninduced effect is inconsistent with the phase of the DAMA annual modulation effect.

No role for the muons direct interaction in the detectors
In addition to the previous arguments, the direct interaction of muons crossing the DAMA set-ups cannot give rise to any quantitatively 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 Montecarlo calculation which takes into account the measured muon intensity distribution, the Gran Sasso overburden map [19], and the geometry of the set-up; three typologies 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 direct muons contributions -for 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 the muons 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 [17]) than those from environmental radioactivity; their typical flux is about 10 −9 neutrons/cm 2 /s [20], 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 ef f 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 [17]: Y = 4.54 × 10 −5 A 0.81 neutrons per muon per g/cm 2 ; alternatively, it can also be expressed as [17]: 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,17] for relatively light nuclei and Y ≃ 4.5×10 −3 neutrons per muon per g/cm 2 [17] for lead.
Consequently, the modulation amplitude of the single-hit 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 muons modulation amplitude. Since: assuming the very cautious values: and taking for M ef f the total mass of the heavy shield, 15 ton, one obtains: m ≪ 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 muons 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. [3,4,5,6,22] will be recalled without entering in details, while in the second one the case of the neutron capture on Iodine will be deeply analysed.

... and no role for environmental neutrons
Environmental neutrons cannot give any significant contribution to the annual modulation measured by the DAMA experiments [3,4,5,6,22]. In fact, the thermal neutron flux surviving the multicomponent DAMA/LIBRA shield has been determined by studying the possible presence of the 24 Na from neutron activation of 23 Na in NaI(Tl). In particular, the 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.) [3], 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. Even assuming cautiously a 10% modulation (of whatever origin 1 ) 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 [3,5]: < 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 limit -limited by the sensitivity of the method -has been found: < 2.2×10 −7 n cm −2 s −1 (90% C.L.) [3], 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 [3,5].
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
Ref. [23] has claimed that environmental neutrons (mainly thermal and/or epithermal), occasionally produced by high energy muons interactions, once captured by Iodine might contribute to the modulation observed by DAMA 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 Sect. 3.1 and 3.5), moreover it has already been rejected in ref. [9,10]; anyhow, in the following we will focus just on its main argument avoiding to comment several other wrong statements present in ref. [23].
The 128 I decay schema is reported in Fig. 4. When 128 I decays in 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. [23] 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 [3,4] 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 [3,4]; 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 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. [23] 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 K-shell EC, the efficiency to detect events in the energy interval around 30 keV for one 128 I decay has been evaluated by the Montecarlo code to be 5.8%. Hence, one can obtain a limit on possible activity of 128 I (a 128 ): a 128 < 15 µBq/kg (90% C.L.). [25] and the text. Bottom -Data points: the DAMA measured modulation amplitude as a function of the energy. Histogram (color online): 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.
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 [2]. 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 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. [26] it is argued that delayed phosphorescent pulses induced by the muons 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. [26] in our ref. [11]. We will just focus on the main argument.
Since the µ flux in DAMA/LIBRA is about 2.5 µ/day (see Sect. 2), the total µ modulation amplitude in DAMA/LIBRA is about: 0.015 × 2.5 µ/day ≃ 0.0375 µ/day (1.5% muon modulation has been adopted, see Sect. 2). The single-hit modulation amplitude measured in DAMA/LIBRA in the (2-6) keV energy interval is instead [3,4]: 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 (hours, days?). But: i) such a hypothesis would imply a 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 details in [2] (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 single-hit 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. [26] (similar considerations are also reported in ref. [27]); 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. [27] 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 the other arguments reported in ref. [27] that are actually already addressed elsewhere in the present paper, which 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 rates -that is the averages on all the detectors (j index) of f lat j (i.e. the single-hit rate of the j-th detector averaged over the annual cycle) -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 (black solid line) above 0.0022 d −1 , obtained when considering the (2-6) keV single-hit residuals of Fig. 1 of ref. [4], is reported for comparison. While a principal mode is present for a frequency equal to 2.735 × 10 −3 d −1 , corresponding to a period of ≃ 1 year, no statisticallysignificant peak is present at lower frequencies and, in particular, at frequency corresponding to 6 years period. This implies that no evidence for a long term modulation in the counting is present. Thus, further again no role for muons, as single-hit events. It is also shown for comparison 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 (black solid line) [4]. The calculation has been performed according to ref. [28], including also the treatment of the experimental errors and of the time binning. As it 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. No statistically-significant peak is present at lower frequencies and, in particular, at frequency corresponding to 6 years period. This implies that no evidence for a long term modulation in the counting rate is present.
it should already be expected on the basis of the many other arguments discussed in this paper (and just one suffices).

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) detectors and/or the surroundings of the set-up cannot give any appreciable contribution to the observed (2-6) keV single-hit 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. 8 (see Appendix) reports the expected σ(A) (solid line) as a function of an effective area, A ef f , 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 ef f , 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 [4], is shown as (red) dashed line in Fig. 8; the σ(A) curve approaches this value just for A ef f ≫ 50 m 2 , while -as we will demonstrate -even in the most cautious case the A ef f value is much smaller. In fact, to have an idea on the possible sizes of A ef f , four interesting cases are shown in Fig. 8: i) the muon direct interaction in the NaI(Tl) DAMA/LIBRA detectors (hypothetically producing either very short range particles or phosphorescence pulses as discussed in section 3.7, ...), corresponding to A ef f equal to the DAMA/LIBRA exposed surface: 0.13 m 2 ; ii) the effective area equal to the one calculated (by a Montecarlo simulation) just considering the 1/r 2 dependence for the flux of the "products", without including any shield effect, corresponding to A ef f ≃ 1.65 m 2 ; iii) the effective area equal to the area of the heavy passive shield, A ef f ≃ 1.75 m 2 ; iv) the effective area equal to the area of the heavy passive shield plus the neutron moderator, A ef f ≃ 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 ef f within the previous considered cases (that is A ef f ≪ 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 gathered 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 [23,26,27,29,30] 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.

Appendix
We give here some elements to properly evaluate the expected standard deviation of the annual modulation amplitude in an experiment studying the DM annual modulation signature, as DAMA does. Two results are compared with the experimental values and the curve of Fig. 8 is justified.
Let us simplify the result of an annual modulation experiment as a set of N i counts per each day i in the interesting energy region. The index i identify the day in the year and, therefore, ranges from 0 to 365. The experimental information about the data taking can be gathered in another quantity that is the set of the daily exposure w i = M ∆t i ∆Eη, where M is the exposed mass, ∆t i is the live time during the day, t i is the time of the bin (center-bin value), ∆E is the energy bin and η is the efficiency. A new variable can be defined: where: i) c i = cosω(t i − t 0 ); ii) ω = 2π T ; iii) 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 . 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 = V ar(X); the variance can be calculated directly from the experimental data hypothesizing that N i is distributed as poissonian variable (that is V ar(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 V ar(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. (2) as: The sensitivity reachable on the modulation amplitude is given by: with . 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 single-hit 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. 8 the experimental value of σ(S m ) and the expected σ(S m ) are reported as a (red) 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 ef f , can be defined as the area around the DAMA/LIBRA set-up 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 ef f muons/day and its relative annual modulation is ∆ ≃ 0.015 [13,14]. Let us assume for a while as true the scenario: i) where each muon, crossing the NaI(Tl) detectors and the A ef f 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 ef f = 0.13 m 2 , as for the case of direct interaction of muons producing e.g. hypothetical delayed phosphorescence pulses [26] (hypothesis already discarded above), each muon should have to give rise tō ε ≃ 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.
To obtain the latter equation one can profit of the fact that N i can be written as the sum of two components, the number of background events and the number of the events hypothetically induced by muons: and the V ar(N µ i ) can be calculated, considering that N µ i has two sources of fluctuation: i) one associated to the number of crossing muons; ii) one associated to the hypothetical production of low-energy events. Without loosing generality, the two terms can be written as: All this justifies eq. (10). One can see that for large number ofε the contribution of the fluctuation of the number of muons is largely dominant. 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 ef f 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.
Therefore, the major contribution to the error for A ef f < 100 m 2 is from the fluctuation of the number of muons crossing the NaI(Tl) detectors and the A ef f area; in addition, this contribution becomes larger for smaller A ef f , 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. (13), is reported in Fig. 8 as function of the effective area A ef f .
For example, if A ef f 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. 8, 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 ef f 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 setup 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.