Rotation Curve of the Milky Way and the Dark Matter Density

We review the~current status of the~study of rotation curve (RC) of the~Milky Way, and~present a~unified RC from the~Galactic Center to the galacto-centric distance of about 100 kpc. The~RC is used to directly calculate the~distribution of the~surface mass density (SMD). We then propose a~method to derive the~distribution of dark matter (DM) density in the~in the~Milky Way using the~SMD distribution. The~best-fit dark halo profile yielded a local DM density of $\rho_\odot = 0.36\pm 0.02$ GeV/cc. We also review the~estimations of the~local DM density in the~last decade, and~show that the~value is converging to a~value at $\rho_\odot=0.39\pm 0.09$ GeV/cc.


Introduction
The rotation curve (RC) of the Milky Way was obtained by observations of galactic objects in the non-MOND (MOdified Newtonian Dynamics)frame work. The existence of the dark halo (DH) has been confirmed by the analysis of the observed RCs, assuming that Newtonian dynamics applies evenly to the result of the observations. In this article, current works of RC observations are briefly reviewed, and a new estimation of the local dark matter (DM) density is presented in the framework of Newtonian dynamics.
An RC is defined as the mean circular velocity V rot around the nucleus plotted as a function of the galacto-centric radius R. Non-circular streaming motion due to the triaxial mass distribution in a bar is crucial for kinematics in the innermost region, though it does not affect the mass determination much in the disk and halo. Spiral arms are another cause for local streaming, which affect the mass determination by several percent, while they do not influence the mass determination of the dark halo much.
There are several reviews on RCs and mass determination of galaxies [1 ; 2 ; 3 ].
In this review, we revisit recent RC studies and determination of the local DM density in our Milky Way. In Section 2, we briefly review the current status of the RC determinations along with the methods. In Sections 3 and 4 we propose a new method to use the surface mass density (SMD) directly calculated from a unified RC to estimate the local DM density, and apply it to the newest RC of the Milky Way up to radius of ∼ 100 kpc. We adopt the galactic constants: (R 0 , V 0 )=(8.0 kpc, 238 km s −1 ) [4 ; 5 ], where R 0 is the distance of the Sun from the galactic center (GC) and V 0 is the circular velocity of the local standard of rest (LSR) at the Sun [6 ].

Progress in the Last Decades
The galactic RC is dependent on the galactic constants. Accordingly, the uncertainty and error in the RC include uncertainties of the constants. Currently recommended, determined, or measured values are summarized in Table 1, where they appear to be converging to around ∼ 8.0 − 8.3 and ∼ 240 km s −1 . In this paper, we adopt R 0 = 8.0 kpc and V 0 = 238 km s −1 from the recent measurements with VERA (VLBI Experiments for Radio Astrometry) [4 ; 5 ].
The RC of the galaxy has been obtained by various methods as described in the next

Methods to Determine the Galactic RC
The particular location of the Sun inside the Milky Way makes it difficult to measure the rotation velocity of the galactic objects. Sophisticated methods have been developed to solve this problem, as briefly described below.

Tangent-Velocity Method
Inside the solar circle (−90 • ≤ l ≤ 90 • ), the galactic gas disk has tangential points, at which the rotation velocity is parallel to the line of sight and attains the maximum radial velocity v r max (terminal or tangent-point velocity). The rotation velocity V (R) at galacto-centric distance R = R 0 sin l is calculated simply correcting for the solar motion.

Radial-Velocity + Distance Method
If the distance r of the object is measured by spectroscopic and/or trigonometric observations, the rotation velocity is obtained by geometric conversion of the radial velocity, distance, and the longitude. The distance has to be measured independently, often using spectroscopic distances of OB stars, and the distances are assumed to be the same as those of associated molecular clouds and HII (ionized hydrogen) regions, whose radial velocities are observed by radio lines. Since the photometric distances have often large errors, obtained RC plots show large scatter.

Trigonometric Method
If the proper motion and radial velocity along with the distance are measured at the same time, or from different observations, the 3D velocity vector, and therefore the rotation velocity, of any source is uniquely determined without being biased by assumption of circular motion as well as the galactic constants. VLBI (very long baseline interferometer) measurements of maser sources [42 ; 4 ; 5 ; 44 ] and optical/IR trigonometry of stars [46 ; 20 ] have given the most accurate RC.

