The Rotation Curve, Mass Distribution, and Dark Matter Content of the Milky Way from Classical Cepheids

With the increasing number of large stellar survey projects, the quality and quantity of excellent tracers for studying the Milky Way are rapidly growing, one of which is the classical Cepheids. Classical Cepheids are high-precision standard candles with very low typical uncertainties (<3%) available via the mid-infrared period–luminosity relation. About 3500 classical Cepheids identified from the Optical Gravitational Lensing Experiment, All-Sky Automated Survey for Supernova, Gaia, Wide-field Infrared Survey Explorer, and Zwicky Transient Facility survey data have been analyzed in this work, and their spatial distributions show a clear signature of Galactic warp. Two kinematical methods are adopted to measure the Galactic rotation curve (RC) in the Galactocentric distance range of . Gently declining RCs are derived by both the proper motion (PM) method and three-dimensional velocity vector (3DV) method. The largest sample of classical Cepheids with the most accurate 6D phase-space coordinates available to date are modeled in the 3DV method, and the resulting RC is found to decline at the relatively smaller gradient of (−1.33 ± 0.1) . Comparing to results from the PM method, a higher rotation velocity ((232.5 ± 0.83) ) is derived at the position of the Sun in the 3DV method. The virial mass and local dark matter density are estimated from the 3DV method, which is the more reliable method, and GeV, respectively.


Introduction
The mass distribution and dark matter density profiles of the Milky Way are not just key probes of its assembly history (e.g., Lake 1989;Read et al. 2008;Deason, Belokurov & Sanders 2019), but also provide crucial clues for the cosmological context of galaxy formation (e.g., Dubinski 1994;Ibata et al. 2001;Lux et al. 2012). The two distributions are usually studied in the frame work of the 'standard' Cold Dark Matter model (ΛCDM for short, where the Λ refers to the density of 'dark energy'). In this cosmological model, the energy density of the Universe comprises approximately 5% of baryons, 27% of dark matter and 68% of dark energy. The rotation (or circular velocity) curve measurement is a classical way to deliver an indirect measurement of these profiles of the Milky Way (Volders 1959;Freeman 1970;Bosma & van der Kruit 1979;van Albada et al. 1985;Sofue et al. 2009).
Specifically, the Galactic rotation curve (RC) is the mean circular velocity around the center of the Galaxy as a function of galactocentric distance measured in the disk-mid plane. The RC is has been derived with various methods and various tracer objects moving in the gravitational potential of the Galaxy (e.g., Wilkinson & Evans 1999;Weber & de Boer 2010;Sofue 2012;Nesti & Salucci 2013). For example, the RC of the Galactic inner region has been derived by the tangent-point method associated with H I regions (Gunn et al. 1979;Levine et al. 2008;Sofue et al. 2009). Comparing to the tangent-point method, methods using stars, dwarf galaxies or globular clusters with distances and at least one of the velocity components (radial velocity and/or proper motions) are considered as the better measurement for the Galactic (inner and outer regions) RC (e.g., Smith et al. 2007;Honma et al. 2007;Bovy et al. 2012;Bovy & Rix 2013;Bhattacharjee et al.2014;Kafle et al 2014;Reid et al. 2014;Bowden et al. 2015;Binney & Wong 2017;Pato & Icco 2017;Russeil et al. 2017;Ablimit & Zhao 2017;Katz et al. 2018;Monari et al. 2018;Sohn et al. et al. 2018). Recently, the measured number of tracers with accurate 6 dimensional (D) phase-space information is increasing rapidly, with the growing numbers of sky surveys, such as SDSS, Gaia, WISE, ZTF, OGLE, ASAS, Gaia-ESO, APOGEE, etc., and these data enable us to measure more precise RC. matics and dynamics of the Galaxy, such as RR lyrae stars (Ablimit & Zhao 2017;Medina et al. 2018;Utkin et al. 2018;Wegg et al. 2018) and Cepheids (e.g., Kawata et al. 2018). Frink, Fuchs & Wielen (1995) derived the Galactic rotation curve from proper motions of 144 Cepheids. Subsequently, Pont et al. (1997) constructed the RC of the Galaxy from radial velocities of 48 classical Cepheids distributed in the outer disc region between the Galactocentric distance 10 kpc and 15 kpc. Gnaciński (2019) obtained the RC by adopting three kinematic approaches by using 160, 228 and 120 classical Cepheids from the catalog of Mel'nik et al.(2015). They showed that the slope of the RC lies between a flat RC and a Keplerian rotation curve. However, Mróz et al. (2019) tracked the RC from the 6D phase-space information of 773 classical Cepheids, and they found a relatively flat rotation curve. They did not estimate mass distribution and dark matter content of the Milky Way.
In this work, we have selected and analyzed about 3500 classical Cepheids which have precise distances and measured the Milky way rotation curve using the proper motion method (Gnaciński 2019) and 3D velocity vector method (Reid et al. 2009). In §2, we introduce the classical Cepheids data collected for this work. Two methods to calculate the rotation velocities of classical Cepheids are introduced, and the resulted rotation curve & its constraint on the mass and dark matter profile of our Galaxy are given and discussed in §3. The concluding remarks are presented in §4.

