Dynamical modelling of the Galactic bulge and bar:
the Milky Way’s bar pattern speed, stellar, and dark matter mass distribution
Abstract
We construct a large set of dynamical models of the galactic bulge, bar and inner disk using the MadetoMeasure method. Our models are constrained to match the red clump giant density from a combination of the VVV, UKIDSS and 2MASS infrared surveys together with stellar kinematics in the bulge from the BRAVA and OGLE surveys, and in the entire bar region from the ARGOS survey. We are able to recover the bar pattern speed and the stellar and dark matter mass distributions in the bar region, thus recovering the entire galactic effective potential. We find a bar pattern speed of 39.0\pm 3.5\,\rm{km\,s^{1}\,kpc^{1}}, placing the bar corotation radius at 6.1\pm 0.5\rm{kpc} and making the Milky Way bar a typical fast rotator. We evaluate the stellar mass of the long bar and bulge structure to be M_{\rm{bar/bulge}}=1.88\pm 0.12\times 10^{10}\,\rm{M}_{\odot}, larger than the mass of disk in the bar region, M_{\rm{inner\ disk}}=1.29\pm 0.12\times 10^{10}\,\rm{M}_{\odot}. The total dynamical mass in the bulge volume is 1.85\pm 0.05\times 10^{10}\,\rm{M}_{\odot}. Thanks to more extended kinematic data sets and recent measurement of the bulge IMF our models have a low dark matter fraction in the bulge of 17\%\pm 2\%. We find a dark matter density profile which flattens to a shallow cusp or core in the bulge region. Finally, we find dynamical evidence for an extra central mass of \sim 0.2\times 10^{10}\,\rm{M}_{\odot}, probably in a nuclear disk or disky pseudobulge.
keywords:
methods: numerical – Galaxy: bulge – Galaxy: kinematics and dynamics – Galaxy: structure – Galaxy: centrepdfpagelabels=false, pageanchor=false, hyperfootnotes=false, plainpages=false,hypertexnames=falsehyperref
1 Introduction
Although it is well established that the Milky Way hosts a central barred bulge which causes nonaxisymmetric gas flow (Peters, III, 1975; Binney et al., 1991) and asymmetries in the nearinfrared light (Blitz & Spergel, 1991; Weiland et al., 1994) and star counts (Nakada et al., 1991; Stanek et al., 1997); our understanding of this structure has dramatically improved in the last decade. The discovery of the socalled split red clump in the galactic bulge (Nataf et al., 2010; McWilliam & Zoccali, 2010) and the later 3D mapping of the bulge density by Wegg & Gerhard (2013) showed that the galactic bulge has a boxy/peanut (B/P) shape, similarly to bulges formed in Nbody models by buckling of a vertically unstable stellar bar (Combes et al., 1990; MartinezValpuesta et al., 2006).
The existence of the bar outside of the bulge initially revealed by Hammersley et al. (1994) has been subject to controversy as to whether it is a separate structure from the bulge. First studies indicated a misalignment between the bulge and the bar, leading to the hypothesis of a double bar system in the inner Milky Way (Benjamin et al. 2005; LópezCorredoira et al. 2005; CabreraLavers et al. 2008; but see also MartinezValpuesta & Gerhard 2011). Recently, Wegg et al. (2015, hereafter W15) demonstrated by combining the VVV, UKIDSS, GLIMPSE and 2MASS catalogues that the galactic bulge smoothly segues into the long bar. Both components appear at a similar angle, showing that the galactic bulge and the long bar in the Milky Way are consistent with being a single structure that became vertically thick in its inner part, similarly to the buckled bars of Nbody models.
In Portail et al. (2015a, hereafter P15), we constructed dynamical models of the galactic bulge by combining the 3D density of Red Clump Giants (RCGs) in the bulge from Wegg & Gerhard (2013) with bulge kinematics from the BRAVA survey (Rich et al., 2007; Kunder et al., 2012) using the MadetoMeasure method. In this paper, we extend the MadetoMeasure modelling to the entire bar region by taking advantage of the recent measurement of W15 on the bar outside the bulge together with stellar kinematics in the bar region from the ARGOS survey (Freeman et al., 2013; Ness et al., 2013). The goal is to combine stellar density and kinematics under the constraint of dynamical equilibrium in order to recover the effective potential in the bar region, i.e. the stellar and dark matter mass distribution together with the bar pattern speed.
Extending the modelling from the bulge to the entire bar region is not a straightforward task. The galactic long bar extends about 5\,\rm{kpc} from the Galactic Centre (GC) as shown by W15 and thus reaches radii where the outer disk contribution to the potential is important. Hence, modelling the bar region also requires modelling the disk potential into which the bar is embedded. In addition, the dark matter contribution to the radial force increases with galactocentric radius reaching about 50\% at the solar radius (Read, 2014). As a consequence, the global mass distribution of Galaxy has to be taken into account in order to produce a good model of the galactic bar region. The building of good initial conditions for the MadetoMeasure modelling is particularly challenging.
The paper is organized as follows. In Section 2, we briefly describe the MadetoMeasure method and the problem posed by the initial conditions. In Section 3, we construct a static density model of the inner 10\,\rm{kpc} for the Galaxy by combining the current knowledge of the bulge, bar, disk and dark matter density. This density model is used in Section 4 to tailor a set of Nbody models with different bar pattern speeds that already broadly match the MilkyWay mass distribution. In Section 5, we discuss the different data sets used in this study to constrain the models and summarize the modelling procedure in Section 6. In Section 7, we show the effect of the main modelling parameters on the bulge dynamics, and in Section 8 we analyse a large number of models, recover the bar pattern speed of the Milky Way and identify our best model. In Section 9, we summarize our constraints on the stellar and dark matter mass distribution that arises from our models, and discuss our results in the light of other works in Section 10. We finally conclude in Section 11. The impatient reader can read first Sections 8 and 9 where we use our modelling to recover the effective potential in the Galaxy. Our best model is the first nonparametric model of the entire bar region and may be made available upon request to the authors.
2 MadetoMeasure modelling of the Galaxy
2.1 M2M modelling
Stellar dynamical equilibria for galaxies can be studied with momentbased methods (Binney et al., 1990; Cappellari et al., 2009), classical distribution functionbased methods (Dejonghe, 1984; Qian et al., 1995), actionsbased methods (Binney, 2010; Sanders & Binney, 2013), orbitbased methods (Schwarzschild, 1979; Thomas et al., 2009) or with the MadetoMeasure (M2M) method (Syer & Tremaine, 1996; De Lorenzi et al., 2007). In these M2M models, an initial selfgravitating Nbody model is used to provide a discrete sample of a distribution function reasonably close to the system of interest. This Nbody model is then slowly adapted by modifying the weights of the Nbody particles such as to make the model reproduce a given set of constraints. The Nbody weights can hence be seen simultaneously as mass elements (Nbody point of view), or as weights for the orbits traced by the particles (Schwarzschild’s method point of view). The method was extended to allow the fitting of observational data and implemented as the NMAGIC code by De Lorenzi et al. (2007). The M2M method has been used in both extragalactic (e.g. De Lorenzi et al., 2008, 2009; Das et al., 2011; Morganti et al., 2013; Zhu et al., 2014) and Galactic context (e.g. Bissantz et al. 2004; Long et al. 2013; Hunt & Kawata 2014; P15). We heavily modified the NMAGIC code for the purpose of modelling barred disk galaxies.
Formally, the M2M method works as follows. Any observable y of a system can be written in terms of the distribution function f(\bm{z}) of the system by
y=\int K(\bm{z})f(\bm{z})\,d^{6}\bm{z}  (1) 
where K is the kernel corresponding to the observable and \bm{z} the phasespace vector. In the M2M method, f(\bm{z}) is discretely sampled via a set of N particles with particle weights w_{i}(t). Equation 1 is then evaluated by
y(t)=\sum_{i=1}^{N}K(\bm{z}_{i}(t))w_{i}(t)  (2) 
where \bm{z}_{i}(t) is the phasespace coordinate of particle i at time t.
The M2M method consists of adjusting the particle weights w_{i} in order to maximize a given profit function F. This is achieved by a simple gradient descent in which the particle weights are evolved with time according to
\frac{dw_{i}}{dt}=\varepsilon w_{i}\frac{\partial F}{\partial w_{i}}  (3) 
and \varepsilon is a numerical factor that sets the typical timescale of the weight evolution. The profit function F usually consists of a chisquare term, which drives the model towards the data, and an entropy term to regularize the particle model. We describe the M2M formalism in more detail in Section 6.1.
Note that the M2M method only weights particles and does not have the ability to create new Nbody orbits. It is thus a very efficient modelling technique provided all the Nbody orbits required to fit the data are already present in the initial model. Building good initial conditions, including its dark matter component, is a major issue when modelling the inner 10\,\rm{kpc} of the Milky Way.
2.2 The problem of the initial conditions
To model the galactic bar, we need an initial Nbody model of a barred stellar disk that provides a firstguess discrete sampling of the final model. The classical way to build such initial conditions is to evolve a nearequilibrium Nbody stellar disk in a live dark matter halo. During the evolution the disk becomes unstable and forms a bar, that later forms a B/P bulge thought the buckling instability mechanism (Combes et al., 1990; MartinezValpuesta et al., 2006). This process has been used widely in order to build models that can then be compared to data (MartinezValpuesta et al., 2006; Athanassoula, 2007; Shen et al., 2010) but suffers from three major limitations:

The evolution process is nonlinear and very sensitive to initial conditions. As noted by (Sellwood & Debattista, 2009) large changes in the evolved bar model can result from simply changing the seed of the random number generator used in building the initial conditions.