Disk-Thickness Method
The errors in the above methods are mainly caused by the uncertainty of the distance measurements. This disadvantage is eased by the HI-disk thickness method [39 ; 40 ]. The angular thickness of the HI disk along an annulus ring is related to can be used to determine the rotation velocity by combining with radial velocity distribution along the longitude.

Pseudo-RC from Non-Disk Objects
Beyond or outside the galactic disk, globular clusters and satellite galaxies are used to estimate the pseudo-circular velocity from their radial velocities based on the Virial theorem, assuming that their motions are at random, or the rotation velocity is calculated by where v g is the galacto-centric radial velocity. On the other hand, Huang et al. (2016) [29 ] have recently employed more sophisticated, probably more reliable, method to solve the Jeans equations for the non-disk stars and clusters.

Unified RC
A RC covering a wide region of the galaxy has been obtained by compiling the existing data by re-scaling the distances and velocities to the common galactic constants representing the central massive black hole. Outer RC beyond R ∼ 30 kpc has been determined from the radial motions of satellite galaxies and globular clusters.
The RC determination has been improved recently by compiling a large amount of data from a variety of spectroscopic as well as trigonometric measurements from radio to optical wavelengths. An extensive compilation of the data of rotation velocities of the galactic disk has been published recently, and is available as an internet data base [25 ; 26 ; 47 ; 48 ; 27 ; 28 ]. ; 30 ]re-scaled to the galactic constants of (8.0 kpc, 238 km s −1 ) following the method described in [34 ].Although individual data points are largely scattered, their averages well coincide with the unified RC. In the figures we also compare the data with the RC by [29 ]up to ∼ 100 kpc without rescaling, which also coincides with the other data within the scatter.
We here comment on the property of the unified RC built by averaging the published data. It must be remembered that the averaging procedure does not satisfy the condition of statistics in the strict meaning, because the data are compiled from different authors using a variety of instruments and analysis methods, which makes it difficult to evaluate common statistical weights for the used data points. So, remembering such a property, in view that the unified RC well approximates the original curves as well as for its convenience for the determination of the mass distribution by the least-squares and/or χ 2 fitting, we shall employ it in our present analysis.

Mass Components
The rotation velocity is related to the gravitational potential, hence to the mass distribution, where Φ i is the gravitational potential of the i-th component and V i is the corresponding circular velocity. The rotation velocity is often represented by superposition of the central black hole (BH), bulge, disk, and the dark halo as Here, the subscript BH represents black hole, b stands for bulge, d for disk, and h for the dark halo. The contribution from the black hole can be neglected in sufficiently high accuracy, when the dark halo is concerned. The mass components are usually assumed to have the following functional forms.

Massive Black Hole
The GC of the Milky Way is known to nest a massive black hole of mass of M BH ∼ 4×10 6 M ⊙ [17 ; 9 ; 10 ].The RC is assumed to be expressed by a curve following the Newtonian potential of a point mass at the nucleus.

De Vaucouleurs Bulge
The commonly used SMD profile to represent the central bulge, which is assumed to be proportional to the empirical optical profile of the surface brightness, is the de Vaucouleurs law [49 ], where Σ be is the value at radius R b enclosing a half of the integrated surface mass [2 ].Note that the de Vaucouleurs surface profile, also the exponential disk, has a finite value at the center. The volume mass density ρ(r) at radius r for a spherical bulge is calculated using the SMD by and the mass inside R is 8 The circular velocity is thus obtained by More general form e −(R/re) n called the the Sérsic law is discussed in relation to its dynamical relation to the galactic structure based on the more general profile [50 ; 51 ].

Exponential Disk
The galactic disk is generally represented by an exponential disk [52 ],where the SMD is expressed as Here, Σ d is the central value, R d is the scale radius. The total mass of the exponential disk is given by M disk = 2πΣ dc R 2 d . The RC for a thin exponential disk is expressed by [53 ] V where y = R/(2R d ), and I i and K i are the modified Bessel functions.
The dark halo is described in the next section

Dark Halo
The existence of dark halos in spiral galaxies has been firmly evidenced from the well estab- indicates that the isothermal model can be ruled out in representing the Milky Way's halo.

