Analysis of the theoretical bias in dark matter direct detection

Fitting the model"A"to dark matter direct detection data, when the model that underlies the data is"B", introduces a theoretical bias in the fit. We perform a quantitative study of the theoretical bias in dark matter direct detection, with a focus on assumptions regarding the dark matter interactions, and velocity distribution. We address this problem within the effective theory of isoscalar dark matter-nucleon interactions mediated by a heavy spin-1 or spin-0 particle. We analyze 24 benchmark points in the parameter space of the theory, using frequentist and Bayesian statistical methods. First, we simulate the data of future direct detection experiments assuming a momentum/velocity dependent dark matter-nucleon interaction, and an anisotropic dark matter velocity distribution. Then, we fit a constant scattering cross section, and an isotropic Maxwell-Boltzmann velocity distribution to the simulated data, thereby introducing a bias in the analysis. The best fit values of the dark matter particle mass differ from their benchmark values up to 2 standard deviations. The best fit values of the dark matter-nucleon coupling constant differ from their benchmark values up to several standard deviations. We conclude that common assumptions in dark matter direct detection are a source of potentially significant bias.


Introduction
Dark matter constitutes about five sixth of the total matter in the observable Universe [1], and it forms large spheroidal halos hosting the majority of the astrophysical structures, including our Milky Way [2,3].
The particles forming the Milky Way dark matter halo have so far escaped detection, and the dark matter mass, as well as its interactions, remain unknown. The detection of dark matter particles in the solar neighborhood through their scattering off the nuclei of a terrestrial detector would provide us with a direct measurement of the dark matter particle mass and of the strength of the dark matter-nucleon interaction [4]. To accomplish this goal is the aim of the dark matter direct detection technique [5].
The direct detection technique has experienced a great advancement in the past few years. Current experiments can probe spin-independent dark matter-nucleon scattering cross sections at the level of 10 −45 cm 2 [6,7], or below [8], for weakly interacting dark matter candidates of mass around 30 GeV, and cross sections of approximately 10 −41 cm 2 , for masses in the 5-10 GeV range [9,10]. The next generation of direct detection experiments will exploit target masses at the ton-scale, improving current exclusion limits on the dark matter-nucleon scattering cross-section by about two order of magnitudes [11][12][13][14][15].
Fitting the dark matter particle mass and coupling constants to direct detection data, assumptions about the dark matter-nucleon interaction and the local dark matter velocity distribution are necessary 1 [22]. Dark matter-nucleon scattering cross sections independent of the momentum transfer and of the relative velocity, and a Maxwell-Boltzmann velocity distribution are common assumptions in this field [23]. The former assumption is motivated by the small velocity of the dark matter particles in the Milky Way, the latter by the simplicity of the Maxwell-Boltzmann distribution (i.e. a self-consistent distribution generated by the density profile of an isothermal sphere [24]). Although both assumptions are wellmotivated, other interaction types and velocity distributions are equally plausible, and could be considered in the data analysis [25].
In the dark matter-nucleon scattering, momentum and velocity independent interaction operators arise when dark matter couples to the nuclear charge density operator, or to the nuclear spin current density operator [26]. The former operator generates the so-called spin-independent interaction, the latter the so-called spin-dependent interaction. Despite the small velocity of the dark matter particles in the Milky Way, these operators do not necessarily generate the leading contributions to the non relativistic limit of the dark matter-nucleus scattering cross section. There are interesting examples where momentum and velocity independent interaction operators are forbidden by symmetry arguments, or parametrically suppressed [27][28][29][30][31][32][33][34].
In studying the velocity distribution of dark matter particles in the solar neighborhood, various approaches have been investigated, including Markov Chain Monte Carlo analyses [35][36][37][38] and N-body simulations [39][40][41]. As shown by the N-body simulations, the Maxwell-Boltzmann distribution only provides an approximate description of the particles forming the Milky Way dark matter halo. For instance, N-body simulations predict that the dark matter velocity distribution is anisotropic in a broad range of galactocentric distances, in contrast to the standard Maxwell-Boltzmann approximation [42]. Interesting polynomials expansions [43], and decompositions in streams [44] to model general dark matter velocity distributions have been recently proposed.
In summary, standard assumptions regarding the dark matter interactions, and velocity distribution might oversimplify the interpretation of future direct detection experiments, or even be incorrect, thereby introducing a bias in the data analysis. In this paper, we perform a quantitative study of the theoretical bias in dark matter direct detection, with a focus on assumptions regarding the dark matter-nucleon interaction, and the dark matter velocity distribution. We address this problem within the effective theory of isoscalar dark matternucleon interactions mediated by a heavy spin-1 or spin-0 particle [26,28].
First, we simulate future dark matter direct detection data assuming momentum and velocity dependent dark matter-nucleon interactions, and an anisotropic dark matter velocity distribution. Then, we analyze the simulated data assuming either a constant dark-matter nucleon scattering cross section, or an isotropic Maxwell-Boltzmann velocity distribution, hence introducing a bias in the analysis. Comparing the best fit points with the original benchmark points, we estimate the bias in dark matter direct detection potentially induced by incorrect theoretical assumptions.
We present the first study of the bias in dark matter direct detection based on the general effective theory of isoscalar dark matter-nucleon interactions. We cover both particle physics and astrophysical aspects of the problem. Previous studies restrict the analysis to the familiar spin-independent and spin-dependent interactions [45], or to a bias of astrophysical nature [46,47].
The paper is organized as follows. In Sec. 2, we briefly review the theoretical framework and the statistical methods adopted in the analyses. Sec. 3 is devoted to a description of the simulated data. In Sec. 4 we present our results. We conclude in Sec. 5. Useful dark matter response functions are listed in the appendix. Table 1. Non-relativistic operators appearing in Eq. (2.1). The operators O i are the same as in Ref. [30]. We adopt the following notation: 1 χ 1 N is the identity, q is the momentum transfer, v ⊥ χN is the χ-nucleon transverse relative velocity operator, S χ 1 N and 1 χ S N are the dark matter and nucleon spin operators, respectively. χ denotes the dark matter particle.

