Determination of the local dark matter density using K-dwarfs from Gaia DR2

Local dark matter density, ρdm, is one of the crucial astrophysical inputs for the estimation of detection rates in dark matter direct search experiments. Knowing the value also helps us to investigate the shape of the Galactic dark halo, which is of importance for indirect dark matter searches, as well as for various studies in astrophysics and cosmology. In this work, we performed kinematics study of stars in the solar neighborhood to determine the local dark matter density. As tracers we used 95,543 K-dwarfs from Gaia DR2 inside a heliocentric cylinder with a radius of 150 pc and height 200 pc above and below the Galactic mid plane. Their positions and motions were analyzed, assuming that the Galaxy is axisymmetric and the tracers are in dynamical equilibrium. We applied Jeans and Poisson equations to relate the observed quantities, i.e. vertical position and velocity, with the local dark matter density. The tilt term in the Jeans equation is considered to be small and is therefore neglected. Galactic disk is modelled to consist of a single exponential stellar disk, a thin gas layer, and dark matter whose density is constant within the volume considered. Marginalization for the free parameters was performed with Bayesian theorem using Markov Chain Monte Carlo (MCMC) method. We find that ρdm = 0.0116 ± 0.0012 M /pc3 or ρdm = 0.439 ± 0.046 GeV/cm3, in agreement within the range of uncertainty with results of several previous studies.


Introduction
Efforts to understand the properties of dark matter are constantly being made through both indirect and direct experiments. The later detect interactions between dark matter particles and nuclei in the detector materials, which are usually assembled underground. For direct experiments, local dark matter density, ρ dm , is an important parameter for estimation of the detection rate of dark matter particles. ρ dm is defined as the value of dark matter density in the solar neighborhood. Although there are no tight restrictions set by the term "local", ρ dm is commonly determined by examining a small volume around the Sun, typically within the height of a few hundred parsecs up to ∼ 10 3 pc from the Galactic mid plane [1].
Aside from its importance to estimate the detection rate of dark matter particles, it is also intriguing to obtain the value of ρ dm , because knowing it can help revealing the actual shape of the Galactic dark halo -which has been assumed to have a spherical symmetry. This is done by comparing ρ dm with the value obtained from global dark matter density profile through decomposition of Milky Way's rotation curve [2,3]. If the value of ρ dm is smaller, then the Galaxy has a prolate-shaped halo. A higher value of ρ dm indicates that Milky Way's halo is oblate with a possibility of a dark-disk structure.
The value of ρ dm from determinations in the past few years lies within the range of 10 −3 M /pc 3 and 0.05 M /pc 3 [4][5][6][7][8][9]. Methods to determine local dark matter density are constantly being updated along with the increase of the quality and amount of observational data, such as obtained by Gaia, a mission initiated by the European Space Agency (ESA) with astrometry as the main focus. In particular Gaia measures positions and velocities of stars in the Galaxy [10]. Data obtained by Gaia has a precision up to miliarcseconds for position and proper motion measurements, and a range of magnitude of 3 ≤ G ≤ 21 for photometric measurement. Gaia's photometrics G passband covers a wavelength range between 330 nm and 1050 nm. At the moment, the most recent data available is Gaia DR2, published on 25th April 2018 [11].

Method
One of the possible approaches to determine the local dark matter density is by observing the motion of stars near the Sun as tracers. Stars move based on their perceived gravitational potential, while gravitational potential is generated by mass, both luminous or dark matter. If we are able to model the distribution of baryonic matter, then the difference between the total density and the density of baryonic matter informs the density of non-baryonic dark matter.
The populations of stars selected as tracers are, where possible, the ones in dynamical equilibrium. Dynamical equilibrium is important as it can simplify the collisionless Boltzmann equation. Generally, early-type stars of class O and B are never in dynamic equilibrium due to their short lifetime. Several similar works by Garbari et al. (2012) [4], Xia et al. (2016) [7], and Sivertsson et al. (2018) [8], for instance, with the consideration that a long lifetime will produce a larger chance for a star population to reach dynamical equilibrium, selected longer-lived stars as tracers, rather than early-type stars.

