Aboulaye Fall
email: abdoulaye.fall@telecom-sudparis.eu

Elodie Le Cren

Kevin Lenglé

Catherine Lepers

Yaneck Gottesman

Monique Thual

Laurent Bramerie

Denis Molin

Pierre Sansonetti

Dennis Van Ras

Versatile Graded-Index Multi-Mode Fiber for High Capacity Single-and Multi-Mode Optical Home Network

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Versatile Graded-Index Multi-Mode Fiber for High Capacity Single-and Multi-Mode Optical Home Network Aboulaye Fall [START_REF] Ph | Multiformat Home Networks using silica fibers[END_REF] , Elodie Le Cren (2) , Kevin Lengle (2) , Catherine Lepers [START_REF] Ph | Multiformat Home Networks using silica fibers[END_REF] , Yaneck Gottesman [START_REF] Ph | Multiformat Home Networks using silica fibers[END_REF] , Monique Thual (2) , Laurent Bramerie (2) , Denis Molin [START_REF] Sim | High-Speed Multi-mode Fiber Transmission by Using Mode-Field Matched Center-Launching Technique[END_REF] , Pierre Sansonetti [START_REF] Sim | High-Speed Multi-mode Fiber Transmission by Using Mode-Field Matched Center-Launching Technique[END_REF] , Dennis Van Ras [START_REF] Sim | High-Speed Multi-mode Fiber Transmission by Using Mode-Field Matched Center-Launching Technique[END_REF] , Michel Gadonna (2) , Charles Populaire [START_REF] Tonobe | Dual-Mode Fiber (DMF) capable of single-and multi-mode transmission[END_REF] , Georges Martin [START_REF] Tonobe | Dual-Mode Fiber (DMF) capable of single-and multi-mode transmission[END_REF] , Laurent Valencia [START_REF] Tonobe | Dual-Mode Fiber (DMF) capable of single-and multi-mode transmission[END_REF] , Philippe Guignard

Introduction

Several major evolutions will contribute to a drastic increase of the bit rate required in the Home Network in a near future. Continuously growing bit rate in the access networks with FTTH deployments will soon boost the bit rate demand up to 10 Gbit/s. Other technologies are also widely used to deliver contents, as terrestrial and satellite television. There is a real interest to transport these latter signals in the home in their native radiofrequency format, to ensure a total compatibility with the customer's devices as TV sets or satellite demodulators.

A unique home network able to carry all these signals in their different formats is of high interest, and silica optical fiber technology appears to be the only media providing the required capacity to support the various network systems generation 1 .

We report in this paper on the design, realization and characterization of a home network optical fiber (HNOF) able to be used like a multi-mode fiber (MMF) when using 850nm multi-mode sources, and ensure singlemode propagation when using single-mode sources operating at 1310nm or 1550nm. This solution allows supporting cost effective multimode systems as well as single-mode systems that would arise in the longer term due to bit rate increase and need of more flexibility 1 . This potential for further evolution is particularly important for a medium which will be installed in the walls of the house for decades.

Fiber Design

Single-mode and multi-mode transceivers, operating at 1310/1550nm and 850nm, are respectively ended with standard single-mode fiber (SMF) and standard 50µm multi-mode fiber (MMF) pigtails. The HNOF is thus excited by SMF or MMF. After propagation through the HNOF, the excited optical waves are coupled into SMF or MMF pigtails of the receivers.

The compatibility with multi-mode terminations constraints HNOF to be highly multi-mode at 850nm like OM3/4 multi-mode fibers in order to be used with high speed Vertical Cavity Surface Emitting Lasers (VCSEL), typically used at this operating wavelength. As a consequence, HNOF is also highly multi-mode at 1310 & 1550nm. The Effective Modal Bandwidth (EMB) requirements of 500MHz.km at 850nm have been calculated with the IEEE Spreadsheet model 2 to support 10Gb/s error-free transmission at 850nm over 100m (2000MHz.km is required to equal the system performance of OM3 fibers). Therefore, a graded-index profile is recommended.

Single-mode transmissions at high bit rate (≥ 10Gb/s) over MMF is limited by the modal dispersion, that can be efficiently, and passively, mitigated by using the mode-matched center launch technique 3 . In that case, only the fundamental mode is excited in the MMF. It can be achieved using a tapered fiber between SMF pigtail and HNOF 3 , or by adapting the refractive index profile of HNOF in order to match the fundamental mode field diameter (FMFD) to that . For sake of cost efficiency and simplicity, we have chosen the matching option. Instead of using a multi-mode profile with a step index at the core center not compatible with the above mentioned EMB, we tune the gradedindex core diameter (CD) and the numerical aperture (NA) of the fiber to minimize the higher order mode (HOM) excitation. The coupled power ratio between HOM and fundamental mode (LP01) has been simulated for CD between 20µm and 80µm and NA between 0.2 and 0.4 using overlap integral techniques [START_REF] Gholami | Physical Modeling of 10GbE Optical Communication Systems[END_REF] (Fig. 1). This ratio appears to be correlated to:

NA CD NA CD C + ⋅ - = 01 . 0 ) , ( (1)
with CD in µm at both 1310 & 1550nm. CD=40µm and NA=0.35 is a good tradeoff between connection losses in multi-mode propagation at 850nm and fundamental mode matching in single-mode propagation at 1310/1550nm. Min EMB beyond 500 MHz.km is achievable despites the large NA with wellknown alpha shaped profiles 6 .

Fiber Characterization

We have characterized three HNOF samples realized with versatile PCVD process. They exhibit C-parameters between -0.07 and -0.05. To achieve large NA, the cladding has been down-doped. Refractive index profiles are reported in Fig. 2. Min EMB has been calculated from the Differential Mode Delay measurements using standardized weight functions rescaled to CD=40µm 7 . At 850nm, min EMB of 250 MHz.km, 660 MHz.km and 745 MHz.km have been obtained for fiber#1, fiber#2 and fiber#3 respectively. It means that higher data bit rate will be obtained with fiber#2 and #3 at 850nm. Fiber attenuation is dominated by Rayleigh scattering, and does not exceed 1.90, 0.38 & 0.23 dB/km at 850nm, 1310nm & 1550nm. At 1550nm, HNOF bending losses are less than 1 dB measured on 20 turns with a bending radius R=5mm (0.05 dB/turn compared to 0.15 dB/turn for G.657.B3 fiber).

Single-mode Propagation Characterization

Single-mode propagation of fiber#1 (100m) excited by a single-mode pigtail is characterized with OFDR [START_REF] Soller | High resolution optical frequency domain reflectometry for characterization of components and assemblies[END_REF][START_REF] Robillart | Spectral and Temporal Phase Measurement by Optical Frequency-Domain Reflectometry[END_REF] and scan mode [START_REF] Nguyen | Modal Decomposition Technique for Multi-mode Fibers[END_REF] techniques at 1310nm and 1550nm. An OM3 fiber is characterized with the same techniques as the fiber reference.

For OFDR measurements, a time-frequency analysis has been developed. It consists in representing the different interferograms recorded by varying the central wavelength of the OFDR probe (Fig. 3) with a controlled temporal resolution. This permits to measure the temporal responses of the fiber, which depends on the different modes actually excited, over a wide spectrum. The diagrams obtained with 2.3fs/m resolution, reported in the insets of Fig. 3, show bright lines that give the fiber group delays of the excited modes as a function of the wavelength. The three lines seen with the OM3 fiber correspond to three excited mode groups, This behavior is confirmed by the scan mode 10 measurements (Fig. 4) performed at 1520nm with 3.9fs/m resolution (HNOF) and 30fs/m (OM3) that allows to identify the excited spatial modes. The measured coupled power ratio between HOM and the fundamental mode is -30dB with HNOF and -18.6dB with OM3 fiber in agreement with modeling (cf. Fig. 1). The LP01 mode field diameters (MFD) calculated from the scan mode measurements confirm that that of fiber #1 is closer (4% less than SMF MFD) to that of the SSMF, than that of the OM3 fiber (15% more than SMF MFD).

Additional experimental analysis has permitted to check that bending does not affect single-mode propagation.

System Assessment

System assessment is performed on 100m fiber length (fiber #1-3 and bend insensitive OM3) at 1550 nm in single-mode regime. The signal issued from a small form-factor pluggable (SFP) chirped managed laser driven by a 2 31 -1 pseudo random binary sequence (PRBS) at 10Gb/s. Bit error rate (BER) measurements are performed as a function of the receiver input power (Fig. 5). Error free propagation is observed for all fibers and power penalty between 0 dB (fiber#2) and 0.3 dB (fiber#1) are measured at BER=10 -9 compared to the back to back configuration (reference). Under bending constraints (10 turns, R=1cm), negligible impact on performance has been measured on fibers #1 and #2 whereas 0.4 dB additional penalty are observed for fiber #3 and OM3. Similar measurements have been performed in multimode transmission with a 10 Gb/s VCSEL SFP at 850 nm without any error for fibers #2 and #3, with a power penalty at 10 -9 respectively of 2.25 dB and 2.75 dB. No impact is observed under bending constraints. Note that fibre #1 was assessed only at 7 Gb/s in multi mode regime due to the its lower min EMB. In our lab, no modal noise power penalties have been observed at 1550nm with the OM3 fiber despites the excitation of several modes as previously explained. This may be different in house and would have to be re-assessed in this environment since modal noise is expected to be sensitive to temperature fluctuations and mechanical stress that are generally not controlled at home.

Conclusions

We have demonstrated the feasibility of a versatile graded-Index multi-mode fiber for high capacity single and multi-mode optical home network. Among its main advantages, it can be either used with single-mode pigtailed transceivers at 1310/1550nm or multi-mode pigtailed transceivers at 850nm. To obtain such a result, a graded-index multi-mode fiber has been optimized so that only the fundamental mode is excited when coupled with a standard single-mode fiber at 1310/1550nm while offering large effective modal bandwidth at 850nm under VCSEL excitations. 10Gb/s multi-mode and single-mode error-free transmissions over 100m have been demonstrated.

Fig. 1 :

 1 Fig. 1: HOM optical power over fundamental mode optical power in HNOF excited by a SMF

Fig. 2 :Fig. 4 :Fig. 3 :

 243 Fig. 2: Refractive index profile of the realized fiber #1, fiber #2 and fiber #3 in comparison with target profile

Fig. 5 :

 5 Fig. 5: BER measurements with 10 Gb/s NRZ signal on 100 m HNOF fibers and impact of bending constraints (R=1cm, 10 turns) : (a) single-mode, (b) multi-mode

Acknowledgements

This work is funded by the FUI12 RLDO project.


Nathanaël Schaeffer

Efficient Spherical Harmonic Transforms aimed at pseudo-spectral numerical simulations

In this paper, we report on very efficient algorithms for the spherical harmonic transform (SHT). Explicitly vectorized variations of the algorithm based on the Gauss-Legendre quadrature are discussed and implemented in the SHTns library which includes scalar and vector transforms. The main breakthrough is to achieve very efficient on-the-fly computations of the Legendre associated functions, even for very high resolutions, by taking advantage of the specific properties of the SHT and the advanced capabilities of current and future computers. This allows us to simultaneously and significantly reduce memory usage and computation time of the SHT. We measure the performance and accuracy of our algorithms. Even though the complexity of the algorithms implemented in SHTns are in O(N 3 ) (where N is the maximum harmonic degree of the transform), they perform much better than any third party implementation, including lower complexity algorithms, even for truncations as high as N = 1023.

Introduction

Spherical harmonics are the eigenfunctions of the Laplace operator on the 2-sphere. They form a basis and are useful and convenient to describe data on a sphere in a consistent way in spectral space. Spherical Harmonic Transforms (SHT) are the spherical counterpart of the Fourier transform, casting spatial data to the spectral domain and vice versa. They are commonly used in various pseudo-spectral direct numerical simulations in spherical geometry, for simulating the Sun or the liquid core of the Earth among others [START_REF] Glatzmaier | Numerical simulations of stellar convective dynamos. i. the model and method[END_REF][START_REF] Sakuraba | Effect of the inner core on the numerical solution of the magnetohydrodynamic dynamo[END_REF][START_REF] Christensen | A numerical dynamo benchmark[END_REF][START_REF] Brun | Large scale flows in the solar convection zone[END_REF][START_REF] Wicht | Theory and modeling of planetary dynamos[END_REF].

All numerical simulations that take advantage of spherical harmonics use the classical Gauss-Legendre algorithm (see section 2) with complexity O(N 3 ) for a truncation at spherical harmonic degree N . As a consequence of this high computational cost when N increases, high resolution spherical codes currently spend most of their time performing SHT. A few years ago, state of the art numerical simulations used N = 255 [START_REF] Sakuraba | Generation of a strong magnetic field using uniform heat flux at the surface of the core[END_REF].

However, there exist several asymptotically fast algorithms [START_REF] Driscoll | Computing fourier transforms and convolutions on the 2sphere[END_REF][START_REF] Potts | Fast algorithms for discrete polynomial transforms[END_REF][START_REF] Mohlenkamp | A fast transform for spherical harmonics[END_REF][START_REF] Suda | A fast spherical harmonics transform algorithm[END_REF][START_REF] Healy | Ffts for the 2-sphereimprovements and variations[END_REF][START_REF] Tygert | Fast algorithms for spherical harmonic expansions[END_REF], but the overhead for these fast algorithms is such that they do not claim to be effectively faster for N < 512. In addition, some of them lack stability (the error becomes too large even for moderate N ) and flexibility (e.g. N + 1 must be a power of 2).

Among the asymptotically fast algorithms, only two have open-source implementations, and the only one which seems to perform reasonably well is SpharmonicKit, based on the algorithms described by [START_REF] Healy | Ffts for the 2-sphereimprovements and variations[END_REF]. Its main drawback is the need of a latitudinal grid of size 2(N +1) while the Gauss-Legendre quadrature allows the use of only N + 1 collocation points. Thus, even if it were as fast as the Gauss-Legendre approach for the same truncation N , the overall numerical simulation would be slower because it would operate on twice as many points. These facts explain why the Gauss-Legendre algorithm is still the most efficient solution for numerical simulations.

A recent paper [START_REF] Dickson | Importance of explicit vectorization for CPU and GPU software performance[END_REF] reports that carefully tuned software could finally run 9 times faster on the same CPU than the initial non-optimized version, and insists on the importance of vectorization and careful optimization of the code. As the goal of this work is to speed-up numerical simulations, we have written a highly optimized and explicitly vectorized version of the Gauss-Legendre SHT algorithm. The next section recalls the basics of spherical harmonic transforms. We then describe the optimizations we use and we compare the performance of our transform to other SHT implementations. We conclude this paper by a short summary and perspectives for future developments.

2 Spherical Harmonic Transform (SHT)

Definitions and properties

The orthonormalized spherical harmonics of degree n and order -n ≤ m ≤ n are functions defined on the sphere as:

Y m n (θ, φ) = P m n (cos θ) exp(imφ) (1)
where θ is the colatitude, φ is the longitude and P m n are the associated Legendre polynomials normalized for spherical harmonics

P m n (x) = (-1) m 2n + 1 4π (n -|m|)! (n + |m|)! (1 -x 2 ) |m|/2 d |m| dx |m| P n (x) (2)
which involve derivatives of Legendre Polynomials P n (x) defined by the following recurrence:

P 0 (x) = 1 P 1 (x) = x nP n (x) = (2n -1) x P n-1 (x) -(n -1)P n-2 (x)
The spherical harmonics Y m n (θ, φ) form an orthonormal basis for functions defined on the sphere:

2π 0 π 0 Y m n (θ, φ)Y k l (θ, φ) sin θ dθ dφ = δ nl δ mk (3)
with δ ij the Kronecker symbol, and z the complex conjugate of z. By construction, they are eigenfunctions of the Laplace operator on the unit sphere:

∆Y m n = -n(n + 1)Y m n (4)
This property is very appealing for solving many physical problems in spherical geometry involving the Laplace operator.

Synthesis or inverse transform

The Spherical Harmonic synthesis is the evaluation of the sum

f (θ, φ) = N n=0 n m=-n f m n Y m n (θ, φ) (5)
up to degree n = N , given the complex coefficients

f m n . If f (θ, φ) is a real-valued function, f -m n = f m n .
The sums can be exchanged, and using the expression of Y m n we can write

f (θ, φ) = N m=-N   N n=|m| f m n P m n (cos θ)   e imφ (6)
From this last expression, it appears that the summation over m is a regular Fourier Transform.

Hence the remaining task is to evaluate

f m (θ) = N n=|m| f m n P m n (cos θ) (7)
or its discrete version at given collocation points θ j .

Analysis or forward transform

The analysis step of the SHT consists in computing the coefficients

f m n = 2π 0 π 0 f (θ, φ)Y m n (θ, φ) sin θ dθ dφ (8)
The integral over φ is obtained using the Fourier Transform:

f m (θ) = 2π 0 f (θ, φ)e -imφ dφ (9)
so the remaining Legendre transform reads

f m n = π 0 f m (θ)P m n (cos θ) sin θ dθ (10)
The discrete problem reduces to the appropriate quadrature rule to evaluate the integral (10) knowing only the values f m (θ j ). In particular, the use of the Gauss-Legendre quadrature replaces the integral of expression 10 by the sum

f m n = N θ j=1 f m (θ j )P m n (cos θ j )w j ( 11
)
where θ j and w j are respectively the Gauss nodes and weights [START_REF] Temme | Gauss quadrature[END_REF]. Note that the sum equals the integral if f m (θ)P m n (cos θ) is a polynomial in cos θ of order 2N θ -1 or less. If f m (θ) is given by expression 7, then f m (θ)P m n (cos θ) is always a polynomial in cos θ, of degree at most 2N . Hence the Gauss-Legendre quadrature is exact for N θ ≥ N + 1.

A discrete spherical harmonic transform using Gauss nodes as latitudinal grid points and a Gauss-Legendre quadrature for the analysis step is referred to as a Gauss-Legendre algorithm.

3 Optimization of the Gauss-Legendre algorithm

Standard optimizations

Let us first recall some standard optimizations found in almost every serious implementation of the Gauss-Legendre algorithm. All the following optimizations are used in the SHTns library.

Use the Fast-Fourier Transform The expressions of section 2 show that part of the SHT is in fact a Fourier transform. The fast Fourier transform (FFT) should be used for this part, as it improves accuracy and speed. SHTns uses the FFTW library [START_REF] Frigo | The design and implementation of FFTW3[END_REF], a portable, flexible and highly efficient FFT implementation.

Take advantage of Hermitian symmetry for real data When dealing with real-valued data, the spectral coefficients fulfill f -m n = f m n , so we only need to store them for m ≥ 0. This also allows the use of faster real-valued FFTs. Take advantage of mirror symmetry Due to the defined symmetry of spherical harmonics with respect to a reflection about the equator

P m n (cos(π -θ)) = (-1) n+m P m n (cos θ)
one can reduce by a factor of 2 the operation count of both forward and inverse transforms.

Precompute values of P m n

The coefficients P m n (cos θ j ) appear in both synthesis and analysis expressions (7 and 10), and can be precomputed and stored for all (n,m,j). When performing multiple transforms, it avoids computing the Legendre polynomial recursion at every transform and saves some computing power, at the expense of memory bandwidth. This may or may not be efficient, as we will discuss later.

Polar optimization High order spherical harmonics have their magnitude decrease exponentially when approaching the poles as shown in Figure 1. Hence, the integral of expression 10 can be reduced to

f m n = π-θ mn 0 θ mn 0 f m (θ)P m n (cos θ) sin θ dθ (12)
where θ mn 0 ≥ 0 is a threshold below which P m n is considered to be zero. Similarly, the synthesis of f m (θ) (eq. 7) is only needed for θ mn 0 ≤ θ ≤ π -θ mn 0 . SHTns uses a threshold θ mn 0 that does not depend on n, which leads to around 5% to 20% speed increase, depending on the desired accuracy and the truncation N .

On-the-fly algorithms and vectorization

It can be shown that P m n (x) can be computed recursively by

P m m (x) = a m m 1 -x 2 |m|/2 (13) P m m+1 (x) = a m m+1 xP m m (x) (14) P m n (x) = a m n xP m n-1 (x) + b m n P m n-2 (x) (15)
with

a m m = 1 4π |m| k=1 2k + 1 2k (16) a m n = 4n 2 -1 n 2 -m 2 (17) b m n = - 2n + 1 2n -3 (n -1) 2 -m 2 n 2 -m 2 (18)
The coefficients a m n and b m n do not depend on x, and can be easily precomputed and stored into an array of (N + 1) 2 values. This has to be compared to the order N 3 values of P m n (x j ), which are usually precomputed and stored in the spherical harmonic transforms implemented in numerical simulations. The amount of memory required to store all P m n (x j ) in double-precision is at least 2(N + 1) 3 bytes, which gives 2Gb for N = 1023. Our on-the-fly algorithm only needs about 8(N + 1) 2 bytes of storage (same size as a spectral representation f m n ), that is 8Mb for N = 1023. When N becomes very large, it is no longer possible to store P m n (x j ) in memory (for N 1024 nowadays) and on-the-fly algorithms (which recompute P m n (x j ) from the recurrence relation when needed) are then the only possibility.

We would like to stress that even far from that storage limit, on-the-fly algorithm can be significantly faster thanks to vector capabilities of modern processors. Most desktop and laptop computers, as well as many high performance computing clusters, have support for Single-Instruction-Multiple-Data (SIMD) operations in double precision. The SSE2 instruction set is available since year 2000 and currently supported by almost every PC, allowing to perform the same double precision arithmetic operations on a vector of 2 double precision numbers, effectively doubling the computing power. The recently introduced AVX instruction set increases the vector size to 4 double precision numbers. This means that P m n (x) can be computed from the recursion relation 15 (which requires 3 multiplications and 1 addition) for 2 or 4 values of x simultaneously, which may be faster than loading pre-computed values from memory. Hence, as already pointed out by [START_REF] Dickson | Importance of explicit vectorization for CPU and GPU software performance[END_REF], it is therefore very important to use the vector capabilities of modern processors to address their full computing power. Furthermore, when running multiple transforms on the different cores of a computer, the performance of on-the-fly transforms (which use less memory bandwidth) scales much better than algorithms with precomputed matrices, because the memory bandwidth is shared between cores. Superscalar architectures that do not have double-precision SIMD instructions but have many computation units per core (like the POWER7 or SPARC64) could also benefit from on-the-fly transforms by saturating the many computation units with independent computations (at different x).

Figure 2 shows the benefit of explicit vectorization of on-the-fly algorithms on an intel Xeon E5-2680 (Sandy Bridge architecture with AVX instruction set running at 2.7GHz) and compares on-the-fly algorithms with algorithms based on precomputed matrices. With the 4-vectors of AVX, the fastest algorithm is always on-the-fly, while for 2-vectors, the fastest algorithm uses precomputed matrices for N 200. In the forthcoming years, wider vector architecture are expected to become widely available, and the benefits of on-the-fly vectorized transforms will become even more important.

Runtime tuning

We have now two different available algorithms: one uses precomputed values for P m n (x) and the other one computes them on-the-fly at each transform. The SHTns library compares the time taken by those algorithms (and variants) at startup and chooses the fastest, similarly to what the FFTW library [START_REF] Frigo | The design and implementation of FFTW3[END_REF] does. The time overhead required by runtime tuning can be several order of magnitude larger than that of a single transform. The observed performance gain varies between 10 and 30%. This is significant for numerical simulations, but runtime tuning can be entirely skipped for applications performing only a few transforms, in which case there is no noticeable overhead. Figure 2: Efficiency (N + 1) 3 /(2tf ) of various algorithms, where t is the execution time and f the frequency of the Xeon E5-2680 CPU (2.7GHz). On-the-fly algorithms with two different vector sizes are compared with the algorithm using precomputed matrices. Note the influence of hardware vector size for on-the-fly algorithms (AVX vectors pack 4 double precision floating point numbers where SSE3 vectors pack only 2). The efficiency of the algorithm based on precomputed matrices drops above N = 127 probably due to cache size limitations.

Multi-threaded transform

Modern computers have several computing cores. We use OpenMP to implement a multi-threaded algorithm for the Legendre transform including the above optimizations and the on-the-fly approach. The lower memory bandwidth requirements for the on-the-fly approach is an asset for a multi-threaded transform because if each thread would read a different portion of a large matrix, it can saturate the memory bus very quickly. The multi-threaded Fourier transform is left to the FFTW library.

We need to decide how to share the work between different threads. Because we compute the P m n on the fly using the recurrence relation 15, we are left with each thread computing different θ, or different m. As the analysis step involve a sum over θ, we choose the latter option.

From equation 7, we see that the number of terms involved in the sum depends on m, so that the computing cost will also depend on m. In order to achieve the best workload balance between a team of p threads, the thread number i (0 ≤ i < p) handles m = i + kp ≤ N , with integer k from 0 to (N + 1)/p.

For different thread number p, we have measured the time T s (p) and T a (p) needed for a scalar spherical harmonic synthesis and analysis respectively (including the FFT).

Figure 3 shows the speedup T (1)/T (p), where T (p) is the largest of T s (p) and T a (p), and T (1) is the time of the fastest single threaded tranform. It shows that there is no point in doing a parallel transform with N below 128. The speedup is good for N = 255 or above, and excellent up to 8 threads for N ≥ 511 or up to 16 threads for very large transform (N ≥ 2047).

Performance comparisons

Table 1 reports the timing measurements of two SHT libraries, compared to the optimized Gauss-Legendre implementation found in the SHTns library (this work). We compare with the Gauss-Legendre implementation of libpsht [START_REF] Reinecke | Libpsht -algorithms for efficient spherical harmonic transforms[END_REF], a parallel spherical harmonic transform library targeting very large N , and with SpharmonicKit 2.7 (DH) which implements one of the Driscoll-Healy fast algorithms [START_REF] Healy | Ffts for the 2-sphereimprovements and variations[END_REF]. All the timings are for a complete SHT, which includes the Fast Fourier Transform. Note that the Gauss-Legendre algorithm is by far (a factor of order 2) the fastest algorithm of the libpsht library. Note also that SpharmonicKit is limited to N + 1 being a power of two, requires 2(N + 1) latitudinal colocation points, and crashed for N = 2047. The software library implementing the fast Legendre transform described by [START_REF] Mohlenkamp | A fast transform for spherical harmonics[END_REF], libftsh, has also been tested, and found to be of comparable performance to that of SpharmonicKit, although the comparison is not straightforward because libftsh did not include the Fourier Transform. Again, that fast library could not operate at N = 2047 because of memory limitations. Note finally that these measurements were performed on a machine that did not support the new AVX instruction set.

In order to ease the comparison, we define the efficiency of the SHT by (N + 1) 3 /(2T f ), where T is the execution time (reported in Table 1) and f the frequency of the CPU. Note that (N +1) 3 /2 reflects the number of computation elements of a Gauss-Legendre algorithm (the number of modes (N +1)(N +2)/2 times the number of latitudinal points N +1). An efficiency that does not depend on N corresponds to an algorithm with an execution time proportional to N 3 .

The efficiency of the tested algorithms are displayed in Figure 4. Not surprisingly, the Driscoll-Healy implementation has the largest slope, which means that its efficiency grows fastest with N , as expected for a fast algorithm. It also performs slightly better than libpsht for N ≥ 511. However, even for N = 1023 (the largest size that it can compute), it is still 2.8 times slower than the Gauss-Legendre algorithm implemented in SHTns. It is remarkable that SHTns achieves an efficiency very close to 1, meaning that almost one element per clock cycle is computed for N =≥ 511. Overall, SHTns is between two and ten times faster than the best alternative.

Accuracy

One cannot write about an SHT implementation without addressing its accuracy. The Gauss-Legendre quadrature ensures very good accuracy, at least on par with other high quality implementations.

The recurrence relation we use (see §3.2) is numerically stable, but for N 1500, the value P m m (x) can become so small that it cannot be represented by a double precision number anymore. To avoid this underflow problem, the code dynamically rescales the values of P m n (x) during the recursion, when they reach a given threshold. The number of rescalings is stored in an integer, which acts as an enhanced exponent. Our implementation of the rescaling does not impact performance negatively, as it is compensated by dynamic polar optimization: these very small values are treated as zero in the transform (eq. 7 and 11), but not in the recurrence. This technique ensures good accuracy up to N = 8191 at least, but partial transforms have been performed successfully up to N = 43600.

To quantify the error we start with random spherical harmonic coefficients Q m n with each real part and imaginary part between -1 and +1. After a backward and forward transform (with orthonormal spherical harmonics), we compare the resulting coefficients R m n with the originals Q m n . We use two different error measurements: the maximum error is defined as

ǫ max = max n,m |R m n -Q m n |
while the root mean square (rms) error is defined as

ǫ rms = 2 (N + 1)(N + 2) n,m |R m n -Q m n | 2
The error measurements for our on-the-fly Gauss-Legendre implementation with the default polar optimization and for various truncation degrees N are shown in Figure 5. The errors steadily increase with N and are comparable to other implementations. For N < 2048 we have ǫ max < 10 -11 , which is negligible compared to other sources of errors in most numerical simulations.

Conclusion and perspectives

Despite the many fast spherical harmonic transform algorithms published, the few with a publicly available implementation are far from the performance of a carefully written Gauss-Legendre algorithm, as implemented in the SHTns library, even for quite large truncation (N = 1023). Explicitly vectorized on-the-fly algorithms seem to be able to unleash the computing power of nowadays and future computers, without suffering too much of memory bandwidth limitations, which is an asset for multi-threaded transforms.

The SHTns library has already been used in various demanding computations (eg. [START_REF] Schaeffer | On the reflection of alfvén waves and its implication for earth's core modelling[END_REF][START_REF] Augier | A new formulation of the spectral energy budget of the atmosphere, with application to two high-resolution general circulation models[END_REF][START_REF] Figueroa | Modes and instabilities in magnetized spherical couette flow[END_REF]. The versatile truncation, the various normalization conventions supported, as well as the scalar and vector transform routines available for C/C++, Fortran or Python, should suit most of the current and future needs in high performance computing involving partial differential equations in spherical geometry.

Thanks to the significant performance gain, as well as the much lower memory requirement of vectorized on-the-fly implementations, we should be able to run spectral geodynamo simulations at N = 1023 in the next few years. Such high resolution simulations will operate in a regime much closer to the dynamics of the Earth's core.

Figure 1 :

 1 Figure 1: Two associated Legendre polynomials of degree n = 40 and order m = 33 (blue) and m = 36 (red), showing the localization near the equator.

Figure 3 :

 3 Figure 3: Speedup obtained with multiple threads using OpenMP (gcc 4.6.3) on a 16 core intel Xeon E5-2680 (Sandy Bridge architecture with AVX instruction set running at 2.7 GHz). N 63 127 255 511 1023 2047 4095 libpsht (1 thread) 1.05 ms 4.7 ms 27 ms 162 ms 850 ms 4.4 s 30.5 s DH (fast) 1.1 ms 5.5 ms 21 ms 110 ms 600 ms NA NA SHTns (1 thread) 0.09 ms 0.60 ms 4.2 ms 28 ms 216 ms 1.6 s 11.8 s Table 1: Comparison of execution time for different SHT implementations. The numbers correspond to the average execution time for forward and backward scalar transform (including the FFT) on an Intel Xeon X5650 (2.67GHz) with 12 cores. The programs were compiled with gcc 4.4.5 and -O3 -march=native -ffast-math compilation options.

Figure 4 :

 4 Figure 4: Efficiency (N +1) 3 /(2T f ) of the implementations from Table1, where T is the execution time and f the frequency of the Xeon X5650 CPU (2.67GHz) with 12 cores.

Figure 5 :

 5 Figure 5: Accuracy of the on-the-fly Gauss-Legendre algorithm with the default polar optimization.

  Table 1, where T is the execution time and f the frequency of the Xeon X5650 CPU (2.67GHz) with 12 cores.

Acknowledgements

The author thanks Alexandre Fournier and Daniel Lemire for their comments that helped to improve the paper. Some computations have been carried out at the Service Commun de Calcul Intensif de l'Observatoire de Grenoble (SCCI) and other were run on the PRACE Research Infrastructure Curie at the TGCC (grant PA1039).


Julien Lavalle
email: lavalle@in2p3.fr

Stefano Magni
email: stefano.magni@univ-montp2.fr

Making sense of the local Galactic escape speed estimates in direct dark matter searches

may come

Making sense of the local Galactic escape speed

estimates in direct dark matter searches

I. INTRODUCTION

Direct detection (DD) of dark matter (DM) [1,2] has reached an impressive sensitivity thanks to developments in background rejection, increase in the target masses, and also boosted by the advent of dual-phase xenon detectors [3][4][5][6][7][8][9]. Weakly interacting massive particle (WIMP) candidates [10][11][12][13] are the main focus of such searches. In the last decade and half, important theoretical and experimental efforts have been invested in the inspection of the signal reported by the DAMA and DAMA/LIBRA Collaborations [14,15], an annual modulation of the event rate a priori consistent with a WIMP signal [16]. Although some experiments have reported signal-like events in the light WIMP mass region ∼ 10 GeV that might be consistent with the DAMA signal [7,[17][18][19], while not statistically significant, other negative results [9,20,21] tend to rule out the simplest interpretation in terms of spin-independent (SI) scatterings of WIMPs on nuclei. It also turns out that this WIMP mass range is not favored by indirect searches in the antiproton channel if the relic abundance is set by an s-wave annihilation to quarks [START_REF] Lavalle | Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics[END_REF][START_REF] Cerdeno | [END_REF]. Nevertheless, although the signal-like features show up very close to the experimental thresholds, and while there are some attempts to provide explanations of the DAMA signal in terms of backgrounds [24][25][26], a careful investigation of the phenomenology of light WIMPs is still worth. We also note that this mass range has also serious theoretical motivations from the model-building point of view (e.g. [27][28][29]).

The experimental sensitivity to light WIMPs strongly depends on both instrumental and astrophysical effects. The former is mostly related to the energy threshold and resolution of the detector, for given target nuclei and assuming perfect background rejection. The latter is related to the tail of the WIMP speed distribution function (DF) in the observer's rest frame, as a detection threshold energy converts into larger speed thresholds for lighter WIMPs. The (annual average of the) maximum speed a WIMP can reach in the laboratory is the sum of the escape speed v esc , i.e. the speed a test particle would need to escape the Galactic gravitational potential at the Sun's location, and of the circular velocity of the Sun. These are the two main astrophysical parameters relevant to assess the experimental sensitivity in the low WIMP mass regime. The latter is continuously investigated, and many studies have been revisiting and improving its estimation since the recommendation of 220 km/s by Kerr & Lynden-Bell [30]. As for v esc , which should be a real physical cut-off in the velocity DF of the gravitationally bound WIMPs as measured in the dark halo rest frame, constraints are more scarce. Better understanding and better estimating the escape speed should therefore strongly benefit direct searches for light WIMP candidates.

The local Galactic escape velocity can be estimated with different approaches. One can go through a global halo modeling from photometric and kinematic data (e.g. [31,32]) and compute the gravitational potential from the resulting mass model. This is an indirect method that mostly relies on rotation curve fitting procedures. A more direct method consists in trying to use those nearby high-velocity stars which are supposed to trace the tail of the corresponding stellar phase-space DF, which, though not necessarily equal to the DM DF, should still also vanish at the escape speed. A few studies have pushed in this direction since the seminal work by Leonard & Tremaine [33]. In particular, two results have been published based on the data of the RAVE survey, a massive spectroscopic survey whose first catalog was released in 2006 [34]. One of these results was based upon the first data release [35] (hereafter S07), and provided an estimate of v esc = 544 +64 -46 at 90% confidence level (CL). This value is now used in the so-called standard/simplest halo model (SHM), a proxy for the local DM phase-space conventionally employed to derive DD exclusion curves (see e.g. [36] for associated uncertainties). The latest RAVE analysis, which relies on the fourth data release and different assumptions, was published in 2014 [37] (hereafter P14) and found v esc = 533 +54 -41 at 90% CL, consistent with the previous result while slightly reducing the central value and the statistical errors. Among improvements with respect to S07, P14 used much better distance and velocity reconstructions, and a more robust star selection. This means that systematic errors should have decreased significantly. We will focus on this result in the following.

Whereas it is straightforward to derive the impact of this new escape speed range on the exclusion curves or signal contours, or to include it as independent constraints or priors in Bayesian analyses, it is worth stressing that this would lead to inconsistent conclusions in most cases, as, in particular, the range itself relies on a series of assumptions. In this paper, we will examine these assumptions and derive self-consistent local DM phase-space models, accounting for correlations in the astrophysical parameters. We will see that these correlations translate into DD uncertainties quite different from those that would be obtained by using the 90% CL range blindly. Although the way astrophysical parameters affect direct detection searches is rather well understood (e.g. [38][39][40]), the related uncertainties are still being refined.

There have been several studies mostly based on rotation curves which tried to bracket these uncertainties (e.g. [41][42][43][44]). Here, we provide a complementary view focusing on observational constraints that rely on different observables, namely non co-rotating high-velocity stars, with emphasis on the escape speed. This is particularly relevant to light WIMP searches. The outline of this paper is the following. We will first review the results obtained in [37] and derive a local DM phase-space DF consistent with their assumptions. Then, thanks to this phase-space modeling, we will convert their 90 and 99% CL ranges into uncertainties in DD exclusion curves. We will eventually discuss these results in light of other complementary and independent astrophysical constraints before concluding.

II. LOCAL DM PHASE-SPACE FROM RAVE-P14 ASSUMPTIONS AND RESULTS

A. Brief description of RAVE-P14 and main parameter assumptions

The RAVE-P14 analysis is based on a sample of stars from the RAVE catalog which are not co-rotating with the Galactic disk (34 or 69 stars after selection cuts, mostly hard -> 300 km/s-or weak -> 200 km/svelocity cut, respectively). It is complemented with data from another catalog [45], which were used to derive the best-fit range for the escape speed recalled above (19 or 17 stars). Therefore, the full sample contains highvelocity halo stars which are assumed to probe the highvelocity tail of the phase-space DF. We will only consider results obtained with hard velocity cut in the following.

Assuming steady state and that the phase-space DFs of DM and non co-rotating stars mostly depends on energy, one may introduce the concept of escape speed v esc consistently with Jeans theorem, by stating that the phasespace DF must vanish for speeds greater than v esc ( x) = -2 φ( x), where φ( x) is the gravitational potential at position x (e.g. [START_REF] Binney | Galactic Dynamics: Second Edition[END_REF]). Then, observational constraints on v esc from stellar velocities can only be derived if the true shape of the underlying velocity DF is known. P14 used the ansatz proposed in [33] for the high-velocity tail of the DF of stars, which reads

f (v) ∝ (v esc -v) k θ(v esc -v),
where v = | v| and k is a free index calibrated from galaxy simulations. This ansatz is general enough and holds provided v is close to v esc . It differs from the one used in S07,

f (v) ∝ (v 2 esc -v 2 ) k θ(v esc -v).
Actually, P14 tested both ansatze with cosmological simulations performed in [START_REF] Scannapieco | [END_REF], which include baryons and star formation, and found the former to be more consistent with the simulation results -note that S07 also employed cosmological simulations, though less recent, to calibrate their reconstruction method. This is a significant systematic difference between S07 and P14, though the quantitative impact is not spectacular: when performing their likelihood analysis from S07 ansatz, P14 find v esc = 537 +37 -33 km/s instead of their nominal result v esc = 533 +54 -41 (90% CL). Another systematic difference between S07 and P14 likelihood analyses, which allows the latter to increase statistics, is that P14 applied a correction to the selected star velocities to "relocate" them at the radial position of the Sun r before running their likelihood. Given a line-of-sight velocity component v || ( x) for a star located at position x, the correction reads

v || ( x ) = v || ( x) × |φ( x )/φ( x)|,
where x is the position of the Sun and φ is the total gravitational potential. These corrected speeds are supposed to describe the local DF more reliably and are expected to reduce the systematic uncertainties in the determination of v esc , while increasing the statistics.

S07 actually dealt with that issue by selecting stars within a small radial range centered around r , but from more recent distance estimates, it turns out that their sample is likely biased (see discussion in P14). We may remark that despite the improvements in the P14 methodology with respect to S07, correcting the star velocities by the gravitational potential automatically introduces a Galactic mass model dependence. Even though the gravitational potential is not expected to vary much for nearby stars, this aspect is important when one wants to use P14 results in the context of DD. We will show how to cope with that in the next (sub)sections.

Before going into the details of the Milky Way (MW) mass modeling, it is worth recalling some key parameters that are used in P14. First, r is fixed to the central value found in [48], r = 8.28 kpc .

(

) 1
For the circular velocity of the Sun, v c , three cases were considered, each associated with different results for v esc : (i) 220 km/s, (ii) 240 km/s, (iii) free v c -we will discuss these cases into more details later on. Finally, the peculiar velocity of the Sun was fixed to the vector elements found in [49]:

U = 11.1 km/s
(2) V = 12.24 km/s W = 7.25 km/s .

B. Milky-Way mass model

Since P14 likelihood relies on star "relocation", a procedure in which measured stellar speeds are corrected by the gravitational potential, the escape speed estimate must be used consistently with the Galactic mass model they assumed. The latter is made of three components: a dark matter halo, a baryonic bulge, and a baryonic disk.

We note that the baryonic content is fixed, while the DM halo parameters are left free for the P14 likelihood analysis. There is therefore a direct correlation between the dark halo parameters and the escape speed as a result of P14 analysis.

The baryonic bulge is described with a spherical Hernquist density profile [50] as follows:

ρ b (r) = M b 2 π r b /r (r + r b ) 3 , (3)
where r denotes the Galactocentric radius, G Newton's constant, r b a scale radius, and M b the total bulge mass. The parameters values are given in Table I. To this profile is associated a gravitational potential:

φ b (r) = -G M b (r+r b ) . ( 4
)
The disk is modeled from an axisymmetric Miyamoto-Nagai profile [51]:

ρ d (R, |z|) = z 2 d M d 4 π      R d R 2 + (R d + 3 z 2 + z 2 d )(R d + z 2 + z 2 d ) 2 R 2 + (R d + z 2 + z 2 d ) 2 5/2 (z 2 + z 2 d ) 3/2      , (5)
where M d is the total disk mass, and R d and z d are radial and vertical scale lengths.

The parameters values are given in Table I. This profile converts into the following potential:

φ d (R, |z|) = -G M d R 2 + (R d + z 2 + z 2 d ) 2 . (6)
Finally, the DM halo was modeled from a spherical NFW profile [52]:

ρ(r) = ρ s x -1 (1 + x) -2 , (7)
where x = r/r s , and ρ s and r s are a scale density and a scale radius, respectively, which are left as free parameters in P14 1 . The associated gravitational potential 1 Note that if a prior on vc is considered, then only one DM halo reads:

φ dm (r) = -4 π G ρ s r 3 s r ln 1 + r r s (8) -→ r→0 -4 π G ρ s r 2 s ,
which is minimal and finite at the Galactic center as shown by the limit above. Note that P14 also considered an adiabatically contracted NFW, but we will not discuss this case here. Actually, P14 mostly used the concentration parameter c and the total Milky Way mass instead of r s and ρ s , parameter remains free. which is strictly equivalent. The former is defined as:

c 340 = R 340 r s , (9)
where R 340 is the radius that encompasses the DM halo such that the average DM density is 340 times the critical density ρ c (taking a Hubble parameter of H 0 = 73 km/s/Mpc):

R 340 = 3 M dm (R 340 ) 4 π × 340 ρ c 1/3 . ( 10
)
Here, M dm is the full dark matter mass enclosed within a radius of R 340 . The total Milky Way mass M 340 is defined as the sum of the DM and baryonic components inside R 340 ,

M 340 = M tot (R 340 ) = M dm (R 340 ) + M b + M d . (11)
C. Circular and escape speeds: converting the original RAVE-P14 results

The MW mass model defined in the previous section is used to connect the RAVE data, i.e. the stellar line-ofsight velocities, to the circular velocity v c in the MW disk (z = 0) at the position of the Sun. Indeed, from classical dynamics, the circular velocity v c in the MW disk (z = 0) is related to the radial gradient of the gravitational potential:

v 2 c (R, z = 0) = R ∂ R φ tot (R, z = 0) , (12)
where the expression is given in its axisymmetric form (r 2 = R 2 + z 2 ), and where φ tot is the sum of all contributions to the gravitational potential (baryonic components and DM).

At distance r , one has therefore:

v 2 c (r , z = 0) = r G M b (r + r b ) 2 + r M d (r 2 + R2 d ) 3/2 + 4 π ρ s r s x 2 ln(1 + x ) - x 1 + x , (13)
where Rd = R d + z d .

In P14, three different choices were made for v c : 220 km/s, 240 km/s, and free parameter for the likelihood. In the last case, the best-fit to the data is obtained for v c = 196 km/s, a value significantly smaller than most recent estimates (e.g. [53,54]). Nevertheless, this fit includes a prior on the concentration parameter based on Ref. [55], which strongly affects the result (the best-fit obtained for the concentration is c 340 = 5). We will discuss this into more details in Sec. III D. We now come to the main observable constrained in P14: the escape speed at the solar location. It is in principle defined from the global gravitational potential at the position of the Sun, as v esc = -2 φ tot (r ). Nevertheless, even though the potential is usually assumed to vanish at infinity, the presence of nearby galaxies such as Andromeda, located at ∼ 790 kpc [56,57], makes it tricky to define the potential difference relevant to the escape. As an educated guess, P14 assumed a boundary

R max = 3 R 340 , (14)
which roughly lies in the range 500-600 kpc. Accordingly, the escape speed is defined from the (positive) potential difference ψ such that

v esc (r ) ≡ 2 ψ tot (r ) , (15)
ψ tot (r ) ≡ -(φ tot (r ) -φ tot (R max )) . ( 16
)
While ψ tot contains all matter components, this definition will also be retained for individual components.

Since the baryonic MW mass model components are fixed, a pair of (v c , v esc ) directly converts into a pair of NFW parameters (ρ s , r s ), or equivalently (M 340 , c 340 ). It is the latest pair of parameters that was used to illustrate the likelihood results in P14, in the case of free v c (see their Figs. 13 and14). We have used Eqs. (13) and (15) to convert P14 likelihood results into the plane v esc -v c , more relevant to derive uncertainties in direct detection, as it will be discussed in Sec. III. This is displayed in Fig. 1. The original frame of the P14 plot (their Fig. 13, left panel) in the plane M 340 -c 340 is reported as the yellowish curves. The 10(1)% of maximum likelihood band is shown as the dark (light) blue area, while the 90 (99)% CL log-normal prior on the concentration parameter appears as a dark (light) orange band, and the concentration model of Ref. [55] is displayed as a dashed red curve. The final contours accounting for the concentration prior is shown as black curves. P14 best-fit points obtained after enforcing v c = 220 (240) km/s (no prior on concentration) are shown as green (magenta) points, with associated 90% CL errors -the former is often referred to as the primary result of P14. We also report the iso-DM-density curves (gray lines) at the so- lar location, ρ = ρ(r ), in units of GeV/cm 3 . Finally, the green band shows additional constraints on the ratio v c /r taken from Ref.

[54] that we will discuss in Sec. III C. We clearly see from Fig. 1 that the P14 results, because of the initial assumptions, induce strong correlations among the local DM parameters relevant to predictions in direct DM searches. Therefore, a blind use of the escape velocity constraints, as very often done (e.g. by taking arbitrary values of local DM density), can already be classified as inconsistent. Fig. 1 will be our primary tool to make a consistent interpretation of P14 results in the frame of direct DM detection.

We may also note that the trend of the blue band is in contrast with the claim in P14 that their v esc estimate anticorrelates with v c as a result of their sample selection (biased to negative longitudes), which explains why v esc decreases from 533 to 511 km/s as v c increases from 220 to 240 km/s -we have illustrated this anticorrelation in Fig. 1 with dotted lines. This actually comes from the method used in P14 to extract the likelihood region in the c 340 -M 340 plane 2 (see their Fig. 13): the authors kept the 2 We warmly thank Tilmann Piffl for clarifying this issue.

posterior PDF of v esc frozen to the shape obtained with v c = 220 km/s, while varying only c 340 and M 340 . This means that they did not recompute the velocities of their stellar sample according to the changes in v c induced by those in c 340 and M 340 . Such an approximation is hard to guess from the text, and even more difficult from their Fig. 13 showing the likelihood band in the c 340 -M 340 plane. Indeed, they also report iso-v c curves on the same plot, so we may expect that the corresponding posterior PDF for v esc was taken accordingly. We remind that in principle the pair v c -v esc is strictly equivalent to the pair c 340 -M 340 , so we may have expected to recover the anticorrelation claimed for the former pair from the contours obtained for the latter. While the goal of P14 in the v cfree case was mostly to investigate how to improve the matching between Galactic models and the primary fit results, with a focus on the Milky Way mass, this somehow breaks the self-consistency of the analysis. This has poor impact on the Galactic mass estimate, which was the main focus in P14, but this affects the true dynamical correlation that v esc should exhibit with the other astrophysical parameters. Unfortunately, improving on this issue would require access to the original data, which is beyond the scope of this paper. Therefore, we will stick to this result in the following, and further comment on potential ways to remedy this limitation at the qualitative level in Sec. III D 2.

To conclude this section, we calculate and provide in Table II the DM halo parameters associated with the P14 best-fit points. We stress that to each value of v esc found in the table corresponds a specific value of the local DM density. This kind of correlations should be taken into account for a proper use of P14 results. It is obvious that the relevance of the values and ranges for these parameters should be questioned in light of complementary constraints. We will discuss this issue in Sec. III D.

III. UNCERTAINTIES IN DIRECT DM SEARCHES

In this section, we sketch a method to make a selfconsistent interpretation of P14 results, and show how they convert into astrophysical uncertainties beyond only uncertainties in the escape speed. We will explicitly show the impact of dynamical correlations between the local circular and escape speeds and the DM parameters.

A. Generalities

We recall here the main equations and assumptions used to make predictions in direct DM searches, focusing on the spin-independent class of signals (reviews can be found in e.g. [11,40,58]). Assuming elastic collisions between WIMPs of mass m χ and atomic targets of atomic mass number A and mass M (and equal effective couplings f n and f p of WIMPs to neutrons and protons), the differential event rate per atomic target mass in some DD experiment reads:

dR dE r = ρ m χ A 2 σ p F 2 (q) 2 µ 2 p η(v min (E r ), t) , (17)
where µ p is the WIMP-proton reduced mass, q is the exchanged momentum, E r = q 2 /(2 M ) is the nuclear recoil energy, σ p = (4/π) µ 2 p f 2 p is the WIMP-proton cross section, F (q) is a nuclear form factor, and η is the truncated mean inverse local WIMP speed beyond a threshold v min :

η(v min (E r ), t) = v>vmin d 3 v f ⊕ ( v, t) v . (18)
Here, v = | v| is the WIMP speed in the observer frame, and f ⊕ is the associated DF (normalized to unity over the full velocity range). The threshold speed

v min (E r , m χ , M ) = M E r /(2 µ 2 ) , ( 19
)
where µ is the WIMP-nucleus reduced mass, is the minimal speed to achieve a recoil energy of E r . A simple Galilean transformation allows to connect the WIMP velocity DF in the local frame f ⊕ to that in the Milky Way frame f :

f ⊕ ( v, t) = f ( v + v ⊕ (t)) . (20)
We will further discuss the form of f in Sec. III B. The time dependence is explicitly shown to arise from the Earth motion around the Sun v ⊕/ (t), since the Earth velocity in the MW frame is given by v ⊕ (t) = v + v ⊕/ (t). In a frame whose first axis points to the Galactic center and the second one in direction of the rotation Galactic flow in the disk, we may parameterize the solar velocity as follows:

v = (U , V + v c , W ) , (21)
where v c is the circular speed and the peculiar components are given in Eq. ( 2). For the Earth motion around the Sun, featured by v ⊕/ (t), we will use the prescription of Ref. [59], whose accuracy has been recently confirmed again in [60].

To compute the event rate above a threshold for a given experiment, or more generally inside an energy bin (or interval) i, one has further to account for the experiment efficiency (E r ) and (normalized) energy resolution G(E r , E r ) such that:

R i (t) = ∞ 0 dE r (E r ) dR(E r ) dE r E max i E min i dE r G(E r , E r ) . (22)
A sum over the fractions of atomic targets must obviously be considered in case of multi-target experiments (or to account for several isotopes). To very good approximation, we can define the average and modulated rates, Ri and Ri as:

Ri = 1 2 {R i (t max ) + R i (t min )} (23) Ri = 1 2 {R i (t max ) -R i (t min )} , (24)
where t max 152 (t min 335) is the day of the year where the rate is maximum (minimum), which can revert in some cases depending on the WIMP mass for a given energy threshold (and vice-versa).

One can already read the trivial scaling in the local density ρ off Eq. ( 17), which is still already quite relevant to interpret P14 results self-consistently. Indeed, we have seen in the previous section that ρ was correlated with the circular speed v c as well as to the escape speed v esc by construction. We investigate the impact of the velocity DF in the next subsection.

B. Toward consistent DM velocity distributions

Here we shortly introduce the usual way to deal with the WIMP velocity DF before specializing to the ergodic assumption. A nice review can be found in Ref. [39]. is conventionally used to compare experimental sensitivities and results. Beside fixing the local DM density ρ , the SHM specifies the WIMP velocity distribution as a truncated Maxwell-Boltzmann (MB) distribution, based on the assumption of an isothermal sphere, and defined in the MW rest frame (also called local standard of rest -LSR) as:

f 1 (v) (km/s) -1 Lavalle & Magni (2014) Maxwell-Boltzmann-like (v c =220 km/s) Ergodic NFW+baryons (v c =220 km/s) Ergodic NFW+bulge (v c <220 km/s) Ergodic NFW only (v c <220
f 1 (v) (km/s) -1 Lavalle & Magni (2014) MB-like (v c =220 km/s) mean MB-like (v c =220 km/s) 90% upper MB-like (v c =220 km/s) 90% lower Ergodic (v c =220 km/s) mean Ergodic (v c =220 km/s) 90% upper Ergodic (v c =220 km/s) 90% lower MB-like (v c =240 km/s) mean MB-like (v c =240 km/s) 90% upper MB-like (v c =240 km/s) 90% lower Ergodic (v c =240 km/s) mean Ergodic (v c =240 km/s) 90% upper Ergodic (v c =240 km/s) 90% lower
f MB ( v) =      N 0 e -v 2 v 2 c -e - v 2 esc v 2 c ∀v ≤ v esc 0 ∀v > v esc (25)
where N 0 allows a normalization to unity over d 3 v, and the right-hand side term enforces the phase-space to vanish at v esc . The SHM is also sometimes defined with a hard velocity cut-off, which corresponds to neglecting the second term of the right-hand side above. As this func-tion is here defined in the LSR, a proper use in DD merely rests upon the Galilean transformation of Eq. ( 20).

The SHM is defined with fixed parameters that we summarize below:

ρ shm = 0.3 GeV/cm 3 ( 26
)
v shm c = 220 km/s v shm esc = 544 km/s .
The escape speed value is actually the central value found in the first RAVE analysis, S07 [35]. It is noteworthy that the velocity dispersion relation

σ 2 MB = 3 v 2 c /2 (27)
predicted from the Jeans equations in the absence of anisotropy for the isothermal sphere is no longer valid as soon as the by-hand truncation at v esc is introduced. While the difference is quantitatively not important, especially for either low circular speed and large escape speed, the consistency of the SHM velocity distribution is already broken at this stage. Anyway, the SHM provides a useful approximation to compare different results, and is particularly convenient as η (see Eq. 18) takes an analytic expression (see e.g. Ref. [36]). The SHM may be (and is actually widely) used to study the impact of astrophysical uncertainties in DD predictions, while keeping in mind that it relies on not fully consistent theoretical grounds. The real dynamical correlations among the different astrophysical parameters can still be implemented by hand in the SHM, but both the shape of the velocity distribution and the velocity dispersion are fixed.

Ergodic velocity distribution

Another approach to build a self-consistent velocity distribution is based on considering functions of integrals of motion, which automatically satisfy the Jeans equations [START_REF] Binney | Galactic Dynamics: Second Edition[END_REF]. The simplest case, that we will consider here, arises for steady-state, non-rotating, and spherical systems where the energy, E(r, v) = m v 2 /2 + m φ(r), is an integral of motion. Then the phase-space distribution can be fully described by positive functions of E: such systems are called ergodic. In that case, it was shown by Eddington [61] that one could relate the mass density profile of the objects population under scrutiny to its phase-space distribution (per phase-space volume and per mass unit) through the so-called Eddington equation:

f ( ) = 1 √ 8 π 2 1 √ dρ dψ ψ=0 + 0 dψ √ -ψ d 2 ρ dψ 2 (28) = 2 √ 8 π 2 1 2 √ dρ dψ ψ=0 + √ d 2 ρ dψ 2 ψ=0 + 0 dψ -ψ d 3 ρ dψ 3 ,
where ρ is, in the case of interest here, the DM density profile and ψ is the full potential difference given in Eq. ( 16) (including DM and baryons) and used to define the escape speed. Note that since both the DM profile and total gravitational potential are monotonous functions of r, each radial slice in ρ can be matched to a slice in ψ, which allows for an unambiguous definition of dρ/dψ. The latest line of the previous equation is simply an integration by parts of the former line, and is more suitable for numerical convergence because of the absence of the 1/ √ -ψ factor which diverges when ψ → . The relative energy (per mass unit) is positive defined and reads

(r, v) ≡ ψ(r) - v 2 2 ( 29
)
such that vanishes at v esc . Since by definition of the phase-space distribution f we have ρ(r) = d 3 v f ( ), we can straightforwardly derive the velocity distribution at radius r (in the MW rest frame):

f erg (v) r = f ( ) ρ(r) . (30)
By construction, this function is positive defined in the range | v| ∈ [0, v esc ], and vanishes outside. This formalism has already been used in the context of DD in several previous studies, for instance in Refs. [62,63] (see also e.g. [42,44,64] for more recent references). In practice, we calculate Eq. ( 28) from a dedicated C++ code. It is important to remark that since the zero of the relative gravitational potential ψ is not met as r → ∞, but instead at r = R max (see Eqs. 14 and 16), the terms of Eq. ( 28) evaluated at ψ = 0 cannot be neglected, as very often done. In the following,

f 1,erg (v) = 4 π v 2 f erg (v) (31)
will denote the DF associated with the velocity modulus v, sufficient to describe the isotropic case.

An important point for this derivation to be consistent is that the system must be spherically symmetric (by construction of the ergodic function, otherwise additional degrees of freedom would be necessary). This is not strictly the case as the gravitational potential associated with the Galactic disk, as given in Eq. ( 6), is axisymmetric. Since the disk does not dominate the gravitational potential, we can still force spherical symmetry while not affecting the circular velocity at the solar position (r , z = 0). We must make sure that the equation

v 2 c (r , z = 0) = R ∂ R φ tot (R, z) R=r z=0 (32)
is verified. As r 2 = R 2 + z 2 , we can safely redefine a spherical potential φtot (r) = φ tot (R = r, z = 0) (33) for which Eq. ( 32) is obeyed. A more refined approach would be to include the angular momentum as a supplementary degree of freedom to allow for anisotropic velocity distributions (see e.g. Refs. [44] and [64] to get some insights), but we leave this potential improvement to further studies. We still note that a few studies have tested ergodic DFs against cosmological simulations, where the angular momentum was shown to play a role essentially far outside the central regions of the dark halo (see e.g.

Ref. [65]). In Fig. 2, we show the differences between the onedimensional truncated MB speed distribution and the ergodic function, in the MW rest frame. For the latter we explicitly plot the impact of the baryonic components of the gravitational potential (neglecting a component leads to lowering v c and v esc ) in the left panel. In this panel, we took the P14 best-fit point v esc = 533 km/s associated with the prior v c = 220 km/s -see Table II. As expected, the peaks and the mean speeds are different, typically with larger values for the ergodic distribution.

In the right panel of Fig. 2 we compare the MB and ergodic DFs for both P14 best-fit points v esc = 533 and 511 km/s (with priors v c = 220 and 240 km/s, respectively), showing also the DFs associated with the 90%CL upper and lower corners. What is interesting to note is that while the peaks of the MB DFs move with v c , the peaks of the ergodic DFs are almost unaffected; only the high-velocity tails are, which merely reflects the changes in v esc .

C. RAVE-P14 results on vesc in terms of direct detection limits

Before converting P14 results on v esc in terms of DD limits, we stress again that not only do they affect the whole WIMP velocity DF (v c , v esc , and the velocity dispersion), but also the local DM density, as already showed in Sec. II C. We first recall the effects of the main parameters on the DD exclusion curves:

• v esc + v c + V : Defines the average WIMP mass threshold m min χ given an atomic target and a recoil energy threshold, by solving v min (E th , m χ , M ) = v esc + v c + V -see Eq. ( 19). This corresponds to the position on the mass axis of the asymptote σ max → ∞ of the upper limit on the SI cross section (the larger v esc and/or v c the lower the mass threshold).

• v c : Impact on the relative position of the maximum sensitivity of a DD experiment (for a given atomic target); a larger v c globally shifts the cross section limit curve to the left while not fully affecting the asymptote at the mass threshold, for which v esc is also relevant.

• Velocity dispersion: The larger the dispersion, the larger the sensitivity peak (in the SHM, it is fixed by v c ).

• Local DM density: Global linear and vertical shift of the exclusion curve.

At this stage, we can already notice a few features from P14 results with the help of Fig. 1. First, we see that there is a strong correlation between the circular velocity v c and the local DM density ρ , while the latter is poorly correlated with v esc . Second, we note that though the two P14 best-fit points with priors on v c exhibits an anticorrelation between v c and v esc , this is no longer the case when v c is left free. While this is somewhat contradictory with P14 comments, this implies that the v c -free case would provide constraints on v esc independent of v c and ρ . Since there are different choices in P14 as for either the priors on v c or the prior on the dark halo concentration, there may be several ways to investigate how the constraints on v esc affect the DD exclusion curves. Here, we first consider the three P14 results independently before comparing them together. For the sake of illustration, we will compute the DD limits by using the LUX results [9] as a reference for Xenon experiments, the SuperCDMS results [20] as a reference for Germanium experiments, and the CRESST-II results [21] as a reference for multi-target experiments. As of today, these experiments provide the most stringent combined bounds on the spin-independent cross section for low WIMP masses (assuming elastic scattering and no isospin violation). More technical details are given in the Appendix.

In Fig. 3, we show the exclusion curves obtained in the σ p -m χ plane with the 90% CL contours associated with the three P14 best-fit points. These curves have been derived by using ergodic DFs that self-consistently correlate the dark halo parameters with v esc and v c . We also report on the same plot the exclusion curves calculated from the SHM, with the parameters given in Eq. ( 26). We may note that the best-fit points associated with priors on v c lead to similar curves and bands, in particular in the low WIMP mass region. This is due to the anticorrelation between v c and v esc , which is such that the sum v c + v esc , relevant to characterize the WIMP mass threshold, remains almost constant. From Fig. 1, we also see that ρ spans only a reduced range of ∼ [0.35, 0.47] GeV/cm 3 for these two points, which leads to less than ±15% of relative difference in the predicted event rate. In contrast, the exclusion curves associated with the v c -free case have (i) a slight offset toward larger WIMP masses while (ii) lying slightly above the others with a larger uncertainty. This is due to (i) the prior on the concentration parameter that forces small values of v c while not affecting v esc , and, as a result, to (ii) lower values of ρ spanning the less favorable, while larger range [0.15, 0.35] GeV/cm 3 . We will further comment on the v c -free case in Sec. III D. In comparison, the SHM curves (aimed at reproducing the limits published by the experimental collaborations) lie in between, and are less constraining than the P14 parameters with priors on v c .

As a reference case to evaluate the uncertainties coming from these P14 results on v esc let us first focus on the best-fit P14 point with the v c = 240 km/s prior, which we may assume to be the most motivated case (see discussion in Sec. III D 1). In Fig. 4, we display the exclusion curves and the associated 90% CL uncertainties for this astrophysical configuration. The LUX (Super-CDMS, CRESST-II) limits are shown in the left (middle, right) panels; the absolute (relative) uncertainties are given in the top (bottom) panels. We also compare the limits obtained from the ergodic DF to those calculated from the SHM model on the one hand, and from the MB DF with the P14 values for the astrophysical parameters on the other hand. The fact that the P14 escape speed estimate v esc = 511 +48 -35 (90% CL) converts into ρ = 0.43 ± 0.05 GeV/cm 3 explains why the ergodic limit is more constraining, by ∼ 40%, than the SHM one (ρ = 0.3 GeV/cm 3 ) over a large part of the depicted WIMP mass range. In contrast, the SHM limit beats the ergodic one at very low WIMP masses because v esc itself affects the effective WIMP mass threshold: it is set to 544 km/s in the SHM (i.e. v c + v esc = 764 km/s), while the P14 best-fit point (for v c = 240 km/s) corresponds to v esc =511 km/s (i.e. v c + v esc = 751 km/s). Interestingly, we also see the differences induced by different DFs while using the same astrophysical inputs, by comparing the MB and ergodic curves. The former is actually more constraining at energy recoils leading to v min larger than the peak velocity of the DFs. This is because the MB DF exhibits a less steep tail at high velocities than the ergodic DF (see Fig. 2, right panel). This illustrates why not only is the escape speed important in the low WIMP mass region, but also the high velocity tail of the DF, and thereby the DF itself (see similar discussion in Ref. [66]). Finally, we may notice that the relative uncertainties in the exclusion curves for this specific P14 point saturate around ∼ ±10% (90% CL) at large WIMP masses, which is mostly set by the allowed range in ρ .

In the low WIMP mass region, the high-velocity tail and v esc come also into play, and the uncertainties strongly degrade, obviously, as the maximum possible recoil energy approaches the experimental energy threshold. This can clearly be observed in the case of LUX, whose efficiency drops for WIMP masses below ∼ 8 GeV (Fig. 4, bottom left panel). However, we may also remark the nice complementarity between the different experiments that employ different target atoms, as one reaches a relative uncertainty ∼ ±30% for the CRESST-2 proxy down to WIMP masses ∼ 3 GeV. The bumps appearing in the bottom middle and right panels of Fig. 4 (corresponding to SuperCDMS and CRESST-II) not only comes from the differences in the atomic targets, but also from the impact of the observed nuclear-recoil-like events on the maximum gap method [67] that we used to derive the exclusion curves (see more details in the Appendix). At this stage, it is still difficult to draw strong conclusions as for the overall uncertainties in the exclusion curves induced by the P14 results without questioning more carefully the initial assumptions or priors -we will proceed so in Sec. III D. Nevertheless, we can already emphasize that a self-consistent use of these estimates of v esc is not straightforward (for instance, one cannot just vary v esc in a given CL range irrespective of the other astrophysical parameters). A proper use must take the correlations between all the relevant astrophysical parameters into account. Indeed, we have just seen that not only is the WIMP mass threshold affected (direct consequence of varying v esc + v c ), but also the global event rate, as ρ FIG. 4: Top panels: experimental limits calculated using the P14 best-fit point with the prior v c = 240 km/s (ρ = 0.43 GeV/cm 3 ; v esc = 511 km/s), and associated 90% CL contour. Bottom panels: relative uncertainties with respect to the central value. Left panels: LUX data. Middle panels: SuperCDMS data. Right panels: CRESST-II data. See text for details.

must be varied accordingly.

D. Discussion

In this section, we wish to examine the previous results in light of independent constraints on the astrophysical inputs. As we saw, P14 provided three best-fit configurations each based on different priors. The most conservative approach would be to relax fixed priors as much as possible, as all astrophysical parameters are affected by uncertainties. We will therefore focus on the v c -free case in the following.

Additional and independent constraints on vc

In Fig. 1, we see that the prior on the concentration model forces the fit to low values of v c , with a central value v c = 196 km/s. Nevertheless, this prior is based upon Ref. [55], which studies the halo massconcentration relation in cosmological simulations that do not contain baryons. Since we now know that baryons may affect the structuration of dark matter significantly on galactic scales (see e.g. [68][69][70]), the motivation to use such a constraint appears to us to be rather weak. In the following, we will therefore relax this prior and rather consider the entire range provided by the P14 likelihood in the v c -v esc plane.

If we relax the prior on the concentration law, we see that although the astrophysical parameters relevant to DD are not much correlated with v esc (see the explanation and warning in Sec. II C), which lies in the range ∼ 500-570 km/s, the local DM density ρ is strongly correlated with the circular velocity v c , which was explored over a large range of ∼ 160-260 km/s. We recall that P14 originally explored the M 340 -c 340 plane that we converted into v c -v esc for the sake of DD interpretations.

The circular velocity v c can actually be constrained independently of any assumption on the dark halo modeling. Indeed, for instance, one can use measures of trigonometric parallaxes and proper motions of stars to reconstruct the local kinematics, where the only model ingredients to consider are the circular and peculiar motions of the Sun in the Milky Way rest frame, and the distance of the Sun to the Galactic center. It turns out that the ratio v c /r is often better constrained that v c and r individually. More precisely, in the absence of prior on the peculiar velocity of the Sun, these observations constrain the ratio (v c + V )/r (see for instance Refs. [32,54,71]). Since both r and the peculiar velocity of the Sun have been fixed in P14, we could still use these independent constraints on v c /r to delineate an observationally motivated range for v c . Fortunately, the recent study in Ref. [54] on local kinematics, based on a large statistics of parallaxes and proper motions of masers associated with young and massive stars, provides results that can be used rather straightforwardly as some of their priors are similar to those of P14. This avoids to rely on the (v c + V )/r ratio, which induces the risk of underestimating the errors for fixed values of V and r . Indeed, in their Bayesian fit model B1 (associated with the cleanest data sample), the authors of Ref. [54] take as a prior on the peculiar velocity the results of Ref. [49], also chosen by P14 (see Eq. 2). They then derive a best fit value for r = 8.33 ± 0.16 kpc, fully consistent with the 8.28 kpc adopted in P14. Accordingly, they get the following range for v c :

v c = 243 ± 6 km/s (1 σ), (34)
that we can use directly (we will actually use the 2-σ range ±12 km/s for a more realistic discussion on the uncertainties -the posterior PDF is found close to Gaussian [54]). We report this range as a green band in Fig. 1, which represents an independent constraint in the v c -v esc plane. We emphasize that this additional constraint depends very little, if not at all, on the dark halo model, as these local kinematics fits. We only need to make sure that our Galactic mass models are consistent with a last feature of model B1 of Ref. [54] regarding the local derivative of the circular velocity, dv c /dR = -0.2 ± 0.4 km/s/kpc, pointing to an almost flat local rotation curve.

We have checked that this is the case almost everywhere in the relevant part of the uncertainty band of Fig. 1 within three standard deviations. The largest tensions with respect to the range aloowed for dv c /dR arise for the lowest values allowed for v esc , which force the scale radius r s of the NFW profile to be smaller than r and thereby affect the local slope of the rotation curve. This is not favored by independent analyses (see e.g. Ref. [32]), and interestingly enough, this disfavored zone also corresponds to the lowest values of the local DM density.

In Fig. 1, we have reported this 2-σ range for v c . It crosses the full blue band of the v c -free case of P14, which allows to define a constrained region in the plane v c -v esc , that we will use to determine the uncertainties implied for DD exclusion curves. Since the blue band corresponds to likelihood values largest than 1% of the maximum, we may interpret the corresponding DD uncertainties as ∼ 2σ uncertainties. We note that this region allows for large values of the local DM density, up to 0.55 GeV/cm 3 . This is actually consistent with the tendency found in recent independent studies (see e.g. Refs. [37,72,73]), and reinforces the potential of DD experiments.

In Fig. 5 we show the uncertainties obtained when considering the combination of the P14 v c -free case with the additional constraints on v c discussed above. These plots are similar to those of Fig. 4 as only the uncertainty contours change. In particular, we have also reported the exclusion curves obtained with the P14 "v c =240 km/s" best-fit configuration, which are shown to lie within the contours of the v c -free case (plus independent and additional constraints on v c ). This can easily be understood from Fig. 1. The changes in the uncertainty contours can again be understood in terms of local DM density, which is found in the range [0.37,0.57] GeV/cm 3 . This sets the relative uncertainties to ∼ ±20% for large WIMP masses, while they further degrade toward low WIMP masses be-cause of the additional effects from v esc and v c .

Speculating beyond P14

Because of the caveats affecting the v c -free case of P14 (see discussion in Sec. II C), and in order to try to recover a better consistency in the Galactic mass modeling, we may try to further speculate about what a self-consistent v c -free band would look like. P14 argument that v c and v esc should linearly anticorrelate is sound as it is based on a purely geometrical reasoning (the fact that the stellar sample is completely biased toward negative longitudes). Therefore, it is quite reasonable to think that the dotted lines in Fig. 1, that relate the two best-fit points and errors associated with the 220 and 240 km/s priors on v c , may represent a more consistent range than the v c -free band. As before, we also investigate how the uncertainties export to DD, using again the additional 2-σ constraints on v c from Ref. [54] -see Eq. (34).

We show our results in Fig. 6, similar to those of Fig. 5, though likely more consistent with the original data used in P14 (because of the supposed anticorrelation between v esc and v c ). As expected, the relative uncertainties improve on the entire WIMP mass range. They improve down to ∼ ±10% in the large WIMP mass region. While not spectacular, this still illustrates the need to account for as many astrophysical correlations as possible in deriving the DD limits.

IV. CONCLUSION

In this article, we have studied the impact of the recent estimate of the Galactic escape speed v esc from Ref. [74] (P14) in deriving the DD exclusion curves. We note that this observable is difficult to reconstruct, and the method used in P14, as recognized by the authors, is potentially subject to systematic errors. Nevertheless, these constraints are independent of those coming from studies of rotation curves, and thereby may provide complementary information on the WIMP phase-space. We have shown that the conversion of these results is highly non-trivial, as the constraints on v esc are obtained from a series of assumptions that relate each value of v esc to a different set of parameters for the dark halo profile. This implies that one cannot use the different v esc ranges provided by P14 blindly and irrespective of these assumptions (as it is often done by implementing them with flat or Gaussian priors in likelihood or Monte Carlo Markov chain calculations). We have assumed spherical symmetry and gone beyond the Maxwell-Boltzmann approximation by considering ergodic WIMP velocity DFs. This method ensures a self-consistent physical connection between the DFs and the underlying Galactic mass model. In particular, all local variables relevant to DD calculations, i.e. the average WIMP speed, the dispersion velocity, the local DM density, are consistently dynamically correlated in this approach.

We have studied the three best-fit configurations provided by P14, two with priors on the circular velocity v c (220 and 240 km/s), and one with v c left free. These configurations are shown in the v c -v esc plane in Fig. 1. As the first two appeared to us too specific, we concentrated on the latest, which was originally optimized for Milky Way-mass estimates (we had to convert the P14 results from the c 340 -M 340 plane to the v c -v esc plane). We discovered that the anticorrelation between v c and v esc arising from the biased locations of the stars of the RAVE-P14 sample toward negative longitudes, and found in the v c = 220-240 km/s cases, was no longer present. This is actually due to the use of the posterior PDF for v esc with the v c = 220 km/s prior as a proxy, which could hardly be divined from Fig. 13 in P14. We therefore further speculated on what a fully self-consistent v c -free case could look like, and considered this guess as an alternative. Finally, we accounted for independent constraints on v c from Ref. [54], which have the advantage of being almost fully independent of the Galactic mass model. This has driven us to favor large v c regions, around 240 km/s, which are in principle associated with lower escape speeds in P14.

We have translated these P14 results in terms of DD exclusion curves focusing on the LUX, SuperCDMS, and CRESST-II experiments as references. We emphasize that a self-consistent use of these P14 results imply rather large values of the local DM density, and thereby more constraining exclusion curves (up to ∼40% more). Interestingly enough, this is consistent with several independent recent results on the local DM density (see e.g. Refs. [37,72,73]). This is good news for direct DM searches as this tends to increase their potential of discovery or exclusion. We have also investigated the associated relative uncertainties, and shown that they are highly non-trivial as P14 values of v esc are correlated with other astrophysical parameters, as already stated above. We have shown that taking P14 results at face value (plus eventually additional independent constraints on v c ) converts into moderate uncertainties, down to ∼ ±10% in the regime where the experiments can trigger on the whole phase-space (large WIMP masses). When getting closer to the high-velocity tail (low WIMP masses), the uncertainties explode as the experimental efficiency drops, but we have illustrated the nice complementarity between the experiments using different target atoms in this regime. This complementarity allows to maintain a moderate uncertainty of ∼ ±20% down to WIMP masses of a few GeV. Nevertheless, since the SHM value for the sum v c + v esc is 764 km/s (based on S07), i.e. slightly less than the 751 km/s found from the more recent P14 estimate (with the prior v c = 240 km/s), the generic outcome is that we find an effective WIMP mass threshold slightly heavier than in the SHM.

There are several limitations in our analysis. First, the Galactic mass model is rather simple and the baryonic content has been fixed in P14. Unfortunately, we cannot FIG. 6: Same as Fig. 5, but for the v c -v esc range supposed to properly take the anticorrelation into account (cite dotted lines in Fig. 1).

go beyond this choice as this would no longer be consistent with P14 results. Second, we made the assumption that the WIMP phase-space was governed by energy in a spherically symmetric system, which led to the derivation of ergodic velocity DFs. While this approach selfconsistently correlates the local velocity features and the local DM density to the full gravitational potential, it remains to be investigated in detail whether it reliably captures the dynamics at stake in spiral galaxies. Some works do indicate that this approach provides a reasonable description of cosmological simulation results (see e.g. Ref. [65,66]), but it is obvious that more studies will be necessary to clarify this issue. We plan to discuss this into more detail in a forthcoming paper. Still, it is less ad-hoc an assumption than using the Maxwell-Boltzmann velocity distribution. Finally, we stress that our study is complementary to those based e.g. on rotation curves, as it relies on different, and independent, observational constraints. Further merging these different sets of constraints would be interesting in the future. Moreover, because of these dynamical correlations inherently arising in the local astrophysical parameters, which was continuously underlined along this paper, several improvements could also be expected in the complementary use of direct and indirect detection constraints in order to exclude or validate some WIMP scenarios. recoils. Of these events, 24 fall into the calibration nuclear recoil band. Nevertheless, as these events are distributed only in the lower half of the band, while WIMP events should spread over the full range, it is unlikely that a significant part of them really comes from WIMP scatterings. In [76], several configurations of event counting are considered to feed the MGM procedure: 0, 1, 3, 5 or all the 24 events. The best matching with the experimental result is obtained for the MGM method run with between 0 and 1 event. As the former better reproduces the experimental curve at low WIMP masses, we adopt this one.

For the experimental efficiency after cuts, we have used Fig. 9 of [9]. As in [76], we set the counting efficiency to 0 below 3 keVnr. We use a Gaussian energy resolution, with a dispersion σ(E nr ) = σ PMT Q(E nr ), with Q(E nr ) = 4.131 Enr/KeV + 0.690 and an S1 single photoelectron resolution of σ PMT = 0.37 photoelectrons.

An indicative conversion from S1 and S2 signals to E nr can be inferred from the contour lines in Fig. 4 of [9]. We have used the relation

E nr = S1 L y L eff S e S n (A.1)
with the light yeld L y = 3 photo-electrons/keVee and the scintillation quenching S e = 0.54 for electron recoils and S n = 0.93 for nuclear recoils from Ref. [77]. Even if the value of the Lindhard factor L eff in liquid xenon is still subject to debate, we simply assumed the value L ef f = 0.14 in this analysis which is not aimed at investigating the compatibility of exclusions with putative signal regions. Finally, we use the Helm Form Factor [58] to model effects of the nuclear shape on the elastic scattering.

CRESST-II

The CRESST experiment is a multi-target detector made of CaWO 4 crystals -calcium (Ca), oxygen (O) and Tungsten (W). We base our limits on the recent results released in Ref. [21], relying on an exposure (before cuts) of 29.35 kg•day. The energy range used in the data analysis is [0. 6,40] keV, for which we consider all the events collected, i.e. around 77 events. The energies of these events can be inferred from Fig. 1 and from the inset of Fig. 2 of Ref. [21]. We obtain the 90% CL upper limits by means of the MGM, while the original analysis employs the Optimum Interval Method. We use the following atomic target fractions: f Ca = 1/6, f O = 4/6, and f W = 1/6. We accounted for a Gaussian energy resolution with σ = 0.107 keV, and took the experimental efficiency after cuts from the blue curve in Fig. 3 of Ref. [21].

SuperCDMS

SuperCDMS is a detector using Ge crystals. A recent data analysis focused on low-mass WIMPs has been released in [20], based on an exposure of 577 kg•day. The energy range used in the data analysis is [1.6, 10] keV nr , in which 11 candidate events have been observed for a background expectation of 6.198 events. Their energies are listed in Table I of Ref. [20]. While the experiment limit is derived from the optimum interval method without background subtraction [67], we still use the MGM which gives similar results. We take into account the experimental efficiency that we obtain from the red curve of Fig. 1 in Ref. [20], and which is an exposure-weighted sum of the measured efficiency for each detector and period. We assumed a Gaussian energy resolution, and for lack of better information we used the same energy dispersion as in CDMSLite [78], σ = 14 eV ee (which corresponds to 87.5 eV nr after conversion by means of the Lindard theory, as done in Ref. [78]).

M b = 1 . 5 ×

 15 10 10 M r b = 0.6 kpc Disk M d = 5 × 10 10 M R d = 4 kpc z d = 0.3 kpc

FIG. 2 :

 2 FIG. 2: Comparison between a truncated MB DF and an ergodic DF. Left panel: Assuming the parameters of P14 best-fit point with v c = 220 km/s, we show the MB curve and several curves for ergodic case to make the impact of the baryonic components explicit. Right panel: Comparison of the MB-like DFs with the ergodic DFs assuming the parameters of the P14 best-fit points with both v c = 220 and 240 km/s. Here, we also show the DFs associated with the upper and lower corners of the 90%CL interval.

TABLE I :

 I Values of the parameters associated with the baryonic components of the mass model.

  10% of Max L RAVE, 1% of Max L Combined contours RAVE best, v c =220 km/s Maccio et al. (2008) Simulations, 10% of Max L Simulations, 1% of Max L Constraints from Reid et al., 2σ RAVE best, v c =240 km/s

		Lavalle & Magni (2014)
	RAVE parameter space RAVE,	2 5 0 1 . 3 0 0

FIG. 1: P14 likelihood results converted in the plane v esc -v c . See text for details.

TABLE II :

 II Values of the DM halo parameters associated with P14 best-fit points. The reported errors correspond to 90% CL and 99% CL, the latter being inferred from the PDFs reported in Fig.11in P14. The 99%CL cannot be inferred for the "v c = 240 km/s" case from P14. For the "v c -free" case, the errors are taken from Fig.13of P14see also Fig.1of this paper; we report the largest uncertainties corresponding to likelihood values larger than 1% of the maximum. See text for further details.

	0.006
	0.005
	0.004
	0.003
	0.002
	0.001
	0 0.000	100	200 speed v in LSR (km/s) 300	400	500

  RAVE best Fitv c =220 km/s RAVE best Fit v c =240km/s RAVE best Fit v c km/s free, prior on c 340 SHM Unc. on LUX, v c =220 km/s Unc. on LUX v c =240 km/s Unc. on LUX, free v c Unc. on CRESST2, v c =220 km/s Unc. CRESST2, v c =240 km/s Unc. on CRESST2, free v c Unc. on SuperCDMS, v c =220 km/s Unc. on SuperCDMS, v c =240 km/s Unc. on SuperCDMS, free v c

	10 -39	Lavalle & Magni (2014)
	10 -40
	)
	σ pSI (cm 2

FIG.

3

: P14 90% CL results on v esc converted into DD exclusion curves.

  .H.M., v c =220 km/s, M.B. RAVE best, v c =240 km/s, M.B. RAVE Best, v c =240 km/s, Ergodic RAVE 1% of Max L, free v c , Ergodic RAVE 10% of Max L, free v c , Ergodic .H.M., v c =220 km/s, M.B. RAVE best, v c =240 km/s, M.B. RAVE Best, v c =240 km/s, Ergodic RAVE 1% of Max L, free v c , Ergodic RAVE 10% of Max L, free v c , Ergodic .H.M., v c =220 km/s, M.B. RAVE best, v c =240 km/s, M.B. RAVE Best, v c =240 km/s, Ergodic RAVE 1% of Max L, free v c , Ergodic RAVE 10% of Max L, free v c , Ergodic .H.M., v c =220 km/s, M.B. RAVE best, v c =240 km/s, M.B. RAVE best, v c =240 km/s, Ergodic RAVE 1% of Max L, free v c , Ergodic RAVE 10% of Max L, free v c , Ergodic .H.M., v c =220 km/s, M.B. RAVE best, v c =240 km/s, M.B. RAVE best, v c =240 km/s, Ergodic RAVE 1% of Max L, free v c , Ergodic RAVE 10% of Max L, free v c , Ergodic v c =240 km/s, Ergodic RAVE 1% of Max L, free v c , Ergodic RAVE 10% of Max L, free v c , Ergodic FIG. 5: Top panels: experimental limits calculated with a combination of the P14 v c -free analysis with additional constraints on v c from Ref. [54]. Bottom panels: relative uncertainties with respect to the v c = 240 km/s best-fit point of P14. Left panels: LUX data. Middle panels: SuperCDMS data. Right panels: CRESST-II data. See text for details.

		10 -42			Lavalle & Magni (2014)	10 -37	Lavalle & Magni (2014)	10 -37		Lavalle & Magni (2014)
						10 -38		10 -38
		10 -43				10 -39		10 -39
	σ pSI (cm 2 )				σ pSI (cm 2 )	10 -42 10 -41 10 -40	σ pSI (cm 2 )	10 -41 10 -40
						10 -43		10 -42
						10 -44	10 1 m χ (GeV/cm 3 )	10 0 10 -43	10 1 m χ (GeV/cm 3 )	10 2
		10 1			Lavalle & Magni (2014)	10 1	Lavalle & Magni (2014)	10 1	Lavalle & Magni (2014) S.H.M., v c =220 km/s, M.B. RAVE best, v c =240 km/s, M.B. RAVE best,
	σ pSI /σ pSI,BestErgodic	10 0			σ pSI /σ pSI,BestErgodic	10 0	σ pSI /σ pSI,BestErgodic	10 0
		10 -1	10 1	m χ (GeV/cm 3 )	10 2	10 -1	10 1 m χ (GeV/cm 3 )	10 -1	10 1 m χ (GeV/cm 3 )	10 2

SSSSS

  .H.M., v c =220 km/s, M.B. RAVE best, v c =240 km/s, M.B. RAVE Best, v c =240 km/s, Ergodic RAVE Forced Correlation v c /v esc , Ergodic .H.M., v c =220 km/s, M.B. RAVE best, v c =240 km/s, M.B. RAVE Best, v c =240 km/s, Ergodic RAVE Forced Correlation v c /v esc , Ergodic .H.M., v c =220 km/s, M.B. RAVE best, v c =240 km/s, M.B. RAVE Best, v c =240 km/s, Ergodic RAVE Forced Correlation v c /v esc , Ergodic .H.M., v c =220 km/s, M.B. RAVE best, v c =240 km/s, M.B. RAVE best, v c =240 km/s, Ergodic RAVE Forced Correlation v c /v esc , Ergodic .H.M., v c =220 km/s, M.B. RAVE best, v c =240 km/s, M.B. RAVE best, v c =240 km/s, Ergodic RAVE Forced Correlation v c /v esc , Ergodic

		10 -42			Lavalle & Magni (2014)	10 -37	Lavalle & Magni (2014)	10 -37		Lavalle & Magni (2014)
						10 -38		10 -38
		10 -43				10 -39		10 -39
	σ pSI (cm 2 )				σ pSI (cm 2 )	10 -42 10 -41 10 -40	σ pSI (cm 2 )	10 -41 10 -40
						10 -43		10 -42
						10 -44	10 1 m χ (GeV/cm 3 )	10 0 10 -43	10 1 m χ (GeV/cm 3 )	10 2
		10 1		Lavalle & Magni (2014) RAVE Forced Correlation v c /v esc , Ergodic RAVE best, v c =240 km/s, Ergodic RAVE best, v c =240 km/s, M.B. S.H.M., v c =220 km/s, M.B.	10 1	Lavalle & Magni (2014)	10 1		Lavalle & Magni (2014)
	σ pSI /σ pSI,BestErgodic	10 0			σ pSI /σ pSI,BestErgodic	10 0	σ pSI /σ pSI,BestErgodic	10 0
		10 -1	10 1	m χ (GeV/cm 3 )	10 2	10 -1	10 1 m χ (GeV/cm 3 )	10 -1	10 1 m χ (GeV/cm 3 )	10 2

SSSSS

[hep-ex].

ACKNOWLEDGMENTS

We thank Tilmann Piffl for very useful exchanges about the details of the data analysis conducted in P14. We are also grateful to Benoit Famaey and Paul McMillan for valuable discussions about kinematics studies. This work was partly supported by the CFT-Théorie-IN2P3 and by the Labex OCEVU. The bibliography was generated thanks to the JabRef software [75].

Appendix: Direct detection data and limits

For every experiment, we briefly summarize the most important quantities that allow us to reproduce the experiment results.

LUX

LUX is a liquid xenon time-projection chamber experiment which currently runs with a 370 kg target. We use the results released by the collaboration in Ref. [9]. The first run described there collected 85.3 days of data, with a detector fiducial xenon mass of 118 kg. The next 300-day run expected for end 2014 should improve the sensitivity by a factor of 5.

As the published results are based on a private likelihood analysis, we cannot reproduce them exactly. Instead, we follow the approach detailed in [76], relying on the maximum gap method (MGM) [67] over an S1 range of 2-30 photoelectrons, the same as the one used in the private likelihood. In this range, 160 events are observed, all consistent with the predicted background of electronic