Determining the Mass of Dark Matter Particles with Direct Detection Experiments

In this article I review two data analysis methods for determining the mass (and eventually the spin-independent cross section on nucleons) of Weakly Interacting Massive Particles with positive signals from direct Dark Matter detection experiments: a maximum likelihood analysis with only one experiment and a model-independent method requiring at least two experiments. Uncertainties and caveats of these methods will also be discussed.


Introduction
There is strong evidence that more than 80% of all matter in the Universe is dark (i.e., interacts at most very weakly with electromagnetic radiation and ordinary matter). The dominant component of this cosmological Dark Matter should be due to some yet to be discovered, non-baryonic particles. Weakly Interacting Massive Particles (WIMPs) χ arising in several extensions of the Standard Model of electroweak interactions are one of the leading candidates for Dark Matter. WIMPs are stable particles with masses roughly between 10 GeV and a few TeV and interact with ordinary matter only weakly (for reviews of WIMPs and some other possible candidates for Dark Matter, see Refs. [1,2,3]).
Currently, the most promising method to detect different WIMP candidates is the direct detection of the recoil energy deposited in a low-background laboratory detector by elastic scattering of ambient WIMPs on the target nuclei [4,5,6] 1 . The basic expression for the differential event rate for elastic WIMP-nucleus scattering is given by [1]: Here R is the direct detection event rate, i.e., the number of events per unit time and unit mass of detector material, Q is the energy deposited in the detector, F (Q) is the elastic nuclear form factor, f 1 (v) is the one-dimensional velocity distribution function of the WIMPs impinging on the detector, v is the absolute value of the WIMP velocity in the laboratory frame. The constant coefficient A is defined as where ρ 0 is the WIMP density near the Earth and σ 0 is the total cross section ignoring the form factor suppression. The reduced mass m r,N is defined by where m χ is the WIMP mass and m N that of the target nucleus. Finally, v min is the minimal incoming velocity of incident WIMPs that can deposit the energy Q in the detector: and v max is related to the escape velocity from our Galaxy at the position of the Solar system, v esc . It was found that, by using a time-averaged recoil spectrum dR/dQ, and assuming that no directional information exists, the normalized one-dimensional velocity distribution function of incident WIMPs, f 1 (v), can be solved from Eq.(1) directly as [7] f 1 (v) = N −2Q · d dQ where the normalization constant N is given by Note that, firstly, because f 1 (v) in Eq. (6) is the normalized velocity distribution, the normalization constant N here is independent of the constant coefficient A defined in Eq. (2). Secondly, the integral in Eq.(7) goes over the entire physically allowed range of recoil energies: starting at Q = 0, and the upper limit of the integral has been written as ∞. However, it is usually assumed that the WIMP flux on the Earth is negligible at velocities exceeding the escape velocity v esc . This leads thus to a kinematic maximum of the recoil energy The velocity distribution function of halo WIMPs reconstructed by Eq. (6) is independent of the local WIMP density ρ 0 as well as of the WIMP-nucleus cross section σ 0 . However, not only the overall normalization constant N given in Eq. (7), but also the shape of the velocity distribution, through the transformation Q = v 2 /α 2 in Eq. (6), depends on the WIMP mass m χ involved in the coefficient α defined in Eq. (5). In fact, any (assumed) value of m χ will lead to a well-defined, normalized distribution function f 1 (v) when one uses Eq. (6). Hence, m χ can be extracted from a single recoil spectrum only if one makes some assumptions about the velocity distribution f 1 (v). In contrast, by comparing two (or more) velocity distributions reconstructed from different recoil spectra with different target nuclei, one could avoid using these assumptions and estimate the WIMP mass model-independently.
The remainder of this article is organized as follows. In Sec. 2 I first review a method for determining the WIMP mass with only one direct detection experiment. In Sec. 3 I present a model-independent method for determining m χ by combining two experimental data sets. Numerical results based on Monte Carlo simulations of future experiments and uncertainties and caveats of these two methods will also be discussed. I conclude in Sec. 4. Some technical details for the data analysis will be given in an appendix.

With one experiment
In this section I review the method for determining the WIMP mass with only one direct detection experiment based on a maximum likelihood analysis [8,9,10,11].

Maximum likelihood analysis
I first describe briefly some (standard) theoretical models/assumptions for fitting the elastic WIMP-nucleus scattering spectrum to experimental data. Then I discuss the determination of the WIMP mass by a maximum likelihood analysis. Note here that only the most commonly used models/assumptions are described as examples to show which information is required for the maximum likelihood analysis; however, it should be understood that other models or assumptions can also be used.