Dark matter-nucleon effective theory
The most general Lagrangian describing the dark matter-nucleon interactions arising from the exchange of a heavy spin-0 or spin-1 particle is given by the linear combination where χ + (χ − ) and N + (N − ) are the positive (negative) frequency parts of the non-relativistic dark matter and nucleon fields, respectively. The operators O i , with i = 1, 3, . . . , 11, are restricted by Galilean invariance, energy and momentum conservation and hermiticity. They are listed in Tab. 1. The interaction operators O 1 and O 4 generate the familiar spinindependent and spin-dependent interactions. In Eq. (2.1), c p i and c n i denote the coupling constants for protons and neutrons, respectively. Their linear combinations define the isoscalar (c 0 i ) and isovector (c 1 i ) coupling constants. In this paper we restrict our investigations to isoscalar interactions, i.e., we set c 1 i = 0. We collectively denote the isoscalar coupling constants by c. Following Ref. [30], we define the coupling constants c 0 i with dimension (mass) −2 .
From the Lagrangian in Eq. (2.1), we calculate the expected number of dark matter scattering events in a target material. The differential rate of scattering events per unit time and per unit detector mass is given by where ρ χ is the local dark matter density, m χ is the dark matter particle mass, and ξ T is the mass fraction of the nucleus T in the target material. In Eq. (2.3), P tot is the square modulus of the transition amplitude in the non-relativistic limit, i.e. |M N R | 2 , averaged over initial spins and summed over final spins. P tot can be written as a combination of nuclear and dark matter response functions: In Eq. (2.4), v and µ T are the dark matter-nucleus relative velocity and reduced mass, respectively, whereas y = (qb/2) 2 , q is the momentum transfer, and b is the oscillator parameter in the independent-particle harmonic oscillator model [30]. j N and j χ are the nucleus spin and the dark matter particle spin, respectively. For definiteness, we assume that the dark matter particle has spin j χ = 1/2.
are defied in Eq. (41) of Ref. [30]. We evaluate them using our FORTRAN version of the Mathematica package described in Ref. [30]. The dark matter response functions and R τ τ ∆Σ are given in the Appendix. In Eq. (2.3), the angle brackets denote the average where v min (q) = q/2µ T is the minimum velocity that a dark matter particle must have in order to transfer a momentum q to the target nucleus, and f is the local dark matter velocity distribution in the galactic rest frame boosted to the detector frame. In Eq. (2.6), v e (t) is the time-dependent Earth velocity in the galactic rest frame.
In the simulations of Sec. 3, we assume the self-consistent anisotropic velocity distribution proposed in Ref. [37]. This velocity distribution has been obtained from a generalized Eddington's inversion formula applied to the galactic model of Refs [35,48]. We set the relevant astrophysical parameters at their mean values, i.e. blue line in the left panel of Fig. 6 in Ref. [37]. In Sec. 4.2, we also use a Maxwell-Boltzmann distribution ) truncated at the local escape velocity v esc = 544 km s −1 , with v 0 = 220 km s −1 and ρ χ = 0.3 GeV cm −3 (i.e. the standard dark matter halo).