Jeans equation
As we can assume the stars in the Galaxy as a collisionless system, to study its dynamics, we can use the collisionless Boltzmann equation where f is a six dimensional phase space distribution (three dimensions for space and three dimensions for velocity) of stars in the Galaxy, v is star velocity, and Φ is gravitational potential.
The collisionless Boltzmann equation fulfills d f /dt = 0. This means that the phase space distribution for a particle does not change with respect to time. Hence, if we follow a star in its orbit, the six-dimensional phase space density surrounding the star is constant. For a system in dynamical equilibrium, the term ∂ f /∂t is 0. This study assumes that the populations of stars used as tracers are in dynamical equilibrium.
If the collisionless Boltzmann equation is integrated with respect to velocity, the Jeans equation will be acquired. In a cylindrical coordinate system, the Jeans equation can be expressed as where ν is the tracer number density, and σ is the tracer vertical velocity dispersion. Respectively, − dΦ dz = K z , 1 R ∂ ∂R (Rνσ Rz ) = T (z), and 1 R ∂ ∂φ (νσ φz ) = A(z) are called vertical acceleration, tilt term, and axial term. We assume that the Milky Way is axisymmetric, meaning the axial term in the Jeans equation can be neglected [4,7].
In the Jeans equation, the motion of radial and vertical components are only coupled through the first term in the right hand side called the tilt term. Read (2014) [1] explains how the tilt term can be neglected with respect to other terms for the region in proximity to the Galactic plane, reducing Jeans equation into only in the z direction

Poisson Equation
The relation between gravitational potential and matter density is explained through the Poisson equation where ρ tot being the total matter density. If Equation 4 is expressed in a cylindrical coordinate system, it reads Under the assumption that the Galaxy is axisymmetric, the potential derivative in the φ direction can be neglected, hence 1 4πG The term 1 ∂R is known as the rotation curve term. Observations show that for the region in proximity to the Sun, the rotation curve of the Galaxy tends to be flat, resulting in a small value of ∂V 2 C /∂R. The value of this term can be calculated using Oort's constants A and B [12], 1 4πGR with A and B as Oort's constants. Bovy (2017) [13] gives the value of A = 15.3 ± 0.4 km s −1 kpc −1 , and B = −11.9 ± 0.4 km s −1 kpc −1 , which we used in our work. Substituting those values, we have the value of rotation curve term to be around ∼ (−10 −22 ). Therefore, that term shall be neglected with respect to the other terms. Thus, Poisson equation is reduced to only in the z direction, The total matter density, ρ total , is assumed to consist of stars, gas, and dark matter.

Vertical Acceleration Modelling
Vertical acceleration is the first derivative of potential with respect to the vertical distance, Substituting to Equation 3, we will obtain To model vertical acceleration, we follow Xia et al. (2016) [7]. Several additional assumptions need to be taken: • the vertical distribution of stellar mass follows an exponential profile, ρ * (z) = ρ * ,0 exp (−z/h).
• gas is distributed in a negligible thickness. Hence, the surface density of the gas, Σ gas , is assumed to be independent of z. In this study, we used the value given by McKee et al.
• the dark matter density is constant within the considered volume, ρ dm (z) = ρ dm .
Using those assumptions, the vertical distribution of matter can be expressed as The expression for vertical acceleration can be acquired by integrating Poisson equation that gives where Σ * is the surface mass density of the stellar disk. The vertical velocity dispersion can be obtained analytically by integrating Equation 11 To compute the integral, the tracer number density is modeled with an exponential profile thus, Equation 15 becomes that gives where Equation 18 is the main equation to be solved in this study to find the local dark matter density ρ dm , based on the observation of tracers' number density ν(z), and the tracers' vertical velocity dispersion σ z .