Simple model distributions
The simplest semi-realistic model halo is a Maxwellian halo. The one-dimensional velocity distribution function in the rest frame of our Galaxy can be expressed as [6,1,7] Here v 0 ≃ 220 km/s is the orbital velocity of the Sun in the Galactic frame, and is the normalization constant which satisfies Note that the second term on the right-hand side of Eq. (9) has been introduced to keep the velocity distribution continuous at v = v esc . Substituting Eq.(9) into Eq.(1), the integral over the velocity distribution function can be calculated as vesc v min where v min = α √ Q in Eq.(4) has been used. Note that, in the v esc → ∞ limit, N Gau → 4/ √ πv 3 0 and the integral approaches to (2/ √ πv 0 ) e −α 2 Q/v 2 0 . On the other hand, when we take into account the orbital motion of the Solar system around the Galaxy as well as that of the Earth around the Sun, the velocity distribution function should be modified to [6,1,7] for v ≤ v esc , with the normalization constant Here v e is the Earth's velocity in the Galactic frame [5,1,2]: v e (t) = v 0 1.05 + 0.07 cos 2π(t − t p ) 1 yr ; (15) t p ≃ June 2nd is the date on which the velocity of the Earth relative to the WIMP halo is maximal. Consequently, an analytic form of the integral over this velocity distribution can be given as vesc v min For practical, numerical uses, an approximate form of the integral over f 1 (v) was introduced as [6] vesc v min where c 0 and c 1 are two fitting parameters of order unity. Not surprisingly, their values depend on the Galactic orbital and escape velocities, the target nucleus, the threshold energy of the experiment, as well as on the mass of incident WIMPs. Note that, the characteristic energy and thus the shape of the recoil spectrum depend highly on the WIMP mass: for light WIMPs (m χ ≪ m N ), Q ch ∝ m 2 χ and the recoil spectrum drops sharply with increasing recoil energy, while for heavy WIMPs (m χ ≫ m N ), Q ch ∼ const. and the spectrum becomes flatter.

Local WIMP density
Currently, the most commonly used value for the local WIMPs density in Eq.(2) is given as [1,2] However, so far it can be estimated only by means of the measurement of the rotational velocity of our Galaxy. Due to our location inside the Milky Way, it is more difficult to measure the accurate rotation curve of our own Galaxy than those of other galaxies. Thus an uncertainty of around a factor of 2 has been usually adopted [1,2]:

Spin-independent WIMP-nucleus cross section
In most theoretical models, the spin-independent (SI) WIMP interaction on a nucleus with an atomic mass number A > ∼ 30 dominates the spin-dependent (SD) interaction [1,2]. Additionally, for the lightest supersymmetric neutralino, which is perhaps the best motivated WIMP candidate [1,2], and for all WIMPs which interact primarily through Higgs exchange, the SI scalar coupling is approximately the same on both protons p and neutrons n. The "pointlike" cross section σ 0 in Eq.(2) can thus be written as where and f p is the effective χχpp four-point coupling, A is the atomic mass number of the target nucleus.

Nuclear form factor
For the SI cross section, an analytic nuclear form factor can be used. The simplest one is the exponential form factor, first introduced by Ahlen et al. [12] and Freese et al. [5]: Here Q is the recoil energy transferred from the incident WIMP to the target nucleus, is the nuclear coherence energy and is the radius of the nucleus. The exponential form factor implies a Gaussian form of the radial density profile of the nucleus. This Gaussian density profile is simple, but not very realistic. Engel has therefore suggested a more accurate form factor [13], inspired by the Woods-Saxon nuclear density profile [1,2], Here j 1 (x) is a spherical Bessel function, is the transferred 3-momentum, is the effective nuclear radius 2 with 3 and is the nuclear skin thickness. 2 In the literature, the form factor given in Eq.(26) is also known as the "Helm" form factor with [14,6] where R A ≃ 1.23 A 1/3 − 0.6 fm, r 0 ≃ 0.52 fm, s ≃ 0.9 fm.