The bar properties such as pattern speed, bar length and bar strength cannot be easily controlled a priori.
An example of such a buckled bar Nbody model is the model M85 of P15 shown in Figure 1. When scaled to a bar halflength of 5\,\rm{kpc}, M85 has a very short disk scalelength of only \sim 1.2\,\rm{kpc}, resulting in a very low contribution of the outer disk to the potential at the end of the bar region together with a lack of mass on long bar orbits. M85 was good enough to model only the bulge in P15 but does not suit our purpose here to model the inner 10\,\rm{kpc} of the Galaxy. To get around the three limitations of pure Nbody evolution shown above, we use a variant the M2M method to tailor initial conditions in a twostep process. In Section 3, we first create a mass density model of the Milky Way by adding together our bestguess densities for the bulge, bar, disk and dark halo. This density is then imprinted on M85 using a variant of the M2M method in Section 4, for different bar pattern speeds. At the end of this process, we obtain a family of Nbody models with different bar pattern speeds that have broadly the right mass distribution and provide suitable initial conditions for modelling the inner 10\,\rm{kpc} of the Milky Way by fitting real data.
3 Density model of the Galaxy for tailoring initial conditions
In this section, we construct a mass density model of the inner 10\,\rm{kpc} of the Galaxy by combining the 3D densities of the B/P bulge, bar, outer disk and dark matter halo. This density model is only used to tailor our initial conditions for the actual M2M modelling which is performed in Section 6. Throughout the paper, we place the Sun in the galactic plane at a distance R_{0}=8.2\,\rm{kpc} (BlandHawthorn & Gerhard, 2016) from the GC. The bar is oriented at an angle of \alpha=28\lx@math@degree with the SunGC line of sight, consistent with the measurement of 27\lx@math@degree\pm 2\lx@math@degree from the bulge RCGs (Wegg & Gerhard, 2013) and the range 28\lx@math@degree33\lx@math@degree measured by W15 from the long bar RCGs. Following BlandHawthorn & Gerhard (2016), we assume that the local standard at rest (LSR) is on a circular orbit at V(R_{0})=238\,\rm{km\,s^{1}}, and a peculiar motion of the Sun in the LSR of (U,V,W)=(11.1,12.24,7.25)\,\rm{km\,s^{1}} (Schönrich et al., 2010). This set of assumption predicts a solar tangential velocity of 250\,\rm{km\,s^{1}}, in good agreement with several recent measurements of 248\pm 6\,\rm{km\,s^{1}} (Schönrich, 2012, from SEGUE data), 242^{+10}_{3}\,\rm{km\,s^{1}} (Bovy et al., 2012, from APOGEE data), 244\pm 5\,\rm{km\,s^{1}} (Sharma et al., 2014, from RAVE data) and 251\pm 5\,\rm{km\,s^{1}} (Reid et al., 2014, from maser velocities). All parameters of the density model are summarized in Table 1.
3.1 Unified bulge and bar structure as traced by RCGs
Recently, Wegg & Gerhard (2013) constructed the first nonparametric measurement of the 3D density of RCGs in the bulge. They took advantage of the narrow luminosity function of RCGs to directly deconvolve the extinction and completenesscorrected magnitude distributions of bulge stars from the VVV survey and obtain line of sight densities of RCGs. Combining the different lines of sight and assuming eightfold symmetry, they produced a 3D density map of RCGs in a box of (\pm 2.2\times\pm 1.4\times\pm 1.2)\,\rm{kpc} around the principal axes of the bulge. This map together with the BRAVA kinematics was later used in P15 to construct a family of dynamical models of the galactic bulge. As bulge component, we adopt here the 3D density of the fitted model M85 of P15 that reproduces the original RCG density very well, with the advantage of being smooth and complete in the plane where the direct measurement of the density was not possible because of extinction and crowding.
Outside the bulge, W15 combined the VVV, UKIDSS and 2MASS surveys and showed that the bulge smoothly segues from its vertically extended B/P shape to the flat long bar. The bulge and long bar are shown in this later work to be consistent with forming a single structure, oriented at \alpha=28\lx@math@degree from the SunGC line of sight. They estimated the long bar halflength to be 5\,\rm{kpc} and found evidence for an extra superthin bar component existing predominantly near the bar end. They finally fit a parametric model of the long bar density that once added to the fitted bulge model M85 of P15 and convolved with the bulge luminosity function fits well the magnitude distribution of stars across the entire bulge and bar region. Consequently, we complement the bulge model described above using their bestfitting parametric densities of the thin long bar and superthin components. Note that due to their analysis method, the long bar density of W15 does not include the inner disk, smooth background of stars filling the bar region. We add the inner disk density in Section 3.3.
Both the bulge and long bar density were measured using RCGs as tracers. Theoretical models by Salaris & Girardi (2002) show that for a 10\,\rm{Gyr} old stellar population, RCGs trace the stellar mass within 10% for metallicities in the range 1.5\leq[\rm{Fe/H}]\leq 0.2. In the particular case of the galactic bulge and bar, Ness et al. (2013) find from the ARGOS sample that 95\% of the stars enter this metallicity range. The age of the bulge is still under debate with pieces of contradictory evidence: photometric studies of the colormagnitude diagram (CMD; Zoccali et al., 2003; Clarkson et al., 2008; Calamida et al., 2015) find that the bulge is older than 10\,\rm{Gyr} while spectroscopic age measurements of microlensed dwarfs find evidence for 45\,\rm{Gyr} old population among stars with [\rm{Fe/H}]\geq0.1 (Bensby et al., 2011, 2013). We assume here that the bulge and bar are 10\,\rm{Gyr} old, implying that the RCG density of the bulge and long bar considered above are proportional to the stellar density with expected variations of less than 10\%. We call the proportionality factor between stellar mass density and RCG density masstoclump ratio ^{1}^{1}1Note that, as in P15 our definition of the masstoclump also includes the red giant branch bump stars, the number of RCGs + red giant branch bump stars is better defined than the number of RCGs only (see Wegg & Gerhard 2013)., denoted as \rm{M/n_{RCG}} by analogy with the mass tolightratio. We make the fiducial assumption that the bulge and the thin long bar have the same masstoclump ratio, as expected if they both formed at the same time. The origin of the superthin bar is still unclear as stated by W15 but its stellar population is expected to be younger given its extremely short scaleheight (45\,\rm{pc}), and therefore it is likely to have a lower masstoclump ratio than the bulge. Assuming a constant star formation rate and a Kroupa initial mass function (IMF) as in W15, the masstoclump ratio of the superthin bar is a factor of 1.6 times smaller than that of a 10\,\rm{Gyr} old bulge.
3.2 Empirical determination of the masstoclump ratio in the bulge
The masstoclump ratio can be predicted by stellar population synthesis models using an IMF, a stellar age distribution and a metallicity distribution as in P15. The most recent measurement of the galactic bulge IMF is from Calamida et al. (2015) who used ultradeep Hubble Space Telescope (HST) photometry to recover the IMF in the mass range 0.15\leq M/M_{\odot}\leq 1. They find an IMF in good agreement with the Kroupa IMF (Kroupa, 2001), for which we computed in P15 a masstoclump ratio of \rm{M/n_{RCG}}=984 for a 10\,\rm{Gyr} old population.
An alternative and more direct approach is to combine stellar mass measurements in some bulge field with the observed number of RCGs in that field, in analogy with the method of Valenti et al. (2016). This approach is advantageous since it does not rely on stellar population models or parametrization of the IMF. Zoccali et al. (2000) used HST photometry in the NICMOS field and after cleaning for disk contamination they find a stellar mass in that field of M_{\rm NICMOS}=570\,\rm{M}_{\odot} (see the revision of the mass in the NICMOS field in Valenti et al. 2016). The NICMOS field has an area of only 408 square arcseconds and does not contain many giant stars. To improve the statistics on the number of RCGs, we follow Valenti et al. (2016) and consider a larger 15\arcmin beam centred on the NICMOS field and rescale the mass and number of RCGs by the ratio of the area of the two fields. We use the completeness and extinctioncorrected VVV catalogue of Wegg & Gerhard (2013) and identify RCGs statistically as the excess above the smooth background of stars in the extinction corrected magnitude distribution. Figure 2 shows the magnitude distributions of VVV stars in our larger field centred on the positions of the NICMOS fields with the identified RCGs above the background of stars.
With this approach, we find a masstoclump ratio of \rm{M/n_{RCG}}=1015. We estimate the error on this figure of about 10%, mostly due to systematic effect arising in defining the smooth background of stars on to which the RCGs sit (see Wegg & Gerhard 2013). Note that the uncertainty on the lowmass end of the IMF due to unseen dwarfs has only a small effect on the masstoclump ratio. Variations from a substellar slope of 0.3 (Kroupa, 2001), as also adopted in the revised mass of the NICMOS field (Valenti et al., 2016), to either a slope of 1.33 (Zoccali et al., 2000) or a lognormal parametrization lead to variations of the masstoclump ratio of only \pm 3\%. Our direct measurement of the masstoclump ratio is in good agreement with the predicted value of 984 for a Kroupa or Calamida IMF. In all the following, we adopt the fiducial value of \rm{M/n_{RCG}}=1000 for the main stellar population in the bulge and bar together with a lower value of \rm{M/n_{RCG}}=600 for the superthin bar component. We show in Section 8 the effect of a 10% smaller or larger masstoclump ratio.
3.3 Stellar disk
Our prime interest in modelling the disk is to obtain a reasonable disk contribution to the potential in the bar region and disk foreground contamination when observing the bulge and bar. Outside the bar region (R\geq 5.5\,\rm{kpc}) we adopt an axisymmetric stellar disk structure with scalelength h_{R,*} and exponential scaleheight h_{Z,*}. From papers based on infrared data BlandHawthorn & Gerhard (2016) concluded that h_{R,*}=2.6\pm 0.5\,\rm{kpc} with the shorter disk scalelengths in the range 2.12.6\,\rm{kpc} usually favored by dynamical studies of stellar kinematics or microlensing optical depth towards the bulge (Wegg et al., 2016). We adopt here a fiducial scalelength of h_{R,*}=2.4\,\rm{kpc} and scaleheight of h_{Z,*}=300\,\rm{pc} (Jurić et al., 2008) and test in Section 9 the effect of a shorter scalelength of 2.15\,\rm{kpc} (Bovy & Rix, 2013) and 2.6\,\rm{kpc} (Jurić et al., 2008).
Following Bovy & Rix (2013), we add the interstellar medium contribution to the potential by modelling it as an additional thin disk with scalelength h_{R,\rm{ism}}=2\times h_{R,*} and scaleheight h_{R,\rm{ism}}=130\,\rm{pc}. The disks are normalized to a baryonic local surface density inside 1.1\,\rm{kpc} above and below the plane of \Sigma_{1.1}(R_{0})=51\,\,\rm{M}_{\odot}.pc^{2} among which 38\,\,\rm{M}_{\odot}.pc^{2} are stars and 13\,\,\rm{M}_{\odot}.pc^{2} are interstellar medium (Bovy & Rix, 2013).
Inside the bar region, very little is known about the structure of the disk component that surrounds the bar and bulge. We can fortunately constrain this inner disk by combining our knowledge of the bulge and of the disk at the boundary of the bar region (i.e. 5.5\,\rm{kpc}). In the bulge region, approximated as the interior of an ellipse reaching 2.2\,\rm{kpc} along the bar major axis and 1.4\,\rm{kpc} along the intermediate axis, the ‘bulge’ model already represents the entire stellar density. Hence the total surface density along the intermediate axis has to smoothly transition between its value at 5.5\,\rm{kpc} and the bulge value at 1.4\,\rm{kpc}. We construct the inner disk in the bar region by first interpolating the total surface density along the intermediate axis between 5.5 and 1.4\,\rm{kpc} using the logarithm of a quadratic Bézier curve interpolation, whose control points are defined to ensure continuity of the derivative. The disk surface density is then constructed assuming a linear decrease of the ellipticity of the disk isocontours between the bulge and the boundary of the bar region. This procedure is shown in Figure 3 where the bar and bulge model (left) and its additional disk component (centre) sum up to a total density that smoothly joins the bulge to the disk at the end of the bar region. Although several assumptions enter in joining the bulge, bar and disk as described above, it results in a reasonable global density model for the Milky Way that suits well our purpose of tailoring the initial conditions for the M2M modelling as already stated at the beginning of this section.
3.4 Dark matter halo
By adding up the bulge, bar and disk as described above, we obtain a 3D density model of the baryonic mass in the Milky Way.
Figure 4 shows the rotation curve V_{c}(r) of these baryonic mass models together with the composite rotation curve from Sofue et al. (2009) rescaled to a distance to the GC and a local circular velocity of (R_{0},V_{0})=(8.2\,\rm{kpc},238\,\rm{km\,s^{1}}) as described in Section 2.2. This baryonic model is insufficient to match the rotation curve of the Milky Way, and we therefore require dark matter within the framework of Newtonian dynamics. Recent studies of dark matter simulations such as Navarro et al. (2010) showed that the innermost regions of dark matter halos were better represented by the Einasto density profile (Einasto, 1965) than by the NFW profile. Inspired by this, we adopt the threeparameter Einasto density profile given by
\rho_{\rm{DM}}(m)=\rho_{0}\,{\rm exp}\left\{\left(\frac{2}{\alpha}\right)% \left[\left(\frac{m}{m_{0}}\right)^{\alpha}1\right]\right\}  (4) 
where m=\sqrt{x^{2}+y^{2}+(z/q)^{2}} is the elliptical radius for an assumed vertical flattening of q=0.8 (Piffl et al., 2014). We can constrain the halo parameters using the rotation curve but caution has to be taken inside the bar region. The data from Sofue et al. (2009) shown in Figure 4 are a combination of different data sets and rely mostly on the tangentpoint method from terminal velocity measurement of CO and \rm{H_{I}} gas inside the solar circle. Because of the influence of the bar on the gas flows, the tangent point method is likely to be flawed inside the bar region (Englmaier & Gerhard, 1999; Chemin et al., 2015) so we exclude data points inside 6\,\rm{kpc} from the GC. We also exclude all data points outside the solar circle, as the rescaling to our assumptions for (R_{0},V_{0}) would require us to take into account the nature of the different data sets entering the work of Sofue et al. (2009). The rotation curve between 6\,\rm{kpc} and R_{0} provides a good constraint on the average value and slope of the dark matter circular velocity at the solar position, but is not sufficient to constrain the dark matter in the inner region. The determination of the dark matter contribution in the inner Galaxy requires proper dynamical modelling and is addressed in detail in Section 7.1. At this stage we assume a dark matter mass inside 2\,\rm{kpc} of 0.5\times 10^{10}\,\rm{M}_{\odot}, resulting in the rotation curve shown in Figure 4, consistent with our bulge models in P15. In Section 7.2, we relax the constraint on the dark matter mass inside 2\,\rm{kpc} during the modelling process and adapt it directly to match the bulge kinematics.
A summary of the diverse parameters of our fiducial model for tailoring the initial conditions is given in Table 1.
Parameter  Fiducial value  Reference  Section  
Geometry  Distance to GC R_{0}  8.2\,\rm{kpc}  BlandHawthorn & Gerhard (2016)  Section 3 
Bar angle  28\lx@math@degree  Average between Wegg & Gerhard (2013) and W15  Section 3  
Bulge  3D density  Fitted M85  P15  Section 3.1 
Masstoclump ratio \rm{M/n_{RCG}}  1000  Direct measurement, Kroupa + Calamida IMF  Section 3.2  
Bar  Thin component  Analytical density  W15  Section 3.1 
Superthin component  Analytical density  W15  Section 3.1  
\rm{M/n_{RCG}} of thin bar  Same as bulge  –  Section 3.1  
\rm{M/n_{RCG}} of superthin bar  600  Stellar population models  Section 3.1  
Stellar disk  Scalelength h_{R,*}  2.4\,\rm{kpc}, middle of range 2.152.6\,\rm{kpc}  BlandHawthorn & Gerhard (2016)  Section 3.3 
Scaleheight h_{Z,*}  300\,\rm{pc}  Jurić et al. (2008)  Section 3.3  
Local surface density \Sigma_{1.1,*}(R_{0})  38\,\rm{M}_{\odot}.\,\rm{pc}^{2}  Bovy & Rix (2013)  Section 3.3  
ISM disk  Scalelength h_{R,\rm{ism}}  2\times h_{R,*}  Bovy & Rix (2013)  Section 3.3 
Scaleheight h_{Z,\rm{ism}}  130\,\rm{pc}  Bovy & Rix (2013)  Section 3.3  
Local surface density \Sigma_{1.1,\rm{ism}}(R_{0})  13\,\rm{M}_{\odot}.\,\rm{pc}^{2}  Bovy & Rix (2013)  Section 3.3  
Dark matter halo  Profile  Bestfitting Einasto  –  Section 3.4 
Flattening  0.8  Piffl et al. (2014)  Section 3.4  
Outer constraints  V_{c}(6\,\rm{kpc}\leq r\leq R_{0}\,\rm{kpc})  Sofue et al. (2009)  Section 3.4  
Inner constraint  M(<2\,\rm{kpc})=0.5\times 10^{10}\,\,\rm{M}_{\odot}  P15  Section 3.4  
Dynamical parameters  LSR circular velocity V_{0}  238\,\rm{km\,s^{1}}  Schönrich (2012), Reid et al. (2014)  Section 3.4 
Solar motion in the LSR (U,V,W)  (11.1,12.24,7.25)\,\rm{km\,s^{1}}  Schönrich et al. (2010)  Section 3  
Bar pattern speed \Omega_{b}  Systematic search in the range 2550\,\rm{km\,s^{1}\,kpc^{1}}  –  Section 3 
4 Tailoring initial models for modelling the Milky Way
In this section, we use a variant of the M2M method to create a family of Milky Way models with a specified mass distribution and different bar pattern speeds. We adiabatically adapt M85 to the density model of the Galaxy described above and change slowly the bar pattern speed, hence gaining full control on the effective potential of the model.
4.1 Adiabatic adaptation of the initial conditions
We first evaluate the initial stellar and dark matter mass distribution of M85 on respectively a Cartesian grid of \pm 12\,\rm{kpc}\times\pm 12\,\rm{kpc}\times\pm 2\,\rm{kpc} and a radial grid extending to 40\,\rm{kpc} with flattening q=0.8. We then integrate our fiducial stellar and dark matter density target from Section 3 in the grid cells and define 25 intermediate targets, logspaced between the initial model and the target mass distribution. We modify the bar pattern speed by defining 25 intermediate corotation radii, linearly spaced between the initial model corotation and the target corotation.
The adiabatic adaptation of M85 to the target model then consists of 25 iterations of the following procedure:

Perform an M2M fit of the stellar and dark matter mass distributions to the current target mass distribution while evolving the model in the current rotating potential at constant pattern speed. We run this M2M fit for one time unit, corresponding to the period of one circular orbit at 4\,\rm{kpc}.

Update the potential to the new particle masses and adapt the pattern speed at which the potential rotates to place the corotation at the next intermediate target value.