Statistical framework
We now introduce the statistical methods that we will use in Sec. 4 in order to analyze the synthetic data introduced in Sec. 3. We adopt frequentist and Bayesian statistical methods based on the likelihood function, L(d|Θ). Here d denotes an array of simulated data.
Θ ≡ (θ 1 , θ 2 ) is a 2-dimensional array composed of m χ and of one of the coupling constants c 0 i . The index i depends on the simulated data that we analyze (see Secs. 3 and 4 for concrete applications). The likelihood function that we use in the analyses is defined in Sec. 3.3. In the frequentist approach, we use the likelihood function to construct an effective chi-square defined as ∆χ 2 eff ≡ −2 ln L/L max , where L max is the absolute maximum of the likelihood function. From ∆χ 2 eff , we derive approximate 2D frequentist confidence intervals, which admit a rigorous statistical interpretation when the regularity conditions underlying Wilk's theorem apply. We maximize the likelihood function to obtain the best fit points of the model parameters.
In the Bayesian approach, we calculate the posterior probability density function (PDF) of the model parameters P(Θ|d), applying Bayes' theorem to the likelihood function L(d|Θ), (2.7) In Eq. (2.7), π(Θ) is the prior PDF. Any prior information on the model parameters must be included in this term. In absence of any empirical guideline, we assume log-priors for m χ and c. Tab. 2 shows the model parameters and the support of the associated prior PDFs. E(d) is the evidence, and in the present analysis it plays the role of a normalization constant. From the posterior PDF, we derive the posterior means of the model parameters. The posterior mean of the parameter θ 1 , for instance, is defined as follows Bayesian analog of the 2D frequentist confidence intervals are the 2D credible regions. 95% 2D credible regions, for instance, contain the 95% of the total posterior probability in the θ 1 -θ 2 plane. By construction, P(θ in 1 , θ in 2 |d) ≥ P(θ out 1 , θ out 2 |d), for (θ in 1 , θ in 2 ) inside the credible region, and (θ out 1 , θ out 2 ) outside the credible region. We define as "bias parameter" β of the variable X, the absolute value of the difference between the best fit value of X, i.e. X bf , and the value of X at the benchmark point, i.e. X P , divided by the error ∆: ∆ is defined as half of the width of the 95% confidence interval of the variable X. In Sec. 3 we provide a detailed description of the likelihood function used in the analyses. We take advantage of the Multinest program [49][50][51] with parameters set at n live = 20000 and tol=10 −4 in order to sample the likelihood function. We use our own routines to calculate direct detection scattering rates and evaluate the likelihood function. Finally, we use the programs Getplots [52] and Matlab to produce the figures.

Analysis
This section is devoted to the simulation of our synthetic direct detection data. First, we describe the ton-scale detectors adopted in the simulations (Sec. 3.1), then we explain how to simulate synthetic data from a given detector and a fiducial point in the parameter space of the theory defined in Sec. 2 (Sec. 3.2). We analyze these data in Sec. 4, using the frequentist and Bayesian statistical methods described in Sec. 2.2.

Type
Prior range Prior type log-prior For the benchmark points P 9 -P 24 we extend the prior range of log 10 (m χ /GeV) to 4.