Extended likelihood function
Now we are ready to put all pieces for predicting the elastic WIMP-nucleus scattering spectrum together and then fit this spectrum to experimental data by maximizing the logarithm of the extended likelihood function [10]: Here is the expected event number with the (assumed) exposure of the experiment, E, N tot is the total number of events recorded in one (simulated) experiment, Q a are measured recoil energies in the data set between the minimal and maximal cut-off energies, Q min and Q max , and is the total event rate. Note that, firstly, the definition of L in Eq.(34) takes into account the fact that the event number N tot and the measured recoil spectrum E(dR/dQ) of each (simulated) experiment are not fixed. Secondly, except c 0 and c 1 in Eq. (17), there are two fitting parameters in the extended likelihood function L, i.e., the WIMP mass m χ (involved in α) and the SI WIMP-proton cross section σ SI χp .

Numerical results
Here I show some numerical results with 10,000 simulated experiments based on Monte Carlo simulations performed by A. Green [10,11]. 76 Ge has been chosen as the target nucleus with a threshold energy of 10 keV. A three-dimensional Maxwellian velocity distribution in the Galactic rest frame for an isotropic isothermal WIMP halo, taking into account the Earth's motion around the Sun with v 0 = 220 km/s and v esc = 540 km/s, and the Helm form factor in Eqs.(26), (27), (29), and (30) have been used. The standard assumption for the local WIMP density of 0.3 GeV/cm 3 has been adopted. Note that the simulations demonstrated here as well as in the next section for the method combining two experimental data sets are based on several simplified assumptions 4 . Firstly, the sample to be analyzed contains only signal events, i.e., is free of background. Active background suppression techniques [16,17,18] 5 should make this condition possible. Secondly, all experimental systematic uncertainties as well as the uncertainty on the measurement of the recoil energy have been ignored. The energy resolution of most existing detectors is so good that its error can be neglected compared to the statistical uncertainty for the foreseeable future.

Statistical uncertainty
Figs. 1 show the distributions of the best-fit WIMP mass m χ and SI WIMP-proton cross section σ SI χp on the cross section versus WIMP mass plane. The input WIMP mass and the cross Figure 1: Distributions of the best-fit WIMP mass and SI WIMP-proton cross section on the cross section versus WIMP mass plane. The input WIMP mass and the cross section are 100 GeV and 10 −7 pb, respectively. The exposures have been assumed to be 3 ×10 3 (left) and 3 ×10 4 (right) kg-day and the corresponding expected event numbers are 78 and 780, respectively. In each frame, the contours contain 68% and 95% of the simulated experiments. See the text for further details (Plots from [10]).
section are 100 GeV and 10 −7 pb, respectively. The exposures have been assumed to be 3 × 10 3 (left) and 3 × 10 4 (right) kg-day and the corresponding expected event numbers are 78 and 780 6 , respectively. It can be seen that, especially for the smaller exposure, the distribution is asymmetric and there are (significantly) more experiments with best-fit masses and cross sections larger than the input values. Quantitatively, for a WIMP mass of 100 GeV with ∼ 80 events, the 1σ and 2σ statistical uncertainties are +40 −35 GeV and +100 −50 GeV, respectively [10]. Fig. 2 shows the 95% (solid) and 68% (dotted) confidence limits on the best-fit WIMP mass as functions of the input WIMP mass. The input SI WIMP-proton cross section has been set here as 10 −8 pb. The assumed exposures are 3 × 10 3 , 3 × 10 4 , and 3 × 10 5 kg-day, respectively. We see here that since, as mentioned above, the shape of the recoil spectrum varies significantly with the WIMP mass for light WIMP masses (m χ < m N ), the WIMP mass (and also the cross section) can be fitted with a higher accuracy: the 1σ and 2σ statistical uncertainties for m χ = 25 GeV are ±4 GeV and +8 −7 GeV, for m χ = 50 GeV are +15 −12 GeV and +22 −19 GeV, respectively [10]. In contrast, the weak dependence of the shape of the recoil spectrum on the WIMP mass for heavy WIMP masses (m χ ≫ m N ) means that it will be more difficult or even impossible to extract the WIMP mass with O(100) events, if WIMPs are (much) heavier than the target nucleus [10]. Note that the dependence of the shape of the recoil spectrum on the WIMP mass as well as on that of the target nucleus suggests that heavy nuclei, e.g., Xe, would be able to measure the mass of heavy WIMPs more accurately; however, the rapid decrease of the nuclear form factor with increasing recoil energy, which occurs for heavy nuclei, means that, due to less expected events, this is in fact not necessarily the case.