Multiply all particle velocities by a factor \sqrt{\bm{r}\cdot\bm{\nabla}\Phi_{\rm new}}/\sqrt{\bm{r}\cdot\bm{\nabla}\Phi_{% \rm old}} where \bm{r} is the particle position and \Phi_{\rm old} and \Phi_{\rm new} are the potential respectively before and after the potential update. This step is necessary as the circular velocity in the new potential is different from that of the original model.
Each M2M fit to an intermediate target mass distribution is performed using Equation 3 with the following profit function
F=\frac{1}{2}\sum_{k,j}(\Delta_{j})^{2}  (5) 
where \Delta_{j} is the difference between the model mass and the target mass in cell j. Since model observables can be noisy in regions of space where the particle density is low, we use temporal smoothing by replacing the instantaneous model observable y(t) at time t by its temporally smoothed value \tilde{y}(t) defined as
\tilde{y}(t)=\int y(t\tau)e^{\alpha\tau}\,d\tau\leavevmode\nobreak\ .  (6) 
where 1/\alpha is the temporal smoothing timescale.
We then continue the M2M fit to the final target mass distribution and corotation radius for another 25 time units. The M2M fit can lead to a broad distribution of weights that has the undesired effect of lowering the model resolution and potentially introducing clumps in the potential due to very massive particles. We added to NMAGIC the particle resampling algorithm described in Dehnen (2009). This algorithm consists of creating a new particle model with equalweights particles by resampling the original set of particles with a probability proportional to their weights. When multiple selections of a given particle occur, we evolve the parent particle for one orbital time (estimated as T\sim 2\pi\sqrt{\frac{r}{f_{r}}} where r and f_{r} are the radius and radial acceleration of the particle) and select multiple particles along the trajectory in phasespace sampled by the parent particle.
Following this procedure, we create a set of Nbody models that broadly matches the static density model of Section 3 with bar pattern speeds in the range 2550\,\rm{km\,s^{1}\,kpc^{1}}. All models have 10^{6} equal weight stellar particles and 10^{6} dark matter particles. Figure 5 shows the surface densities of three Nbody models with pattern speeds of 25, 35 and 45\,\rm{km\,s^{1}\,kpc^{1}}. As expected, different pattern speeds lead to slightly different bar shapes but the global mass distribution of the target density is anyway well reproduced in these adiabatically adapted Nbody models.
4.2 Integration and potential solver
The particle model is integrated using a driftkickdrift adaptive leapfrog algorithm in the full gravitational potential rotating at a constant patten speed. The gravitational potential is computed directly from the particle mass distributions using the hybrid grids method described in the appendix of Sellwood (2003). This hybrid method combines a grid based potential solver on a flat cylindrical grid to evaluate the disk potential with a spherical harmonics potential solver on a spherical grid to evaluate the potential of the dark matter halo. For the cylindrical potential solver, we use the 3D polar grid code from Sellwood & Valluri (1997), in a cylindrical grid extending to 12\,\rm{kpc} in radius and \pm 2\,\rm{kpc} in the vertical direction. As a spherical solver, we use the spherical harmonics solver of De Lorenzi et al. (2007) up to order 8 on a spherical grid extending to 40\,\rm{kpc}.
In addition to the hybrid grid method, we modified the 3D polar grid code in order to allow the resolution of strong vertical gradients. In the original code from Sellwood & Valluri (1997), the particle mass distribution is softened using a spherical cubic spline density kernel where the softening scale should not be smaller than the planar scale of the grid cells in order to give accurate results. To keep the number of planar cells under control and still resolve strong vertical gradient we replace the spherical softening by an oblate softening with vertical axial ratios of 0.2. In the end, we define the grid parameters in order to obtain a planar resolution of 100\,\rm{pc} and a vertical resolution of \sim 20\,\rm{pc}.
5 NMAGIC data constraints for the galactic bulge, bar and disk
In the previous section, we built a set of initial Nbody models with different pattern speeds that already broadly matches the MilkyWay bulge, bar, and disk density. In this section, we describe the different data sets to which we fit our Nbody models in Section 6. For each data set k and observable j, we describe the NMAGIC kernels K_{j}^{k} to be used in Equation 2 for applying observational selection bias to the particles when observing the Nbody models. An overview of the spatial coverage of the different datasets described in this section is plotted in Figure 6.
5.1 Density and kinematics of the inner Galaxy from the ARGOS survey
The Abundance and Radial velocity Galactic Origin Survey (ARGOS) is a large spectroscopic survey of about 28000 stars of the galactic bulge and inner disc. It was designed to sample RCGs all the way from the near disk (\sim 4.5\,\rm{kpc} from the Sun) to the far side of the GC (\sim 13\,\rm{kpc} from the Sun). From the mediumresolution spectra, Ness et al. (2013) estimated various stellar parameters including radial velocity, stellar temperature and surface gravity. These parameters together with the intrinsically narrow luminosity function of RCGs allow the determination of relatively accurate distances. Hence, the ARGOS survey provides structural information about the inner disk, bulge and bar together with the radial velocity field in a wide spatial range extending up to l=20\lx@math@degree. In Section 5.1.1, we briefly review the selection strategy of the ARGOS survey and compute its selection function. In Section 5.1.2, we determine distances for all ARGOS stars and compute the mean velocity and velocity dispersion in several distance modulus bins in each field. Finally in Section 5.1.3 we describe the NMAGIC observable kernels that map the survey selection strategy.
5.1.1 The ARGOS selection function
The ARGOS stars were selected from the 2MASS pointsource catalogue (Skrutskie et al., 2006), according to a selection procedure fully described in Freeman et al. (2013). When computing the ARGOS selection function, the three following points need to be considered.

About 1000 stars are selected for each field. The sampling of a large distance range from \sim 4.5 to \sim 13\,\rm{kpc} is achieved by selecting randomly \sim 330 stars in three magnitude bins defined in the Iband (I\sim 1314, 1415 and 1516).

The 2MASS stars from which the ARGOS stars are selected in (i) do not correspond to the full pointsource 2MASS catalogue but only to a highquality subsample of it. This subsample is defined by a blue color cut (JK_{s})_{0}\geq 0.38 to remove foreground disk contamination, two magnitude cuts 11.5\leq K_{s}\leq 14, and additional criteria to exclude stars with large photometric errors, contamination, blends and low photometric quality. The subsample is biased towards the bright stars for which highquality imaging is easier to achieve.