Dark Halo Models
There have been various proposed DH models, which may be categorized into two types: Beta model with β = 1 [58 ] : Burkert model [55 ; 56 ] : Brownstein model [57 ] : On the other hand, the central cusp models are often represented by the following functions. NFW model [62 ; 59 ] : Moore model [60 ; 61 ] with α = 1.5: Figure 2 shows schematic density profiles for various DH models with h = 10 kpc combined with the de Vaucouleurs bulge and exponential disk, where the halo density is normalized at R = 20 kpc..

Cusp vs Cored Halo
The density profiles for the NFW (Navarro, Frenk and White), Moore, Burkert, β, and Brownstein models are almost identical beyond the core radius h, where they tend to ∝ R −3 . Differences among the models appear within the Solar circle. The cusp models (NFW and Moore models) predict steep increase of density toward the center with a singularity. The cored halo models predict a mild and low density plateau in the center with the peak densities not much differing from each other within a factor of two. However, the Burkert model has a singularity with the density gradient being not continuous across the nucleus.
Most of the DH models predict lower density in the innermost galaxy by two to several orders of magnitudes than the bulge's density. This implies that the DH does not much influence the kinematics in the inner galaxy. Namely, it is practically impossible to detect the DM cusp by analyzing the RC. Only the Moore model predicts cusp density exceeding the bulge's density in the very center at R <∼ 0.1 pc, whereas the applicability of the model to such small sized region is not obvious [61 ].

Central DM Density
If we assume that the functional form of the NFW model is valid in the very central region, the SMD at R ∼ 100 pc could be estimated to be about Σ ∼ 2. Interestingly, the column density of DM, hence brightness (flux/steradian) of selfannihilation emission (γ-ray) stays almost constant against the radius and is therefore constant regardless the resolution of the detector. On the other hand, the emission measure ∼ ρ 2 R varies as ∝ R −1 , hence, the brightness of collision-origin emission (γ or microwave haze) increases toward the center [e.g., 65 ], so that the detection rate will increase with the detector's resolution.
Another concern about the DM cusp is the kinetic energy of individual particles. In order for the cusp to be stationary, the particles must be bound to the gravitational potential, so that the particle's speed must be lower than the escaping velocity v ∼ √ 2V rot ∼ 300 km s −1 . This will give a constraint on the cross section σ A of the DM annihilation, if the collision rate σ A v is fixed by the detection of DM-origin emissions.
The cored halo models (isothermal, Burkert, Brownstein, and the β models) predict a mild and finite-density plateau with scale radius of h (∼ 10 kpc). Their central densities are also several orders of magnitude less than the bulge's density, hence do not contribute to the kinematics of the gas and stars in the GC.
Finally, it is emphasized that the discussion of a central DM cusp has been obtained based on the theoretical predictions, but neither on any observed phenomenon nor on measured quantity from the RC or kinematics of GC objects. This makes a strong contrast to the observed fact of the dark halo in the outer galaxy, firmly evidenced by the analyses of the RC.

SMD from RC
In the decomposition method of the RC, the resulting mass distribution depend on the assumed functional forms of the model profiles. In order to avoid this inconvenience, the RC can be used to directly calculate the surface mass distribution without employing any functional form. Only an assumption has to be made, either if the galaxy's shape is a sphere or a flat disk.
On the assumption of spherical distribution, the mass inside radius R is given by Then the surface-mass density (SMD) Σ S (R) at R is calculated by where If the galaxy is assumed to be a flat thin disk, the SMD Σ d (R) is calculated by solving Poisson's equation (Freeman 1970; Binney and Tremaine 1987) by Here, K is the complete elliptic integral, which becomes very large when x ≃ R.
The SMD distributions in the galaxy for the sphere and flat-disk cases have been calculated for the recent RCs [2 ].In this paper we apply the same method to the here obtained unified RC (Figure 1). Since we aim at studying the dark halo, which is postulated to be rather spherical than a flat disk, we assume spherical mass distribution. The calculated SMD distribution is shown in Figure 3.
The SMD is strongly concentrated toward the center, reaching a value as high as ∼ 10 5 M ⊙ pc −2 within R ∼ 10 pc, representing the core of the central bulge with the extent of several hundred pc. It is followed by a straightly declining profile from R ∼ 2 to 8 kpc in the semi-logarithmic plot, representing the exponential nature of the galactic disk. In the outer galaxy beyond ∼ 8 kpc, the SMD profile tends to be displaced from the straight disk profile, and is followed by an extended outskirt with a slowly declining profile, representing a massive halo extending to the end of the RC measurement at ∼ 100 kpc.