Data Selection
We collected our sample from several classical Cepheid catalogs as follows: the All-Sky Automated Survey for Supernovae (ASAS-SN) Variable stars catalog (Shappee et al. 2014;Jayasinghe et al. 2018), the classical Cepheids sample by Skowron et al. (2019a, b) basically from the Optical Gravitational Lensing Experiment (OGLE) (Udalski et al. 2015;Udalski et al. 2018), classical Cepheids from the European Space Agency (ESA) mission Gaia (Gaia Collaboration 2016, 2018Ripepi et al. 2019), and the classical Cepheids catalog by  from the Wide-field Infrared Survey Explorer (WISE) (Wright et al. 2010). We added new classical Cepheids identified from the Zwicky Transient Facility (ZTF) catalog (Bellm et al. 2019) by Chen et al.(2020). We made a cross-match of all the Cepheids from different catalogs in order to remove multiple entries. Then, we selected Cepheids which have mid-infrared (W 1, W 2, W 3 and W 4 bands) magnitudes from the All WISE catalogue. We calculated heliocentric distances (D h ) based on the relations given in Wang et al. (2018) with W 1, W 2, W 3 and W 4 bands, and took average values for each Cepheid (also see Skowron et al. (2019a) for same calculation method). Recently, it has been discussed that distances derived from mid-infrared period-luminosity relations are more accurate than distances obtained from parallaxes (e.g., Mróz et al. 2019). After deriving distances, we keep classical Cepheids with |z| ≤ 4 kpc, and we have 3483 classical cepheids (Galactic longitude (l) and latitude (b) distributions are shown in the upper left panel of The spatial distributions are shown in the Figure 1, and all distributions show the clear Galactic warp which is reported by Skowron et al. (2019a, b) and . The 3D positions of Cepheids and galactocentric distances (r) in the Cartesian coordinate system are calculated as x = R − D h cos b cos l, y = D h cos b sin l, z = D h sin b and r = x 2 + y 2 + z 2 , where R is the distance from the Sun to the Galactic center, and the recent most accurate value, 8.122 ± 0.031 kpc (GRAVITY collaboration et al. 2018), is adopted. The projection of galactocentric distance on the Galactic plane (R) is as follows, (1) 3. Modeling the rotation curve

The halo model
The rotation velocity at a radius R from the center of an axisymmetric mass distribution is related to the total gravitational potential within R and mass M (< R) (at z ∼ 0), where Φ and G are the gravitational potential and gravitational constant, respectively. If we consider the bulge, thin disk, thick disk and dark matter halo for the Galactic potential, which for the respective contributions are Φ bulge , Φ thin , Φ thick and Φ halo , and velocity contributions to the RC from different components are given by, The Navarro-Frenk-White (NFW) model (Navarro et al. 1996(Navarro et al. , 1997 which is derived from the simulations in the CMD scenario of galaxy formation has been widely used for modeling the dark matter halo (e.g., Sofue 2012; Wang et al. 2018). We assume that density distributions of all stellar components are well-known, and the velocity contribution of the dark matter halo is fitted by searching for the best-parameters by using the Markov Chain Monte Carlo method. For the fitting model, the Miyamoto-Nagai potential model (Miyamoto & Nagai 1975) and a spherical Plummer potential (Plummer 1911) are used for the thin/thick disks and bulge, respectively. We take the parameter values of the enclosed mass, the scale length, and the scale height from the model I of Pouliasis et al.(2017).
The NFW dark matter density profile is described as, where ρ crit = 3H 2 /8πG, and H = 70.6 km s −1 Mpc −1 is taken for the Hubble constant. The quantity of δ c is the characteristic overdensity of the halo. Here, r s = R vir /c is the scale radius, where c is so-called concentration parameter, and R vir is the virial radius. R vir is related to the virial mass as M vir = 200ρ crit 4π 3 R 3 vir (see Navarro et al. (1996Navarro et al. ( , 1997 for more details). In the next subsections, the rotation curves from different kinematical models and fitting results are discussed.

The rotation curve from proper motions
After measuring the proper motion of the star and setting the solar rotation speed as V c, = 233.6 ± 2.8 km s −1 (Mróz et al. 2019), then assuming a circular orbit for the Cepheid, the following formula gives the rotation velocity (Gnaciński 2019), where the transverse velocity V t = Dµ l , and µ l is the proper motion in the Galactic longitude (multiplied by cosb). The stars with |z| > 0.5 kpc are excluded, and unphysical velocities caused by small or negative denominators are removed (see Gnaciński (2019) for the same selection criterion), so only 591 classical Cepheids are left from whole classical Cepheids for this kinematical modeling. Among our sample, there are 168, 324 and 411 classical Cepheids distributed in the Galactocentric range of R > 12 kpc, R > 10 kpc and R > 8 kpc, respectively. Figure 2 shows µ l and the calculated rotation velocities of 591 classical Cepheids. More than 98% of µ l have uncertainties less than 0.2 mas yr −1 , and this leads to small uncertainties in the rotation velocity calculation. The number of analyzed classical Cepheids in this work is about twice that used in Gnaciński (2019), and we have more stars in the outer disc which is helpful to improve the accuracy of the RC measurement. Figure 3 shows the rotation velocity distribution from R = 4 kpc and 19 kpc (see Table  1 for the values), and the linear function fitted from it is, This yields a gently declining rotation curve with a small gradient of (−1.45 ± 0.16) km s −1 kpc −1 , and indicates the rotation velocity at the position of the Sun as V c (R ) = 222.91 ± 2.08 km s −1 . By fixing the contributions of baryonic components of the Galaxy (see model I of Pouliasis et al.(2017)), we estimated the mass and the properties of the Milky Way's dark matter halo with the NFW profile (fitted results in Figure 3), and we derived M vir = (6.63 ± 0.67) × 10 11 M , corresponding to a viral radius R vir = 178.57 ± 5.42 kpc. We obtained the concentration of c = 12.36 ± 0.42 and a scale radius of r s = 14.45 ± 0.46 kpc. The indicated characteristic density is ρ 0 = (1.05 ± 0.12) × 10 7 M kpc −3 , and dark matter density at the location of the Sun is ρ DM, = 0.28 ± 0.04 GeV cm −3 .

The rotation curve from 3D velocity vector
The rotation velocity can also be determined from the 3D velocity vector if the three quantities of radial velocity and proper motions are available. Reid et al. (2009) described the calculation formulas of stellar motions by using radial velocity (V h ) and proper motions, which we adopt here: U -velocity component toward the Galactic center, V -velocity component along with the Galactic rotation, W -toward the North Galactic pole. The optimizing , where V c, and dVc dR are the Sun's rotation velocity and fitted parameter, is adopted for deriving rotation velocities (see Reid et al. (2009) for more details). For the peculiar (noncircular) solar motions with respect to the local standard of rest, the values of U = 11.1 ± 1.3 km s −1 , V = 12.24 ± 2.1 km s −1 and W = 7.3 ± 0.7 km s −1 are taken from Schönrich et al. (2010).
The proper motions of the sample are obtained from the Gaia DR2, and the radial velocities are derived by cross-matching with Gaia DR2 and LAMOST DR6 data (e.g., Zhao et al. 2006Zhao et al. , 2012. We excluded five Cepheids known in the binary systems, and we put extra constraints of |z| ≤ 2.0 kpc and radial velocity uncertainty < 20 km s −1 to remove 11 objects in order to reduce uncertainties. It is well known that the radial velocity uncertainty may be larger for a single star when it's measured near the pulsation phase (Stibbs 1955). However, the uncertainties of variable stars caused by the pulsation need further investigations, and it may not clearly affect the statistical result (see Ablimit & Zhao 2017). For the 3D  Cleaned Sample. There are likely some objects in 1078 star sample, which may be members of binary systems (and unrecognized with incorrect astrometric solutions) or categorized erroneously as classical Cepheids which as such and may actually be another type of variable. There are also some classical Cepheids with observed velocity components about 4σ (σ is dispersion of residuals) larger than the mean. Considering these possibilities and uncertainties, we selected 963 classical Cepheids from 1078 stars as the cleaned sample, and derived rotation velocities of the cleaned sample are shown in Table 1. The measured RCs from the cleaned sample and all 1078 sample can be fitted by the same linear function.
The rotation curve from this method is gently decreasing with a derivative of (−1.33 ± 0.1) km s −1 kpc −1 . The slope of the curve and rotation velocity at the location of the Sun (V c (R ) = 232.5±0.83 km s −1 ) are in a good agreement with the results of Mróz et al. (2019), as about 70% of data points in the sample overlap with that of Mróz et al. (2019). However, there are more than 300 different stars in this work. In particular, our sample has more stars in the outer disc, which improves the accuracy of the RC, and is helpful to put more accurate constraint on the distribution of dark matter in the Milky Way. Comparing to the virial mass from the proper motion method, we derived a higher viral mass in this method, M vir = (8.22 ± 0.52) × 10 11 M with a corresponding viral radius R vir = 191.84 ± 4.12 kpc. The resulted concentration and scale radius are c = 13.04 ± 0.34 and r s = 14.71 ± 0.42 kpc, respectively. The estimated characteristic density and dark matter density at the location of the Sun are ρ 0 = (1.20±0.1)×10 7 M kpc −3 and ρ DM, = 0.33±0.03 GeV cm −3 , respectively.

Comparison and discussion
There are 366 common classical Cepheids modeled in two methods, and different tracers are selected due to different criterion for different methods. The discrepancy of two methods' results are within 10%. The most important advantage of our sample is the accuracy of the distances which have uncertainties at a level of 2-3%. We have small uncertainties in our results (see values of uncertainties in Table 1), however, only bootstrapping uncertainties without the systematic uncertainties are considered in this work (see Eilers et al. (2019) for analysis of the possible systematic uncertainties). The effect of the asymmetric drift is not considered in the calculation of this work due to the very small systematic uncertainty it causes (e.g., estimated as ±0.28 km s −1 by Kawata et al. 2018). Within 19 kpc, all systematic uncertainties added up (i.e. caused by uncertainties of distances, uncertainty of R and the asymmetric drift etc.) only affect the RC measurement at a 5% level. It is well known that the motions of stars are affected by Galactic substructures (e.g. Grand, Kawata & Cropper 2014;Bovy et al.2015;Kawata et al. 2018;Martinez-Medina et al. 2019). We did not use stars located at R < 4.0 kpc in order to reduce the influence of other structures like the Galactic bar.
The slopes of the rotation curves from two methods are gently decreasing, as favored by the recent discoveries (e.g., Mróz et al. 2019;Eilers et al. 2019). They are not as flat as demonstrated in Sofue et al. (2009) and Reid et al. (2014), and it is not as steep as showed in Gnaciński (2019). This indicates that the dark matter content would not possibly so high or so low claimed in those previous works. The result (see the cross-point between the rotation curve of all stellar components and dark matter halo in Figure 3) from the proper motion method suggests that the dark matter halo dominates the Galactic rotation when R 14.5 kpc, and this is in good agreement with recent finding by Eilers et al. (2019). However, based on the 3D velocity method (as shown in Figure 5), the dark matter halo dominates the rotation velocity if R 12.5 kpc. The comparison of two velocity distributions from two methods give a same dip-like feature, there are a clear declining at the distance around ∼10 kpc, and this is consistent with the similar dip claimed by previous works (Sofue et al 2009;Kafle et al. 2012;McGauph 2018). However, there is no dip in the results of the cleaned sample with the 3D velocity method (Eilers et al. 2019).
Effect of uncertainties of baryonic mass components. The estimation of the dark matter halo profiles rely on observational results of the baryonic mass components. Recently, de Salas et al.(2019) discussed that the dark matter density estimation is more sensitive to the uncertainties of the baryonic components rather than the uncertainties of the rotation velocities. They found a different uncertainty (±0.149 GeV cm −3 ) of the dark matter density with the same velocities, and it is about 3 times of what Eilers et al. (2019) find. They also show that using different model such as the NFW dark matter profile and Einasto dark matter profile also gives a uncertainty of ±0.036 GeV cm −3 . Comparing to the Galactic disc mass in some previous works (e.g., Smith et al. 2007), we took relatively higher masse for the Galactic (thin + thick) disc from Pouliasis et al.(2017). Thus, we may underestimate the halo profiles. We examine it by taking a very simple example, and we run the model by using 5.0 × 10 10 M for the whole disc mass instead of 7.888 × 10 10 M (thin + thick discs) as in this work from Model I of Pouliasis et al.(2017). We found that the dark matter density goes up to 0.408 GeV cm −3 when we reduce the baryonic mass of Galactic disc in the estimation modeling, and it is 0.078 GeV cm −3 higher than the value (0.33 GeV cm −3 ) derived from Model I of Pouliasis et al.(2017). This supports the statement given by de Salas et al.(2019). The future observational data may provide better constraints on the baryonic components 2 .
Interestingly, the density derived from our 3D velocity method is basically consistent with the estimated dark matter density by de Salas et al.(2019) which is in a range of (0.3 − 0.4) GeV cm −3 . Our local dark matter density estimated from 3D velocity method is in very good agreement with the local dark matter density (0.32 − 0.36 GeV cm −3 ) inferred from fitting models to the Gaia DR2 Galactic rotation curve and other data .

Conclusion
We have analyzed 3483 classical Cepheids selected from thousands of classical Cepheids identified by several survey projects (e.g., OGLE, ASAS-SN, Gaia, WISE and ZTF), and constructed the rotation velocity distribution of the Milky Way between the Galactocentric distance 4 kpc and 19 kpc by using two different methods. The distances of these classical Cepheids have the typical uncertainties of < 3% (which is crucial in the analysis of the rotation curve), and 3D spatial distributions show a vary clear Galactic warp feature claimed by previous works (see the section 2). 591 and 1078 classical Cepheids have been analyzed by using the proper motion and 3D velocity methods, and most of observed uncertainties of proper motions and radial velocities are less than 0.2 mas yr −1 and 20 km s −1 , respectively. This represents the largest classical Cepheid sample analyzed to date. We apply the NFW profile approach to simulate the dark matter content of the Milky Way. Our main findings are, 1. The different methods or/and different sample would give different results in some extent. The uncertainties of baryonic components also have important role in the estimation of dark matter profiles. The result of the proper motion method shows that the dark matter halo is main contributor to the Galactic rotation when the distance R 14.5 kpc, while the 3D velocity modeling demonstrates that the Galactic rotation curve is dominated by the dark matter halo at R 12.5 kpc. The rotation curve constructed by both method are gently declining. The rotation curve from 3D velocity method is decreasing more gently with a derivative of (−1.33 ± 0.1) km s −1 kpc −1 . The rotation velocity at the position of the Sun ((232.5 ± 0.83) km s −1 ) obtained from the 3D velocity method is about 10 km s −1 faster than the rotation velocity of the Sun derived from the proper motion method.
2. The best-estimation with the NFW profile based on the rotation curve of the 3D velocity method generates a higher viral mass (M vir = (0.822 ± 0.052) × 10 12 M ) with the corresponding radius of R vir = 191.84 ± 4.12 kpc and concentration of c = 13.04 ± 0.34. At the same time, the predicted local dark matter density (ρ DM, = 0.33 ± 0.03 GeV cm −3 ) is also higher than the estimated value from the proper motion modeling.
We are grateful to Anna-Christina Eilers for useful discussions and comments to improve the manuscript. We thank Xiaodian Chen for providing the new classical Cepheids catalog identified from the ZTF data. This work was supported by National Natural Science This publication made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/ gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/ dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This work has used the data products from the Wide field Infrared Survey Explorer (WISE), which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. The work also have used the data from the Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST) which is a National Major Scientific Project built by the Chinese Academy of Sciences. Funding for the project has been provided by the National Development and Reform Commission. LAMOST is operated and managed by the National Astronomical Observatories, Chinese Academy of Sciences.