The 2MASS catalogue is rapidly incomplete in crowded fields. This affects mostly the three fields at (l,b)=(0\lx@math@degree,5\lx@math@degree) and (\pm 5\lx@math@degree,5\lx@math@degree) where 2MASS is only \sim 45\% complete at K_{s}\sim 14.
We evaluate the incompleteness of the 2MASS survey by comparing with the deeper VVV survey where available (Saito et al., 2012), completeness corrected by Wegg & Gerhard (2013). In all the ARGOS fields at the edges of the VVV coverage, 2MASS and VVV do not deviate significantly from each other over the magnitude range 11.5\leq K_{s}\leq 14. Given that the effect of crowding decreases with increasing l and b, we consider that 2MASS is complete in all the ARGOS fields out of the VVV coverage area.
Extinction is evaluated on a starbystar basis using the RayleighJeans Color Excess method (Majewski et al. 2011; W15) given by:
A_{K_{s}}=\frac{A_{K_{s}}}{E(JK_{s})}[(JK_{s})(JK_{s})_{\rm RCG})]  (7) 
where (JK_{s})_{\rm RCG} is the colors of RCGs and \frac{A_{K_{s}}}{E(JK_{s})} is a constant that depends on the extinction law. For consistency with Wegg & Gerhard (2013), we adopt (JK_{s})_{\rm RCG}=0.674 (Gonzalez et al., 2011) and \frac{A_{K_{s}}}{E(JK_{s})}=0.528 (Nishiyama et al., 2006).
To take points (ii) and (iii) into account, we construct the selection function C({K_{s}}_{0}) that gives the probability for a star of extinctioncorrected magnitude {K_{s}}_{0} to belong to the highquality photometric 2MASS subsample from which the ARGOS stars were selected. This is evaluated empirically as shown in Figure 7. We finally define C({K_{s}}_{0}) as the fraction of stars with magnitude {K_{s}}_{0} from the completenesscorrected 2MASS catalogue that are also in the highquality 2MASS subsample.
To correct for the selection bias (i) we assign a weight w_{k} to each of the ARGOS stars, corresponding to the fraction of selected stars out of the number of stars from the 2MASS subsample of good photometric quality present in the considered Iband bin.
Figure 7 illustrates this selection procedure for the field at (l,b)=(0\lx@math@degree,5\lx@math@degree), the field most affected by selection effects.
5.1.2 Density and kinematics as a function of distance
Assuming that all ARGOS stars are RCGs and using the starbystar extinction of Equation 7, we estimate the distance moduli \mu_{K_{s}} of all ARGOS stars as
\mu_{K_{s}}=K_{s}A_{K_{s}}M_{K_{s},\rm{RCG}}  (8) 
where we adopt an absolute magnitude of RCGs of M_{K_{s},\rm{RCG}}=1.72 for consistency with Wegg & Gerhard (2013).
Stars in the ARGOS sample are either real RCGs, for which distances are accurately determined by Equation 8, or red giants that happen to be at the right distance to appear with similar color and magnitude as the real RCGs. For these giants, the distance estimate given by Equation 8 can be wrong by several magnitudes. To minimize their effect, we follow Ness et al. (2013) and evaluate the absolute magnitude M_{K_{s}} of every star using the surface gravity \log g measurement obtained from fitting the spectra together with the PARSEC isochrones (Bressan et al., 2012; Chen et al., 2014; Tang et al., 2014) and assuming a 10\,\rm{Gyr} old population, a Kroupa IMF and the overall metallicity distribution of all the ARGOS stars. We then statistically remove nonRCG stars by replacing the weights w_{k} of all ARGOS stars with w_{k}\times\omega(M_{K_{s}}), where \omega is a weighting function depending on the inferred absolute magnitude. The uncertainty in \log g is \sim 0.3\,\rm{mag}, which is equivalent to an uncertainty in M_{K_{s}} of 0.7\,\rm{mag}. Original work from the ARGOS team chose for \omega a tophat function around M_{K_{s},\rm{RCG}}, identifying in this way what they call the ‘probable RCGs’. In order to take advantage of the full sample of stars, we adopt instead the weighting \omega(M_{K_{s}})=G_{[M_{K_{s},\rm{RCG}},0.7]}(M_{K_{s}}) where G_{[\mu,\sigma]} is the Gaussian function of mean \mu and standard deviation \sigma.
Finally, we bin the distance modulus space in bins of 0.25\,\rm{mag} and compute for each field the number of stars n_{j,m}, the mean radial velocity \textrm{v}_{j,m} and radial velocity dispersion \sigma_{j,m} of the stars in field j and distance modulus bin m, taking into account each star weight w_{k}. Errors in those quantities are computed by 1000 bootstrap resamplings.
5.1.3 NMAGIC observable for the ARGOS data
In order to map the observational criteria to observing an Nbody model, we first need to study the stellar population that falls into the ARGOS sample in more detail. Using the PARSEC isochrones for a 10\,\rm{Gyr} old population with a Kroupa IMF and the overall metallicity distribution of all the ARGOS stars, we predict the luminosity function \Phi of the stars that matches the ARGOS color and magnitude cuts. The RCGs then appear as a sharp peak over a large background of red giants, located at M_{K_{s},\rm{RCG}}=1.47, slightly fainter than the assumed value of M_{K_{s},\rm{RCG}}=1.72 used by Wegg & Gerhard (2013) for measuring the 3D density of the galactic bulge. For internal consistency between our data sets we shift the isochrone luminosity function to agree with the maximum of the RCG peak at M_{K_{s},\rm{RCG}}=1.72.
Let us now observe an Nbody model by effectively turning the particles into stars and applying the selection procedure of the ARGOS survey. We consider a particle at distance modulus \mu_{i} and note f_{i}(\mu_{K_{s}}) the distribution of observed distances of mock stars drawn from \Phi(K_{s}) at the position of the considered particle when assuming that they are all RCGs. A mock star with an absolute magnitude M_{K_{s}} will be inferred to lie at a distance \mu_{K_{s}}=M_{K_{s}}+\mu_{i}M_{K_{s},\rm{RCG}}. f_{i}(\mu_{K_{s}}) can be written as
\begin{split}\displaystyle f_{i}(\mu_{K_{s}})=&\displaystyle\Phi(\mu_{K_{s}}% \mu_{i}+M_{K_{s},\rm{RCG}})\times\overline{\omega}(\mu_{K_{s}}\mu_{i}+M_{K_{s% },\rm{RCG}})\\ &\displaystyle\times C(\mu_{K_{s}}+M_{K_{s},\rm{RCG}})\end{split}  (9) 
where C(K_{s0}) is the selection function computed in Section 5.1.1. The mean weight \overline{\omega}(M_{K_{s}}), average of the weights obtained after statistical removal of nonRCGs from a mock measurement of M_{K_{s}}, with accuracy 0.7\,\rm{mag} is given by
\overline{\omega}(M_{K_{s}})=\int\omega(M)\times G_{[M_{K_{s}},0.7]}(M)\,dM.  (10) 
The distribution of distances f_{i}(\mu_{K_{s}}) resulting from Equation 9 is shown in Figure 8 for three particles at at \mu_{i}=13.5, 14.5 and 15.5. The distribution is narrow in all cases, showing that the ARGOS survey provides accurate distances. The spreading in distances due to the background of giants is minimized thanks to the selection bias towards nearby stars and also the extra information provided by the measurement of \log g.
As NMAGIC observables we adopt the number count and first and second massweighted velocity moments. The corresponding kernels for a field j and distance modulus m to be used in Equation 2 are given by
K_{j,m}^{\rm{ARGOS},0}=\delta_{j,m}^{\rm{ARGOS}}(\bm{z}_{i})  (11) 
K_{j,m}^{\rm{ARGOS},1}=\delta_{j,m}^{\rm{ARGOS}}(\bm{z}_{i})\times\frac{% \textrm{v}_{i}}{W_{j,m}^{\rm{ARGOS}}}  (12) 
and
K_{j,m}^{\rm{ARGOS},2}=\delta_{j,m}^{\rm{ARGOS}}(\bm{z}_{i})\times\frac{% \textrm{v}_{i}^{2}}{W_{j,m}^{\rm{ARGOS}}}  (13) 
where \textrm{v}_{i} is the radial velocity of particle i, W_{j,m}^{\rm{ARGOS}} is given by
W_{j,m}^{\rm{ARGOS}}=\sum_{i}w_{i}\delta_{j,m}^{\rm{ARGOS}}(\bm{z}_{i})  (14) 
and \delta_{j,m}^{\rm{ARGOS}} by
\delta_{j,m}^{\rm{ARGOS}}(\bm{z}_{i})=\begin{cases}\int_{\mu(m)}^{\mu(m+1)}f_{% i}(\mu)\,d\mu&\text{if $i\in$ field $j$, }\\ 0&\text{otherwise}\end{cases}  (15) 
with \mu(m) and \mu(m+1) the boundaries of the distance modulus bin m.
5.2 Magnitude distribution of the bulge and bar
In the region delimited by l\leq 40\lx@math@degree and b\leq 9\lx@math@degree, we use the combined catalogue of the VVV, UKIDSS and 2MASS surveys from W15. For each line of sight, this catalogue consists of histograms of distance moduli of stars \mu_{K_{s}} defined as in Equation 8. Each of these histograms shows an exponential background distribution of stars plus an overdensity of stars due to the RCGs located in the bulge or bar, as shown in Figure 2 for one bulge field. Information on the density is very hard to extract from the exponential background distribution of stars as it arises from the convolution of the line of sight density with the luminosity function of giant and dwarf stars that is very broad, poorly known and likely to vary from field to field and along the line of sight. Hence, we restrict our use of the histograms to the fields that show a significant excess of stars above the exponential background. To do so, we follow W15 and first fit a Gaussian plus an exponential to the distribution of distance moduli, separately in each field. We remove from further considerations all lines of sight where either the RCG bump is not detected to at least three sigma or the exponential background slope is too small, indicating incompleteness (see W15). The spatial coverage of the remaining fields is plotted in Figure 6.
The kernels of our model observables in a field j and distance modulus bin m bounded by \mu(m) and \mu(m+1) are given by
K_{j,m}^{\rm{hist}}(\bm{z}_{i})=\begin{cases}\frac{1}{\rm{M/n_{RCG}}}\times% \int_{\mu(m)}^{\mu(m+1)}\Phi(\mu\mu_{i})\,d\mu&\text{if $i\in$ field $j$,}\\ 0&\text{otherwise}\end{cases}  (16) 
where \rm{M/n_{RCG}} is the masstoclump ratio described in Section 3.2 and \Phi the luminosity function for RCGs only, expressed in distance modulus (see equation 17 in W15). The exponential background of stars, absent from the model observables, needs to be introduced before comparing model to data. To do so, we follow W15 and first fit an exponential distribution y_{j}^{e} to the difference between the model observable y_{j}^{\rm{hist}} and full data histogram Y_{j}^{\rm{hist}}. This exponential background of stars is then included in the datamodel comparison by replacing y_{j}^{\rm{hist}} by y_{j}^{\rm{hist}}+y_{j}^{e} in Equation 24.
For all fields with b\geq 1.35\lx@math@degree, we assume a 10\,\rm{Gyr} old stellar population with a Kroupa IMF and the metallicity distribution of Zoccali et al. (2000) in Baade’s window and compute \Phi from the PARSEC isochrones. After removing the exponential background we scale \Phi to our fiducial masstoclump ratio of 1000 (see Section 3.2).
For b<1.35\lx@math@degree, as already discussed in Section 3.1, W15 found evidence for a superthin bar component, present mostly near the bar end. This superthin bar is likely to be formed by younger stars as indicated by its extremely small scaleheight of only 45\,\rm{pc}. Detailed modelling of the superthin component is beyond the scope of this paper, but as our goal is to constrain the gravitational potential, its mass has to be included in the modelling. We follow W15 and assume a constant star formation rate for the superthin component to compute its luminosity function. The superthin bar population has a masstoclump ratio of 600, lower than the old bar population. For fields with b\leq 1.35\lx@math@degree we use a superposition of the thin bar and superthin bar populations, using the ratio of the densities given by the parametric models of W15.
For efficiency we combine the original data of W15 in cells of 2\lx@math@degree\times 1.2\lx@math@degree in galactic coordinates for b\geq 1.35\lx@math@degree and cells of 2\lx@math@degree\times 0.6\lx@math@degree for b<1.35\lx@math@degree and symmetrize with respect to the b=0\lx@math@degree.
5.3 3D density of the bulge and outer disk
In the bulge we constrain the stellar density using the 3D density of RCGs measured by Wegg & Gerhard (2013), scaled to stellar mass density using our fiducial masstoclump ratio. The map covers a box of (\pm 2.2\times\pm 1.4\times\pm 1.2)\,\rm{kpc} along the bulge principal axes but is incomplete within \pm 150\,\rm{pc} above and below the galactic plane because of large extinction and crowding. We found in P15 that the vertical density profile of our Nbody models of B/P bulges was well represented by a vertical \rm{sech}^{2} profile. Thus, to fill the midplane gap in the RCG density map we use the fiducial extrapolation of P15, obtained by fitting a \rm{sech}^{2} profile to each vertical slice through the bulge. This extrapolation provides us with the full 3D density of what we call the smooth bulge. In Section 7.3, we consider an extra inplane disk component and show that indeed an inplane over density is required to match our bulge kinematic data. We finally integrate the RCGs map on a grid of (30\times 28\times 32) cells, and the corresponding observable kernels K_{j}^{\rm{RCG}} are given by
K_{j}^{\rm{RCG}}(\bm{z}_{i})=\begin{cases}\rm{M/n_{RCG}}^{1}&\text{if $i\in$ % cell $j$,}\\ 0&\text{otherwise.}\end{cases}  (17) 
Outside the bar region, for cylindrical radius larger than 5\,\rm{kpc}, we use the 3D density of the disk of Section 3.3 evaluated on a large 3D density grid using a massincell kernel similarly to Equation 17.
5.4 Bulge kinematics from the BRAVA survey
The BRAVA survey is a large spectroscopic survey of about 10000 M giant stars, mostly towards the bulge (Rich et al., 2007; Howard et al., 2008; Kunder et al., 2012). We use only the fields with l\leq 10\lx@math@degree as we found in P15 that the disk contamination could become significant outside the bulge. The selected BRAVA fields provide 82 measurements of the mean radial velocity and velocity dispersion through the bulge.
As NMAGIC observables we use here the first and second weighted velocity moments whose kernels are given by:
K_{j}^{\rm{BRAVA},1}(\bm{z}_{i})=\delta_{j}^{\rm{BRAVA}}(\bm{z}_{i})\times% \frac{\textrm{v}_{i}}{W_{j}^{\rm{BRAVA}}}  (18) 
and
K_{j}^{\rm{BRAVA},2}(\bm{z}_{i})=\delta_{j}^{\rm{BRAVA}}(\bm{z}_{i})\times% \frac{\textrm{v}_{i}^{2}}{W_{j}^{\rm{BRAVA}}}  (19) 
where \textrm{v}_{i} is the radial velocity of particle i and \delta_{j}^{\rm{BRAVA}} is the selection function of the BRAVA survey and W_{j}^{\rm{BRAVA}} is given by
W_{j}^{\rm{BRAVA}}=\sum_{i}w_{i}\delta_{j}^{\rm{BRAVA}}(\bm{z}_{i}).  (20) 
As shown in P15 the BRAVA survey is biased towards nearby stars and the selection function is given by:
\delta_{j}^{\rm{BRAVA}}(\bm{z}_{i})=\begin{cases}r_{i}^{1.4}&\text{if $i\in$ % field $j$ with $y_{i}<3.5\,\rm{kpc}$,}\\ 0&\text{otherwise.}\end{cases}  (21) 
5.5 Bulge proper motions from the OGLEII survey
Proper motion data in the bulge for more than half a million stars have been measured by the OGLE survey. In this paper, we chose to use the OGLE proper motions to compare with model predictions rather than using them in the model fitting. We use the proper motion dispersions of RCG in the bulge from the OGLEII survey as computed by Rattenbury et al. (2007) from the proper motion catalogue of Sumi et al. (2004). RCGs are selected in an ellipse of the dereddened I_{0} versus (VI)_{0} CMD centred on the expected locus of the red clump at I_{0}=14.6 and (VI)_{0}=1.0. Stars with proper motion error larger than 1\,\rm{mas\,yr^{1}} in either the l or the b direction were excluded from the sample, as well as stars with total proper motion larger than 10\,\rm{mas\,yr^{1}} which are likely to belong to the foreground disk. We compute the selection function by using again the PARSEC isochrones for a 10\,\rm{Gyr} population, a Kroupa IMF and the metallicity distribution of ARGOS stars in the bulge. The cut in proper motion error introduces a bias towards nearby stars as the error is more likely to be large for faint stars. To model this effect we first compute in each field j the fraction of stars that pass the error cut threshold as a function of their extinction corrected I_{0} magnitude, denoted by C_{\sigma_{\mu},j}(I_{0}). Then from the isochrones we compute for each distance \mu the distribution of I_{0} magnitudes C_{j}(\mu,I_{0}) of stars that end up inside the ellipse selection region of the I_{0} versus (VI)_{0} diagram. The final selection function f^{\rm OGLE}(\mu) is then given by
f^{\rm{OGLE}}(\mu)=\int C_{j}(\mu,I_{0})\times C_{\sigma_{\mu},j}(I_{0})\,dI_{0}  (22) 
and is shown in Figure 9 for the field at (l,b)=(1.0\lx@math@degree,3.7\lx@math@degree). The large theoretical contribution of low distance modulus is due to faint mainsequence stars in the disk that are close enough to fall into the ellipse selection of the CMD. As the stellar population of the nearby disk is very different from the old bulge population, the selection function should not be trusted at small distances. We adopt a simple distance cut and discard contribution from any particle at distances less than 5\,\rm{kpc}. Disk contamination is also likely to be removed in the data thanks to the cut in total proper motion.
Observable kernels are defined similarly to Equations 18 and Equations 19 by replacing the radial velocity by the proper motion in the heliocentric frame, with the selection function \delta_{j}^{\rm{OGLE}} given by:
\delta_{j}^{\rm{OGLE}}(\bm{z}_{i})=\begin{cases}f^{\rm OGLE}(\mu_{i})&\text{if% $i\in$ field $j$ with $r_{i}>5\,\rm{kpc}$}\\ &\text{ and $\sqrt{{\mu^{*}_{l}}_{i}^{2}+{\mu^{*}_{b}}_{i}^{2}}<10\,\rm{mas\,% yr^{1}}$,}\\ 0&\text{otherwise.}\end{cases}  (23) 
The errors in proper motion dispersion quoted by Rattenbury et al. (2007) are statistical errors, very small due to the large number of observed stars. However, Rattenbury et al. (2007) noted that adjacent fields could show variations in proper motion dispersion of up to 0.2\,\rm{mas\,yr^{1}}. To take those systematics into account, we replace the error bars by adding in quadrature to statistical error a systematic error of 0.1\,\rm{mas\,yr^{1}}.
6 Dynamical modelling of the bar region
In this section, we use the M2M method to fit the Nbody models constructed in Section 4 to the data described in Section 5.
6.1 M2M formalism
Indexing again an observable by j and a data set by k, the difference between the model observable y_{j}^{k}(t) and the real data Y_{j}^{k} is evaluated through the residual \Delta_{j}^{k}(t) defined as
\Delta_{j}^{k}(t)=\frac{y_{j}^{k}(t)Y_{j}^{k}}{\sigma(Y_{j}^{k})}  (24) 
where \sigma(Y_{j}^{k}) is the error on Y_{j}^{k}(t).
Following De Lorenzi et al. (2007), we match our data adopting the profit function
F=\frac{1}{2}\sum_{k,j}\lambda_{k}(\Delta_{j}^{k})^{2}+\mu S.  (25) 
The first term in Equation 25 is a weighted chisquare term where the \lambda_{k} are numerical weights for the different data sets (see Long & Mao 2010). The second term is an entropy term, forcing the particle weight distribution to remain narrow around some predetermined prior values, improving hence the convergence of the individual particle weights. We use the pseudoentropy S of Morganti & Gerhard (2012) given by
S=\sum_{i=1}^{N}w_{i}\,\left[\log\left(\frac{w_{i}}{\hat{w_{i}}}\right)1\right]  (26) 
where \hat{w_{i}} are prior weights, chosen to be the mean stellar weight and mean dark matter weight for respectively the stellar and dark matter particles.
Using the observable kernels K_{j}^{k} described in Section 5, Equation 3 becomes
\begin{split}\displaystyle\frac{dw_{i}}{dt}=\varepsilon w_{i}\bigg{[}&% \displaystyle\mu\log\left(\frac{w_{i}}{\hat{w_{i}}}\right)\\ &\displaystyle+\sum_{k}\lambda_{k}\sum_{j}\left(K_{j}(\bm{z}_{i})+w_{i}\frac{% \partial K_{j}(\bm{z}_{i})}{\partial w_{i}}\right)\frac{\Delta_{j}^{k}(t)}{% \sigma(Y_{j}^{k})}\bigg{]}.\end{split}  (27) 
When the observable kernels do not depend on the weights, Syer & Tremaine (1996) and De Lorenzi et al. (2007) showed that the observables converge exponentially on a timescale O(\varepsilon^{1}) provided the initial model is sufficiently close to the solution. We find in practice that weightdependent kernels, as chosen for the ARGOS and BRAVA surveys, can be introduced without affecting the convergence provided (i) the derivative of the kernel with respect to the weight is taken into account in Equation 27, (ii) some other observables constrain the absolute value of the weights on the spatial domain where the weightdependent kernel is nonnull. We numerically evaluate the convergence of the weights by following the method of Long & Mao (2010) considering that a given particle weight has converged if its maximum relative deviation from its mean value on the previous orbit is smaller than 10\%.
6.2 Fitting procedure and parametrization
After the adiabatic evolution of the initial model, a typical NMAGIC fit consists of three phases:

Temporal smoothing: evolve and observe the model for a time \rm{T_{smooth}} to initialize the observables.

M2M fit: evolve the model and modify the particle weights according to Equation 27 for a time \rm{T_{M2M}}.