Systematic uncertainties
Different sources of the systematic uncertainties in this model-dependent analysis have been considered [10,11]. Figs. 3 show the distributions of the best-fit WIMP mass and cross section with different input orbital velocity of the Solar system: v 0 = 200 (left) and 240 (right) km/s, while the standard value of v 0 = 220 km/s has been used for the data analysis. As shown here, for an input WIMP mass of 100 GeV, there could be an ∼ ±20 GeV shift in the best-fit WIMP mass combined with an ∼ ±10 −8 pb (∼ 10%) shift in the SI WIMP-proton cross section caused by the ± 20 km/s difference between the real and the assumed orbital velocities [10]. Moreover, the larger the real orbital velocity, the less the expected event number (with a fixed exposure), and thus the larger the statistical uncertainties on both the WIMP mass and SI WIMP-proton cross section one could obtain.
More detailed illustrations and discussions about the effects of varying the underlying WIMP mass and cross section, the detector target nucleus, the exposure, the minimal and maximal cutoff energies, the orbital velocity of the Solar system, as well as the background event rate and its spectrum can be found in Refs. [10,11,19].

Combining two experiments
In this section I first review the model-independent method for reconstructing the WIMP mass by using two experimental data sets with different target nuclei 7 . Then I also describe an extension of this method for estimating (or at least constraining) the SI WIMP-proton cross section.

Model-independent determination
As mentioned in the introduction, the normalized one-dimensional velocity distribution function of incident WIMPs can be solved from Eq.(1) directly and, consequently, its generalized moments can be estimated by [20] v Here v(Q) = α √ Q, Q (min,max) are the experimental minimal and maximal cut-off energies, is an estimated value of the measured recoil spectrum (dR/dQ) expt (before the normalization by the exposure E) at Q = Q min , and I n (Q min , Q max ) can be estimated through the sum: where the sum runs over all events in the data set that satisfy Q a ∈ [Q min , Q max ]. Note that, firstly, by using the second Eq.(37) v n (v(Q min ), v(Q max )) can be determined independently of the local WIMP density ρ 0 , of the velocity distribution function of incident WIMPs, f 1 (v), as well as of the WIMP-nucleus cross section σ 0 . Secondly, as shown later, r(Q min ) and I n (Q min , Q max ) are two key quantities for this model-independent method, which can be estimated either from a functional form of the recoil spectrum or from experimental data (i.e., the measured recoil energies) directly 8 . However, r(Q min ) and I n (Q min , Q max ) estimated from a scattering spectrum fitted to experimental data are not model-independent any more.

Basic expressions for determining m χ
By requiring that the values of a given moment of f 1 (v) estimated by Eq.(37) from two detectors with different target nuclei, X and Y , agree, m χ appearing in the prefactor α n on the right-hand side of Eq.(37) can be solved as [21]: where and R n,Y can be defined analogously. Here n = 0, m (X,Y ) and F (X,Y ) (Q) are the masses and the form factors of the nucleus X and Y , respectively, and r (X,Y ) (Q min,(X,Y ) ) refer to the counting rates for detectors X and Y at the respective lowest recoil energies included in the analysis. Note that, firstly, the general expression (40) can be used either for spin-independent or for spindependent scattering, one only needs to choose different form factors under different assumptions. Secondly, the form factors in the estimate of I n,X and I n,Y using Eq.(39) are also different.
On the other hand, by using the theoretical prediction that the SI WIMP-nucleus cross section given in Eq. (21) dominates, and the fact that the integral over the one-dimensional WIMP velocity distribution on the right-hand side of Eq.(1) is the minus-first moment of this distribution, which can be estimated by Eq.(37) with n = −1, one can easily find that [20] Note that the exposure of the experiment, E, appears in the denominator. Since the unknown factor ρ 0 |f p | 2 on the left-hand side above is identical for different targets, it leads to a second expression for determining m χ [20] Here m (X,Y ) ∝ A (X,Y ) has been assumed, min,X r X (Q min,X ) F 2 X (Q min,X ) and similarly for R σ,Y .

χ 2 -fitting
In order to yield the best-fit WIMP mass as well as to minimize its statistical error by combining the estimators for different n in Eq.(40) with each other and with the estimator in Eq.(43), a χ 2 function has been introduced [20] where for i = −1, 1, 2, . . . , n max , and the other n max + 2 functions f i,Y can be defined analogously. Here n max determines the highest moment of f 1 (v) that is included in the fit. The f i are normalized such that they are dimensionless and very roughly of order unity in order to alleviate numerical problems associated with the inversion of their covariance matrix. Note that the first n max + 1 fit functions depend on m χ only through the overall factor α and that m χ in Eqs.(46a) and (46b) is now a fit parameter, which may differ from the true value of the WIMP mass. Finally, C in Eq.(45) is the total covariance matrix. Since the X and Y quantities are statistically completely independent, C can be written as a sum of two terms 9 :

