Investigating Earth shadowing effect with DAMA/LIBRA-phase1

In the present paper the results obtained in the investigation of possible diurnal effects for low-energy single-hit scintillation events of DAMA/LIBRA-phase1 (1.04 ton $\times$ yr exposure) have been analysed in terms of an effect expected in case of Dark Matter (DM) candidates inducing nuclear recoils and having high cross-section with ordinary matter, which implies low DM local density in order to fulfill the DAMA/LIBRA DM annual modulation results. This effect is due to the different Earth depths crossed by those DM candidates during the sidereal day.


Introduction
The present DAMA/LIBRA experiment [1,2,3,4,5,6,7,8,9,10,11,12,13], as the former DAMA/NaI [14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43] has the main aim to investigate the presence of DM particles in the galactic halo by exploiting the model-independent DM annual modulation signature (originally suggested in Refs. [44,45]). In particular, they have cumulatively reached a model independent evidence at 9.3σ C.L. for the presence of DM particles in the galactic halo by exploiting the DM annual modulation signature [4]. Recently the results obtained by investigating the presence of possible diurnal variation in the low-energy single-hit scintillation events collected by DAMA/LIBRA-phase1 (1.04 ton × yr exposure) have been released and analysed in terms of a DM second order modelindependent effect due to the Earth diurnal rotation around its axis [12]. In particular, the data were analysed using the sidereal time referred to Greenwich, often called GMST. No diurnal variation with sidereal time has been observed at the reached level of sensitivity, which was not yet adequate to point out the effect searched for there. In the present paper those experimental data are analysed in terms of an effect -named "Earth Shadow Effect" -which could be expected for DM candidate particles inducing nuclear recoils; this effect would be induced by the variationduring the day -of the Earth thickness crossed by the DM particle in order to reach the experimental set-up. It is worth noting that a similar effect can be pointed out only for candidates with high cross-section with ordinary matter, which implies low DM local density in order to fulfill the DAMA/LIBRA DM annual modulation results. Such DM candidates could get trapped in substantial quantities in the Earth's core; in this case they could annihilate and produce secondary particles (e.g. neutrinos) and/or they could carry thermal energy away from the core, giving potentiality to further investigate them.

The Earth Shadow Effect
During a sidereal day the Earth shields a terrestrial detector with a varying thickness, and this induces a variation of the flux of the DM candidates impinging the detector, mainly because of the modification of their veloc-ity distribution, f ( v). It is worth noting that this Earth Shadow Effect is very small and could be detectable only in case of candidates with high cross-section with ordinary matter (i.e. present in the galactic halo with small abundance).
The detector (and the hosting laboratory) velocity in the Galactic frame can be written as: where: (i) v LSR is the velocity of the Local Standard of Rest (LSR) because of the rotation of the Galaxy; (ii) v ⊙ is the Sun peculiar velocity with respect to LSR; (iii) v rev (t) is the velocity of the orbital motion of the Earth around the Sun and (iv) v rot (t) is the velocity of the rotation of the Earth around its axis. The two latter terms change as function of the sidereal time, t. Using the galactic coordinate frame (that is x axis towards the galactic center, y axis following the rotation of the Galaxy and the z axis towards the galactic North pole), one gets: v LSR = (0, v 0 , 0), where v 0 = (220 ± 50) km/s (uncertainty at 90% C.L.) [22,49] is the local velocity, and v ⊙ = (9, 12, 7) km/s [50]. The DM particles in the galactic halo have a velocity distribution g( w), which depends on the considered galactic halo model. Ref. [22] has shown many possible scenarios for the galactic halo; in the following we consider the isothermal halo model just because of its simplicity: with A normalization constant and v esc escape velocity, assumed in the following equal to 650 km/s, as often considered in literature; however, it is also affected by uncertainty. However, no sizeable differences are observed in the outcome when a different value of v esc = 550 km/s is considered, more closer to the 90% C.L. range of the RAVE Survey results [51]. In the laboratory frame the DM velocity distribution f ( v) is obtained from eq. 2 straightforward since w = v + v lab . To evaluate the expected daily variation of the DM particles velocity distribution due to the Earth Shadow Effect, it is necessary to estimate the time dependence of the θ angle, the "zenith distance" of v lab (i.e. the distance between v lab and the zenith, see Fig. 1). This can be determined by astrophysical considerations studying the Earth's rotation around its axis.
The simplest way to calculate θ(t) is in the equatorial coordinate system [12] where theê ecs 1 axis is directed towards the vernal equinox, and Figure 1: Schematic view of the DM particles impinging direction on a detector; the x ′ , y ′ , z ′ represent the laboratory frame coordinate system. Left: schematic representation of the correlation between the thickness, d, crossed by the considered DM candidates to reach a laboratory (hypotetically placed at the geographic North pole) and the DM impinging angle, θ in . Right: schematic representation of the experimental condition considered in the text: detector placed at the Gran Sasso National Laboratory (LNGS) with the z ′ axis in the vertical direction and the x ′ axis pointing to the vernal equinox. We define v s = v LSR + v ⊙ . In this section, when a numerical calculation is employed, we assume v 0 = 220 km/s; hence, v s = 232.28 km/s, and in the equatorial coordinate system: θ ecs = 42 • .18 (the "zenith distance") and ϕ ecs = −46 • .14 (azimuth angle fromê ecs 1 ). For simplicity in the following of this section in order to offer an estimate of the Earth Shadow Effect diurnal behaviour we consider v rev (t) equal to its annual mean value, i.e. zero. To introduce the Earth motion around its axis, firstly we consider the simplest case of a laboratory at North Pole and we define the horizontal coordinate system with z ′ axis directed as shown in Fig.1, lef t and x ′ axis directed towards a given longitude λ 0 . In this coordinate system (N label) the velocity of the Earth can be written as: v N lab = v s , θ N = θ ecs , and ϕ N = ϕ ecs − ω rot (t + t 0 ) = −(ϕ 0 + ω rot (t + t 0 )), where ϕ 0 = −ϕ ecs , ω rot = 2π/T d with T d = 1 sidereal day, t sidereal time referred to Greenwich, and t 0 = 24λ 0 /2π sidereal hours.