Data
The considered volume in this study is a cylinder with a radius of 150 pc (0.9% of the radius of the Milky Way disk) and a half-height of z = 200 pc with the Sun on the cylinder's axis and the Galactic plane dividing the cylinder into two equal parts. There is no strict rules for choosing the considered volume size and shape. However, some considerations need to be taken into account. First, the considered space should be relatively small in comparison to the Milky Way's disk, which means that the Galaxy rotation velocity in the considered volume is relatively constant. Second, the considered volume has to be relatively small compared to the dark matter scale-length, which is 19 kpc [1], so that the dark matter density in the considered volume can be assumed to be constant. Finally, a cylindrical shaped was chosen to ease the data processing steps, which involve two equations which are functions of z. Stellar population that is used as tracers must satisfy these following criteria: 1. the population must be old enough that it is possible for it to have reached dynamical equilibrium. Hence, collisionless Boltzmann equation can be simplified by omitting ∂ f /∂t.
2. the population must contain a sufficiently large number of stars to represent statistically good sample size.
3. the members of the population must be spread throughout the considered volume.
Unfortunately, Gaia DR2 has several limitations that narrow the options for tracers. The shortfall of Gaia relevant to our case are: the range of effective temperatures for stars with T e f f data lies within 3, 000K < T e f f < 10, 000K and only stars within the magnitude range of 4 < G < 13 and effective temperature range of 3, 550 < T e f f < 6, 900K have radial velocity data. Realizing Gaia DR2 restrictions and concerning our criteria for tracers, we select Kdwarfs as tracers in this study. Our selection process gives 95,543 K-dwarfs with positional data, so that they can be used in fitting of tracer number density. Among them, 46,466 having data of position and velocity and can be used to calculate vertical velocity dispersion.

The likelihood
The free parameter marginalization was performed based on Bayesian theorem using Markov Chain Monte Carlo (MCMC) method. Bayesian theorem can be expressed as The likelihood function used in this study is, where L is the likelihood, σ data, j (z) is the vertical velocity dispersions as a function of height obtained from observational data, σ model, j (z) is the vertical velocity dispersion as a function of height based on the model explained in section 2, and σ j (z) is the vertical velocity dispersion errors from observational data.

Tracer Number Density Fitting
The considered cylindrical volume was divided in 10 pc bins in the z direction. The tracer number density was obtained by calculating the number of tracers in a bin divided by with the volume of the bin. The Galaxy was assumed to be symmetrical to the Galactic plane. The fitting error of ν(z 0 ) and h t was approximated by equations, and where N,z, ν(z 0 ) , h t , and m are the number of tracers in a bin, the average height of tracers in a bin, the tracer number density error on the Galactic plane, the tracer altitude scale error, and the number of parameters, respectively. We have two parameters fitting here, namely ln[ν(z 0 )] and h t . Meanwhile, b 0 and b 1 are errors of parameters b 0 and b 1 . From the fitting, the value of tracer number density on Galaxy plane, ν(z 0 ) = 0.0024 ± 0.0001 pc −3 , and tracer altitude scale, h t = 427.25 ± 1.54 pc, were obtained.

Tracer Vertical Velocity Dispersion Calculation
Tracer vertical velocity, v z , was obtained using data of equatorial coordinates, proper motion, parallax, and radial velocity. Velocity transformations from µ α * , µ δ and v r into v x , v y , and v z were done using library astropy.coordinates. Then, v z error was calculated respecting error propagation principle using these following equations, where In the equations above, µ b is proper motion in the galactic latitude direction, v t,b is tangential velocity in the galactic latitude direction, δ G = 27.12825 • is the declination of the north galactic pole, and α G = 192.85948 • is the right ascension of the North Galactic Pole [15]. The tracer vertical velocity dispersion, σ data (z), for every bin was calculated using equation wherev is the average vertical velocity for each bin. Velocity dispersion error for each was calculated following steps described by Pryor & Meylan (1993) [16]. Data processing was then continued by running the MCMC program, using velocity dispersion, velocity dispersion error, tracer height-scale, and tracer number density previously obtained (see sub section 3.2). The solution to Equation 18 was determined using MCMC program. Local dark matter density (ρ dm ), stellar surface density (Σ * ), stellar scale-height (h * ), and tracer scale-height (h t ) were set as free parameters, whose values were acquired based on the likelihood distribution within a certain range of priors.