Matching the cut-off energies
The basic requirement of the expressions for determining m χ given in Eqs. (40) and (43) Note that α defined in Eq.(5) is a function of the true WIMP mass. Thus this relation for matching optimal cut-off energies can be used only if m χ is already known. One possibility to overcome this problem is to fix the cut-off energy of the experiment with the heavier target, minimize the χ 2 (m χ ) function defined in Eq.(45), and estimate the cut-off energy for the lighter nucleus by Eq.(48) algorithmically [20].

Numerical results
Here I show some numerical results for the reconstructed WIMP mass based on Monte Carlo simulations. The upper and lower bounds on the reconstructed WIMP mass are estimated from the requirement that χ 2 exceeds its minimum by 1 11 . 28 Si and 76 Ge have been chosen as two target nuclei. The scattering cross section has been assumed to be dominated by spin-independent interactions. The shifted Maxwellian velocity distribution given in Eq.(13) (the second term involving v esc has been neglected) with v 0 = 220 km/s, v e = 1.05 v 0 12 , and v esc = 700 km/s and the Woods-Saxon form factor in Eq.(26) have been used. The threshold energies of two experiments have been assumed to be negligible and the maximal experimental cut-off energies are set as 100 keV. 2 × 5,000 experiments have been simulated. In order to avoid large contributions from very few events in the high energy range to the higher moments [7], only the moments up to n max = 2 were included in the χ 2 fit.

Statistical uncertainty
In Figs. 4 the dotted (green) curves show the median reconstructed WIMP mass and its 1σ upper and lower bounds for the case that both Q max,Si and Q max,Ge have been fixed to 100 keV. As argued earlier, the values of a given moment of the WIMP velocity distribution estimated by Eq.(37) do not agree when the same maximal cut-off energy for both experimental data sets is used. This causes a systematic underestimate of the reconstructed WIMP mass [21] which can be seen obviously here.
The solid (black) curves were obtained by using Eq.(48) for matching the cut-off energy Q max,Si perfectly with Q max,Ge = 100 keV and the true (input) WIMP mass, whereas the dashed (red) curves show the case that Q max,Ge = 100 keV, and Q max,Si has been determined by minimizing χ 2 (m χ ; Q max,Si ). As shown here, with only 50 events on average before cuts (upper frame) from each experiment, the algorithmic process seems already to work pretty well for WIMP masses up to ∼ 500 GeV. For m χ < ∼ 100 GeV the median WIMP mass determined in this way overestimates its true value by 15 to 20%; however, the true WIMP mass always lies within the median limits of the 1σ statistical error interval estimated by the algorithmic Q max matching procedure up to even m χ = 1 TeV [20].