Deformation of the DM velocity distribution due to the Earth Shadow Effect
To study the experimental data in terms of possible Earth Shadow Effect, a Montecarlo code has been developed to simulate the propagation of the DM candidates elastically scattering off Earth nuclei in their travel in the Earth towards the underground experimental site. For such a purpose useful information has been gathered about the Earth composition and density. The Montecarlo code numerically estimates the velocity distribution -in the laboratory coordinate system -of the impinging DM particles after having crossed the Earth; such velocity distribution depends on the mass of the DM candidate, on its cross-section on nucleons, on the initial unperturbed velocity distribution, on the sidereal time, and on the latitude and longitude of the laboratory: f lab (v, t|m DM , σ n ). Then, this velocity distribution has been used to evaluate -in an assumed framework -the expected counting rate as a function of the sidereal time in order to be compared with the experimental data.
In this section details are given about the assumptions adopted in the simulation, as in particular: the Earth model, the mean free path and path reconstruction of such DM candidate, the adopted interaction model (nuclear form factor, scaling law, etc.) and the f lab (v, t|m DM , σ n ) estimation.

The Earth model
An Earth model has been assumed in order to estimate the signal variation due to the Earth Shadow Effect; in particular, the matter density and composition of the Earth have to be considered. The simulation adopts a simplified Earth model starting by the Preliminary Reference Earth Model (PREM) [52]. Just three main Earth layers with a constant density and homogeneous distribution in each one (these values are averaged over the PREM density distribution behaviour) are considered: i) the Inner Core; ii) the External Core; iii) the Mantle. The densities and mass percentage for each layer are given in Table 1; in this simplified model the rare Table 1: Density values, ρ L , and i-th nucleus mass percentage, δ i , adopted in the present calculations for the layers of the considered Earth model [52].
isotopes (mass percentage lower than 0.1%) have been neglected.