Relaxation: evolve and observe the model for a time \rm{T_{relax}} to test the stability of the fit.
We use 4, 40 and 16 internal time units for respectively \rm{T_{smooth}}, \rm{T_{M2M}} and \rm{T_{relax}}, with one internal unit of time corresponding to the period of one circular orbit at 4\,\rm{kpc} (i.e. about 125\,\rm{Myr}). We adopt a smoothing timescale of 1/\alpha=1 model time units. The force of change parameter \varepsilon is fixed to \varepsilon=10^{1}\times w_{0} where w_{0} is the mean particle mass. We determined the \lambda_{k} of Equation 25 by analyzing the contribution to the bracket term of Equation 27 of each set of observables separately. We found that the values of \lambda_{\rm density}=1 for all density observables and \lambda_{\rm kinematics}=10 for all kinematic observables give approximately the same median strength of the force of change to each set of observables. We thus adopted these values for the \lambda_{k} and checked that only minimal differences are obtained for values of \lambda_{\rm kinematics}=5 or 20. The priors of the entropy term of Equation 26 are fixed to the mean particle weight and the magnitude of the entropy is fixed to \mu=10^{4}. This value allows more than 97% of the particle to converge while resulting in a weight distribution that extends only to \pm 1\rm{dex} from the priors. With these settings and the observables described above, each NMAGIC run typically requires \sim 190 CPU corehour to be completed.
7 Understanding the bulge dynamics
In this section, we show how the bar pattern speed, dark matter density and inplane stellar bulge density influence the bulge kinematics and how we can recover and constrain the effective potential by fitting our models to the data described in Section 5.
7.1 Signature of the pattern speed and inner dark matter halo in the bulge kinematics
The first column of Figure 10 shows the BRAVA kinematics for three models with different pattern speeds between 25 and 45\,\rm{km\,s^{1}\,kpc^{1}}, fitted to our bulge and long bar data in the same dark matter halo. The effect of the bar pattern speed is clearly visible both in the mean velocity and velocity dispersion. Increasing pattern speed leads to an increase in mean velocity but a decrease in the velocity dispersion. The virial theorem provides intuitive explanation of this: for a given mass distribution and hence a given potential energy, a larger pattern speed places more kinetic energy in pattern rotation, leaving less energy available to build up random motions. The second panel of Figure 10 shows the BRAVA kinematics for three models with the same pattern speed and stellar mass density but different dark matter masses in the bulge. As already found in P15, large masses increase the dispersion but leave the mean velocity essentially unchanged. Hence, for a given pattern speed and stellar density model, we can recover the dark matter mass in the bulge from the BRAVA velocity dispersions.
7.2 Recovering the best dark matter halo
In Section 3.4, we constructed a firstguess dark matter density in the Galaxy by fitting an Einasto profile to the rotation curve between 6\,\rm{kpc} and R_{0} under the constraint of an assumed dark matter mass inside 2\,\rm{kpc}, M_{\rm DM}(<2\,\rm{kpc}). Now from our dynamical modelling, we can adapt the value of M_{\rm DM}(<2\,\rm{kpc}) to what would be required in order to best match the BRAVA velocity dispersions. This is done iteratively during the NMAGIC fit: every eight time units we evaluate the value of a factor denoted by \mathcal{F} to be applied on the model velocity dispersions in the BRAVA fields in order to best fit the data. Heuristically, \mathcal{F}^{2} estimates the multiplicative factor to be applied to the total dynamical mass in the bulge in order to best fit the BRAVA dispersions. Thus, \mathcal{F}^{2}>1 (<1) indicates that more (less) mass in the bulge is required to get a better agreement with the data. Therefore, every eight time units we increase M(<2\,\rm{kpc}) by \Delta M(<2\,\rm{kpc})=1.0\times 10^{10}\,\,\rm{M}_{\odot}\times(\mathcal{F}^{% 2}1) and redefine the dark matter density in the entire galaxy by fitting again the Einasto density to the rotation curve data and the new value for M_{\rm DM}(<2\,\rm{kpc}). The prefactor of 10^{10}\,\,\rm{M}_{\odot} sets the rate at which the inner dark matter is adapted, such that at end of the fit M(<2\,\rm{kpc}) has converged to the BRAVA dispersion for the considered pattern speed and stellar mass distribution. We assume that the halo is an oblate spheroid with a vertical flattening of 0.8 throughout the Galaxy. There is not enough power in our data to see a significant effect of the dark matter shape that is not degenerate with the dark matter profile and therefore more complicated 3D shapes of the dark matter density are not justified. In practice, we expand the target halo density in spherical harmonics up to order eight and fit the dark matter particles to the expansion with the M2M method using the spherical harmonics kernels extensively described in previous uses of NMAGIC (see De Lorenzi et al. 2007 and subsequent work for more detail).
7.3 The missing central mass
Using the modelling procedure defined above and the best dark matter halo that matches the overall BRAVA dispersion, we find evidence for a lack of central mass concentration in the models. In Figure 11 the purple line and points show the model predictions of the BRAVA dispersions along the minor axis and the OGLE proper motion dispersions along the b direction for our fiducial bulge density and a pattern speed of 40\,\rm{km\,s^{1}\,kpc^{1}}. The underestimation of the central dispersion is very clear, and very little freedom is available either in the masstoclump ratio or in the pattern speed to correct for this. As NMAGIC adjusts the orbital structure in order to best match the kinematics, the only way remaining to increase the central dispersion is to deepen the gravitational potential by adding an extra central mass component. Motivated by the massive nuclear disk found by Launhardt et al. (2002) and considering the fact that the vertical structure of the bulge density could be steeper than our fiducial \rm{sech}^{2} extrapolation, we model the missing mass by assuming that it is distributed as an elongated exponential disk following the bar, whose density is given by
\rho_{\rm c}\propto\exp\left(\frac{\sqrt{x^{2}+(y/0.5)^{2}}}{h_{\rm r}}\right% )\times\exp\left(\frac{z}{h_{\rm z}}\right)  (28) 
where x, y, and z are coordinates along the principal axes of the bar in \,\rm{kpc}. When h_{\rm r} and h_{\rm z} are too small, the large central concentration leads to the death of the peanut shape as already noted by Athanassoula et al. (2005). When h_{\rm r} and h_{\rm z} are too large, the central concentration is not sufficient to provide the central potential we need to match our kinematic data. After experimenting with several combinations of h_{\rm r} and h_{\rm z}, we adopt the values of 250\,\rm{pc} and 50\,\rm{pc}, respectively.
The other three lines in Figure 11 show the predictions of models containing an extra barlike component with the density of Equation 28 normalized to masses of M_{c}=0.1, 0.2 and 0.3\times 10^{10}\,\,\rm{M}_{\odot}. We see that we need an additional mass of about 0.2\times 10^{10}\,\,\rm{M}_{\odot} in order to reproduce both the BRAVA minor axis dispersions and the OGLE proper motions. In the lower plot of Figure 11 we show the sideon projection of the RCG density in the bulge assuming that this extra mass is stellar and traced by red clump stars. With our parametrization, most of the extra central component is located in the inner midplane \pm 150\,\rm{pc} strip where the RCG density has not been directly measured. The interpretation of this extra mass is discussed in more detail in Section 10.2.
8 Key parameters of the Milky Way’s effective potential
In this section we explore systematically the key parameters of the Milky Way’s effective potential. We focus on the three parameters that have the largest impact on our data sets: (i) the bar pattern speed \Omega_{b}, (ii) the masstoclump ratio \rm{M/n_{RCG}} and (iii) the extra central mass noted M_{\rm c}. We explore nine values of \Omega_{b} between 25\,\rm{km\,s^{1}\,kpc^{1}} and 50\,\rm{km\,s^{1}\,kpc^{1}} in steps of 2.5\,\rm{km\,s^{1}\,kpc^{1}}, three values for \rm{M/n_{RCG}}, 1000 (our fiducial) and 900 and 1100 corresponding to the range of values consistent with the bulge stellar population (see Section 3.2), and five values for M_{c} between 0.1\times 10^{10}\,\,\rm{M}_{\odot} and 0.3\times 10^{10}\,\,\rm{M}_{\odot} in steps of 0.05\times 10^{10}\,\,\rm{M}_{\odot}. This results in a 3D grid of 135 NMAGIC runs for which we will evaluate how well they reproduce each individual data set.
In Figure 12 we give an overview of our 135 simulations. In this figure, the area of the different wedges for each model shows how well that simulation performs in reproducing the various data set, larger area meaning better agreement. The area of each wedge is proportional to 1\tilde{\chi}^{2} where \tilde{\chi}^{2} is the \chi^{2} rescaled so that the best simulation has \tilde{\chi}^{2} of 0 and the median simulation has \tilde{\chi}^{2} of 1, separately for each observable. Some interesting features are directly visible in Figure 12. From the ARGOS and BRAVA kinematics, we see that the masstoclump ratio has some degeneracy with both the pattern speed and the central mass: low masstoclump ratio tends to prefer lower pattern speed and higher central mass. The model at (\rm{M/n_{RCG}},M_{c},\Omega_{b})=(1000,0.2\times 10^{10}\,\,\rm{M}_{\odot},40% \,\rm{km\,s^{1}\,kpc^{1}}) as indicated by the black square in Figure 12 is able to reproduce well all data sets simultaneously and therefore constitutes our best model.
In order to recover the range of \rm{M/n_{RCG}}, M_{c} and \Omega_{b} around the best model that is allowed by the data, we would ideally construct a global likelihood for all our data sets together and search for the 3D region of the grid that contains the 68\% most likely models. However, this approach is not applicable here for two reasons:

Comparing different data sets in a purely statistical way is dubious when the fitted data sets, such as the RCG density in the bulge and the OGLE proper motions, have dominant errors that are not statistical. We are more interested in reproducing all data sets fairly well at once than in maximizing the formal total likelihood, which may be dominated by one particular aspect because of greatly different number of observables or unaccounted systematic errors.

The evaluation of the range of models consistent with some data within some confidence limit is not a straightforward task in M2M modelling. Morganti et al. (2013) showed that the magnitude of the \Delta\chi^{2} that allows a fair estimation of the range of models that is compatible with the data can be much larger than that expected, due to modelling noise.
Instead, we adopt a more phenomenological approach where we use our knowledge of the bar dynamics gained from the experiments of Section 7 to explore the 3D grid of models in more detail. In Section 8.1, we study the masstoclump/pattern speed degeneracy and recover the pattern speed of the galactic bar. In Section 8.2, we study the masstoclump/central mass degeneracy and recover the best value for the central mass. We finally show how our best model compares to the different data sets in Section 8.3.
8.1 Bar pattern speed and corotation
We showed in Section 7.1 that the pattern speed had a clear signature in the mean radial velocity observables. Hence, we focus on the mean velocity of both the BRAVA and ARGOS surveys. For each combination of pattern speed and masstoclump ratio, we first search for the best value of M_{c}, as described in the next subsection. In Figure 13, we show the \chi^{2} per datapoint as a function of the pattern speed \Omega_{b} for different masstoclump ratios and the corresponding best value of M_{c}. Good fits to the data are found in the range \Omega_{b}\sim 3542.5\,\rm{km\,s^{1}\,kpc^{1}} depending on the masstoclump ratio. An increase in masstoclump ratio of 10\% requires an increase in \Omega_{b} of \sim 2.5\,\rm{km\,s^{1}\,kpc^{1}}. We also notice that the ARGOS data systematically prefers pattern speeds \sim 2.5\,\rm{km\,s^{1}\,kpc^{1}} larger than the BRAVA data. Forming a joined \chi^{2} between ARGOS and BRAVA and assuming flat priors on the masstoclump ratio in the range 9001100, we find a mean of the best pattern speeds to be \Omega_{b}=39\,\rm{km\,s^{1}\,kpc^{1}}. As already discussed, the evaluation of statistical errors on measurements and parameters from M2M modelling is usually problematic since the classical \Delta\chi^{2}=1 method tends to underestimate the real error. Morganti et al. (2013) developed a method to better estimate the statistical error using a value of \Delta\chi^{2} corresponding to the scatter of the \chi^{2} surface from the models around the minimum. Using this method we find statistical errors lower than 1\,\rm{km\,s^{1}\,kpc^{1}}, smaller than the systematics arising from both the degeneracy with the masstoclump ratio and the systematic offset between the ARGOS and BRAVA data sets. Adding in quadrature an error of 2.5\,\rm{km\,s^{1}\,kpc^{1}} from both these sources of systematics, we conclude that the pattern speed of the Milky Way bar is \Omega_{b}=39\pm 3.5\,\rm{km\,s^{1}\,kpc^{1}}. Using the composite rotation curve of Sofue et al. (2009) rescaled to (R_{0},V_{0})=(8.2\,\rm{kpc},238\,\rm{km\,s^{1}}), the corotation radius of the bar is found at R_{\rm{cr}}=6.1\pm 0.5\,\rm{kpc}.
8.2 Central mass distribution
We showed in Section 7.3 the necessity of an additional central mass component M_{c} for matching the inner BRAVA dispersions and OGLE proper motions in the b direction. In Figure 14, we show the \chi^{2} per datapoint of the best pattern speed models as a function of M_{c} for all the BRAVA and OGLE fields within 3.5\lx@math@degree from the centre (\sim 500\,\rm{pc} at the distance of the GC). We see that the BRAVA radial velocity prefers an additional mass of 0.2\times 10^{10}\,\,\rm{M}_{\odot} with a slight degeneracy with the masstoclump ratio where a 10\% increase in \rm{M/n_{RCG}} leads to a decrease of M_{c} of about 0.05\times 10^{10}\,\,\rm{M}_{\odot}. The signature of M_{c} in the \sigma_{b} proper motions appears to be in slight tension with the BRAVA dispersion, having a systematic preference for higher values of M_{c}. Note however that these proper motions are simply predictions and thus are not expected to perfectly fit the data. We show in Figure 16 that for a central mass of 0.2\times 10^{10}\,\,\rm{M}_{\odot} our best model provides a very good fit to the proper motion dispersions in both directions, staying on average within 5% of the data, even though it systematically underestimates it. In addition, by looking at the \chi^{2} of the RCG density as shown in the bottom plot of Figure 14, we find that large values of M_{c} tend to weaken or destroy the peanut shape of the bulge, as was already shown by previous studies (e.g Shen & Sellwood 2004; Athanassoula et al. 2005).
Since no value of M_{c} simultaneously gives a best fit to the OGLE \sigma_{b} predictions and the RCG density fit, we estimate M_{c} based on the BRAVA central dispersions to be M_{c}=0.200.05\times(\rm{M/n_{RCG}}1000)/100\,\times 10^{10}\,\,\rm{M}_{\odot} and recognize the presence of possible unaccounted systematic effects.
Since a 10\% increase in masstoclump ratio requires a 2.5\,\rm{km\,s^{1}\,kpc^{1}} increase in pattern speed and a 0.05\times 10^{10}\,\,\rm{M}_{\odot} decrease in central mass we define two boundary models around our best model at each end of this three dimensional degeneracy, (\rm{M/n_{RCG}},M_{c},\Omega_{b})=(900,0.25\times 10^{10}\,\,\rm{M}_{\odot},37% .5\,\rm{km\,s^{1}\,kpc^{1}}) and (\rm{M/n_{RCG}},M_{c},\Omega_{b})=(1100,0.15\times 10^{10}\,\,\rm{M}_{\odot},4% 2.5\,\rm{km\,s^{1}\,kpc^{1}}). Those two models are indicated by the grey squares in Figure 12 and are used in Section 9 to quantify uncertainties in measuring the stellar and dark matter mass distributions.
8.3 Best fitting dynamical model of the Galaxy
Starting with the bulge kinematics, Figure 15 shows how our best model compares to the BRAVA mean velocities and velocity dispersions. The agreement is overall very good and is an improvement over P15 where the dispersions at 6\lx@math@degree were always underestimated.
The largest improvement over P15 is seen in the proper motion dispersions in the l and b directions shown in Figure 16. In nearly all fields, the model is within 10% of the OGLE proper motions even though the proper motion data were not included in the fitting procedure. We notice however that the model tends to systematically underestimate the data by about 5\%. An increase in \sigma_{b} can be obtained by increasing the central mass but at the cost of losing the peanut shape; or by lowering the pattern speed at the cost of a worse agreement with the mean radial velocity. In the l direction, an increase in \sigma_{l} can be obtained by increasing the bar angle, at the cost of a worse fit to the RCGs histograms in the bar. Since a 10\% variation, i.e. 0.20.3\,\rm{mas\,yr^{1}} is only slightly larger than the fieldtofield variation of the data (\sim 0.2\,\rm{mas\,yr^{1}}, Rattenbury et al. 2007), we consider our model as already consistent with the OGLE proper motions.
Extending to the bar region, Figure 17 shows the ARGOS mean velocity as a function of distance. The streaming motion along the bar is very visible at latitude b=5\lx@math@degree from the twist in the mean velocity as a function of distance. Here again the model performs very well in reproducing the data.
The transition between the bulge and the bar, together with the inplane structure of the bar, is shown in Figure 18.
The model does a very good job at reproducing the RCG distribution for b\geq 2\lx@math@degree. In the plane, and mostly in the region 12\lx@math@degree\leq l\leq 22\lx@math@degree, the model cannot find a good fit to the very narrow distribution of RCGs along the line of sight. This is probably due to the superthin bar component found by W15. The investigation of the detailed structure of the superthin bar is beyond the scope of this paper but could be addressed in the future using the APOGEE data. The APOGEE survey, as an infrared spectroscopic survey, can penetrate the high extinction in the plane and provides stellar kinematics and chemical abundances in the superthin bar region. By looking at clumps in chemical space, Hogg et al. (2016) already found evidence for a very young stellar component that may correspond to the superthin bar of W15.
Anticipating the future use of APOGEE to constrain the dynamical models further, we show in Figure 19 a comparison between our best model and the latest APOGEE kinematics from Ness et al. (2016). Our model is already in very good qualitative agreement with APOGEE even though a more quantitative comparison would require the modelling of the APOGEE selection function which is not included here.
9 Stellar and dark matter mass distribution in the Milky Way
In this section we study the stellar and dark matter mass distribution of our best model and evaluate the effect of variations in the different modelling assumptions. As a range of reasonable variations around our best model, we will consider the following:

The two boundary models of Section 8 found at each end of the three dimensional degeneracy valley between masstoclump ratio, bar pattern speed and additional central mass;

Varying the bar angle from our fiducial 28\lx@math@degree to either 23\lx@math@degree or 33\lx@math@degree;

Varying the dark matter flattening from our fiducial 0.8 to either 0.6 or 0.4;
These variations will be used to evaluate errors (systematic) on the mass parameters found from our best model. All errors quoted in the next two sections refer to half the range of values found in the fiducial and variation models. A summary of the mass parameters of our best model is given in Table 2.
9.1 Bulge, bar, inner and outer disk
In Figure 20 and Figure 21, we show the surface density map and profiles of our best model with the range of the model variations, assuming that the additional mass component is stellar as discussed in Section 10.2. The entire Galaxy is to some level nonaxisymmetric. This is very clear at the edge of the bulge where the surface density along the minor axis at 2\,\rm{kpc} is a factor of about 4 smaller than along the major axis. Both profiles cross at 6.3\,\rm{kpc} from the centre, close to the corotation radius. Beyond corotation, the surface density becomes larger along the minor axis than along the major axis, as one would expect based on the linear perturbation of near circular orbit (Dehnen, 2000; Binney & Tremaine, 2008). Along the major axis, the density of the bar is close to exponential. This is a characteristic of latetype barred galaxies as revealed by the \rm{S^{4}G} survey (Elmegreen et al., 2011). Traditionally, the complex structure of the stellar density of the Galaxy has been divided into a small number of discrete components, mainly bulge, bar and disk for which we would wish to measure mass and shape. We will focus here on three definitions:

Bar, bulge and disk structure from stellar mass ‘photometric’ profiles: From the major and minor axis profiles in Figure 21, we see clearly three regimes: (a) the outer disk, nearly axisymmetric and exponential outside \sim 5.3\,\rm{kpc} (b) the inner disk, axisymmetric with a nearly constant surface density inside 5.3\,\rm{kpc} (c) the bar/bulge, i.e. the bar which formed a bulge in its inner part. This photometric definition of the bar and inner disk components has the advantage to be easily applicable to external galaxies. By integrating the surface density associated with the ‘photometric’ bulge and bar, we find a ‘photometric’ bar/bulge stellar mass of M_{\rm{bar/bulge}}=1.88\pm 0.12\times 10^{10}\,\,\rm{M}_{\odot}, among which 1.34\pm 0.04\times 10^{10}\,\,\rm{M}_{\odot} is located in the bulge and 0.54\pm 0.04\times 10^{10}\,\,\rm{M}_{\odot} in the long bar. The ‘photometric’ stellar mass associated with the inner disk within the bar region is found to be M_{\rm{inner\ disk}}=1.29\pm 0.12\times 10^{10}\,\,\rm{M}_{\odot}, lower than the mass of the bar/bulge structure. The errors in these quantities are systematics derived from the variation models presented previously.

Nonaxisymmetric long bar from 2D ‘photometry’: An alternative way to define the bar component is to search for non axisymmetries in the faceon surface density. We first define a maximum axisymmetric model using the minor axis profile and remove this maximum axisymmetric model from the original surface density map. By integrating the residuals for radii inside 5.3\,\rm{kpc} and outside the bulge, we find a nonaxisymmetric long bar mass of M_{\rm{nonaxi}}=0.46\pm 0.03\times 10^{10}\,\,\rm{M}_{\odot}. This measure is similar to but slightly lower than the previous ‘photometric’ long bar mass since the inner disk has a similar but slightly lower surface density than our maximum axisymmetric model. From both estimates combined, the bar outside the bulge has slightly less than half the mass of the disk in the same radial range.

Barfollowing orbits: Unlike some spiral structures, the bar is not a density wave: stars that form the bar stay in the bar and orbit mostly along elongated orbits of the x_{1} family and its descendants (Contopoulos & Grosbøl, 1989). Since we have access to individual orbits in our dynamical model, we can directly identify orbits that compose the bar using the method of Portail et al. (2015b). We integrate all particles and compute the dominant frequencies of the time variation of the cylindrical radius f_{r} and bar major axis position f_{x} in the corotating frame. Barsupporting orbits are found in the vicinity of f_{r}/f_{x}=2, i.e have two radial oscillations for one period along the bar. In the bar region, particles that do not follow the bar have Rosettalike orbits in the bar frame for which f_{r}/f_{x}\neq 2; they build the inner disk. This orbitbased definition of the bar is more elegant than the ‘photometric’ definition and is also closer to what makes the bar a separate component but is in general not observable in external galaxies for which individual orbits are usually unknown. We find a stellar mass on barsupporting orbit of 1.04\pm 0.06\times 10^{10}\,\,\rm{M}_{\odot}. This estimate misses the nonbar following orbits in the bulge.
W15 determined from a combination of VVV, UKIDSS and 2MASS data the length of the bar and found a halflength of 5.0\pm 0.2\,\rm{kpc}. Since our model is the first nonparametric fit of the galactic bulge and bar it is important to see how it compares to direct determination from the data. We follow W15 and focus on the following three methods to measure the bar halflength:

L_{\rm{drop}}: radius at which the ellipticity of the bar drops the fastest;

L_{\rm{prof}}: radius at which the major and minor axes agree within 30%;

L_{\rm{m=2}}: radius at which the relative m=2 component of the Fourier decomposition of the surface density drops to 20% of its maximum value.
We measure L_{\rm{drop}}=5.77\,\rm{kpc}, L_{\rm{prof}}=5.12\,\rm{kpc} and L_{\rm{m=2}}=5.02\,\rm{kpc}. By taking the mean of those three measurements and the standard deviation of the three measurements applied on all our variation models, we find a bar halflength of the galactic bar of 5.30\pm 0.36\,\rm{kpc}, in good agreement with the measurement of 5.0\pm 0.2\,\rm{kpc} found by W15 from the their component fit of the RCG magnitude distributions.
9.2 The dark matter mass distribution in the Milky Way
In Section 7.2, we showed how we recognize and evaluate the need of dark matter in the bulge from the BRAVA kinematics once the stellar density and bar pattern speed are fixed. Since we have a \sim 10\% accurate measurement of the masstoclump ratio in the bulge, we have access to the dark matter mass distribution in the bulge. In the volume of the box in which the RCG density was measured, i.e. a box of (\pm 2.2\times\pm 1.4\times\pm 1.2)\,\rm{kpc} along the bar principal axes, the mass budget of the bulge is found to be as follows: stars as traced by RCGs account for 1.32\pm 0.08\times 10^{10}\,\,\rm{M}_{\odot}, dark matter accounts for 0.32\pm 0.05\times 10^{10}\,\,\rm{M}_{\odot} and 0.2\times 10^{10}\,\,\rm{M}_{\odot} of additional mass are required in the centre, probably stars in the nuclear disk (see Section 10.2). The resulting total dynamical bulge mass is 1.85\pm 0.05\times 10^{10}\,\rm{M}_{\odot}, in excellent agreement with our estimation of 1.84\pm 0.07\times 10^{10}\,\,\rm{M}_{\odot} found in P15 from modelling the bulge only. Altogether our best model has a low dark matter fraction in the bulge of only 17\%\pm 2\%. The small error quoted here represents half the range of dark matter fractions found in the set of the fiducial and variation models considered and could underestimate the true error on the determination of the dark matter fraction in the galactic bulge. It arises because variations in the total mass, stellar mass and central mass of the bulge in our model variations approximately compensate leaving the dark matter mass in the bulge nearly unchanged.
Going away from the bulge, the dark matter halo profile of our best model is shown in Figure 22 together with the range of profiles from the variation models. We find a local dark matter density of \rho_{\rm{DM}}(R_{0})=0.0132\pm 0.0014\,\,\rm{M}_{\odot}.\,\rm{pc}^{3}, in very good agreement with the recent measurement of \rho_{\rm{DM}}(R_{0})=0.0154\pm 0.0023\,\,\rm{M}_{\odot}.\,\rm{pc}^{3} from Piffl et al. (2014) for a halo flattening of 0.8 (see Read (2014) for a review). Interestingly, our dark matter profile is cored. For the fitted Einasto density profiles, the presence of a core is a consequence of accounting simultaneously for a low dark matter mass in the bulge, a significant dark matter mass enclosed within the solar radius and a gently rising halo rotation curve.
Would an NFW halo profile be allowed by our data? In the Milky Way, an NFW halo density would be expected to have a scale radius of 2040\,\rm{kpc} (BlandHawthorn & Gerhard, 2016). Thus, in the bar region, the density would be well approximated by a simple power law \rho_{\rm{DM}}(r)\propto r^{1}. Given the baryonic mass distribution of our best model and scaling the dark matter density profile such that the total circular velocity at the solar radius is matched, we find a dark matter mass in the bulge that coincidentally is in very good agreement with our best model value. However, such an NFW halo fails to reproduce the nearly flat total rotation curve observed between 6 and 8\,\rm{kpc}, with a halo circular velocity that is about 17\,\rm{km\,s^{1}} lower than the data at 6\,\rm{kpc} shown in Figure 23. The presence of a core in our best model halo density thus appears as a consequence of the constraint on the flat shape of the total circular velocity in the 68\,\rm{kpc} range, which for our baryonic mass distribution requires the dark matter density to fall off more steeply than \rho_{\rm{DM}}(r)\propto r^{1}. In order to then not overpredict the dark matter mass in the bulge, the dark halo density is forced to become shallower further in.
For an Einsato profile, this results in a central core. However, we note that our constraint on the dark matter density in the bulge arises from a constraint on the dynamical mass of the bulge. Thus, we would also expect good agreement with the data for a steeper halo density profile in the bulge, provided it has the same dark matter mass in the bulge. In order to evaluate the possibility of a steeper density in the inner halo, we perform the following experiment. Assuming for the dark matter inside 2\,\rm{kpc} a power law density \rho\propto r^{\alpha}, we determine for each model the powerlaw index \alpha that would keep constant the dark matter mass enclosed within 2\,\rm{kpc}, while ensuring continuity of the density profile at 2\,\rm{kpc}. Doing so we find from our models power law slopes of \sim0.6, whose range is shown as the dark grey span in Figure 22. We thus conclude that a central powerlaw cusp shallower than \alpha\sim0.6 would also provide a good fit to the data.
A different view at the dark matter contribution to the gravitational potential can be seen in Figure 23 where we plot the resulting rotation curve of our best model and the range of rotation curves provided by the model variations, again considering the additional central mass as entirely baryonic. The dark matter support to the rotation is often expressed as the degree of maximality of the disk (Sackett, 1997), representing the ratio of the stellar rotation velocity to the total rotation velocity at some particular radius. In the case of an axisymmetric disk galaxy, this radius is traditionally chosen to be 2.2 disk scalelength, corresponding to the position of the peak of the disk rotation curve. Since the Milky Way is not axisymmetric and hosts a central bulge, the stellar rotation curve does not peak at 2.2 outer disk scalelengths but instead at only about 1\,\rm{kpc} where a stellar rotation of 185\,\rm{km\,s^{1}} provides 94% of the rotational support. At 2.2 outer disk scalelengths, the stellar contribution to the rotation curve drops to 75\% still within the range of what would be called a maximum disk (Sackett, 1997). This result is in agreement with the microlensing analysis of Wegg et al. (2016).
Mass inside a box of (\pm 2.2\times\pm 1.4\times\pm 1.2)\,\rm{kpc} along the bar axes  