Priors for the Free Parameters
The range of priors for the MCMC method were selected based on preceding knowledge available from literature. The prior for h * refers to the typical scale-height value of the thin disk, where a range of 200 pc to 400 pc was used. The typical scale-height value of the thin disk was taken as a benchmark, because the height range considered in this study was between 0 pc to 200 pc from the Galactic plane and within this range the stellar thin disks should still dominate stellar the distribution [17]. The range allowed for the Σ * sample retrieval was between 20 M /pc 2 to 50 M /pc 2 , refering to the value obtained by Xia et al. (2016) [7] and Zhang et al. (2013) [5], which is 42 ± 5 M /pc 2 dan 39.1 ± 7.1 M /pc 2 . The prior for ρ dm were taken based on previous works such as by Garbari  (2018) [8]. We decided to give it a range of 0 M /pc 3 to 0.1 M /pc 3 . Even though the value of h t has been obtained from number density fitting, it was still set as a free parameter, albeit with a relatively narrower prior value span. This was done to accommodate errors from the fitting results. The prior range used was 422.25 pc -432.25 pc. The value of ρ dm obtained in this study is 0.0116 ± 0.0012 M /pc 3 or ρ dm = 0.439 ± 0.046 GeV/cm 3 . As shown by PDFs of the free parameters on Figure 1, Σ * and h * are not tightly constrained, indicating that the model we used is insensitive to both quantities. This is mainly because the terms from the main equation (Equation 18) containing Σ * and h * are more complex than the linear term containing ρ dm .

Result and Comparison with other Determinations
A number of local dark matter density determinations carried out by several researchers in the past few years are listed in Table 1 and plotted in Figure 2. As shown in Table 1, the value we obtained is smaller than the results obtained by Garbari [9], who also used Gaia DR2 to determine the local dark matter density (but using different tracers from this study), also obtained higher results than is obtained in this study. But within the error ranges, the result of this study agrees with that of Garbari et al. (2012) [4], Xia et al. (2016) [7], Sivertsson et al. (2018) [8], and Buch et al. (2018) [9] for the G-dwarfs tracer.
The difference in the value of ρ dm we obtained with the value of ρ dm in the aforementioned studies can be attributed to several factors as follows: The selection function of each observation can also be different. Selection function can be related to the sensitivity of the instrument used. If the instrument sensitivity differs, it is possible that the objects observed are also different. Generally, instruments with higher sensitivity can observe more objects.
Missions with different locations can also generate different data because the location determines the area of the sky that can be observed. RAVE missions, for example, The differences in the shape of the volume and the height considered cause the total mass in the volume being considered to also be different, and it can make the value of ρ dm obtained to be different as well. On the other hand, all of these determinations assumed a constant value of ρ dm in the considered volume. As ρ dm is in fact not homogeneously distributed [1], the larger the space being reviewed, the less accurate the assumption of a constant ρ dm is.  [6] did not make the stellar disk scale-height, h * , as a free parameter but rather made it as a constant value at 300 pc.

Different applied assumptions
The differences in these assumptions can result in different applied main equation, or different quantities used as free parameters, thus causing differences in the results of ρ dm obtained.

Conclusion
We examined kinematics and positions of 95,543 K-dwarfs from Gaia DR2 to determine the local dark matter density. The value we acquire is 0.0116 ± 0.0012 M /pc 3 or ρ dm = 0.439 ±  [3] who get values of ρ dm 0.3 − 0.4 GeV/cm 3 and ρ dm = 0.385 ± 0.027 GeV/cm 3 , respectively. This gives us a hint that the shape of our Galaxy's dark halo is oblate.