Interactions of the considered DM candidates and path reconstruction
We assume that the considered DM candidates loss their energy elastically scattering off nuclei with spin-independent coupling. The mean free path in the L-th Earth's layer is: where: i) n i = ρ L m i δ i is the number density of the i-th nuclei in the Earth's layer composed by N nuclear species; ii) σ DM,i is the cross-section on the i-th nucleus and m i is the mass of the i-th nucleus; iii) ρ L is the layer's density; iv) δ i is the i-th nucleus mass percentage in the layer. We considered a coherent scattering and a scaling law 2 : where: i) A i is the mass number of nucleus i; ii) σ n is the DM candidatenucleon cross-section; iii) µ i (µ n ) is the DM candidate -nucleus(nucleon) reduced mass. Thus, assuming m i ≃ A i m n , one can write: For simplicity, the deflection of the DM particles crossing the Earth is neglected (i.e. a linear trajectory is considered); in this assumption the Earth thickness crossed by those DM candidates, d, depends only on the impinging angle with respect to the detector, θ in , according to the relation: with R ⊕ Earth radius. This d value is the sum of the distances passed through by the DM candidate in each layer: d = L d L . Defining the maximum radii of the Earth layers R ic , R ec and R m = R ⊕ for Inner Core, External Core and Mantle respectively (see Table 1), the number of layers crossed by DM particles in the considered schema is: 0 for θ in ≥ 90 • ; 1 for 33 • .11 ≤ θ in < 90 • ; 3 for 11 • .05 ≤ θ in < 33 • .11; 5 for 0 • ≤ θ in < 11 • .05 (considering that arcsin (R ec /R ⊕ ) = 33 • .11 and arcsin (R ic /R ⊕ ) = 11 • .05). Thus, in this scenario, the DM particles move in each Earth's layer with a mean free path λ L , that mainly depends on the interaction cross section σ n (see eq. 8). The number of interactions in each layer, n hit , has been estimated as: -Case 1, high interaction cross sections, d L ≥ 50λ L : n hit is relatively high and follows a gaussian distribution with mean value and variance equal to d L /λ L ; -Case 2, small interaction cross sections, d L < 50λ L : a step-by-step approach has been adopted in the simulation. The path between two consecutive interactions, can be used to propagate the particle within the layer as long as In the considered scenario the DM candidate particles interact via SI elastic scattering on nuclei; thus, their energy-loss for each interaction is given by the induced nuclear recoil energy: where E in = m DM v 2 /2 and v are the DM energy and the velocity before the interaction, θ * is the angle of diffusion in the center of mass and r is a kinematic factor: The interaction cross-section is given by [18]: with F 2 i (E R ) nuclear form factor. In addition, for a given velocity v: with (also see eq. 10): Hence, assuming the cross section scaling law given in eq. 7, the interaction cross-section can be written as: Thus, after one interaction, such a DM particle loses an energy in the range 0 ≤ E R ≤ rE in with distribution given in eq. 15. In the following we assume for the nuclear form factor the Helm formula 3 [53,54]: where r 0 = r 2 i − 5s 2 , r i = 1.2A In the simulation, once n hit has been evaluated, the energy loss of the considered DM particle, E R , is estimated for each interaction and the output velocity is: Therefore, the net effect is a modification of the velocity distribution, f lab (v, t|m DM , σ n ); Fig. 3 shows some examples of the obtained velocity distribution for a detector located at LNGS.

The expected interaction rate
In the SI coupling scenario, considered here, the DM candidates scatter off the nuclei in the detector. Their expected interaction rate as a function of the nuclear recoil energy, E R , for a mono-atomic (i nucleus) detector is: where N T is the number density of the target nuclei in the detector and ρ = ξρ 0 is the DM particles density in the galactic halo (ξ is the relative abundance and ρ 0 the overall DM density in the galactic halo). The integral is calculated over all the possible DM particle velocities in the laboratory frame considering the distribution f lab (v) = f lab (v, t|m DM , σ n ).

