Tight Constraint on Local Dark Matter Density from a Joint Analysis of Local Kinematics and Rotation Curves

Tight Constraint on Local Dark Matter Density from a Joint Analysis of Local Kinematics and Rotation Curves Sutirtha Mukherjee1,∗Aditya Singh2,†Viraj Karambelkar3,‡Sayan Mandal3,§and Subhabrata Majumdar3¶ 1Department of Physical Sciences (DPS), Indian Institute of Science Education and Research (IISER), Mohali, India 2Department of Astronomy and Astrophysics, Tata Institute of Fundamental Research, Mumbai 400005, India Cahill Center for Astrophysics, California Institute of Technology, Pasadena, CA 91125, USA and 3Department of Theoretical Physics, Tata Institute of Fundamental Research, Mumbai 400005, India ABSTRACT Dark Matter (DM) particles, by their definition, are very weakly interacting and cannot be seen. They make their presence felt primarily through gravity. In particular, in the Milky Way (MW), we estimate the distribution of DM (and its velocities) by using visible matter tracers (like luminous compact old stars, extended mass like globular clusters, pulsating variable stars etc.) to construct rotation curves and mass of Halo. These rotation curves are also dependent on absolute normalisation given by the Local Standard of Rest (LSR) whose values are again dependent on visible matter tracers. Moreover, the gravitational potential in which the tracers move depend on both the distributions of dark and visible matter which are difficult to disentangle. We show that choice of the LSR, the final rotation curve made out of binning and averaging of individual tracer velocities as well as any prior on the visible matter distribution crucially influences the estimate of the DM density. This results in the wide range of values for the local DM density (as diverse as 0.1−0.9GeVcm−3, but with small error bars) quoted in the literature. Precise knowledge of the distribution of ‘visible’ matter by the ongoing GAIA DR3 survey and the futuristic THEIA mission is, therefore, of paramount importance in our understanding of DM. We have also analyzed Mass distribution of our galaxy following NFW profile. This is especially important when using the astrophysical inferences about DM as inputs in particle DM experiments. I. INTRODUCTION One of the most popular candidates for dark matter particles is the hypothetical weakly interacting massive particle (WIMP). The determination of the density of these particles in the halo of the Milky Way galaxy in general and the solar neighborhood in particular is crucial for many direct detection experiments which measure the rate of nuclear recoil events when the WIMPs scatter off the nuclei of the detector materials, since this rate has a direct dependence on the local density. One approach of determining the densities is to use the rotation curve (RC) data to find the likelihood of the parameters characterizing the density distributions of the various mass components of the galaxy. In general, the visible matter (VM) parameters are fixed (from observational data) and the dark matter (DM) parameters are minimised. A full likelihood analysis (DM and VM) was first done [1] by taking the Navarro-FrenkWhite (NFW) profile [2] for the DM halo, a spheroidal bulge and an axisymmetric disk [3-7]. The halo and disk profiles are normalized to the respective densities at the solar neighborhood. The NFW profile is given by ρDM(r) = ρDM,⊙R⊙ rrs+R⊙ rs+r2where ρDM,⊙is the lo- ∗[email protected] †adit[email protected] ‡vira[email protected] §sa[email protected] ¶Email: [email protected] cal DM density, rsis the scale radius of the halo and R⊙is the solar distance from the galactic center. The bulge and disk density profiles are given respectively by ρb(r) = ρb01 + r2 r2 b−3 2and (in cylindrical coordinates) ρd(R, z) = Σ⊙ 2zdexp −R−R⊙ Rd −|z| zd, where ρb0is the central bulge density, Σ⊙is the local disk surface density, and rband Rdare the scale radii of the bulge and the disk respectively. The parameter zdis the scale height of the disk. From the densities, we can obtain the total gravitational potential Φ(R) at a given radius on the galactic plane, and thus we get the circular rotation speed at R by v2 c(R) = R∂ ∂R Φ(R, z = 0) Now, a Markov Chain Monte Carlo (MCMC) analysis is done using the RC data, and the most likely values of the above parameters (along with their 1σ ranges) is determined. The χ2test statistic used is χ2=PN i=1 vobs,i−vtheo,i σi2, where vobs,i is the observed circular velocity value, vtheo,i is the one theoretically calculated, σithe error in the observed velocity value, and Nis the total number of data points. The likelihood estimation of the parameters have been performed with suitable prior ranges on the parameters. The value of zd is immaterial to our analysis, and it has been fixed [8] at 340 pc. 2 II. MILKY WAY MASS MODEL & LOCAL DARK MATTER We model the Galaxy in terms of three components: •A spherical dark matter (DM) Halo, •A spheroidal visible matter (VM) Bulge, and •An axisymmetric VM disk. We assume a DM halo of the Navarro-Frenk-White (NFW) form [2], where the density ρDM of the halo as a function of galactocentric radius ris given by ρDM(r) = ρDM,⊙R⊙ rrs+R⊙ rs+r2 (1) The VM bulge profile density is given by ρb=ρb01 + r2 r2 b −3 2 (2) while the disk density profile is given by ρd(R, z) = Σ⊙ 2zd exp −R−R⊙ Rd −|z| zd(3) Given a spherically symmetric density distribution ρ(r), the mass enclosed inside a sphere of galactocentric radius rcan be written as M(r) = Z∞ 0 4πr′2ρ(r′)dr′(4) For the NFW profile, it is given by MDM(r)=4πρDM,⊙R0(rs+R0)2log rs+r rs −r r+rs (5) while for the central bulge, it is given by Mb(r)=4πρb0    sinh−1r rbr3 b− rr4 br1 + r rb2 r2+r2 b    (6) The rotational speed v2 cis then given by v2 c=GM(r) r(7) The disk is not spherically symmetric, and the calculation of v2 c(r) is not very straightforward. The steps for this are outlined by Duric, and we show the result here, v2 c,disk(r) = 4πGΣ⊙Rdr 2Rd2 exp R⊙ Rd× I0r 2RdK0r 2Rd−I1r 2RdK1r 2Rd (8) where In(x) and Kn(x) are the Modified Bessel functions of order nof the first and second kinds respectively. We perform the analysis imposing the following priors on the VM disk - •No prior. •Gausssian Prior 1 : Σ⊙= 48 ±8M⊙pc−2, Rd= 2.3±0.6kpc •Gaussian Prior 2 : Σ⊙= 67 ±8M⊙pc−2, Rd= 2.3±0.6kpc III. CONNECTING TRACERS TO DM DISTRIBUTION WITH ROTATION CURVE OF GALAXY : Bhattacharjee et al.(2014) uses several tracers like Cepheids, Planetary Nebulae, Clusters etc. to map the rotation curve of the Milky Way galaxy. Prior to the Gaia mission and the availability of accurate proper motions for several stars, the only available data was spectroscopically determined Line-of-Sight(LoS) velocities. These LoS velocities can be used to reconstruct the rotation curve with the assumption that the tracers all orbit the Milky Way in circular orbits. One of the tracers used are Classical Cepheids, which have accurately determined distances from the Period-Luminosity (PL) relations. Data from Pont et al.(1994) which provide gamma velocities (radial velocities accounting for the pulsation of the surface) and use Period-Luminosity-Colour (PLC) relations and reddening data to determine the distance to the stars. We have knowledge now of the LoS velocities (vh), and the distance to the stars (rh). Since we observe from the Sun’s frame, we also need information about the Sun’s distance from the Galactic Centre, and the velocity of the Sun from the frame of the Galactic Center. The radial distance (R0) from the Galactic Centre (GC), and the Circular Velocity (V0) of the Sun i.e. [R0, V0] defines the Local Standard of Rest (LSR) which follows the mean motion of material near the Sun in the Milky Way. The Sun’s motion is also characterized by its Peculiar Velocity (U⊙, V⊙, W⊙) with respect to the LSR. Shifting from the Sun’s frame to the frame of an observer in the LSR, we can write vLSR =vh+U⊙cos bcos l+V⊙cos bsin l+W⊙sin b(9) This accounts for the Sun’s peculiar motion, which is known from measuring the motion of the stars in the 3 vicinity of the Sun. Now, we need to move from the LSR to the Galactocentric Frame, noting that the LoS velocity measured is the projection of the actual (assumed) circular velocity along the loine of sight. With knowledge of the galactic coordinates (l, b), their relation can be derived as Vc(R) = R R0hvLSR sin lcos b+v0i(10) where Ris the distance of the source from the GC, which is given by R=qR2 0+r2 hcos2b−2R0rhcos bcos l(11) To properly see the features of the Rotation Curve, we bin the data and average the samples in each bin to get a mean Vcand mean Rfor each bin. These values are then plotted. The coverage of the data in the galactic plane can be visualized by a scatterplot between the X and Y coordinates. x=rhcos bsin l y=R0−rhcos bcos l(12) z=rhsin b Using Jean’s Equation for Halo Tracers I Using Jean’s equation (Binney & Tremaine, 2008), we can write V2 c=−σ2 Rdln ntr dln r+dln σ2 R dln r+ 2β(13) where, σR(σT) is the velocity dispersion in the radial(transverse) direction. σ2 T=σ2 ϕ+σ2 θ 2(14) and βis the anisotropy defined as β= 1 −σ2 T σ2 R (15) ntr is the tracer density given by ntr =Number of tracers in a bin 4πR2dR To calculate the velocities in the (R, θ, ϕ) coordinates, we note that the velocities (U2,V2,W2) of the tracer in the Galactocentric frame are along the ( X,Y,Z ) coordinate system defined with the Sun at (0,R0,0) (upto a change in sign). Hence we can convert between them. The 6-D phase parameters provided by the Gaia-DR3 catalog for each star contains the parallax π, the Galactic coordinates (ℓ, b), the radial velocity vrand the two proper motions in equatorial coordinates µαcos δand µδ (which can be transformed into Galactic coordinates as µℓcos band µb). These six variables allow for each star to be placed in the 6D phase space: 3D spatial +3 D velocity. The Galactocentric position in cylindrical coordinates has the following three components: the Galactocentric distance R, the Galactocentric azimuth ϕ1, and the vertical distance Z. Analogously, the Galactocentric velocity in cylindrical coordinates has the following three components: the radial velocity VR, the azimuthal velocity Vϕ, and the vertical velocity VZ. r=1AU π X(r, ℓ, b) = R⊙−rcos ℓcos b, Y(r, ℓ, b) = rsin ℓcos b, Z(r, b) = Z⊙+rsin b, R(r, ℓ, b) = pX2+Y2 ϕ(r, ℓ, b) = tan−1Y X between π 2and 3π 2for X < 0, and V(r, ℓ, b, vr, µℓ, µb) =   VR Vϕ Vz =  −cos ϕsin ϕ0 sin ϕcos ϕ0 0 0 1    U∗+U⊙ V∗+Vg,⊙ W∗+W⊙    U∗ V∗ W∗ =  cos ℓcos b−sin ℓ−cos ℓsin b sin ℓcos bcos ℓ−sin ℓsin b sin b0 cos b   vr vℓ vb  vℓ=rµℓcos b vb=rµb The Galactic parameters used by the previous kinematic Gaia-DR3 analysis (G18):R⊙= 8.34kpc (Reid et al. 2014),Z⊙= 0.027kpc (Chen et al. 2001); (U⊙, V⊙, W⊙) = (11.1,12.24,7.25)km/s(Sch¨onrich et al. 2010);Vg,⊙=Vc,⊙+V⊙, with the rotation speed Vc,⊙= 240 km/s We have knowledge of proper motions in RA (µα) and DEC (µδ) and their errors and covariances from Gaia DR3, i.e. the covariance matrix, Σαδ =σ2 µασµαµδ σµαµδσ2 µδ(16) 4 Transformed to galactic coordinates (Poleski, 2013), we get µl µb=1 cos bC1C2 −C2C1µα µδ(17) and the covariance matrix transforms as Σlb =CΣαδCT(18) where C1 = sin δGcos δ−cos δGsin δcos (α−αG) C2 = cos δGsin (α−αG) cos b=pC12+C22 C=1 cos bC1C2 −C2C1 Now, we find the 3-D velocities (in kms−1), which requires LoS distance rh, of the tracers as Vl= 4.74rhµl Vb= 4.74rhµb The new covariance matrix is given by ΣVlb =FΣlbFT=σ2 VbσVlVb σVlVbσ2 Vl F=4.74µl4.74rh0 4.74µb0 4.74rh where σ2 Vl=V2 lσ2 r r2 h +σ2 l µ2 l σ2 Vb=V2 bσ2 r r2 h +σ2 b µ2 b σVlVb=VlVbσ2 r r2 h +σlb µlµb Now, we can rotate to UVW coordinates at the Sun’s position   U1 V1 W1 =  cos bcos l−sin bcos l−sin l cos bsin l−sin bsin lcos l sin bcos b0   Vr Vb Vl  (20) and the covariance matrix is ΣUV W1=F1ΣVrblFT 1 ΣVrbl =  σ2 Vr0 0 0σ2 VbσVlVb 0σVlVbσ2 Vl   and F1is the transformation matrix in Eq.20. Now, we can add the Sun’s velocity to get the galactocentric velocites. U2 = U1 + U⊙ V2 = V1 + V⊙+V0 W2 = W1 + W⊙ One final rotation to the UVW coordinates at each tracer’s position will give the final values we want.   Us Vs Ws =  cos β−sin β0 sin βcos β0 0 0 1    Vr Vb Vl  sin β=rhcos bsin l R cos β=[R0−rhcos b] cos l R and the final covariance matrix, using the above transformation matrix F2, is Σ = F2ΣUV W1FT 2 LOCAL KINEMATICS For our kinematic analysis, we use data from [18].[18] analyze the kinematics of 412 red giants, located at distances ranging from 1.5−4.5kpc from the Galactic midplane. These 412 stars are selected so as to exclude the thin disk and galactic halo stars and ensure that they belong to the Galactic thick disk. For each star, they calculate the (U,V,W ) velocity components in the Galactic cylindrical reference frame, ( U is positive towards the Galactic center, V is directed in the sense of Galactic rotation and points towards the North Galactic pole). They adopt a solar peculiar velocity of (U,V,W) = (10.0,5.1,7.2)kms−1. In their analysis, they adopt a R⊙= 8.0kpc. The choice of solar motion, however does not affect the velocity dispersions. They divide the sample into 31 bins, with bin centers ranging from 1.5 to 4.5 kpc separated by 0.1 kpc , and calculate the mean (U,V,W) velocities, the dispersion velocities (σU, σV,σW) and the non diagonal terms of the dispersion matrix σ2 UW , σ2 UV , σ2 V W . The errors on the velocities are derived from propagating uncertainties from observed quantities, and the errors on the dispersion velocities are derived by analyzing artificial data through Monte Carlo simulations. For our analysis, we use the σWand σ2 UW data (denoted by σ2 zand σrz ). This data has been previously analyzed by [18] and[19]. We follow the approach described in [18] and [19]. Our analysis is based on the vertical Jean’s equation 5 1 Rν ∂ ∂R (RνσRz) + 1 Rν ∂ ∂ϕ (νσϕz) + 1 ν d dz νσ2 z=−dΦ dz (21) Here Φ is the gravitational potential and νis the stellar number density. On assuming axial symmetry for the Galaxy, the second term on the left hand side drops out, and we obtain 1 Rν ∂ ∂R (RνσRz)+ 1 ν d dz νσ2 z=−dΦ dz =FZ(R, Z) (22) We then write the Poisson equation in cylindrical coordinates Σ(R, z) = −1 2πG Zz 0 dz 1 R ∂(RFR) ∂R +FZ(R, Z)(23) FR=−∂Φ ∂R can be linked to the circular velocity as V2 C=−RFR. As ∂VC ∂R = 0, the integral term can be dropped from the Poisson equation. However, observations show that VCis constant only in the galactic plane. [18] show that approximating the integral in Eq. 23 by its value in the galactic plane (Z = 0) leads to an underestimation of the dark matter density. We adopt this approximation, and integrating 22 gives σ2 z(z) = 1 ν(z)Rz 0ν(z′)Fz(z′)−1 Rν ∂ ∂R (RνσRz)dz′+ C ν(z) where 2πGΣ(z) = dΦ dz =−Fz(24) and C is a normalization parameter. Substituting z = 0 , we have σ2 z(0)ν(0) = C. Following [18], we evaluate 1 Rν ∂ ∂R (RνσRz) (the Tilt term τ) by assuming that radial profiles of the tracer density and σRz can be separated and are exponentials, given by ν(R, z) = ν(z) exp −R hR=ν(R⊙) exp −z hzexp −R hR (25) σRz(R, z) = σRz(z) exp (−R/hσ) (26) We assume hR=hσ= 3.8kpc same as [18], and hz= 0.9kpc. Under these approximations, the tilt term reduces to τ(R, z) = σRz(R, z)1 R−1 hR −1 hσ(27) We assume a linear profile for σRz(z)[18] and as we calculate the velocity dispersions at the solar radius, we have σRz (R⊙, z) = A+B(z−2.5) (28) We fit the Moni Bidin’s data with this profile and obtain best fit values of A = 1033.88 ±64.92 and B = 537.53 ±73.53. Combining all this, we have σ2 z(z) = 1 ν(z)Zz 0 ν(z′) [−2πGΣ (z′)−σRz (R⊙, z′) 1 R⊙ −1 hR −1 hσdz′+σ2 z(0) (29) This is the equation central to the analysis described by [18], and is identical to equation (25) of [18] in their analysis . We split Σ(z) as Σ(z) = ΣDM (z)+ΣBaryon (z) (30) Following the method described by [18], we use Multinest to perform a Bayesian analysis on the [18] σ2 zdata and estimate the parameters in our models for ΣDM and ΣBaryon . We model the local mass density using a single exponential disk, and a NFW dark matter halo. LOCAL KINEMATICS & ROTATION CURVE INDIVIDUAL -VSJOINT ANALYSIS INDIVIDUAL CONSTRAINTS DATASETS, BINNING AND AVERAGING A. Bhattacharjee’s Data Selection •Two independent classes of non-disk stellar tracers1. 4985 BHB stars (SDSS DR8, Xue et al.), 4781 K Giants (SDSS DR9, Xue et al.)- for galactic halo to 100 kpc 2. Heterogeneous sample of 430 objects143 Globular Clusters (Harris, 2010), 118 Red halo giants (Carney et al. 2003, 2008), 108 field blue horizontal branch (Clewley et al. 2004), 38 RR Lyrae (Kinman et al. 2012), 23 dwarf spheroidals (McConnachie 2012) 3. Cut on z and R of tracers, leaving out r <25 kpc. 4. Final sample: 1457 BHB, 2227 K Giants; 16 Globular Clusters, 28 FHB stars, 21 dwarf spheroidals (grouped as Hg) •BHBuniform bin widths, 2 kpc over 25 to 55 kpc. 6 •KGbin width = 2 kpc for 25 to 55 kpc, and 4 kpc for 55 to 103 kpc, anything greater than 103 kpc is one bin. •Hgobject-wise binning, first 6 radial bins with 8 objects each, next 2 bins with 6 objects each, and remaining 5 objects in 1 bin. B. Eiler’s Data Selection •Stars are chosen from upper red giant branch with low surface-gravity, i.e. with a surface gravity of 0≤log g ≤2.2, which selects stars that are more luminous than the red clump. •These need to have existing spectral data from APOGEE and complete G-band information from Gaia, near-infrared from + 2MASS, and W1/W2 infrared bands from WISE (at 3.6 and 4.5 micrometer respectively). •Data quality cuts are made to give 44,784 stars. •For circular velocity curve: 1. Stars within a wedge of 60◦from the Galactic center towards the direction of the Sun, i.e. ±30◦ 2. |z| ≤ 0.5 kpc or lie within 6◦from the Galactic plane, i.e. |z|/R ≤tan(6◦) 3. Exclude stars with a vertical velocity component of |vz|>100 km s−1 4. Stars with low α-element abundances, i.e.[α/Fe]<0.12 5. Gives a total of 23,189 stars for analysis. C. Mroz Data Selection Data from Skowron et al. (2018)- 2177 Galactic Cepheids.Cross-match with Gaia DR2 which gives full information for 832 objects. Known cepheids in binary systems, objects with residual velocities at least 4σtimes larger than mean, where is the dispersion of residuals. Cepheids observed from 2 to 44 times with the median number of seven visits were kept.Gives a final set of 773 Cepheids. To determine the RC of the galaxy, one has to use the kinematical and positional information for some tracer objects moving in the gravitational potential of the galaxy. In general, one measures the line-of-sight (LOS) quantities (positional and kinematical data) and the RC is derived from that. To determine the RC for the disk region, we have to adopt a value for the the local standard of rest (LSR) which corresponds to the position (R0) and velocity (vc0) of the sun with respect to the galactic center, and make the assumption that the tracer objects follow a circular orbit arround the galactic center. From this, the positions and velocities about the galactic center can be obtained. The various tracer objects used for this include HI and HII regions (CO emissions from the latter), Cepheids, Planetary nebulae, etc. For regions extending beyond the visible disk of the galaxy, we look at tracers distributed in the halo of the Milky Way (like dwarf spheroids, globular clusters, K-giant stars, etc.). These tracers are of the non-disk kind, and their motion around the galactic center is typically unsystematic. Under the assumption that these objects are isotropically distributed in the halo, one can use the Jeans to relate their velocity vcand their distance rfrom the center to their observed number densities and velocity dispersions. One disadvantage of this approach is that the velocity anisotropy is completely unknown (since observationally, only the radial velocity dispersion is seen), and can significantly affect results [9]. After the raw RC is generated as above, it is suitably binned and averaged. Different strategies for the above procedure can lead to different results. Eg., The data of Sofue 2013 [10] employs a running average method, where all the raw data points are taken, and the midpoints of the bins are (1 + ϵ) times (ϵ= 0.1 or 0.3) the previous value. The weight in each bin is Gaussian, with σbeing ηtimes the bin-center value. The averaging is done only once. While calculating the mean, the observables ( ror v) at higher radii are given more weight. The formula provided for calculating the error bars shows that errors at higher radii have a higher value. Bhattacharjee et.al. [9] have averaged twice, once over the individual data sets, and then over the first-averaged data, with ”optimal” bin sizes in different radial ranges to capture the small-scale features of the raw RC, as well as to get the closest value to the true RC. Figures 2a and 2bshow the difference in the final RC due to these strategies with the underlying data almost same [9, 10]. The inset of Fig.2 shows the 2D posterior probability contours for the DM parameters, namely ρDM,⊙and rs, obtained from the datasets of Sofue (2013) (henceforth referred to asS13) and Bhattacharjee et.al., both for the LSR of ( 8.0kpc,200 km s−1). The larger error bars in the S13 data, especially around the outer regions of the galaxy, where the effect of the DM is strongly predominant, loosens the posterior constraints on the DM parameters. Also, in Fig.1 we see the rotation curves recreated using the obtained mean values obtained from the MCMC analysis. The blue lines (dased and solid) in either subfigure denote to the curves obtained using flat priors on the disk parameters. Fig.2 shows the comparison between the 1σand 2σ2 D posterior probability contours of the DM parameters, computed for the following cases (with the mean values tabulated in Table IV): 1. Bhattacharjee et.al. (8.0,200) data (with flat priors) with (a) the full dataset, (b) data points above 3 kpc , (c) data points above 5 kpc , and (d) data points above 5 kpc and bulge neglected in the MCMC analysis. 7 FIG. 1. Recreated RC’s from the obtained values of the parameters obtained with weak (flat) and strong (Gaussian) priors on the disk parameters. The dashed (dotted) line represents the RC due to the DM (VM) components, and the color blue (red) signifies weak (strong) priors. The solid line is the combined RC. 2. Sofue. (2013) data (with flat priors) with (a) the full dataset, (b) data points above 3 kpc , (c) data points above 5 kpc . Dataset ρDM⊙GeVcm−3rs(kpc) Bhatt (8.0,200) 0.16+0.01 −0.01 64.83+17.03 −12.94 Bhatt (8.3,244) 0.05+0.01 −0.01 213.15+164.09 −102.76 Bhatt (8.5,220) 0.09+0.02 −0.01 74.58+38.24 −27.56 Sofue (2013) 0.04+0.07 −0.02 106.14+230.52 −86.18 TABLE I: Mean values of the DM parameters obtained with flat priors on the disk parameters. Choice of the Local Standard of Rest We have used three different datasets[9] corresponding to three different choices of the LSR - respectively 8.0kpc,200 km s−1,8.3kpc,244 km s−1and 8.5kpc,220 km s−1. We refer to these henceforth as B80, B83, and B85 respectively. The dataset of S13 [10] uses a LSR of 8.0kpc,200 km s−1. Fig.2 shows the 1σand 2σ2D posterior probability contours of the DM parameters ( ρDM,⊙and rs) for the three choices of the LSR’s after marginalizing over the other (VM) parameters. The corresponding mean values are given in Table I. It is seen that the choice of the LSR has a significant effect on the posterior PDF of the parameters, with nearly exclusive zones in the ρDM,⊙−rs space. In Fig.2, we can see the comparison between the posterior probality contours obtained using different choice of LSR’s for the datasets. The blue and the green contours differentiate between the choices of LSR of ( 8.5kpc,220 km s−1and 8.0kpc,200 km s−1for flat priors (on the disk parameters), while the respective lighter contours do so for the Gaussian priors (as specified in more details in the next section). Parameter Flat Priors Gaussian Priors ρDM⊙GeVcm−30.16+0.01 −0.01 0.18+0.02 −0.01 rs(kpc) 64.83+17.03 −12.94 57.21+140 −10.95 ρDM⊙GeVcm−30.04+0.07 −0.02 0.16+0.07 −0.08 rs(kpc) 106.29+230.52 −86.18 24.16+41.22 −10.77 TABLE II: DM parameters obtained with and without priors on the disk parameters, respectively Bhattacahrjee et.al. (8.0),Sofue 2013. EFFECT OF VISIBLE MATTER PRIORS ON DM CONSTRAINTS The addition of contraints on some of the visible matter parameters causes a significant change in the posterior probabilities of the parameters. In particular, large shifts are seen in the obtained mean values of the DM parameters. In our analysis, we have put gaussian priors [5,11] of Σ0= 48 ±8M⊙pc−2and Rd= 2.3±0.6kpc. Table II lists the obtained mean values and 1σranges of the DM parameters obtained from the Bhattacharjee et.al data set (8,200), and the Sofue 2013 data sets, and Figure 1 shows the theoretically recreated RC’s using these values of the parameters.We have obtained the 1σand 2σ2D posterior probability contours of (i) the DM parameters, (ii) the ρDM,⊙−Rdparameters, and (iii) the ρDM,⊙−Σ⊙ parameters obtained in each case with and without the use of priors on the disk densities, for the Bhattacharjee et.al dataset for three LSR’s - (a) 8.3kpc,244 km s−1, (b) 8.5kpc,220 km s−1, and (c) 8.0kpc,200 km s−1. The obtained mean values for these parameters are tabulated in Table III. 8 FIG. 2. Comparison of posterior probability density with different LSRs with prior i.e. flat prior (left) and with prior(right). Dataset ρDM⊙GeVcm−3rs(kpc) Σ⊙M⊙pc−2Rd(kpc) Bhatt. et.al. (8.0, 200) 0.16(0.18) 64.83(57.21) 53.37(60.72) 2.75(2.72) Bhatt. et.al. (8.3, 244) 0.05(0.37) 213.15(7.53) 213.33(119.88) 4.73(3.91) Bhatt. et.al. (8.5, 220) 0.09(0.29) 74.58(8.81) 138.88(77.56) 4.26(3.60) TABLE III: DM and VM parameters obtained with flat (Gaussian) priors on the disk parameters. Dataset ρDM⊙GeVcm−3rs(kpc) Bhatt 8.0 (Full) 0.16 64.83 Bhatt 8.0 (Above 3 kpc) 0.18 52.82 Bhatt 8.0 (Above 5 kpc, No Bulge) 0.19 46.5 Sofue 2013 (Full) 0.04 106.14 Sofue 2013 (Above 3 kpc) 0.06 65.76 TABLE IV: DM parameters obtained with flat priors on the disk parameters on different modified versions of the datasets.for comparison see FIG. 2 Introducing Bulge Density and Kinematics with RC: All the analysis before this section was done without introducing bulge density and corresponding correction in our model. We performed Monte Carlo Markov Chain using Bhattacharjee datasets for 3 different LSRs but without any prior (flat prior) but introducing bulge density, again by doing that we have seen a significant shift in 6 parameters tabulated in Table V. Further, within the same model using Bhattacharjee 8.0 data, if we combine Eilers data, the combined fitted rotation curve is shown in Fig 3. The mean VM parameter reduced drastically, but again, within the constraint of a 10 percent error bar of the combined Eilers+Bhatt dataset, we recover the similar (Rho DM, Sigma) parameters.The MCMC of 6-D parameter space is shown in Fig. 3 Similarly, we have also combined Mroz rotation curve for cepheids with Bhattacharjee(8.0,200),noting that bulge information is only coming from bhattacharjee’s dataset . The circular velocity of a test particle in the Galaxy is, of course, not a directly measured quantity. The RC of the Parameter (8.0,200)(8.3,244)(8.5,220)EilersMroz ρDM,0[GeV cm−3] 0.32 0.26 0.37 0.38 0.432 rs[kpc] 83.80 93.35 207.73 95.20 100.10 ρb,0[105M⊙pc−3] 0.121 0.533 0.906 0.12000.1250 rb[kpc] 0.03 0.03 0.03 0.03 0.03 Σ0[M⊙pc−2] 160.56 213.65 79.81 35.30 35.9 Rd[kpc] 4.41 4.92 2.94 3.14 2.014 TABLE V : Fitted 6-D Parameters after introducing Bulge density) for comparison of VM-DM parameter see FIG.4 Galaxy has to be derived from the kinematical as well as positional data for an appropriate set of tracer objects moving in the gravitational field of the Galaxy. Except in few cases, the full three-dimensional velocity information of the tracers is not available, and the RC has to be reconstructed from only the measured line-of-sight (LOS) velocity and positional information of various tracer objects in the Galaxy.Since currently not much reliable observational information on βis available, in this article we calculate the circular velocities using the Jeans equation with the velocity anisotropy βof the tracers taken as (1) a radially constant free parameter varying over a possible range of values from β= 0 (corresponding to complete isotropy of the tracers’ orbits) to β= 1 (corresponding to completely radial orbits of the tracers), (2) a radially varying βof the Osipkov–Merritt (OM) form (see Binney & Tremaine 2008, pp. 297–298) given by β=1 + r2 a r2−1 , where rais the “anisotropy radius,” and (3) a radial profile of βobtained from recent large high-resolution hy- 9 2.75 3.00 3.25 3.50 Rd , 1 15 30 45 60 [ M pc 3] 4 8 12 16 20 rs 30000 90000 b , 0 0.04 0.06 0.08 0.1 0.12 rb 0.30 0.35 0.40 0.45 DM ,[ GeVcm 3] 2.75 3 3.25 3.5 Rd , 1 15 30 45 60 [ M pc 3] 4 8 12 16 20 rs 30000 90000 b , 0 0.04 0.06 0.08 0.10 0.12 rb run1 run2 a_1_ a_1_ 2.5 5.0 7.5 10.0 Rd , 1 0 60 120 [ M pc 3] 8 16 24 32 rs 0 200000 b , 0 0.05 0.1 0.15 0.2 rb 0.1 0.2 0.3 0.4 0.5 DM ,[ GeVcm 3] 2.5 5 7.5 10 Rd , 1 0 60 120 [ M pc 3] 8 16 24 32 rs 0 2 b , 0 1e5 0.05 0.10 0.15 0.20 rb run1 run2 run3 FIG. 3. Posterior probability contours of combined Eilers+Bhatt(β13 dataset(left) and with additional 10% constraint on errorbar(right) drodynamical simulations of the formation of late-type spirals like our Galaxy (Rashkov et al. 2013) or β13 model . From Fig. 2 and Table IV,V if we combine Bhattacharjee kinematics ( >5 kpc ) with Eilers TRGB and Mroz Cepheids, both the combined version has degeneracy withing the 2 σuncertainity with their respective VM-DM parameter space but combining Bhattacharjee’s kinematics with >3 kpc Mroz Cepheids shows significant dispacement from its original RC in its parameter space with broken degeneracy as shown in figure below: 0.18 0.24 0.30 0.36 0.42 0.48 DM ,[ GeVcm 3] 16 24 32 40 48 56 64 [ M pc 3] Eilers RC Cepheids RC+Kinematics Bhatt RC+Kinematics Eilers RC+Kinematics Cepheids RC 0.15 0.20 0.25 0.30 0.35 0.40 0.45 DM ,[ GeVcm 3] 16 24 32 40 48 56 64 [ M pc 3] New_run3 Cep_Combined Bhatt80 Eilers_Combined Cep FIG. 4. VM-DM posterior contour comparison after combining different kinematics z > 3 (left ) and z > 5(right) of Bhatt. et al Joint Constraint with GAIA DR3 6-D phase plot: - Finding data from the Gaia archives became the next step to get rotation curves beyond ∼20kpc. But Gaia parallaxes are reliable (p/σp>10, i.e. 10% error in the parallax) upto ∼5kpc from the Sun. So using Gaia parallaxes alone would only result in a rotation curve upto ∼13kpc. - Gaia has detected Cepheids in the Milky Way, in a separate catalogue. The periods can be used to infer parallax from the PL relations. (Metallicity measurements are scarce). The sources can then be cross-matched to the Gaia sources (which has the astrometric and photometric data). RR Lyrae stars may also be used - Gaia parallaxes have been noted to be systematically lower than accurate non-Gaia parallaxes by ∼0.05 mas (Groenewegen, 2018; Riess et al., 2018). (See also Ripepi et al. (2019), which confirms this under-estimated zero point for parallaxes using Cepheids) .Ripepi et al. (2019) also finds PL and PLZ relations which may be used.Complications with using PL relations includes: Non-fundamental mode oscillators;Possible Binaries; Possible errors in the period determination by the Gaia pipeline (Groenewegen, 2018). If we limit ourselves to the Gaia parallaxes, we can then use the sources with measured radical velocities. We can now construct a rotation curve. Cuts were made on the parallax errors (within 10% ), distance from the Galactic Plane ( <0.5kpc ) and vertical velocity component Ws<50 km s−1 This was also compared with a rotation curve constructed using TRGB stars using Jean’s equation for axisymmetric systems (Eilers et al., 2019).Xue et al. (2011) gives LOS velocities and distances to 4985 Blue Horizontal Branch (BHB) stars, of which 4982 were crossmatched with Gaia DR2 (for accurate positions and for proper motions). Similarly, from Xue et al. (2014), we get LOS velocities and distances to 6036 K Giants (all crossmatched) and proper motions from Gaia. Using the Sun’s Galactocentric Distance ( R0) as 8.122 kpc (Eilers et al., 2019), we calculate the radial distance from the 16 where ρcis the critical density, and the overdensity δc=200 3 c3 200 ln(1+c200)−c200 1+c200 (not to be confused with ∆c) and the scale radius rsare unique to each halo, and the concentration parameter is given by c200 =r200 rsIn place of δcρc, ρsis often used, where ρsis a parameter unique to each halo. The total mass of the dark matter halo can then be computed by integrating over the volume of the density out to the virial radius r200 : M=Zr200 0 4πr2ρ(r)dr = 4πρsr3 sln r200 +rs rs−r200 r200 +rs =4 πρsr3 shln (1 + c200)−c200 1+c200 i From the definition of the circular velocity, Vc(r) = qGM(r) r, we can find the circular velocity at the virial radius r200 : V200 =rGM200 r200 . Then the circular velocity for the dark matter halo is given by V2 c(r) = V2 200 1 x ln(1 + cx)−(cx)/(1 + cx) ln(1 + c)−c/(1 + c) where x=r/r200. Although the NFW profile is commonly used, other profiles like the Einasto profile and profiles that take into account the adiabatic contraction of the dark matter due to the baryonic content are also used to characterize dark matter halos. The MW has several baryonic components, including a central nucleus, a bulge, and a disk. While many of their properties remain topics of debate, their masses are reasonably well determined (Bland-Hawthorn and Gerhard 2016). From kinematical and dynamical studies, we know that the mass of the MW must be larger than the sum of the baryonic components; indeed, to maintain the system in stable equilibrium with the observed amplitude of the circular velocity, a large fraction of the MW mass must be invisible, i.e., we cannot measure it directly, but we can infer its presence by its gravitational influence. Despite decades of intense efforts, the estimates of the mass of the MW still show significant scatter. These estimates are very sensitive to assumptions made in the modeling and, in particular, to the shape of the halo in which the Galaxy is embedded. Most mass estimators are limited to the region explored by the available tracer population, whose spatial distribution and kinematics are used to estimate the enclosed mass. FIG. 12. Comparison of Φ V2 cfor Different LSR of Bhatt data . In the FIG 12 the ratio Φ V2 c, where Φ(r) is the potential of the disc and Vcis the circular velocity (rotation speed),where if Φ(r) has contributions from other components (like the halo or bulge), this ratio can also help in understanding the relative dominance of the disk potential compared to other components influencing the rotation.The dimensionless ratio Φ V2 ccan then be expressed as: Φ(r) V2 c(r)=Φ(r) rdΦ(r) dr . . Estimates of the MW’s mass have been obtained based on the kinematics of halo stars, the kinematics of satellite galaxies and globular clusters, the evaluation of the local escape velocity, and the modeling of satellite galaxy tidal streams. Estimates typically range from as low as 0.5× 1012M⊙to as high as 4×1012M⊙(Bland-Hawthorn and Gerhard 2016). These estimates assume that dark matter (DM) is in a quasi-spherical virialized halo around the Galaxy (Navarro et al. 1997; Sanders 2010). We have used the standard NFW profile for dark matter: ρDM(r) = ρs r rs1 + r rs2, where ρsand rsare the characteristic density and scale radius, respectively. Additional components (e.g., bulge or disk) can be included as needed. the critical density is given by ρcrit =3H2 8πG, where His the Hubble constant. Define the virial radius r∆as the radius where the mean density within a sphere is ∆ ·ρcrit: ⟨ρ(r∆)⟩=M(r∆) 4 3πr3 ∆ = ∆ ·ρcrit, where ∆ = 500 or 200. we have numerically Computed the mass enclosed within radius rby integrating the total density profile: M(r) = Zr 0 4πr′2ρtotal(r′)dr′. 17 = Zr 0 4πr′2[ρDM(r′) + ρbulge(r′) + ρdisk(r′)] dr′. Iteratively solve for r500 and r200 by satisfying the mean density condition. The virial mass is: M∆=M(r∆). The concentration parameter is c∆=r∆ rs. Table VII shows the average values of some of the quantities that physically characterize the galaxy, derived using the mean values of the density parameters for the various datasets.Fig. 4 shows the functional dependence of the mean value of the total mass Mtotal(r) enclosed within a ra dius r for all the datasets and with chioices of at and Gaussian priors. FIG. 13. Comparison of the average value of the Mtotal(r) (halo, disk and bulge) from different references(up) ;for different anisotropy obtained from Bhatt(bottom left) and for all the datasets in our work(bottom right), with at falt (solid) and Gaussian (dashed) priors.Blue Red Green and Brown represents B8.0, B8.3, B8.5 and S8.0 respectively from the Table VII.Triangles represent (M500,c500) pairs and circles represnt (Mvir,cvir) pair.Filled and empty symbols represents Flat and Gaussian priors respectively,where as all gray filled circles represent all datasets from upper Table VII. Included are the masses obtained in various references with symbols - [12], [13], [14], [15], and [16]. In Fig 13 the 1st Row image describes milky Way mass profile for the maximum posterior density parameters (black dashed curve) and the corresponding 68%/95% credible intervals (dark/light gray shade), conditioned on the radius and model-averaged over baryonic morphologies[22]. Also plotted are results from several other studies of the Milky Way’s cumulative mass distribution [23-44]. The markers around 50 kpc and 100 kpc are artificially dispersed horizontally so that they are distinguishable. The black arrow denotes the latest lower bound for Mtot from [46]. As with figure 8, note that different choices for R0 and/or V0 among studies can induce apparent discrepancies. DISCUSSIONS FIG. 14. DM VM parameter space for different References in comparison with our work. Fig 4 shows the comparison of some of our results with that of the results obtained in some of other earlier works [12-17] regarding the local DM and VM (disk) densities.This figure underscores the effect of not only the choice of the LSR, the binning and averaging strategies used, and the priors on the visible parameters, but also on the technique used to obtain the best-fit values of the parameters. For instance,Zhang et.al. [14] have analyzed the distribution of the 9000 K-dwarfs in position and velocity space in the vicinity of the Sun, as opposed to the MCMC method used in this work. It is interesting to note that Catena and Ullio, considering an LSR of 8.33 ±0.35kpc,245 ±10.4 km s−1, and a local disk density prior of 48±8M⊙pc−2have obtained the value of ρDM⊙= 0.389 ±0.025GeVcm−3, while in our work we have obtained ρDM⊙= 0.37 ±0.02GeVcm−3with the Bhattacharjee et.al dataset corresponding to the LSR of 8.3kpc,244 km s−1. The ongoing GAIA mission is working to build a very precise 3D map of our galaxy by surveying over a billion stars, studying their motions and other properties.By combining both horizontal(Coherent) and vertical kinematics this will give us more tighter constraint on to break degeneracy on VM-DM parameters without any prior assumed. This is to significantly reduced the uncertainty in the determination of the LSR, and give a much better estimation of the distribution of visible matter in the galaxy. This, in turn gave tighter constraints on the local DM parameters. NOTES For LSR80 - Comparing contour plots of rs−ρDM , we see that ((i) with external priors, the error bars are thicker for density, since the prior likelihood gives weigtage at point far away from the best fit without prior from data, (ii) as prior goes to smaller disk density, the best fit DM moves to higher values. However, the kinematic data is more affected due to priors. (iii) comparing with fixed zd, 18 Datasets ρDM,⊙Σ⊙rsRD,1M200 c200 M500 c500 [GeV/cm3] [M⊙/pc2] [kpc] [kpc] [M⊙] [M⊙] Cepheids-β13-nobulge 0.469 ±0.006 31.1±1.6 16.6±1.1 1.903 ±0.026 1.34 ×1012 13.67 1.10 ×1012 9.43 Cepheids-β13-Bhatt 0.426 ±0.007 37.0±1.5 19.2±1.3 2.156 ±0.027 1.39 ×1012 11.97 1.13 ×1012 8.23 Eiler-nobulge-freevisible 0.525+0.007 −0.004 2.78+0.98 −2.73.2+0.08 −0.17 5.36 ×1011 52.23 4.77 ×1011 37.02 Eiler-freevisible 0.091+0.031 −0.043 98.2+11 −996+40 −70 5.79+0.53 −0.53 Eiler-run2-Bhatt 0.25+0.033 −0.013 66.3+3.1 −8.076+38 −56 3.478+0.006 −0.008 Eiler-β13-Bhatt-run1 0.38+0.03 −0.026 35+6 −79.5+1.4 −2.33.14+0.14 −0.14 6.39 ×1011 18.65 5.38 ×1011 12.98 Eiler-β13-Bhatt-run2 0.332+0.014 −0.01 48.2+3.0 −3.413.8+1.3 −1.63.194+0.061 −0.094 7.43 ×1011 13.50 6.09 ×1011 9.31 Eiler-β13-Bhatt-combined 0.432 ±0.0063 24.3±1.5 6.23 ±0.51 3.26 ±0.38 5.64 ×1011 27.29 4.87 ×1011 19.14 Cepheids-β13-Bhatt-combined-R0R01-free 0.4321+0.0061 −0.0051 35.9+1.3 −1.518.9+1.0 −1.22.17+0.028 −0.028 1.39 ×1012 12.16 1.13 ×1012 8.36 Cepheids-β13-Bhatt-combined-R0R01-fixed 0.4669 ±0.0049 26.3±1.1 13.28+0.84 −0.95 2.014+0.023 −0.026 1.08 ×1012 15.88 8.95 ×1011 11.00 GaiaDR3−Bhatt −combined (r∼200kpc,with bulge,falt prior:) 0.450.04 −0.04 49.304.03 −5.31 62.4031.97 −8.79 7.661.20 −2.60 4.22 ×1012 7.51 3.25 ×1012 5.07 GaiaDR3−Bhatt −combined all star crossmatched (r∼200 kpc with bulge) 0.200.12 −0.05 6.833.41 −3.14 47.5312.12 −18.24 7.821.08 −0.65 4.99 ×1010 7.71 3.86 ×1010 5.22 Dataset M500 (1011M⊙)M(rmax) (1011M⊙)Mvir (1011M⊙)c500 cvir Bhatt. et al. (8.0, 200) 13.30 (13.60) 14.63 (14.79) 25.33 (23.12) 2.54 (2.90) 5.25 (5.91) Bhatt. et al. (8.3, 244) 9.80 (6.12) 12.46 (6.98) 24.75 (7.55) 0.70 (16.86) 1.62 (30.93) Bhatt. et al. (8.5, 220) 9.05 (5.17) 11.08 (6.11) 16.22 (6.54) 1.94 (13.62) 4.03 (25.19) Sofue 2013 (8.0, 200) 3.99 (4.82) 6.60 (6.42) 8.03 (7.07) 1.04 (4.85) 2.24 (9.43) TABLE VII : Summary of mass and concentration parameters for different datasets. again affects the kinematic, and pushed kinematic and combined rhoD to higher values Comparing rhodm-sigma0 plot, it isobvious that DM-VM degenracies dominate kinematic data. Fixing zd, reduces error for kinematic but moves away from RC. By fixing zd or putting priors, kinematic and RC contours move apart, and overlap is only 3 sigma region mostly. With unifirm prior, their 2sigma regions overlap, naturally.Disk with z > 5 kpc is better as there is more overlap between kinematic and RC. z > 5 kpc disk pushes Kinematic rhoDM to lower values, RC values remain almost same. For Kinematics - thetrajectory from Bovy Tremaine , LSR 80 to LSR 85, z > 3kpc vs z > 5kpc disk, and fixed vs free zd analysis can be done in future.It is worth to note that LSR 85 has larger Error bars.Also in FIG 10 it is worth noting that the halo density and halo radius highly dependent upon the star we are selecting and the number of tracers and the comparatively larger errors in ρDM occurring due to the consideration of larger radius r ∼200 kpc unlike others having r ∼25 to 100 kpc Flat galactic rotation curves can be explained either by introducing a dark matter component, like in the ΛCDM model, or by modifying Newtonian dynamics, like in MOND theories. However, it has been noticed that the acceleration systematically goes down at large distances, i.e. low accelerations. This systematic downward trend in the weak gravity segment of the galactic rotation curves was interpreted as an external field effect within the framework of the MOND theories. However, this feature was observed in the ΛCDM model without invoking any modification of Newtonian gravity, as was first found by one of the authors in the EAGLE simulation .They showed that dark halos in ΛCDM generate an acceleration feature like MOND predicts, but that in the even-larger-radius, even-lower-acceleration regions, the dark halos are saturated and the acceleration then decreases as Kepler’s Laws, i.e. Newtonian gravity, but not MOND predict.analyzing Milky ways RC(also M31De-Chang Dai) At very large distances, these data are easier to accommodate in the ΛCDM model. ACKNOWLEDGEMENT We thank the Department of Theoretical Physics of TIFR for allowing us the use of its computational facilities and resources. We also thank Dr. Basudeb Dasgupta, Dr. Rishi Khatri, Prof. Bhuvnesh Jain, Vivek Chaurasia and Simon Birrer; Sravya Kalachaveedu from IISER Mohali and Amogh Srivastava from IIT Bombay for the various useful and fruitful discussions and helps. Appendix A ADQL Quesries for fetched Objects FGKM Spectral Type stars : SELECT vrr.*, gs.pmra, gs.pmdec, gs.radial_velocity, gs.parallax, gs.ra, gs.dec, gs.l, gs.b, gs.parallax_error, gs.radial_velocity_error FROM gaiadr3.gold_sample_fgkm_stars AS vrr JOIN gaiadr3.gaia_source AS gs ON vrr.source_id = gs.source_id OAB Spectral Type stars: SELECT vrr.*, gs.pmra, gs.pmdec, gs.radial_velocity, gs.parallax, gs.ra, gs.dec, gs.l, gs.b, gs.parallax_error, gs.radial_velocity_error 19 FROM gaiadr3.gold_sample_oba_stars AS vrr JOIN gaiadr3.gaia_source AS gs ON vrr.source_id = gs.source_id TRGB: SELECT source_id, phot_g_mean_mag, bp_rp, l, b, ra, dec, pmra, pmdec, radial_velocity, radial_velocity_error, ra_error, dec_error, parallax, parallax_error FROM gaiadr3.gaia_source WHERE phot_g_mean_mag <= 13 AND phot_g_mean_mag >= 8 AND bp_rp >= 1.3 Classical Cepheids: SELECT vrr.*, gs.pmra, gs.pmdec, gs.radial_velocity, gs.parallax, gs.ra, gs.dec, gs.l, gs.b, gs.radial_velocity_error FROM gaiadr3.vari_cepheid AS vrr JOIN gaiadr3.gaia_source AS gs ON vrr.source_id = gs.source_id Blue Horizonntal Branch(BHB) stars: SELECT source_id, phot_g_mean_mag, bp_rp, l, b, ra, dec, pmra, pmdec, radial_velocity, radial_velocity_error, ra_error, dec_error, parallax, parallax_error FROM gaiadr3.gaia_source WHERE parallax_over_error >= 5 AND parallax > 0 AND (4.74 / parallax * pm) >= 145 AND phot_g_mean_mag + 5 + 5 * LOG10(parallax / 1000) < 138.07 * POWER(bp_rp, 6) - 153.85 * POWER(bp_rp, 5) - 40.727 * POWER(bp_rp, 4) + 73.368 * POWER(bp_rp, 3) - 7.4054 * POWER(bp_rp, 2) - 9.5575 * bp_rp + 3.8459 AND phot_g_mean_mag + 5 + 5 * LOG10(parallax / 1000) > -3.2382 * POWER(bp_rp, 3) + 7.1259 * POWER(bp_rp, 2) - 3.583 * bp_rp - 0.2 AND bp_rp < 0.5 AND bp_rp > -0.4; RR Lyrae Stars : SELECT vrr.*, gs.pmra, gs.pmdec, gs.radial_velocity, gs.ra, gs.dec, gs.parallax, gs.l, gs.b, gs.parallax_error, gs.radial_velocity_error FROM gaiadr3.vari_rrlyrae AS vrr JOIN gaiadr3.gaia_source AS gs ON vrr.source_id = gs.source_id; Long Period Variable (LPV) stars or Red Supergiants(RS): SELECT vrr.*, gs.pmra, gs.pmdec, gs.radial_velocity, gs.parallax, gs.ra, gs.dec, gs.l, gs.b, gs.radial_velocity_error, gs.parallax_error, gs.ra_error, gs.dec_error, gs.pmra_error, gs.pmdec_error FROM gaiadr3.vari_classifier_result AS vrr JOIN gaiadr3.gaia_source AS gs ON vrr.source_id = gs.source_id WHERE vrr.best_class_name = ’LPV’ --or ’RS’; Appendix B Catalog query for Crossmatched sources: Following tables from the Catalog has been fetched that are already being crosshatched according to Best nearest Neighbor Gaia Pipeline that uses K-Tree Algorithm to crossmatch considering different EPOCH of different sources : gaiadr2.sdssdr9_best_neighbour gaiadr3.allwise_best_neighbour gaiadr3.apassdr9_best_neighbour gaiadr3.gsc23_best_neighbour gaiadr3.panstarrs1_best_neighbour gaiadr3.ravedr5_best_neighbour gaiadr3.ravedr6_best_neighbour gaiadr3.sdssdr13_best_neighbour gaiadr3.skymapperdr2_best_neighbour gaiadr3.tmass_psc_xsc_best_neighbour gaiadr3.urat1_best_neighbour Gaia Selected Stars Gaia Crossmatched r(kpc) vcΣvcr(kpc) vcΣvc 3.103 200.20 0.129 6.537 232.45 0.001 6.401 234.71 0.002 11.885 223.30 0.003 10.515 215.72 0.001 17.867 211.40 0.640 14.107 212.75 0.151 22.453 193.55 2.565 19.226 201.66 1.042 28.017 189.05 2.588 22.058 192.73 2.808 34.052 193.47 3.241 26.910 206.00 1.484 40.086 178.17 3.707 30.948 214.67 5.071 46.121 202.44 4.429 34.052 191.69 1.991 52.155 192.79 5.513 37.328 211.89 6.420 58.190 193.68 5.812 40.810 195.71 2.247 64.224 180.39 5.988 46.121 207.07 2.672 70.259 160.80 6.992 50.086 213.96 9.114 76.293 221.01 7.572 52.155 198.43 3.017 82.328 179.46 7.263 57.845 199.79 3.199 88.362 181.29 6.463 62.845 212.76 11.790 94.397 105.10 0.045 64.224 210.36 3.592 100.431 122.15 7.087 70.052 181.43 3.805 106.466 102.88 7.317 76.155 218.37 4.059 112.500 110.14 5.379 82.259 181.48 4.535 118.534 142.05 10.738 88.362 186.89 4.804 124.569 153.60 6.172 94.466 159.56 5.246 130.603 224.52 13.065 100.603 169.66 5.344 136.638 106.22 12.929 106.466 184.78 6.068 142.672 144.31 12.071 107.500 209.78 22.625 148.707 105.11 3.510 112.845 156.96 1.665 154.741 139.43 13.745 118.879 145.75 5.640 160.776 134.10 8.724 124.569 170.50 6.723 166.810 144.67 13.357 126.638 220.53 27.217 172.845 200.03 16.294 131.207 151.12 6.225 178.879 73.86 8.220 136.638 175.71 7.856 184.914 110.57 0.836 139.397 160.43 22.826 190.948 95.10 1.706 143.448 151.25 7.029 148.707 176.70 8.204 152.155 195.19 27.771 154.741 263.83 9.002 160.029 197.65 8.745 164.914 214.58 39.791 166.810 166.74 9.461 172.328 193.40 9.618 178.578 146.47 4.904 184.698 162.82 3.169 190.690 213.40 10.783 196.810 162.66 33.932 TABLE VIII :Binned Averaged RC for gaia selected stars and Gaia all stars crossmatched stars. 20 Object Number of Samples Binning Criteria Polynomial Fit Parameters (n0, γ) ntr =n0r 50 kpc −γ(σ0, α)σGSR =σ0r 50 kpc −α 2MASS 3,000,000 0–25 kpc: 20 bins 25–50 kpc: 10 bins 50–75 kpc: 5 bins 75–100 kpc: 4 bins 100–200 kpc: 2 bins each 8.3018 ×10−5,−1.6089 113.137,0.046 AllWISE 3,000,000 0–25 kpc: 20 bins 25–50 kpc: 10 bins 50–100 kpc: 4 bins 100–200 kpc: 2 bins each 7.0648 ×10−5,−1.5722 116.124,0.079 APASS 3,000,000 0–25 kpc: 20 bins 25–50 kpc: 10 bins 50–75 kpc: 5 bins 75–100 kpc: 4 bins 100–200 kpc: 2 bins each 1.9925 ×10−4,−0.6533 99.637,0.030 GSC 3,000,000 0–25 kpc: 20 bins 25–50 kpc: 10 bins 50–100 kpc: 4 bins 100–200 kpc: 2 bins each 1.0812 ×10−4,−1.2691 90.764,0.275 PanSTARRS 3,000,000 0–25 kpc: 20 bins 25–50 kpc: 10 bins 50–100 kpc: 4 bins 100–200 kpc: 2 bins each 1.2640 ×10−4,−1.2737 91.187,0.279 RaVE (DR5) 450,659 0–25 kpc: 20 bins 25–50 kpc: 10 bins 50–100 kpc: 4 bins 100–200 kpc: 2 bins each 4.1230 ×10−5,−0.3575 110.581,0.049 RaVE (DR6) 450,597 0–25 kpc: 20 bins 25–50 kpc: 10 bins 50–100 kpc: 4 bins 100–200 kpc: 2 bins each 4.1680 ×10−5,−0.3294 110.277,0.047 SkyMapper 3,000,000 0–25 kpc: 20 bins 25–50 kpc: 10 bins 50–100 kpc: 4 bins 100–200 kpc: 2 bins each 1.1610 ×10−4,−0.5780 116.882,0.085 Tycho 2,518,330 0–25 kpc: 20 bins 25–50 kpc: 10 bins 50–100 kpc: 4 bins 100–200 kpc: 2 bins each 8.5110 ×10−6,−0.1162 142.653,0.070 URAT 3,000,000 0–25 kpc: 20 bins 25–50 kpc: 10 bins 50–100 kpc: 4 bins 100–200 kpc: 2 bins each 2.7730 ×10−4,−1.2265 104.644,0.119 SDSS (DR9) (with Gaia DR2) 3,000,000 0–25 kpc: 15 bins 25–50 kpc: 5 bins 50–100 kpc: 4 bins 100–200 kpc: 2 bins 2.3200 ×10−4,−1.2289 110.114,0.190 SDSS (DR13) (with Gaia DR3) 3,000,000 0–25 kpc: 15 bins 25–50 kpc: 5 bins 50–100 kpc: 4 bins 100–200 kpc: 2 bins 2.4470 ×10−4,−1.6078 92.567,0.177 TABLE:IX Summary of catalogs with binning criteria and polynomial fit parameters. The polynomial fit parameters are given as (n0, γ) and (σ0, α), following the equations for ntr and σGSR. 21 0 25 50 75 100 125 150 175 r[kpc] 10 5 10 4 10 3 10 2 10 1 V[km/s] Eilers 8.09_233.6_rash_BHB 8.09_233.6_rash_KG Gaia RV Beta13 Bhatt8.09 2.75 3.00 3.25 3.50 3.75 Rd , 1 15 30 45 60 [ M pc 3] 4 8 12 16 20 rs 30000 90000 b , 0 0.04 0.06 0.08 0.1 0.12 rb 0.30 0.35 0.40 0.45 DM ,[ GeVcm 3] 2.75 3 3.25 3.5 3.75 Rd , 1 15 30 45 60 [ M pc 3] 4 8 12 16 20 rs 30000 90000 b , 0 0.04 0.06 0.08 0.10 0.12 rb New_run3 New_run4 New_run5 FIG. 15. (Left) RC all fractional error. (Right) triangle plot rot eilers 8.09 233.6 bhatt beta13 New run 3 4 5. 20 40 60 80 100 120 140 160 180 r[kpc] 0.1 0.2 0.3 0.4 0.5 V[km/s] Beta13_8.3_244 Rash_KG_8.3_244 10 1100101102 r [kpc] 100 150 200 250 300 Velocity [Km/s] Run_1 Run_2 New_Run_3 New_Run_4 New_Run_5 Eilers FIG. 16. (Left) RC beta13 rash KG fractional error. (Right) RC fit eilers bhatt beta13 8.09 233.6 5 runs log scale. 0 20 40 60 80 100 120 r[kpc] 0 50 100 150 200 250 300 V[km/s] GAIA Eilers 8.09_233.6_rash_BHB 8.09_233.6_rash_KG Gaia RV Beta13 0.28 0.32 0.36 0.40 0.44 0.48 DM ,[ GeVcm 3] 16 24 32 40 48 56 64 [ M pc 3] RC_Eilers_run1 RC_Eilers_run2 new_run1 new_run2 new_run2 FIG. 17. (Left) RC all last point removed beta13. (Right) rho sigma0 rot eilers 8.09 233.6 bhatt beta13 5 runs version2. 0 5 10 15 20 25 r[kpc] 0 100 200 300 400 500 V[km/s] 8.09_233.6 8.3_244 8.5_220 8.0_200 FIG. 18. RC comparision 8 8.09 8.3 8.5. 0 5 10 15 20 25 r[kpc] 0 100 200 300 400 500 V[km/s] 8.09_233.6 8.3_244 8.5_220 8.0_200 FIG. 19. RC comparision 8 8.09 8.3 8.5. 22 FIG. 20. All RCs and their contour plot with out constraints FIG. 21. VDF caompared with standard Halo model(SHM) FIG. 22. Left : VDF compared to SHM Right :rho(dm) vs r plot Bhattacharjee for 3 different LSR, why there are two bumps? 23 FIG. 23. Grand RC with log scale(left) and Grand RC(right) FIG. 24. Kinematics+RC for all datasets(Up) , Vertical and radial joint probe (Down) 24 [1] Pijushpani Bhattacharjee, Soumini Chaudhury, Susmita Kundu, and Subhabrata Majumdar. Deriving the velocity distribution of Galactic dark matter particles from the rotation curve data. Phys. Rev. D, 87:083525, Apr 2013. [2] Julio F. Navarro, Carlos S. Frenk, and Simon D. M. White. A universal density profile from hierarchical clustering. The Astrophysical Journal, 490(2):493, 1997. [3] J. A. R. Caldwell and J. P. Ostriker. The mass distribution within our Galaxy - A three component model. The Astrophysical Journal, 251:61-87, December 1981. [4] K. Kuijken and G. Gilmore. The mass distribution in the galactic disc. I - A technique to determine the integral surface mass density of the disc near the sun. Monthly Notices of the Royal Astronomical Society, 239:571-603, August 1989. [5] K. Kuijken and G. Gilmore. The Mass Distribution in the Galactic Disc - II - Determination of the Surface Mass Density of the Galactic Disc Near the Sun. Monthly Notices of the Royal Astronomical Society, 239:605-649, August 1989. [6] K. Kuijken and G. Gilmore. The Mass Distribution in the Galactic Disc - Part III - the Local Volume Mass Density. Monthly Notices of the Royal Astronomical Society, 239:651-664, August 1989. [7] K. Kuijken and G. Gilmore. The galactic disk surface mass density and the Galactic force K(z) at Z = 1.1 kiloparsecs. The Astrophysical Journal Letters, 367:L9-L13, January 1991. [8] H. T. Freudenreich. A cobe model of the galactic bar and disk. The Astrophysical Journal, 492(2):495, 1998. [9] Pijushpani Bhattacharjee, Soumini Chaudhury, and Susmita Kundu. Rotation curve of the milky way out to 200 kpc. The Astrophysical Journal, 785(1):63, 2014. [10] Yoshiaki Sofue. Rotation curve and mass distribution in the galactic center – from black hole to entire galaxy. Publications of the Astronomical Society of Japan, 65(6), 2013. [11] Weber, M. and de Boer, W. Determination of the local dark matter density in our galaxy. Astronomy and Astrophysics, 509:A25, 2010. [12] Johan Holmberg and Chris Flynn. The local surface density of disc matter mapped by hipparcos. Monthly Notices of the Royal Astronomical Society, 352(2):440-446, 2004. [13] Silvia Garbari, Chao Liu, Justin I. Read, and George Lake. A new determination of the local dark matter density from the kinematics of k dwarfs. Monthly Notices of the Royal Astronomical Society, 425(2):1445-1458, 2012. [14] Lan Zhang, Hans-Walter Rix, Glenn van de Ven, Jo Bovy, Chao Liu, and Gang Zhao. The gravitational potential near the sun from segue k-dwarf kinematics. The Astrophysical Journal, 772(2):108, 2013. [15] Weber, M. and de Boer, W. Determination of the local dark matter density in our galaxy. Astronomy and Astrophysics, 509:A25, 2010. [16] Riccardo Catena and Piero Ullio. A novel determination of the local dark matter density. Journal of Cosmology and Astroparticle Physics, 2010(08):004, 2010. [17] Jo Bovy and Scott Tremaine. On the local dark matter density. The Astrophysical Journal, 756(1):89, 2012. [18] Moni Bidin, C., Carraro, G., M´endez, R. A., and Smith, R. Kinematical and Chemical Vertical Structure of the Galactic Thick Disk: Implications for the Thick Disk Formation. The Astrophysical Journal, vol. 751, no. 1, 2012, p. 30, doi:10.1088/0004-637X/751/1/30. [19] Ding, Ping-Jie, et al. Vertical Structure of Galactic Disk Kinematics from LAMOST K Giants. The Astronomical Journal, vol. 162, no. 3, 2021, p. 112, doi:10.3847/1538-3881/ac0892. [20] Sylos Labini, Francesco, et al. Mass Models of the Milky Way and Estimation of Its Mass from the Gaia DR3 Data Set. The Astrophysical Journal, vol. 944, no. 1, 2023, p. 8, doi:10.3847/1538-4357/acb6bd. [21] Bhattacharjee, Pijushpani, et al. Rotation Curve of the Milky Way out to 200 kpc. The Astrophysical Journal, vol. 785, no. 1, 2014, p. 63, doi:10.1088/0004-637X/785/1/63. 25 [22]Karukes, Ekaterina V., et al. ”A robust estimate of the Milky Way mass from rotation curve data.” arXiv, 2019, https://arxiv.org/abs/1912.04296. doi:10.48550/arXiv.1912.04296. [23] M. C. Smith et al., The RAVE Survey: Constraining the Local Galactic Escape Speed, Mon. Not. Roy. Astron. Soc. 379 (2007) 755–772, [astro-ph/0611671]. [24] G. Eadie and M. Juri´c, The Cumulative Mass Profile of the Milky Way as Determined by Globular Cluster Kinematics from Gaia DR2, The Astrophysical Journal 875 (Apr, 2019) 159, [1810.10036]. [25] L. Posti and A. Helmi, Mass and shape of the Milky Way’s dark matter halo with globular clusters from Gaia and Hubble, Astron. Astrophys. 621 (Jan, 2019) A56, [1805.01408]. [26] Z.-Z. Li, Y.-Z. Qian, J. Han, T. S. Li, W. Wang and Y. P. Jing, Constraining the Milky Way Mass Profile with Phase-Space Distribution of Satellite Galaxies, arXiv e-prints (Dec, 2019) arXiv:1912.02086, [1912.02086]. [27] A. A. Williams and N. W. Evans, Haloes light and dark: dynamical models of the stellar halo and constraints on the mass of the Galaxy, Monthly Notices of the Royal Astronomical Society 454 (09, 2015) 698–707, [http://oup.prod.sis.lan/ mnras/article-pdf/454/1/698/3923253/stv1967.pdf]. [28] W. Dehnen, D. McLaughlin and J. Sachania, The velocity dispersion and mass profile of the milky way, Mon. Not. Roy. Astron. Soc. 369 (2006) 1688–1692, [astro-ph/0603825]. [29] E. Vasiliev, Proper motions and dynamics of the Milky Way globular cluster system from Gaia DR2, MNRAS 484 (Apr, 2019) 2832–2850, [1807.09775]. [30] L. L. Watkins, R. P. van der Marel, S. T. Sohn and N. W. Evans, Evidence for an Intermediate-mass Milky Way from Gaia DR2 Halo Globular Cluster Motions, The Astrophysical Journal 873 (Mar, 2019) 118, [1804.11348]. [31] S. T. Sohn, L. L. Watkins, M. A. Fardal, R. P. van der Marel, A. J. Deason, G. Besla et al., Absolute hubble space telescope proper motion (HSTPROMO) of distant milky way globular clusters: Galactocentric space velocities and the milky way mass, The Astrophysical Journal 862 (Jul, 2018) 52. [32] A. A. Williams, V. Belokurov, A. R. Casey and N. W. Evans, On the run: mapping the escape speed across the Galaxy with SDSS, MNRAS 468 (Jun, 2017) 2359–2371, [1701.01444]. [33] K. Malhan and R. A. Ibata, Constraining the Milky Way halo potential with the GD-1 stellar stream, MNRAS 486 (Jul, 2019) 2995–3005, [1807.05994]. [34] C. S. Kochanek, The Mass of the Milky Way galaxy, Astrophys. J. 457 (1996) 228, [astro-ph/9505068]. [35] M. I. Wilkinson and N. W. Evans, The present and future mass of the Milky Way halo, Mon. Not. Roy. Astron. Soc. 310 (1999) 645, [astro-ph/9906197]. [36] T. K. Fritz, A. Di Cintio, G. Battaglia, C. Brook and S. Taibi, The mass of our Galaxy from satellite proper motions in the Gaia era, arXiv e-prints (Jan., 2020) arXiv:2001.02651, [2001.02651]. [37] D. Zaritsky et al., Globular Cluster Streams as Galactic High-Precision Scales—The Poster Child Palomar 5, The Astrophysical Journal 803 (Apr, 2015) 80. [38] S. L. J. Gibbons, V. Belokurov and N. W. Evans, ‘Skinny Milky Way please,’ says Sagittarius, Monthly Notices of the Royal Astronomical Society 445 (10, 2014) 3788–3802, [http://oup.prod.sis.lan/mnras/article-pdf/445/4/3788/6087753/ stu1986.pdf]. [39] F. Nesti and P. Salucci, The Dark Matter halo of the Milky Way, AD 2013, JCAP 1307 (2013) 016, [1304.5127]. [40] P. J. McMillan, Mass models of the Milky Way, Mon. Not. Roy. Astron. Soc. 414 (2011) 2446–2457, [1102.4340]. [41] G. M. Eadie and W. E. Harris, Bayesian Mass Estimates of the Milky Way: The Dark and Light Sides of Parameter Assumptions, The Astrophysical Journal 829 (Oct, 2016) 108, [1608.04757]. [42] P. R. Kafle et al., Kinematics of the Stellar Halo and the Mass Distribution of the Milky Way Using Blue Horizontal Branch Stars, The Astrophysical Journal 761 (Dec., 2012) 98, [1210.7527]. [43] A. J. Deason et al., Broken Degeneracies: The Rotation Curve and Velocity Anisotropy of the Milky Way Halo, Mon. Not. Roy. Astron. Soc. 424 (2012) L44–L48, [1204.5189]. [44] O. Y. Gnedin et al., The Mass Profile of the Galaxy to 80 kpc, Astrophys. J. 720 (2010) L108, [1005.2619]. 45D. Zaritsky, C. Conroy, H. Zhang, A. Bonaca, N. Caldwell, P. A. Cargile et al., A Lower Limit on the Mass of Our Galaxy from the H3 Survey, arXiv e-prints (Sep, 2019) arXiv:1909.02025, [1909.02025].