Fitting by Bulge, Disk, and Dark Halo
In order to separate the dark halo from the disk and bulge components, the well established RC decomposition method has been extensively applied to the RCs [2 ; 3 ].Besides this traditional method, we here propose to use the SMD distribution. For this, we assume three mass components of de Vaucouleurs bulge, exponential disk, and dark halo. In order to represent the, we employ the NFW profile as a 'tool' for its popularity and for the dynamics background based on the extensive numerical simulations.
We employ the least χ 2 fitting method, where χ 2 is defined by with i denoting the value at the i-th data point, and σ i is the standard deviation around each data point in the running averaging procedure of the SMD distribution.
Fitting parameters are the scale radius a d and central SMD Σ 0 d for the disk, and the scale (core) radius h and representative DM density ρ 0 model for the halo. The bulge SMD is fixed to an assumed de Vaucouleurs profile, which is negligible in the present fitting range at R ≥ 1 kpc.
The fitting was obtained between R = 1 and 100 kpc. The fitting result for the NFW halo model is shown in Figure 3. The solid line is the χ 2 fit to SMD, and red, blue, and dashed lines represent the halo, disk, and bulge components, respectively. Note that the semilogarithmic plot makes it visually easier to recognize the dark halo significantly displaced from the exponential disk, which appears as a straight line.
The present method to fit the SMD distribution is essentially the same as that to fit the RC in the sense that both the methods search for a set of four parameters (a d , Σ d , h, and ρ 0 model ) that minimize χ 2 either of SMD or of RC. However, an advantage in using SMD is to make it visually easier to discriminate the disk from the halo in the semi-logarithmic plots ( figure 3), where the disk appears as a straight line.
The best-fit parameters for the disk are determined to be a d = 4.38 ± 0.35 kpc and Σ 0 = (1.28 ± 0.09) × 10 3 M ⊙ pc −2 . Table 3 lists the fitted result along with the minimized χ 2 value.
We also obtained χ 2 fitting using the Burkert, Brownstein, and β profiles, and listed the local DM density and minimized χ 2 in Table 3. In these three models, the χ 2 ∼ 17 − 18 were found to be systematically greater than that for the NFW model (χ 2 = 11.9). The reason for the difference is due to the systematic difference in the functional behavior between NFW and the other three models: NSF has a cusp steeply increasing toward the center with sharpening scale radius, which results in the possibility of finer fitting to the slightly curved SMD profile at R <∼ 10 kpc in the semi-log plot. On the contrary, the other three models predict almost negligible SMD there, so that halo parameters contribute less intensively to the fitting in the innermost region, or the fitting must be done only by the disk's two parameters there, resulting in worse fitting.  The local DM density is a key quantity in laboratory experiments by the direct detection of DM, and has been estimated by a number of authors with a variety of methods.
In Table 4 we list the local DM densities from the literature along with the present value for NFW profile. They are also plotted in Figure 4 against publication years. The ρ ⊙ values seem to be nearly constant in the decade. Averaging all the listed values with an equal weighting yields ρ ⊙ = 0.39 ± 0.09 GeV cm −3 , which may be taken as a 'canonical' value.

Summary
We reviewed the current status of determination of the RC of the Milky Way, and presented a unified RC from the GC to outer halo at R ∼ 100 kpc. The RC was used to directly calculate the SMD without assuming any functional form. The disk appears as a straight line on the semi-logarithmic plot of SMD against R, and is visually well discriminated from the DH having an extended outskirt.
The SMD distribution was fitted by a bulge, disk, and NFW dark halo using the χ 2 method. The best-fit DH profile yielded the local DM density of 0.359 ± 0.017 GeV cm −3 .
We also reviewed the current estimations from the literature in the last decade, which appear to be converging to a mean value of ρ ⊙ = 0.39 ± 0.09 GeV cm −3 .    Figure  1.