The minimal velocity providing
The galactic escape velocity is included in the f (v) definition. Using eq. 15, the expected rate can be rewritten as: Generalizing to detectors with more than one kind of target nuclei (as e.g. in the case of the NaI(Tl) considered here), the expected experimental rate is: where: the G(E det , E ′ ) kernel takes into account the detector's energy resolution (generally through a gaussian convolution) and the K i (E ′ |E R ) kernel takes into account the energy transformation of the nuclear recoil energy in keV electron equivalent. For example the latter kernel can be written in the simplest case of a constant quenching factor, q i , as: . For a discussion about the quenching factors see Refs. [5,8].
The expected differential rate -as well as f lab (v, t|m DM , σ n ) -depends on the time through three different effects: i) the time dependence of the Earth's orbital motion velocity, v rev (t) (see eq. 1); ii) the time dependence of the Earth's rotation velocity around its axis, v rot (t) (see eq. 1); iii) the possible Earth Shadow Effect which depends on σ n . Details about the first two effects, responsible of the model-independent annual and diurnal modulation of the DM signal rate, respectively, are reported in Ref. [12]. Following the same approach the expected rate in an energy interval ∆E k can be written as 4 : where: i) S ′ 0,k (m DM , σ n ) is the time independent component of the expected signal; ii) S ′ m,k (m DM , σ n ) is the annual modulation amplitude, ω = 2π/T with T = 1 yr, t 0 ≃ June 2 nd ; iii) S ′ d,k (m DM , σ n ) is the diurnal modulation amplitude, ω rot = 2π/T d with T d = 1 sidereal day, t d for the case of a detector at the Gran Sasso longitude ranges from 13.94 h to 14.07 h depending on the v 0 value (see Ref. [12]); iv) S ′ d,sh,k (m DM , σ n , t) (whose average value over T d is null) takes into account the signal variation as a function of the sidereal time due to a possible Earth Shadow Effect.
The ratio R dy = S ′ d,k (m DM , σ n )/S ′ m,k (m DM , σ n ) is model independent and it is R dy ≃ 0.016 at LNGS latitude; thus, considering the DAMA/LIBRA-phase1 experimental result on the DM annual modulation, the expected ξσ n S ′ d,k (m DM , σ n ) is order of 10 −4 counts per sidereal day per kg per keV (cpd sid /kg/keV, hereafter) [12]. The reached experimental sensitivity of DAMA/LIBRA-phase1 [12] is not yet enough to observe such a diurnal modulation amplitude; in fact, in the (2-4) keV energy interval considered here, the experimental diurnal modulation amplitude from DAMA/LIBRA-phase1 data is (2.0 ± 2.1) × 10 −3 cpd sid /kg/keV [12]  that ξσ n = 1.1 × 10 −7 pb is compatible with the DAMA/LIBRA-phase1 DM annual modulation result, the obtained amplitude of ξσ n S ′ d,sh,k (t) is of order of 3 × 10 −2 cpd sid /kg/keV. Such value can be studied in DAMA/LIBRA-phase1 (see later).