Smooth bulge traced by RCGs  1.32\pm 0.08\times 10^{10}\,\,\rm{M}_{\odot} 
Nuclear disk  0.20\times 10^{10}\,\,\rm{M}_{\odot} 
Dark matter  0.32\pm 0.05\times 10^{10}\,\,\rm{M}_{\odot} 
Dynamical mass  1.85\pm 0.05\times 10^{10}\,\rm{M}_{\odot} 
Mass of the bulge, bar and inner disk  
‘Photometric” bulge+bar  1.88\pm 0.12\times 10^{10}\,\,\rm{M}_{\odot} 
Inner disk  1.29\pm 0.12\times 10^{10}\,\,\rm{M}_{\odot} 
‘Photometric’ bulge  1.34\pm 0.04\times 10^{10}\,\,\rm{M}_{\odot} 
‘Photometric’ long bar  0.54\pm 0.04\times 10^{10}\,\,\rm{M}_{\odot} 
Nonaxisymmetric long bar  0.46\pm 0.03\times 10^{10}\,\,\rm{M}_{\odot} 
Bar supporting orbits  1.04\pm 0.06\times 10^{10}\,\,\rm{M}_{\odot} 
Dark matter distribution  
Bulge dark matter fraction  17\pm 2\% 
Local dark matter density  0.0132\pm 0.0014\,\,\rm{M}_{\odot}.\,\rm{pc}^{3} 
Disk maximality at 2.2 disk scalelength  75\% 
10 Discussion
10.1 An intermediate bar pattern speed
The bar pattern speed is a fundamental quantity that sets the position of resonances in the Milky Way. It has been the focus of many studies in the last two decades but the scatter between different determination methods remains surprisingly large. The new value determined in this paper from the dynamics of the bulge and long bar is 39\pm 3.5\,\rm{km\,s^{1}\,kpc^{1}}, based on several large photometric and kinematic survey data sets.
A number of previous determinations concluded on rather large pattern speeds in the range of 5065\,\rm{km\,s^{1}\,kpc^{1}} from gas dynamics (Englmaier & Gerhard, 1999; Fux, 1999; Bissantz et al., 2003), continuity of a tracer stellar population (Debattista et al., 2002) or interpretation of stellar stream in the solar neighborhood as a resonance effect (Dehnen, 2000; Antoja et al., 2014). Other studies found lower values in the range 2540\,\rm{km\,s^{1}\,kpc^{1}} from gas dynamics (Weiner & Sellwood, 1999; RodriguezFernandez & Combes, 2008; Sormani et al., 2015; Li et al., 2016) or stellar dynamics (Long et al. 2013; P15).
Our new value for the bar pattern speed agrees very well with the extensive recent analysis of gasdynamical models by Sormani et al. (2015), but it is not consistent with the most precise determination so far based on recent analysis of the Hercules stream; this gives 53\pm 0.5\,\rm{km\,s^{1}} when rescaled to (R_{0},V_{0})=(8.2\,\rm{kpc},238\,\rm{km\,s^{1}}) (Antoja et al., 2014). However, this measurement is model dependent, assuming that the Hercules stream originates from the bar’s outer Lindblad resonance. Future work must show whether an alternative interpretation of the Hercules stream can be found.
Lower or intermediate pattern speeds are supported by the recent measurement of the long bar halflength of 5.0\pm 0.2\,\rm{kpc} by W15. Since the bar cannot extend beyond corotation (Contopoulos, 1980), the corotation radius has to be greater than \sim 5\,\rm{kpc}, putting an upper limit on the pattern speed of about \sim 48\,\rm{km\,s^{1}\,kpc^{1}}. Bar pattern speeds are often expressed in terms of the dimensionless ratio \mathcal{R} between the corotation radius R_{\rm{cr}} and the bar halflength R_{\rm{bar}}. Typical values for external disk galaxies are \mathcal{R}=1.2\pm 0.2 (Elmegreen et al., 1996) with an indication of a correlation with bar strength (Aguerri et al., 1998) and morphological type (Rautiainen et al., 2008). Using the measured R_{\rm{bar}}=5.0\pm 0.2\,\rm{kpc} from W15 and our determination of R_{\rm{cr}}=6.1\pm 0.5\,\rm{kpc}, we find for the Milky Way a ratio of \mathcal{R}=1.22\pm 0.11. This places the Milky Way together with the bulk of external barred spiral galaxies, with a bar that can be classified as a fast rotator (\mathcal{R}\leq 1.4, Debattista & Sellwood (2000)).
10.2 The extra central mass
We found in Section 7.3 kinematic evidence for an additional central concentration of mass that was not included in our fiducial RCG bulge model. We can find good agreement with the kinematics by including an additional central mass component, mostly located in the plane where the RCG density has not been directly measured. Therefore, assuming that the required mass is stellar mass does not violate the 3D density measured from RCGs. In fact, there are several pieces of evidence for a stellar overdensity in the plane and near the GC but no consensus yet on the precise shape and mass of such overdensity. Launhardt et al. (2002) found a very concentrated nuclear bulge component in the central 220\,\rm{pc} from decomposition of the IRAS and COBE DIRBE data and estimate its mass to be 0.14\times 10^{10}\,\,\rm{M}_{\odot}. From star counts, this component has a nearexponential vertical profile with scaleheight 45\,\rm{pc} corresponding to a nuclear stellar disk with axial ratio \sim 35:1 (Nishiyama et al., 2013). Schönrich et al. (2015) found in the APOGEE data kinematic evidence of a nuclear disk extending to \sim 150\,\rm{pc} with a vertical height of 50\,\rm{pc}. They measure a rotation velocity of 120\,\rm{km\,s^{1}} at 150\,\rm{pc} which, assuming an exponential density profile, leads to a nuclear bulge mass in reasonable agreement with the mass estimate from Launhardt et al. (2002). In addition, Debattista et al. (2015) postulated the presence of a \,\rm{kpc} scale nuclear disk to explain the highvelocity peaks in the line of sight velocity distributions of the APOGEE commissioning data (Nidever et al., 2012). Such a largescale disk would not be concentrated enough to sufficiently increase the central velocity dispersion required by the modelling, but could account for part of the mass in the galactic plane in the bulge region.
Hence, the most likely interpretation for our extra central mass component is simply a stellar over density near the centre, not included in the largescale RCG bulge. This could well be a similar component as in the SBb galaxy NGC 4565 where HST and Spitzer photometry revealed an additional disky pseudobulge hidden inside the B/P bulge (Kormendy & Barentine, 2010). More data closer to the centre and extending into the plane such as APOGEE (Majewski, 2012) or GIBS (Zoccali et al., 2014) would be required to investigate the structure of this extra mass further.
Some contribution to the required central mass could also come from metal poor stars not traced by RCGs such as old RR Lyrae stars (Dékány et al., 2013; Pietrukowicz et al., 2015). In the ARGOS fields, metalpoor stars with \rm{[Fe/H]}\leq0.9 are only 7\% of the total sample but since they are more concentrated than the RCG bulge stars they could still play a role in the centre. However, there is still no clear picture about what component these stars belong to and what mass is associated with it. Favorite hypotheses are either an old thick disk (Di Matteo et al., 2015), a small classical bulge (Kunder et al., 2016) or the inner part of the stellar halo (PérezVillegas et al., 2017).
10.3 Comparison with the bulge models of P15
In P15 we made a series of five dynamical models of the galactic bulge with different dark matter fractions called M80 to M90 by combining the 3D density of RCGs in the bulge from Wegg & Gerhard (2013) and the BRAVA kinematics. Our main findings from these models were a measurement of the bulge total mass (stellar + dark matter) of 1.84\pm 0.07\times 10^{10}\,\,\rm{M}_{\odot} and a rather low bulge pattern speed in the range 2530\,\rm{km\,s^{1}\,kpc^{1}}. In comparison, our new best model presented here has a total bulge mass of 1.85\times 10^{10}\,\,\rm{M}_{\odot}, in very good agreement with our previous determination. However, from our more advanced modelling, we find a higher pattern speed of 40\,\rm{km\,s^{1}\,kpc^{1}}. We attribute the difference to the lack of mass surrounding the bulge in the models of P15. As already stated in Section 2.2, bars formed in Nbody models are generally several disk scalelengths long, which consequently underestimate the impact of inner disk and long bar orbits on the bulge dynamics. Hence, caution must be taken when interpreting the dynamics of Nbody B/P bulges, unless the bulge, bar and disk scalelengths are all consistent with each other.
The present best model has a smooth bulge stellar mass of 1.32\times 10^{10}\,\,\rm{M}_{\odot}, an additional central mass of 0.2\times 10^{10}\,\,\rm{M}_{\odot} and a dark matter mass in the bulge of 0.32\times 10^{10}\,\,\rm{M}_{\odot}, summing up to the 1.85\times 10^{10}\,\,\rm{M}_{\odot} stated above. Its closest equivalent in the P15 series of models is the fitted model M82.5 which has a masstoclump ratio of 1014, close to the value of 1000 we measured directly^{2}^{2}2This corrects Fig. 16 of P15 where all model values of the masstoclump ratios are underestimated by \sim 20\% due to an overestimate of the red giant branch bump fraction.. The main differences are the larger bar pattern speed in the new model and the fact that the new model has 0.2\times 10^{10}\,\,\rm{M}_{\odot} less dark matter in the bulge but instead a similar additional mass close to the centre, as required by the central velocity dispersion for a larger pattern speed.
10.4 Dark matter in the Galaxy
Early Nbody dark matter simulations predicted a universal cuspy dark matter density profile (Navarro et al., 1997). This was later confirmed by new simulations with higher resolution able to resolve halo densities on scales smaller than a percent of r_{200} (Navarro et al., 2010), i.e about 2\,\rm{kpc} for the Milky Way halo (BlandHawthorn & Gerhard, 2016). In addition, since the inner part of dark matter halos is actually populated by baryons that collapse through dissipative processes, dark halos should contract, thus exacerbating any preexisting cusp (Blumenthal et al., 1986; Gnedin et al., 2004; Abadi et al., 2010). However, the degree of contraction of the halo varies greatly between authors (see Wegg et al. 2016 for a recent comparison).
In disk galaxies, processes related to the dominant role of baryons in the central parts have been found to be able to transform a primordial cusp into a core. Such processes are the resonance effects of a large primordial stellar bar (Weinberg & Katz 2002; but see Dubinski et al. 2009), supernova feedback (Pontzen & Governato, 2012) and stellar feedback (Chan et al. 2015; Schaller et al. 2015; but see Marinacci et al. 2014). How important those processes are is not yet settled, resulting in a cusp/core controversy similar to that in dwarf galaxies (Moore, 1994; Burkert, 1995; De Blok, 2010).
From our modelling, we find a best model with a low dark matter fraction in the galactic bulge of 17\pm 2\% where the 2\% comes from the variation models of Section 9. In order to match simultaneously a low dark matter fraction in the bulge with the flat rotation curve close to the solar radius R_{0}, the dark matter density of our model has a powerlaw slope that is steeper than \propto r^{1} immediately inside R_{0} and then flattens to a shallow cusp or a core in the bulge region. This result is consistent with the recent work of Wegg et al. (2016) who find that a high baryonic fraction is required to account for the high optical depths towards the bulge measured in the EROSII and MOAII microlensing data. Given the similarity between the best model here and the fiducial model of Wegg et al. (2016), we expect our best model to be also consistent with microlensing constraints towards the bulge; this will be discussed in a later paper.
11 Conclusion
We build a large number of dynamical models of the bar region in the Milky Way using the MadetoMeasure method. We first create a set of Nbody models of barred disks that broadly matches the bulge, bar and outer disk density by adiabatic adaptation of a initial Nbody model. This adiabatic procedure allows us to adapt the dark matter distribution to the rotation curve of the Galaxy inside 10 kpc and also to modify the pattern speed of the galactic bar. We then constrain those models with the stellar density of the bulge and bar as traced by red clump giants from a combination of the VVV, UKIDSS and 2MASS surveys, together with stellar kinematics from the BRAVA, OGLE and ARGOS surveys. We explore a threedimensional parameter space given by the stellar mass fraction in the bulge, the bulge and bar pattern speed and the nuclear disk mass, and provide constraints on the galactic effective potential. Our main conclusions are the following:

Modelling the stellar dynamics in the bar region requires a bar pattern speed of \Omega_{b}=39\pm 3.5\,\rm{km\,s^{1}\,kpc^{1}}, placing corotation at 6.1\pm 0.5\,\rm{kpc} from the Galactic Centre. The ratio of corotation radius to bar halflength (R_{\rm{bar}}=5.0\pm 0.2\,\rm{kpc} from Wegg et al. 2015) of the Galaxy is found to be \mathcal{R}=1.22\pm 0.11, in good agreement with what is seen in external disk galaxies.

We find a total dynamical mass of the bulge of 1.85\pm 0.05\times 10^{10}\,\rm{M}_{\odot}, in excellent agreement with the value found in Portail et al. (2015a) from modelling the bulge only.

We find dynamical evidence for an extra central mass component, not included in our previous bulge models, of about 0.2\times 10^{10}\,\,\rm{M}_{\odot} and probably related to a nuclear disk or disky pseudobulge.

We evaluate from our model the mass of the long bar and bulge structure and find M_{\rm{bar/bulge}}=1.88\pm 0.12\times 10^{10}\,\,\rm{M}_{\odot}, larger than the mass of disk in the bar region, M_{\rm{inner\ disk}}=1.29\pm 0.12\times 10^{10}\,\,\rm{M}_{\odot}. The mass of the long bar is slightly less than half the disk mass in the same radial range. Our models predict a nonexponential surface density for the disk in the bar region and illustrate the transition between the bar region and the outer disk.