Statistical fluctuation
In order to study the statistical fluctuation of the reconstructed WIMP mass by algorithmic Q max matching in the simulated experiments, an estimator δm has been introduced as [20  Here m χ,in is the true (input) WIMP mass, m χ,rec its reconstructed value, m χ,lo1(2) are the 1 (2) σ lower bounds satisfying χ 2 (m χ,lo(1,2) ) = χ 2 (m χ,rec ) + 1 (4), and m χ,hi1 (2) are the corresponding 1 (2) σ upper bounds. It has been found that the error intervals of the median reconstructed WIMP mass are quite asymmetric; similarly, the distance between the 2σ and 1σ limits can be quite different from the distance between the 1σ limit and the central value [20] 13 . The definition of δm in Eq.(49) takes these differences into account, and also keeps track of the sign of the deviation: if the reconstructed WIMP mass is larger (smaller) than the true one, δm is positive (negative). Moreover, |δm| ≤ 1 (2) if and only if the true WIMP mass lies between the experimental 1 (2) σ limits. Fig. 5 shows the distribution of δm calculated from 5,000 simulated experiments with 50 events on average before cuts for a rather light WIMP mass of 50 GeV. In this case simply fixing both Q max values to 100 keV still works fine (see the upper frame of Figs. 4). However, the distributions for both fixed Q max and optimal Q max matching look somewhat lopsided, since the error interval is already asymmetric, with m χ,hi1 − m χ,rec > m χ,rec − m χ,lo1 . The overestimate of light WIMP masses reconstructed by algorithmic Q max matching shown in Figs. 4 is reflected by the dashed (red) histogram here, which has significantly more entries at positive values than at negative values. These distributions also indicate that the statistical uncertainties estimated by minimizing χ 2 (m χ ) are indeed overestimated, since nearly 90% of the simulated experiments have |δm| ≤ 1 [20], much more than ∼ 68% of the experiments, that a usual 1σ error interval should contain. 13 Recall that the same asymmetry has also been observed by the maximum likelihood analysis.  Fig. 5, except that the input WIMP mass has been increased to 200 GeV. In the right frame the average event number (before cuts) in each experiment have, in addition, been increased from 50 to 500. Note that the bins at δm = ±5 are overflow bins, i.e., they also contain all experiments with |δm| ≥ 5. (Plots from Ref. [20]).
Unfortunately, as shown in Figs. 6, when the true (input) WIMP mass increases to 200 GeV and the expected event number (before cuts) increases to 500 (right frame), the situations become less favorable. While optimal Q max matching seems to approach very slowly to be Gaussian and the overestimated statistical errors become a little bit more reliable for larger event numbers [20], the errors estimated by the algorithmic procedure for determining Q max,Si are not very reliable in the simulations.
More detailed illustrations and discussions about algorithmic Q max matching with different detector materials or with data sets generated in different halo models, as well as about the statistical fluctuation in the analysis can be found in Ref. [20].

Estimating the SI WIMP-proton coupling
In the maximum likelihood analysis discussed in Sec. 2, the SI WIMP-proton cross section σ SI χp is the second fitting parameter that, combined with the WIMP mass m χ , maximizes the extended likelihood function L calculated from an assumed WIMP velocity distribution.
In contrast, as shown above, by combining two experimental data sets, one can estimate the WIMP mass m χ without knowing the WIMP-nucleus cross section σ 0 . Conversely, by means of Eq.(42), one can also estimate or at least constrain the SI WIMP-proton coupling, |f p | 2 , from experimental data directly without knowing the WIMP mass [22].

Making an assumption for the local WIMP density
In Eq.(42) the WIMP mass m χ on the right-hand side can be determined by the method described above, r(Q min ) and I 0 can also be estimated from one of the two data sets used for determining m χ or from a third experiment. Nevertheless, due to the degeneracy between the local WIMP density ρ 0 and the coupling |f p | 2 , one cannot estimate both of them independently. The simplest way is making an assumption for the local WIMP density ρ 0 14 .

Numerical results
The left frame of Figs. 7 shows the reconstructed SI WIMP-proton coupling |f p | 2 rec as a function of the input WIMP mass m χ,in . Following simulations for the reconstruction of the WIMP mass, 28 Si and 76 Ge were chosen as two target nuclei for estimating m χ in Eq.(42). In order to avoid complicated calculations of the correlation between the error on the reconstructed m χ and that on the estimator of I 0 , a second, independent data set with 76 Ge was chosen as the third target for estimating I 0 . The SI WIMP-proton cross section was set as 10 −8 pb. Each experimental data set has 50 events on average under the common experimental cut-off energy Q max chosen as 100 GeV.
It can be seen that the reconstructed |f p | 2 values are underestimated for WIMP masses > ∼ 100 GeV. This systematic deviation is caused mainly by the underestimate of I 0 . However, in spite of this systematic deviation (and in fact due to the fairly large statistical uncertainty), the true value of |f p | 2 always lies within the 1σ statistical error interval. Moreover, for a WIMP mass of 100 GeV, one could in principle already estimate the SI WIMP-proton coupling with a statistical uncertainty of only ∼ 15% with just 50 events from each experiment. Recall that this is much smaller than the systematic uncertainty of the local Dark Matter density (of a factor of 2 or even larger).
Combining the estimate for the SI WIMP-proton coupling with the estimate for the WIMP mass, the right frame of Figs. 7 shows the reconstructed coupling |f p | 2 rec and the reconstructed WIMP mass m χ,rec on the cross section (coupling) versus WIMP mass plane 15 . It is important to note that, as shown here, |f p | 2 and m χ can be estimated separately and from experimental data directly with neither prior knowledge of each other nor an assumption for the WIMP velocity distribution.

Summary and conclusions
In this article I reviewed the methods for the determination(s) of the mass (and eventually the spin-independent cross section on nucleons) of Weakly Interacting Massive Particles with positive signals of their elastic scattering off target nuclei in direct Dark Matter detection experiments.
With only one experiment, the WIMP mass combined with its SI cross section on nucleons could be estimated by the maximum likelihood analysis using a theoretically predicted scattering spectrum fitted to the measured recoil energies. If WIMPs are light (m χ < m N ), the shape of the recoil spectrum is sensitive to their mass, then the WIMP mass (and also the cross section) can be estimated with a higher accuracy; however, in case WIMPs are (much) heavier than the target nucleus (m χ > ∼ 200 GeV), the recoil spectrum becomes nearly independent on m χ and it is then more difficult or even impossible to estimate the WIMP mass reasonably with O(100) events.
The maximum likelihood analysis depends on the prior assumption for the velocity distribution of halo WIMPs as well as on the local WIMP density. For a WIMP mass of 100 GeV, an ∼ 10% measurement uncertainty on the orbital velocity of the Solar system could cause an ∼ 20% systematic error on the best-fit WIMP mass combined with an ∼ 10% error on the SI WIMP-proton cross section.
In order to determine the WIMP mass without making any assumption for the WIMP velocity distribution, I described a second method based on the reconstruction of (the moments of) the WIMP velocity distribution function from two experiments with different target nuclei. This method can be used without knowing the WIMP-nucleus cross section. The only information needed is the measured recoil energies. By matching the maximal cut-off energies of two experiments one could in principle estimate the WIMP mass up to ∼ 500 GeV with O(50) events from each experiment.
Nevertheless, the algorithmic procedure for determining the maximal cut-off energy of the experiment with the lighter target nucleus by minimizing χ 2 could overestimate the WIMP mass by 15 to 20% if WIMPs are light, or lead to unreliable error estimates if WIMPs are heavy. The latter could become worse with larger event samples. However, the fact that optimal Q max matching works well in all cases, for both the median reconstructed WIMP mass and its statistical error, gives us hope that a better algorithm for Q max matching can be found which only relies on the data.
Additionally, by combining two (or three) experimental data sets one could also estimate the spin-independent WIMP-proton coupling without knowing the WIMP mass. Although, due to the degeneracy between the local WIMP density and the WIMP-nucleus cross section, one needs to adopt the local Dark Matter density (as the unique assumption), at least an upper bound on this coupling could be given. In fact, for a WIMP mass of 100 GeV, with O(50) events from each experiment, a statistical uncertainty of ∼ 15% could be reached. This is much smaller than the systematic uncertainty on the local Dark Matter density (of a factor of 2 or even larger).
In summary, by means of currently running and projected experiments using detectors with 10 −9 to 10 −11 pb sensitivities [16,17,18] (see footnote 5), we stand a good chance of detecting Dark Matter particles, if Dark Matter indeed consists (mainly) of WIMPs. Then the methods presented here can be used to estimate the mass (and eventually the cross section on nucleons) of Dark Matter particles. This information (perhaps combined with information from indirect detection experiments [19]) will allow us not only to constrain the parameter space in different extensions of the Standard Model of particle physics, but also to identify WIMPs among new particles produced at colliders (hopefully in the near future). Once one is confident of this identification, one can use further collider measurements of the mass and couplings of WIMPs. Together with the reconstruction of the velocity distribution of halo WIMPs [7], this will then yield a new determination of the local WIMP density. On the other hand, knowledge of the WIMP couplings will also permit prediction of the WIMP annihilation cross section. Together with information on the WIMP density, this will allow one to predict the event rate in the indirect Dark Matter detection [1,2] as well as to test our understanding of the early Universe. since Finally, since all I n are determined from the same data, they are correlated with where the sum again runs over all events with recoil energy between Q min and Q max . And the correlation between the errors on r(Q min ), which is calculated entirely from the events in the first bin, and on I n is given by cov(r(Q min ), I n ) = r(Q min ) I n (Q min , Q min + b 1 ) note that the sums I i here only count in the first bin, which ends at Q = Q min + b 1 .
On the other hand, with a functional form of the recoil spectrum (e.g., fitted to experimental data), (dR/dQ) expt , one can use the following integral forms to replace the summations given above. Firstly, the average Q−value in the nth bin defined in Eq.(A5) can be calculated by For I n (Q min , Q max ) given in Eq.(39), we have and similarly for the covariance matrix for I n in Eq.(A12), Remind that (dR/dQ) epxt is the measured recoil spectrum before the normalization by the exposure. Finally, I i (Q min , Q min + b 1 ) needed in Eq.(A13) can be calculated by Note that r(Q min ) and I n (Q min , Q min + b 1 ) should be estimated by Eqs.(A9) and (A17) with r 1 , k 1 and Q s,1 estimated by Eqs.(A3), (A4), and (A8) in order to use the other formulae for estimating the (correlations between the) statistical errors without any modification.
A.2 Statistical errors on m χ given in Eqs. (40) and (43) The expression for m χ | v n given in Eq.(40) leads to a lengthy expression for its statistical error: Here a short-hand notation for the six quantities on which the estimate of m χ depends has been introduced: c 1,X = I n,X , c 2,X = I 0,X , c 3,X = r X (Q min,X ) ; and similarly for the c i,Y . Estimators for cov(c i , c j ) have been given in Eqs.(A12) and (A13). Explicit expressions for the derivatives of R n,X with respect to c i,X are: ∂R n,X ∂I n,X = n + 1 n   F 2 X (Q min,X ) 2Q (n+1)/2 min,X r X (Q min,X ) + (n + 1)I n,X F 2 X (Q min,X )   R n,X , and ∂R n,X ∂r X (Q min,X ) = 2 n   Q (n+1)/2 min,X I 0,X − (n + 1)Q 1/2 min,X I n,X 2Q (n+1)/2 min,X r X (Q min,X ) + (n + 1)I n,X F 2 X (Q min,X )   ×   F 2 X (Q min,X ) 2Q 1/2 min,X r X (Q min,X ) + I 0,X F 2 X (Q min,X )   R n,X ; (A20c) explicit expressions for the derivatives of R n,Y with respect to c i,Y can be given analogously. Note that, firstly, factors R n,(X,Y ) appear in all these expressions, which can practically be cancelled by the prefactors in the bracket in Eq.(A18). Secondly, all the I 0,X , I 0,Y , I n,X , I n,Y should be understood to be computed according to Eqs.(39) or (A15) with integration limits Q min and Q max specific for that target. Similar to the analogy between Eqs. (40) and (43), the statistical error on m χ | σ given in Eq.(43) can be expressed as where we have again used the short-hand notation in Eq.(A19); note that c 1,(X,Y ) = I n,(X,Y ) do not appear here. Expressions for the derivatives of R σ,X can be computed from Eq.(44) as ∂R σ,X ∂I 0,X =   F 2 X (Q min,X ) 2Q 1/2 min,X r X (Q min,X ) + I 0,X F 2 X (Q min,X )   R σ,X , ∂R σ,X ∂r X (Q min,X ) =   2Q 1/2 min,X 2Q 1/2 min,X r X (Q min,X ) + I 0,X F 2 X (Q min,X )   R σ,X ; and similarly for the derivatives of R σ,Y .

A.3 Covariance of f i defined in Eqs.(46a) and (46b)
The entries of the C matrix in Eq.(47) involving basically only the moments of the WIMP velocity distribution can be read off Eq.(82) of Ref. [7], with an slight modification due to the normalization factor in Eq.(46a) 16 : cov (f i , f j ) = N 2 m f i f j cov(I 0 , I 0 ) + α i+j (i + 1)(j + 1)cov(I i , I j ) − α j (j + 1)f i cov(I 0 , I j ) − α i (i + 1)f j cov(I 0 , I i ) + D i D j σ 2 (r(Q min )) − (D i f j + D j f i ) cov(r(Q min ), I 0 ) + α j (j + 1)D i cov(r(Q min ), I j ) + α i (i + 1)D j cov(r(Q min ), I i ) .
(A23) 16 Since the last f i defined in Eq.(46b) can be computed from the same basic quantities, i.e., the counting rates at Q min and the integrals I 0 , it can directly be included in the covariance matrix.
Here N m ≡ 1 2Q 1/2 min r(Q min )/F 2 (Q min ) + I 0 , (A24) and for i = −1, 1, 2, . . . , n max ; and A.4 Statistical error on |f p | 2 given in Eq. (42) From Eq.(42), it can easily be found that where N m is defined in Eq.(A24), and The correlation between the error on the reconstructed m χ and that on the estimator of 1/N m , the third term in Eq.(A27), can be neglected in case one uses three independent data sets.