Data analysis
The results, obtained by analysing in the framework of the Earth Shadow Effect the DAMA/LIBRA-phase1 (total exposure 1.04 ton×yr) data, essentially depend on the most sensitive (2-4) keV interval; thus, this is the energy region considered here.
In the present analysis, as in Refs. [5,32], three possibilities for the Na and I quenching factors have been considered: Q I ) the quenching factors of Na and I "constants" with respect to the recoil energy E R : q N a ≃ 0.3 and q I ≃ 0.09 as measured by DAMA with neutron source integrated over the 6.5 − 97 keV and the 22 − 330 keV recoil energy range, respectively [15]; Q II ) the quenching factors evaluated as in Ref. [55] varying as a function of E R ; Q III ) the quenching factors with the same behaviour of Ref. [55], but normalized in order to have their mean values consistent with Q I in the energy range considered there.
Another important effect is the channeling of low energy ions along axes and planes of the NaI(Tl) DAMA crystals. This effect can lead to an important deviation, in addition to the other uncertainties discussed above. In fact, the channeling effect in crystals implies that a fraction of nuclear recoils are channeled and experience much larger quenching factors than those derived from neutron calibration (see [5,30] for a discussion of these aspects). Since the channeling effect cannot be generally pointed out with neutron measurements as already discussed in details in Ref. [30], only modeling has been produced up to now. In particular, the modeling of the channeling effect described by DAMA in Ref. [30] is able to reproduce the recoil spectrum measured at neutron beam by some other groups (see Ref. [30] for details). For completeness, we mention an alternative channeling model, as that of Ref. [56], where larger probabilities of the planar channeling are expected. Moreover, we mention the analytic calculation claiming that the channeling effect holds for recoils coming from outside a crystal and not from recoils produced inside it, due to the blocking effect [57]. Nevertheless, although some amount of blocking effect could be present, the precise description of the crystal lattice with dopant and trace contaminants is quite difficult and analytical calculations require some simplifications which can affect the result. Because of the difficulties of experimental measurements and of theoretical estimate of this channeling effect, in the following it will be either included or not in order to give idea on the related uncertainty.
Thus, the data analysis has been repeated in some discrete cases which allow us to account for the uncertainties on the quenching factors and on the parameters used in the nuclear form factors. The first case (set A) is obtained considering the mean values of the parameters of the used nuclear form factors (see above and Ref. [25]) and of the quenching factors. The set B adopts the same procedure as in Refs. [20,21], by varying (i) the mean values of the measured 23 Na and 127 I quenching factors up to +2 times the errors; (ii) the nuclear radius, r i , and the nuclear surface thickness parameter, s, in the form factor from their central values down to -20%. In the last case (set C) the Iodine nucleus parameters are fixed at the values of case B, while for the Sodium nucleus one considers: (i) 23 Na quenching factor at the lowest value measured in literature; (ii) the nuclear radius, r i , and the nuclear surface thickness parameter, s, in the SI form factor from their central values up to +20%. Finally, three values of v 0 have been considered: i) the mean value: 220 km/s, and ii) two extreme cases: 170 and 270 km/s. Because of the large number of the needed simulations, the mass of the DM candidate and of the cross section on nucleon have been discretized as in the following: six m DM (5 GeV, 10 GeV, 30 GeV, 60 GeV, 100 GeV and 150 GeV) and eight σ n (10 pb, 5 pb, 1 pb, 0.5 pb, 0.1 pb, 0.05 pb, 0.01 pb and 0.005 pb). The expectations are compared with the experimental modelindependent diurnal residual rate of the single-hit scintillation events, measured by DAMA/LIBRA-phase1 in the (2-4) keV energy interval, as function of the sidereal time (see in Fig. 2 of Ref. [12]). Two examples of expected signals are reported in Fig. 5. Here, the used sidereal time bin is 1 hour (24 time bins in the sidereal day) and the experimental residuals are: S exp d (t i ) ± σ d,i . We compute the χ 2 quantity: where The S ′ d,sh,k (m DM , σ n , t i ) have been evaluated for each set of parameters described above.
The χ 2 of eq. 23 is function of only one parameter, ξ. Since the data do not show the presence of significant diurnal variation in the counting rate as already described in Ref. [12], only upper limits for ξ are allowed, once given m DM and σ n . We can define: The ∆χ 2 is a χ 2 with one degree of freedom and is used to determine the upper limit of ξ parameter at 2σ C.L.. Two examples to describe the followed procedure are reported in Fig.  6, where the excluded regions (above the dotted lines) in the ξ vs σ n plane for the cases of m DM = 10 and 60 GeV are shown as obtained on the basis of the Earth Shadow Effect in the given model framework.
The upper limits on ξ can be compared with the positive results from the DM annual modulation signature achieved by DAMA 5 . In particular, DAMA/LIBRA-phase1 reports an annual modulation amplitude in the (2-4) keV energy interval: S exp m = (0.0167 ± 0.0022) cpd/kg/keV, corresponding to 7.6σ C.L. [4]. Here for each set of parameters described above, one can evaluate (see e.g. eq. 22) the ξσ n allowed values as: This corresponds, once including the experimental uncertainties on S exp m , to a band in the ξ vs σ n plane (within the continuous solid line). In Fig. 6 such bands at 2σ C.L. are reported. One can see that for the scenario considered there and for m DM = 10 GeV the upper limits on ξ do not constrain the results of the DM annual modulation. On the contrary, for m DM = 60 GeV the upper limits on ξ do exclude the band with σ n > 0.05 pb and ξ > 10 −3 . The shaded bands in Fig. 6 corresponds to the allowed regions in the ξ vs σ n plane, for the given m DM , from the combined analyses of the DM annual modulation result and of the Earth Shadow Effect in the considered framework.
Finally, for each considered set of parameters the three-dimensional allowed region -calculated as described above -in the parameter's space: ξ, σ n , m DM , is depicted as a surface in Fig. 7, 8 and 9, for v 0 equal to 170, 220, 270 km/s, respectively. We note that the "thickness" of the allowed regions around the shown surfaces is ≤ ±30%; therefore, for simplicity it is not represented in these Figures.
Finally, we recall that other uncertainties not considered in the present paper are present. For example, including other possible halo models sizeable differences are expected in the results as shown e.g. in Refs. [5,22].

Conclusions
The Earth Shadow Effect has been investigated in a given framework considering the model independent results on possible diurnal variation of the low-energy rate of the single-hit scintillation events in the DAMA/LIBRA-phase1 data (exposure: 1.04 ton×yr) reported in Ref. [12]. For the considered DM candidates having high interaction crosssections and very small halo fraction the obtained results constrain at 2σ C.L., in the considered scenario, the ξ, σ n and m DM parameters (see Figs. 7, 8 and 9) when including the positive results from the DM annual modulation analysis of the DAMA/LIBRA-phase1 data [4]. For example, in the considered scenario for quenching factors Q I with channeling effect, B parameters set, v 0 = 220 km/s and m DM = 60 GeV, the obtained upper limits on ξ do exclude σ n > 0.05 pb and ξ > 10 −3 . When also including other uncertainties as other halo models etc. the results would be extended.