We also evaluate the need of dark matter in the inner Milky Way. Using recent measurements of the bulge IMF and more extended data, we now better constrain the stellartodark matter fraction in the bulge and find a preference for a masstoclump ratio of 1000 and a low dark matter fraction of 17\pm 2\% in the bulge. In order to match simultaneously a low dark matter fraction in the bulge with the flat rotation curve close to the solar radius R_{0}, the dark matter density of our model has a powerlaw slope that is steeper than \propto r^{1} immediately inside R_{0} and then flattens to a shallow cusp or a core in the bulge region.
Our bestfitting model is the first nonparametric model of the entire bar region of the Milky Way and can be of significant use for several on going and future Milky Way studies including gas dynamics in realistic galactic bar potentials and chemodynamics of the different stellar populations. The model may be made available upon request to the authors.
Acknowledgements
We thank Manuela Zoccali and Elena Valenti for helpful discussions, and are grateful to Jerry Sellwood and Monica Valluri for making their potential solver code available to us. We acknowledge the great work of the VVV, ARGOS, BRAVA, UKIDSS, OGLE and 2MASS survey teams upon which this paper is built. We also thank the anonymous referee for their careful reading and constructive comments.
References
 Abadi et al. (2010) Abadi M. G., Navarro J. F., Fardal M., Babul A., Steinmetz M., 2010, MNRAS, 407, 435
 Aguerri et al. (1998) Aguerri J. A. L., Beckman J. E., Prieto M., 1998, AJ, 116, 2136
 Antoja et al. (2014) Antoja T., et al., 2014, A&A, 563, A60
 Athanassoula (2002) Athanassoula E., 2002, ApJ, 569, L83
 Athanassoula (2007) Athanassoula E., 2007, MNRAS, 377, 1569
 Athanassoula et al. (2005) Athanassoula E., Lambert J. C., Dehnen W., 2005, MNRAS, 363, 496
 Benjamin et al. (2005) Benjamin R. A., et al., 2005, ApJ, 630, L149
 Bensby et al. (2011) Bensby T., et al., 2011, A&A, 533, A134
 Bensby et al. (2013) Bensby T., et al., 2013, A&A, 549, A147
 Binney (2010) Binney J., 2010, MNRAS, 401, 2318
 Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics, 2nd edn. Princeton
 Binney et al. (1990) Binney J. J., Davies R. L., Illingworth G. D., 1990, ApJ, 361, 78
 Binney et al. (1991) Binney J., Gerhard O. E., Stark A. A., Bally J., Uchida K. I., 1991, MNRAS, 252, 210
 Bissantz et al. (2003) Bissantz N., Englmaier P., Gerhard O., 2003, MNRAS, 340, 949
 Bissantz et al. (2004) Bissantz N., Debattista V. P., Gerhard O., 2004, ApJ, 601, L155
 BlandHawthorn & Gerhard (2016) BlandHawthorn J., Gerhard O., 2016, ARA&A, p. 69
 Blitz & Spergel (1991) Blitz L., Spergel D. N., 1991, ApJ, 379, 631
 Blumenthal et al. (1986) Blumenthal G. R., Faber S. M., Flores R., Primack J. R., 1986, ApJ, 301, 27
 Bovy & Rix (2013) Bovy J., Rix H.W., 2013, ApJ, 779, 115
 Bovy et al. (2012) Bovy J., et al., 2012, ApJ, 759, 131
 Bressan et al. (2012) Bressan A., Marigo P., Girardi L., Salasnich B., Dal Cero C., Rubele S., Nanni A., 2012, MNRAS, 427, 127
 Burkert (1995) Burkert A., 1995, ApJ, 447
 CabreraLavers et al. (2008) CabreraLavers A., GonzálezFernández C., Garzón F., Hammersley P. L., LópezCorredoira M., 2008, A&A, 491, 781
 Calamida et al. (2015) Calamida A., et al., 2015, ApJ, 810, 8
 Cappellari et al. (2009) Cappellari M., et al., 2009, ApJ, 704, L34
 Chan et al. (2015) Chan T. K., Kereš D., Oñorbe J., Hopkins P. F., Muratov A. L., FaucherGiguère C.A., Quataert E., 2015, MNRAS, 454, 2981
 Chemin et al. (2015) Chemin L., Renaud F., Soubiran C., 2015, A&A, 578, A14
 Chen et al. (2014) Chen Y., Girardi L., Bressan A., Marigo P., Barbieri M., Kong X., 2014, MNRAS, 444, 2525
 Clarkson et al. (2008) Clarkson W., et al., 2008, ApJ, 684, 1110
 Combes et al. (1990) Combes F., Debbasch F., Friedli D., Pfenniger D., 1990, A&A, 233, 82
 Contopoulos (1980) Contopoulos G., 1980, A&A, 81, 198
 Contopoulos & Grosbøl (1989) Contopoulos G., Grosbøl P., 1989, A&ARv, 1, 261
 Das et al. (2011) Das P., Gerhard O., Mendez R. H., Teodorescu A. M., de Lorenzi F., 2011, MNRAS, 415, 1244
 De Blok (2010) De Blok W. J. G., 2010, Adv. Astron., 2010, 789293
 De Lorenzi et al. (2007) De Lorenzi F., Debattista V. P., Gerhard O., Sambhus N., 2007, MNRAS, 376, 71
 De Lorenzi et al. (2008) De Lorenzi F., Gerhard O., Saglia R. P., Sambhus N., Debattista V. P., Pannella M., Méndez R. H., 2008, MNRAS, 385, 1729
 De Lorenzi et al. (2009) De Lorenzi F., et al., 2009, MNRAS, 395, 76
 Debattista & Sellwood (2000) Debattista V. P., Sellwood J. A., 2000, ApJ, 543, 704
 Debattista et al. (2002) Debattista V. P., Gerhard O., Sevenster M. N., 2002, MNRAS, 334, 355
 Debattista et al. (2015) Debattista V. P., Ness M., Earp S. W. F., Cole D. R., 2015, ApJ, 812, L16
 Dehnen (2000) Dehnen W., 2000, AJ, 119, 800
 Dehnen (2009) Dehnen W., 2009, MNRAS, 395, 1079
 Dejonghe (1984) Dejonghe H., 1984, A&A, 133, 225
 Dékány et al. (2013) Dékány I., Minniti D., Catelan M., Zoccali M., Saito R. K., Hempel M., Gonzalez O. A., 2013, ApJ, 776, L19
 Di Matteo et al. (2015) Di Matteo P., et al., 2015, A&A, 577, A1
 Dubinski et al. (2009) Dubinski J., Berentzen I., Shlosman I., 2009, ApJ, 697, 293
 Einasto (1965) Einasto J., 1965, Tr. Astrofiz. Instituta AlmaAta, 5, 87
 Elmegreen et al. (1996) Elmegreen B. G., Elmegreen D. M., Chromey F. R., Hasselbacher D. A., Bissell B. A., 1996, AJ, 111, 2233
 Elmegreen et al. (2011) Elmegreen D. M., et al., 2011, ApJ, 737, 32
 Englmaier & Gerhard (1999) Englmaier P., Gerhard O., 1999, MNRAS, 304, 512
 Freeman et al. (2013) Freeman K., et al., 2013, MNRAS, 428, 3660
 Fux (1999) Fux R., 1999, A&A, 345, 787
 Gnedin et al. (2004) Gnedin O. Y., Kravtsov A. V., Klypin A. A., Nagai D., 2004, ApJ, 616, 16
 Gonzalez et al. (2011) Gonzalez O. A., Rejkuba M., Zoccali M., Valenti E., Minniti D., 2011, A&A, 534, A3
 Hammersley et al. (1994) Hammersley P. L., Garzon F., Mahoney T., Calbet X., 1994, MNRAS, 269
 Hogg et al. (2016) Hogg D. W., Casey A. R., Ness M., Rix H.W., ForemanMackey D., 2016, preprint (arXiv:1601.05413)
 Howard et al. (2008) Howard C. D., Rich R. M., Reitzel D. B., Koch A., De Propris R., Zhao H., 2008, ApJ, 688, 1060
 Hunt & Kawata (2014) Hunt J. A. S., Kawata D., 2014, MNRAS, 443, 2112
 Jurić et al. (2008) Jurić M., et al., 2008, ApJ, 673, 864
 Kormendy & Barentine (2010) Kormendy J., Barentine J. C., 2010, ApJ, 715, L176
 Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
 Kunder et al. (2012) Kunder A., et al., 2012, AJ, 143, 57
 Kunder et al. (2016) Kunder A., et al., 2016, ApJ, 821, L25
 Launhardt et al. (2002) Launhardt R., Zylka R., Mezger P. G., 2002, A&A, 384, 112
 Li et al. (2016) Li Z., Gerhard O., Shen J., Portail M., Wegg C., 2016, ApJ, 824, 11
 Long & Mao (2010) Long R. J., Mao S., 2010, MNRAS, 405, 301
 Long et al. (2013) Long R. J., Mao S., Shen J., Wang Y., 2013, MNRAS, 428, 3478
 LópezCorredoira et al. (2005) LópezCorredoira M., CabreraLavers A., Gerhard O. E., 2005, A&A, 439, 107
 Majewski (2012) Majewski S. R., 2012, in American Astronomical Society Meeting Abstracts, 219, 20506
 Majewski et al. (2011) Majewski S., Zasowski G., Nidever D., 2011, ApJ, 739, 25
 Marinacci et al. (2014) Marinacci F., Pakmor R., Springel V., 2014, MNRAS, 437, 1750
 MartinezValpuesta & Gerhard (2011) MartinezValpuesta I., Gerhard O., 2011, ApJ, 734, L20
 MartinezValpuesta et al. (2006) MartinezValpuesta I., Shlosman I., Heller C., 2006, ApJ, 637, 214
 McWilliam & Zoccali (2010) McWilliam A., Zoccali M., 2010, ApJ, 724, 1491
 Moore (1994) Moore B., 1994, Nature, 370, 629
 Morganti & Gerhard (2012) Morganti L., Gerhard O., 2012, MNRAS, 422, 1571
 Morganti et al. (2013) Morganti L., Gerhard O., Coccato L., MartinezValpuesta I., Arnaboldi M., 2013, MNRAS, 431, 3570
 Nakada et al. (1991) Nakada Y., Onaka T., Yamamura I., Deguchi S., Hashimoto O., Izumiura H., Sekiguchi K., 1991, Nature, 353, 140
 Nataf et al. (2010) Nataf D. M., Udalski A., Gould A., Fouqué P., Stanek K. Z., 2010, ApJ, 721, L28
 Navarro et al. (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
 Navarro et al. (2010) Navarro J. F., et al., 2010, MNRAS, 402, 21
 Ness et al. (2013) Ness M., et al., 2013, MNRAS, 430, 836
 Ness et al. (2016) Ness M., et al., 2016, ApJ, 819, 2
 Nidever et al. (2012) Nidever D. L., et al., 2012, ApJ, 755, L25
 Nishiyama et al. (2006) Nishiyama S., et al., 2006, ApJ, 638, 839
 Nishiyama et al. (2013) Nishiyama S., et al., 2013, ApJ, 769, L28
 PérezVillegas et al. (2017) PérezVillegas A., Portail M., Gerhard O., 2017, MNRAS, 464, L80
 Peters, III (1975) Peters, III W., 1975, ApJ, 195, 617
 Pietrukowicz et al. (2015) Pietrukowicz P., et al., 2015, ApJ, 811, 113
 Piffl et al. (2014) Piffl T., et al., 2014, MNRAS, 445, 3133
 Pontzen & Governato (2012) Pontzen A., Governato F., 2012, MNRAS, 421, 3464
 Portail et al. (2015a) Portail M., Wegg C., Gerhard O., MartinezValpuesta I., 2015a, MNRAS, 448, 713 (P15)
 Portail et al. (2015b) Portail M., Wegg C., Gerhard O., 2015b, MNRAS, 450, L66
 Qian et al. (1995) Qian E., De Zeeuw P., van der Marel R., Hunter C., 1995, MNRAS, 274, 602
 Rattenbury et al. (2007) Rattenbury N. J., Mao S., Debattista V. P., Sumi T., Gerhard O., De Lorenzi F., 2007, MNRAS, 378, 1165
 Rautiainen et al. (2008) Rautiainen P., Salo H., Laurikainen E., 2008, MNRAS, 388, 1803
 Read (2014) Read J. I., 2014, J. Phys. G Nucl. Part. Phys., 41, 063101
 Reid et al. (2014) Reid M. J., et al., 2014, ApJ, 783, 130
 Rich et al. (2007) Rich R. M., Reitzel D. B., Howard C. D., Zhao H., 2007, ApJ, 658, L29
 RodriguezFernandez & Combes (2008) RodriguezFernandez N. J., Combes F., 2008, A&A, 489, 115
 Sackett (1997) Sackett P. D., 1997, ApJ, 483, 103
 Saito et al. (2012) Saito R. K., et al., 2012, A&A, 537, A107
 Salaris & Girardi (2002) Salaris M., Girardi L., 2002, MNRAS, 337, 332
 Sanders & Binney (2013) Sanders J. L., Binney J., 2013, MNRAS, 433, 1813
 Schaller et al. (2015) Schaller M., et al., 2015, MNRAS, 451, 1247
 Schönrich (2012) Schönrich R., 2012, MNRAS, 427, 274
 Schönrich et al. (2010) Schönrich R., Binney J., Dehnen W., 2010, MNRAS, 403, 1829
 Schönrich et al. (2015) Schönrich R., Aumer M., Sale S. E., 2015, ApJ, 812, L21
 Schwarzschild (1979) Schwarzschild M., 1979, ApJ, 232, 236
 Sellwood (2003) Sellwood J. A., 2003, ApJ, 587, 638
 Sellwood & Debattista (2009) Sellwood J. A., Debattista V. P., 2009, MNRAS, 398, 1279
 Sellwood & Valluri (1997) Sellwood J. A., Valluri M., 1997, MNRAS, 287, 124
 Sharma et al. (2014) Sharma S., et al., 2014, ApJ, 793, 51
 Shen & Sellwood (2004) Shen J., Sellwood J. A., 2004, ApJ, 604, 614
 Shen et al. (2010) Shen J., Rich R. M., Kormendy J., Howard C. D., De Propris R., Kunder A., 2010, ApJ, 720, L72
 Skrutskie et al. (2006) Skrutskie M. F., et al., 2006, AJ, 131, 1163
 Sofue et al. (2009) Sofue Y., Honma M., Omodaka T., 2009, PASJ, 61, 227
 Sormani et al. (2015) Sormani M. C., Binney J., Magorrian J., 2015, MNRAS, 454, 1818
 Stanek et al. (1997) Stanek K. Z., Udalski A., Szymanski M., Kaluzny J., Kubiak Z. M., Mateo M., Krzeminski W., 1997, ApJ, 477, 163
 Sumi et al. (2004) Sumi T., et al., 2004, MNRAS, 348, 1439
 Syer & Tremaine (1996) Syer D., Tremaine S., 1996, MNRAS, 282, 223
 Tang et al. (2014) Tang J., Bressan A., Rosenfield P., Slemer A., Marigo P., Girardi L., Bianchi L., 2014, MNRAS, 445, 4287
 Thomas et al. (2009) Thomas J., et al., 2009, MNRAS, 393, 641
 Valenti et al. (2016) Valenti E., et al., 2016, A&A, 587, L6
 Wegg & Gerhard (2013) Wegg C., Gerhard O., 2013, MNRAS, 435, 1874
 Wegg et al. (2015) Wegg C., Gerhard O., Portail M., 2015, MNRAS, 450, 4050 (W15 )
 Wegg et al. (2016) Wegg C., Gerhard O., Portail M., 2016, MNRAS, 463, 557
 Weiland et al. (1994) Weiland J. L., et al., 1994, ApJ, 425, L81
 Weinberg & Katz (2002) Weinberg M. D., Katz N., 2002, ApJ, 580, 627
 Weiner & Sellwood (1999) Weiner B. J., Sellwood J. A., 1999, ApJ, 524, 112
 Zhu et al. (2014) Zhu L., et al., 2014, ApJ, 792, 59
 Zoccali et al. (2000) Zoccali M., Cassisi S., Frogel J. A., Gould A., Ortolani S., Renzini A., Rich R. M., Stephens A. W., 2000, ApJ, 530, 418
 Zoccali et al. (2003) Zoccali M., et al., 2003, A&A, 399, 931
 Zoccali et al. (2014) Zoccali M., et al., 2014, A&A, 562, A66