Ton-scale detectors
In the analyses we assume two types of ton-scale detectors. A Germanium detector and a Xenon detector. For the ton-scale Germanium detector, we calculate the differential rate of dark matter scattering events per unit time and per unit detector mass as follows where [11] σ = (0.293) 2 + (0.056 is the energy dependent energy resolution of the detector, and E = 0.3 is the constant experimental efficiency. In Eq. (3.1), E O is the observed energy, E R is the true nuclear recoil energy, and dR/dE R is the the time average of the differential rate defined in Eq. (2.3). We calculate the total number of scattering events in the signal region (E inf , E sup ) as where M T = 1000 × 365 kg-day is the raw exposure of the detector, E inf = 10 keV and E inf = 100 keV. For the ton-scale Xenon detector, we calculate the differential spectrum of dark matter induced photoelectrons (PE), S 1 , as follows In Eq. (3.4), the Gaussian of standard deviation √ nσ PMT , with σ PMT = 0.37, and mean n gives the probability of observing S 1 PE, when n PE have been actually produced. The Poisson distribution of mean ν(E R ) gives the probability of producing n PE from a recoil energy E R . For the Xenon detector we assume the same ν(E R ), energy resolution, nuclear recoil acceptance, and efficiency of LUX [8]. In addition, we consider an exposure of M T = 1000 × 365 kg-day. In Eq. (3.4) dR/dE R denotes the the time average of the differential rate defined in Eq. (2.3). We calculate the total number of scattering events in the signal region (S inf 1 , S sup 1 ) as follows where S inf 1 = 2 PE and S sup 1 = 30 PE. For both detector types we assume the spectrum of irreducible background events, which includes a flat component and an exponentially decreasing component [11]. In Eq. (3.6), the index j identifies the detector type. j = 1 refers to the Germanium detector, whereas j = 2 to the Xenon detector. Accordingly,Ê 1 = E O ,Ê 2 = S 1 , a 1 = 10 keV, b 1 = 100 keV, 1 = 10 keV, a 2 = 2 PE, b 2 = 30 PE and, finally, 2 = 2 PE. For both detectors, we assume one background event in the signal region, which implies η = 0.5 in Eq. (3.6). Accordingly, for the detector of type j, we define the total number of expected events in the signal region as µ In the calculations, we include the most abundant isotopes for each target material, namely: 70 Ge, 72 Ge, 73 Ge, 74 Ge and 76 Ge, for Germanium and 128 Xe, 129 Xe, 130 Xe, 131 Xe, 132 Xe, 134 Xe, and 136 Xe for Xenon. We calculate the nuclear response functions of the different isotopes using our FORTRAN version of the Mathematica package introduced in Ref. [30].

Synthetic data
Given a detector of type j, j = 1, 2, and a benchmark point in the parameter space of the theory defined in Sec. 2, we simulate our synthetic direct detection data as follows. First, we randomly sample the number N j of observed scattering events from a Poisson distribution of mean µ (j) tot (m χ ,ĉ), where (m χ ,ĉ) is the benchmark point. Then, we randomly sample N j recoil energies, or PE, {Ê j } i=1,...,N j from the spectral function We apply this procedure to the detectors described in the previous subsection, and to the 24 benchmark points listed in Tab. 3. In the analyses, we repeat the procedure twice for each benchmark point; once for each detector type. From each benchmark point, we therefore produce two datasets d j , j = 1, 2, one for the Germanium detector and one for the Xenon detector.  Table 3. List of the benchmark points studied in this paper. For each benchmark point, we report the dark matter particle mass, and the value of the coupling constant different from zero at the benchmark point. For the Germanium (j = 1) and Xenon (j = 2) detectors, we also report µ (j) tot (m χ ,ĉ) and N j , respectively, the expected and the "observed" number of scattering events.
At the benchmark points of Tab. 3, the dark matter particle mass varies from 10 GeV to 250 GeV. The 24 benchmark points in Tab. 3 allow to explore the 8 momentum and velocity dependent interaction operators O 3 , O 5 , . . . O 11 . These benchmark points are compatible with present direct detection experiments, and have been proposed in Ref. [54], studying the prospects for direct detection of dark matter. Fig. 1 shows the 24 benchmark points in the 8 planes m χ -c 0 i , with i = 3, 5, . . . , 11.

The likelihood function
For each benchmark point, we generate two set of synthetic data, containing (N j +1) datapoints. One set of data for the ton-scale Germanium detector (j = 1) and one for the ton-scale Xenon detector (j = 2). We assume a Likelihood function for each dataset d j , j = 1, 2. The total Likelihood function of the independent datasets d j , j = 1, 2, is therefore given by The sum appearing in the second line of Eq. (3.10) contains the information on the spectrum of observed events, through the dependence on the functions f (j) , defined in Eq. (3.9).

Fitting procedures
In the analysis of direct detection experiments, it is commonly assumed that: 1. The dark matter particle interacts with the detector nuclei through the interaction operators O 1 and O 4 only (corresponding to the familiar spin-independent and spindependent interactions, respectively).
2. The local population of Milky Way dark matter particles is characterized by a Maxwell-Boltzmann velocity distribution.
Relying on these assumptions, when the actual theory of dark matter violates them, bias the interpretation of direct detection experiments. In order to quantify this bias, we proceed as follows. First, from each benchmark point in Tab. 3, we simulate synthetic direct detection data as explained in Sec. 3.2. Then, we analyze the synthetic data using the methods introduced in Sec. 2.2. We adopt four fitting procedures: 1. Fitting procedure A. Analyzing the data, we assume as dark matter-nucleon interaction operator the same operator O i used in the simulation, and as dark matter velocity distribution our benchmark anisotropic distribution.
2. Fitting procedure B. Analyzing the data, we assume O 1 as the correct dark matternucleon interaction operator, and our benchmark anisotropic distribution as dark matter velocity distribution.
3. Fitting procedure C. Analyzing the data, we assume O 4 as the correct dark matternucleon interaction operator, and our benchmark anisotropic distribution as dark matter velocity distribution.
4. Fitting procedure D. Analyzing the data, we assume as dark matter-nucleon interaction operator the same operator O i used in the simulation, and a Maxwell-Boltzmann distribution with parameters set as in Sec. 2 as dark matter velocity distribution.
Comparing the fitting procedures B, C and D with the fitting procedure A (and the original benchmark points), we estimate the bias in the interpretation of future direct detection experiments potentially induced by incorrect theoretical assumptions.

Results
In this section we analyze our synthetic direct detection data following the four fitting procedures described in Sec

Data analysis with a particle physics bias
In this section, we analyze the synthetic data generated from the benchmark points P 1 -P 24 adopting the fitting procedures A, B and C defined in Sec. 3.4.
We start with a detailed analysis of the benchmark points P 9 -P 16 . We report the results of this analysis in Fig. 2. The benchmark points P 9 -P 16 are characterized by a dark matter particle mass m χ = 50 GeV. Their properties are summarized in Tab. 3. Different panels in Fig. 2 correspond to distinct benchmark points. In the panels, the benchmark points are represented by cyan dots. For each benchmark point, we present our results in terms of 2D 95% confidence intervals in the m χ -f p plane. Red contours, black contours and blue contours correspond to the 2D 95% confidence intervals that we obtain applying, respectively, the fitting procedures A, B and C to the benchmark points P 9 -P 16 . In order to compare the Bayesian and the frequentist approach, in each panel we report the 2D posterior PDF that we derive analyzing the synthetic data with the fitting procedure B. Finally, green (magenta) crosses and stars in the panels identify the best fit points and the posterior means resulting from the fitting procedure B (C).  Table 4. List of bias parameters. For each benchmark point, we report the bias parameters of log 10 (m χ /GeV) and log 10 (f p ) resulting from four different analyses. We denote by b k (x) the bias parameter of log 10 (x), as obtained from the fitting procedure k, with k = A, B, C, D. Sec. 3.4 describes the four different fitting procedures.
First, we comment on the results that we obtain applying the fitting procedure A to the synthetic data generated from the benchmark points P 9 -P 16 . In this analysis, we find that the bias parameter of log 10 (m χ /GeV) is always less than 1, ranging from β = 0.02 at P 9 to β = 0.81 at P 10 . The bias parameter of log 10 (f p ) is larger and it varies from β = 0.25 at P 9 to β = 1.47 at P 14 (see Tab. 4). We conclude that the fitting procedure A leads to a good reconstruction of the benchmark points P 9 -P 16 .
We now apply the fitting procedure B to the synthetic data generated from the benchmark points P 9 -P 16 . The fitting procedure B allows us to study the bias induced by assuming O 1 as the correct dark matter nucleon interaction operator in the determination of the dark matter particle properties. Interestingly, the benchmark points P 9 -P 16 can be divided in two groups, depending on the accuracy within which m χ is determined in this analysis.
A first group includes the benchmark points P 9 , P 11 , P 15 and P 16 . At these benchmark points, the confidence intervals around the best fit values of m χ are large, if compared with the confidence intervals found with the fitting procedure A. In order to further investigate this difference between the fitting procedures A and B, we perform a chi-square goodness of fit test, using the best fit points found within the two approaches. First, we bin the   Figure 3. Synthetic data generated from the benchmark points P 11 (left panel) and P 13 (right panel). We bin the events measured by the Germanium detector in 10 energy bins, and the events measured by the Xenon detector in 10 PE bins. We compare the binned data to the expected number of events at each bin and for each detector. The expected number of events shown in the figures has been obtained integrating Eqs. (3.1) and (3.4). We set the model parameters at the best fit points found using the unbinned likelihood (3.11), and the fitting procedures A, B and C. The fitting procedures A, B and C assume that at P 11 (P 13 ) dark matter interacts with the detector nuclei through the interaction operators O 6 (O 8 ), O 1 and O 4 , respectively. To help the comparison, we have connected points representing different theoretical predictions with colored lines. The notation is the one explained in the legends. For each bin, error bars corresponds to ± the square root of the observed number of events. recoil energies observed by the Germanium detector in 10 energy bins, and the PE observed by the Xenon detector in 10 PE bins. Then, we calculate the expected number of events in the 10 energy bins, and in the 10 PE bins, setting the model parameters at the best fit points found with the fitting procedures A and B. Combining the binned data with the theoretical expectations, we construct a reduced chi-square χ 2 red with 18 degrees of freedom for the fitting procedures A and B. Using the synthetic data generated from the benchmark point P 11 , for instance, we find χ 2 red = 1.34 and χ 2 red = 1.88 for the fitting procedures A and B, respectively. We obtain similar results for the benchmark points P 9 , P 15 and P 16 . This study indicates that for m χ = 50 GeV, the operator O 1 = 1 χ 1 N cannot provide a good fit to the synthetic data generated from the interaction operators Here, we use the same notation introduced in the caption of Tab. 1. Fig. 3 (left panel) shows the binned data, and the expected number of events in each bin for the benchmark point P 11 .
A second group includes the benchmark points P 10 , P 12 , P 13 and P 14 . At these benchmark points, we find confidence intervals based on the fitting procedure B, comparable with those that we obtain using the fitting procedure A. A chi-square goodness of fit test applied to the synthetic data generated from these benchmark points gives as a result similar values  Figure 4. Same as for Fig. 2, but for the benchmark points P 1 -P 8 . In the panels, we do not report the best fit points, since they are essentially superimposed to the posterior means.
of χ 2 red for the fitting procedures A and B. For instance, at the benchmark point P 13 , we find χ 2 red = 1.13 and χ 2 red = 1.15 for the fitting procedures A and B, respectively. We obtain comparable values of χ 2 red for the benchmark points P 10 , P 12 and P 14 . Fig. 3 (right panel) reports the binned data, and the theoretical expectations for the benchmark point P 13 . In summary, for m χ = 50 GeV the operator O 1 provides a good fit to the synthetic data generated from the interaction operators The interaction operator O 1 provides a good fit to these data for the following reason. For m χ = 50 GeV, features in the energy recoil spectrum depending on the velocity or on the momentum transfer, which could distinguish the operators O 5 , O 7 , O 8 and O 9 from the operator O 1 , move below the experimental thresholds (see Fig. 3, right panel). We now comment on the bias parameters resulting from the fitting procedure B applied to the benchmark points P 9 -P 16 . At these benchmark points, the bias parameter of log 10 (m χ /GeV) ranges from β = 0.14 at P 13 to β = 2.49 at P 11 (see Tab. 4). Therefore, the bias resulting from the fitting procedure B is generically larger than the bias found with the fitting procedure A.
We interpret these large values of β as follows. Let us consider for a moment the realist case in which the correct dark matter-nucleon interaction is unknown. In presence of a candidate dark matter signal, a reasonable approach to data analysis consists in fitting the operators O i , i = 1, 3, . . . , 11 to the data independently. A posteriori, a goodness of fit test would identify the correct dark matter-nucleon interaction operator: one would interpret as the correct dark matter-nucleon interaction operator, the operator leading to the value of χ 2 red unambiguously closer to one. From this reasoning, we conclude that a large value of β alone does not necessarily indicate the presence of a real bias in the interpretation of the data. Only when β is large, say larger than ∼ 1, and the values of χ 2 red resulting from the fitting procedures A and B are comparable, there is a significant bias in the results, which cannot be identified from the data directly.
Interestingly, at the benchmark point P 12 , where dark matter interacts with the detector nuclei through the operator O 7 = S N · v ⊥ χN , the bias parameter of log 10 (m χ /GeV) is equal to 1.09, and the bias parameter of log 10 (f p ) is equal to 79.01. At this benchmark point, from a chi-square goodness of fit test we find χ 2 red = 1.15 and χ 2 red = 1.38, binning the synthetic data and using the best fit points of the fitting procedures A and B, respectively.
At the benchmark points P 9 -P 16 , the bias parameter of log 10 (f p ) is always much larger than 1, with the exception of the benchmark point P 16 , where β = 2.63 (see Tab. 4). A bias in log 10 (f p ) is clearly related to the fact that interaction operators other than O 1 and O 4 are momentum and velocity suppressed. Hence, a large value of f p is required in order to compensate for this suppression (see below, for an exception in the case of O 4 ).
The fitting procedure C applied to the synthetic data generated from the benchmark points P 9 -P 16 allows us to study the bias induced by assuming O 4 as the correct dark matter-nucleon interaction operator in the determination of m χ and f p . The results of this analysis are similar to those previously discussed for the case of the interaction operator O 1 . The only significant difference between the two cases is in the best fit value of the dark matter-nucleon coupling constant f p at the benchmark point P 16 . At P 16 , the best fit value of f p is larger than its benchmark value, contrary to what we observe at the other benchmark points. The reason is the following. At the benchmark point P 16 , the dark matter-nucleon interaction is described by the operator O 11 . This interaction operator generates the nuclear response function W τ τ M in the dark matter-nucleus scattering. W τ τ M is significantly larger than the nuclear response function associated with the operator O 4 , for Germanium and Xenon targets [26]. Therefore, a large value of f p = m 2 v c 0 4 /2 compensates for the lower nuclear response function generated by O 4 . Figs. 2 and 3 show, respectively, the 95% 2D confidence intervals and the results of a chi-square goodness of fit analysis based on the fitting procedure C. Tab. 4 contains the values of the bias parameters of log 10 (m χ /GeV) and log 10 (f p ) extracted from the fitting procedure C.
We conclude this subsection with an analysis of the remaining benchmark points. Figs. 4 and 5 show the 2D 95% confidence intervals that we obtain from an analysis of the synthetic data generated from the benchmark points P 1 -P 8 , and P 17 -P 24 , respectively. The benchmark points P 1 -P 8 are characterized by m χ = 10 GeV, whereas at the benchmark points P 17 -P 24 m χ = 250 GeV. For the benchmark points P 1 -P 8 , we find a small bias induced by the fitting procedures B and C in the determination of the dark matter particle mass. The corresponding bias parameter β varies from ∼ 0.15, in the case of the benchmark point P 5 to ∼ 0.75, in the case of the benchmark point P 7 , independently of whether we fit the operator O 1 or the operator O 4 to the synthetic data. This small bias is an expected result, since for light dark matter candidates features in the recoil energy spectrum distinguishing a momentum dependent operator from O 1 and O 4 move below the experimental thresholds.  Figure 6. Study of the benchmark points P 1 -P 8 in a compared analysis of the fitting procedures A and D (see Sec. 3.4 for a definition of the two fitting procedures). The former approach assumes an anisotropic dark matter velocity distribution in the interpretation of the synthetic data, the latter the standard dark matter halo approximation based on a Maxwell-Boltzmann distribution. Blue contours and red contours correspond to the 2D 95% confidence intervals that we obtain from the fitting procedures A and D, respectively. Stars denote the best fit points, whereas dots correspond to the benchmark points.
Finally, at the benchmark points P 17 -P 24 , we reconstruct the dark matter particle mass with large errors in all cases. In addition, the best fit values of log 10 (m χ /GeV) move towards very large dark matter masses, when we analyze the synthetic data adopting the fitting procedures B and C, with the exception of the benchmark points P 12 and P 13 . Tab. 4 shows a list of bias parameters for the benchmark points P 1 -P 8 , and P 17 -P 24 .

Data analysis with an astrophysical bias
We conclude Sec. 4 studying the benchmark points P 1 -P 8 in a compared analysis of the fitting procedures A and D. Our aim is to quantify the bias induced by the Maxwell-Boltzmann approximation to our benchmark anisotropic velocity distribution in the determination of the dark matter particle mass and coupling constants.
We focus on the benchmark points P 1 -P 8 , since at these points the dark matter particle is light. A light dark matter particle can transfer a detectable amount of energy to a target nucleus, only if it has a sufficiently high speed. Since the Maxwell-Boltzmann distribution mostly differs from our benchmark distribution in the high velocity tail (see Fig. 7 in Ref. [37]), we expect that the benchmark points P 1 -P 8 are the most affected by "an astrophysical bias" of the type mentioned above. Fig. 6 shows the 2D 95% confidence intervals that we obtain applying the fitting procedures A (blue contours) and D (red contours) to the synthetic data generated from the benchmark points P 1 -P 8 . Green and magenta stars represent the best fit points of the fitting procedures A and D, respectively. Distinct panels refer to different benchmark points. In all panels, the best fit values of m χ that we obtain with the fitting procedure A are larger, compared to those found with the fitting procedure D. The reason for this general trend is the following. The anisotropic velocity distribution of Ref. [37] predicts less dark matter particles with high speed in the galactic halo, as compared to the Maxwell Boltzmann approximation (again, see Fig. 7 in Ref. [37]). Hence, in order to produce a signal above the experimental thresholds, a larger value of the dark matter particle mass is necessary, when using the fitting procedure A.
The bias parameter of log 10 (m χ /GeV) that we obtain from the fitting procedure D ranges from 0.11 at P 5 , to 1.05 at P 1 . This result is in agreement with the conclusions of Ref. [46], where studying the interaction operator O 1 , the authors find a relative error in the reconstruction of the dark matter particle mass of the order of 1 standard deviation, for certain astrophysical configurations, and with a slightly different definition of bias. Similarly, the bias parameter of log 10 (f p ) varies from 0.07 at P 7 , to 0.87 at P 6 .
Comparing the results of this section with those of Sec. 4.1, we conclude that at present the most serious source of theoretical bias in dark matter direct detection is our poor knowledge of the dark matter-nucleon interaction.

Conclusions
Incorrect assumptions about the dark matter-nucleon interaction, or the local dark matter velocity distribution will bias the interpretation of future direct detection experiments. We have performed a quantitative study of the theoretical bias in dark matter direct detection, focusing on the assumption of momentum/velocity independent dark matter-nucleon interactions, and on the Maxwell-Boltzmann approximation.
We have addressed this problem within the effective theory of isoscalar dark matternucleon interactions mediated by a heavy spin-1 or spin-0 particle. From the 8 momentum and velocity dependent interaction operators predicted by the theory, we have simulated the data of future ton-scale Germanium and Xenon detectors, assuming an anisotropic dark matter velocity distribution. Interpreting the simulated data, we have instead made different assumptions, thereby introducing a theoretical bias in the analysis. We have either assumed a momentum and velocity independent dark matter-nucleon interaction, or a Maxwell-Boltzmann dark matter velocity distribution. Comparing the best fit points with the original benchmark points, we have hence estimated the theoretical bias in dark matter direct detection, using frequentist and Bayesian statistical methods.
We have found examples where the operator O i , i = 1, 4 generating the observed signal, the operator O 1 and O 4 fit the synthetic data equally well, in a chi-square goodness of fit test. However, fitting O 1 and O 4 to the synthetic data, we find a biased reconstruction of m χ and f p .
For instance, we have found that the signal produced in a ton-scale detector by a dark matter particle of mass 50 GeV, coupling to the nucleons through the operator O 7 = S N · v ⊥ χN with a strength f p = 300, can be confused with the signal produced by a dark matter particle of mass 33.8 GeV, coupling to the nucleons through the operator O 1 = 1 χ 1 N with a strength f p = 3 × 10 −4 (see Fig. 2, central-left panel). In this example, the familiar spin-independent interaction (i.e. O 1 ) provides a good fit to the synthetic data generated from O 7 , producing however significantly biased results.
In the last part of the paper, we have found that for a 10 GeV dark matter candidate, incorrect astrophysical assumptions induce an error in the reconstruction of m χ and f p of the order of 1 standard deviation.
We conclude that a potentially important source of theoretical bias inheres within common assumptions in the field of dark matter direct detection. Extracting the correct dark matter-nucleon interaction from a multi-dimensional analysis of the data directly will be important to fully exploit the next generation of direct detection data [54,55]. Marginalizing over (in a Bayesian approach) or profiling out (in a frequentist approach) the astrophysical uncertainties pertaining the dark matter direct detection could mitigate the astrophysical bias identified in this study [46,47].

A Dark matter response functions
In the following, we list the dark matter response functions appearing in Eq. (2.4). In the calculations, we assume for definiteness that the dark matter particle has spin j χ = 1/2.