Dynamical modelling of the Galactic bar

# Dynamical modelling of the Galactic bulge and bar: the Milky Way’s bar pattern speed, stellar, and dark matter mass distribution

Matthieu Portail{}^{1}, Ortwin Gerhard{}^{1}, Christopher Wegg{}^{1} and Melissa Ness{}^{2}
{}^{1} Max-Planck-Institut für Extraterrestrische Physik, Gießenbachstraße, D-85741 Garching, Germany
{}^{2} Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany
E-mail:portail@mpe.mpg.deE-mail:gerhard@mpe.mpg.de
Accepted 2016 October 31. Received 2016 October 28; in original form 2016 August 28
###### Abstract

We construct a large set of dynamical models of the galactic bulge, bar and inner disk using the Made-to-Measure 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: centre
pubyear: 2016pagerange: Dynamical modelling of the Galactic bulge and bar: the Milky Way’s bar pattern speed, stellar, and dark matter mass distributionReferences\PassOptionsToPackage

pdfpagelabels=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 non-axisymmetric gas flow (Peters, III, 1975; Binney et al., 1991) and asymmetries in the near-infrared 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 so-called 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 N-body models by buckling of a vertically unstable stellar bar (Combes et al., 1990; Martinez-Valpuesta 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ópez-Corredoira et al. 2005; Cabrera-Lavers et al. 2008; but see also Martinez-Valpuesta & 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 N-body 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 Made-to-Measure method. In this paper, we extend the Made-to-Measure 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 Made-to-Measure modelling is particularly challenging.

The paper is organized as follows. In Section 2, we briefly describe the Made-to-Measure 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 N-body models with different bar pattern speeds that already broadly match the Milky-Way 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 non-parametric model of the entire bar region and may be made available upon request to the authors.

## 2 Made-to-Measure modelling of the Galaxy

### 2.1 M2M modelling

Stellar dynamical equilibria for galaxies can be studied with moment-based methods (Binney et al., 1990; Cappellari et al., 2009), classical distribution function-based methods (Dejonghe, 1984; Qian et al., 1995), actions-based methods (Binney, 2010; Sanders & Binney, 2013), orbit-based methods (Schwarzschild, 1979; Thomas et al., 2009) or with the Made-to-Measure (M2M) method (Syer & Tremaine, 1996; De Lorenzi et al., 2007). In these M2M models, an initial self-gravitating N-body model is used to provide a discrete sample of a distribution function reasonably close to the system of interest. This N-body model is then slowly adapted by modifying the weights of the N-body particles such as to make the model reproduce a given set of constraints. The N-body weights can hence be seen simultaneously as mass elements (N-body 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 phase-space 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 phase-space 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 time-scale of the weight evolution. The profit function F usually consists of a chi-square 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 N-body orbits. It is thus a very efficient modelling technique provided all the N-body 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 N-body model of a barred stellar disk that provides a first-guess discrete sampling of the final model. The classical way to build such initial conditions is to evolve a near-equilibrium N-body 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; Martinez-Valpuesta et al., 2006). This process has been used widely in order to build models that can then be compared to data (Martinez-Valpuesta et al., 2006; Athanassoula, 2007; Shen et al., 2010) but suffers from three major limitations:

• The evolution process is non-linear 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.

• The bar tends to be 3-5 disk scale-lengths long as already noted by Debattista & Sellwood (2000) and Athanassoula (2002). This is much larger than what would be required to model the Milky Way where this ratio lies within 1.9\pm 0.4 (Bland-Hawthorn & Gerhard, 2016).

An example of such a buckled bar N-body model is the model M85 of P15 shown in Figure 1. When scaled to a bar half-length 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 N-body evolution shown above, we use a variant the M2M method to tailor initial conditions in a two-step process. In Section 3, we first create a mass density model of the Milky Way by adding together our best-guess 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 N-body 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} (Bland-Hawthorn & Gerhard, 2016) from the GC. The bar is oriented at an angle of \alpha=28\lx@math@degree with the Sun-GC 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@degree-33\lx@math@degree measured by W15 from the long bar RCGs. Following Bland-Hawthorn & 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 non-parametric 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 completeness-corrected 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 eight-fold 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 Sun-GC line of sight. They estimated the long bar half-length 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 best-fitting 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 color-magnitude 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 4-5\,\rm{Gyr} old population among stars with [\rm{Fe/H}]\geq-0.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 mass-to-clump ratio 111Note that, as in P15 our definition of the mass-to-clump 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 to-light-ratio. We make the fiducial assumption that the bulge and the thin long bar have the same mass-to-clump 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 mass-to-clump ratio than the bulge. Assuming a constant star formation rate and a Kroupa initial mass function (IMF) as in W15, the mass-to-clump 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 mass-to-clump ratio in the bulge

The mass-to-clump 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 ultra-deep 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 mass-to-clump 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 extinction-corrected 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 mass-to-clump 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 low-mass end of the IMF due to unseen dwarfs has only a small effect on the mass-to-clump ratio. Variations from a sub-stellar 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 mass-to-clump ratio of only \pm 3\%. Our direct measurement of the mass-to-clump 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 mass-to-clump 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 Bland-Hawthorn & Gerhard (2016) concluded that h_{R,*}=2.6\pm 0.5\,\rm{kpc} with the shorter disk scalelengths in the range 2.1-2.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 three-parameter 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 tangent-point 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.

## 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.

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, log-spaced 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 time-scale.

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 equal-weights 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 phase-space sampled by the parent particle.

Following this procedure, we create a set of N-body models that broadly matches the static density model of Section 3 with bar pattern speeds in the range 25-50\,\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 N-body 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 N-body models.

### 4.2 Integration and potential solver

The particle model is integrated using a drift-kick-drift adaptive leap-frog 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 N-body models with different pattern speeds that already broadly matches the Milky-Way bulge, bar, and disk density. In this section, we describe the different data sets to which we fit our N-body 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 N-body 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 medium-resolution 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 point-source 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.

1. 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 I-band (I\sim 13-14, 14-15 and 15-16).

2. The 2MASS stars from which the ARGOS stars are selected in (i) do not correspond to the full point-source 2MASS catalogue but only to a high-quality subsample of it. This subsample is defined by a blue color cut (J-K_{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 high-quality imaging is easier to achieve.

3. 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 star-by-star basis using the Rayleigh-Jeans Color Excess method (Majewski et al. 2011; W15) given by:

 A_{K_{s}}=\frac{A_{K_{s}}}{E(J-K_{s})}[(J-K_{s})-(J-K_{s})_{\rm RCG})] (7)

where (J-K_{s})_{\rm RCG} is the colors of RCGs and \frac{A_{K_{s}}}{E(J-K_{s})} is a constant that depends on the extinction law. For consistency with Wegg & Gerhard (2013), we adopt (J-K_{s})_{\rm RCG}=0.674 (Gonzalez et al., 2011) and \frac{A_{K_{s}}}{E(J-K_{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 extinction-corrected magnitude {K_{s}}_{0} to belong to the high-quality 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 completeness-corrected 2MASS catalogue that are also in the high-quality 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 I-band 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 star-by-star 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 non-RCG 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 top-hat 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 N-body 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 N-body 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 non-RCGs 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 mass-weighted 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 mass-to-clump 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 data-model 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 mass-to-clump 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 mass-to-clump 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 mass-to-clump 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 N-body models of B/P bulges was well represented by a vertical \rm{sech}^{2} profile. Thus, to fill the mid-plane 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 in-plane disk component and show that indeed an in-plane 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 mass-in-cell 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 OGLE-II 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 OGLE-II 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 (V-I)_{0} CMD centred on the expected locus of the red clump at I_{0}=14.6 and (V-I)_{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 (V-I)_{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 main-sequence 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 N-body 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 chi-square 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 pre-determined prior values, improving hence the convergence of the individual particle weights. We use the pseudo-entropy 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 time-scale O(\varepsilon^{-1}) provided the initial model is sufficiently close to the solution. We find in practice that weight-dependent 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 weight-dependent kernel is non-null. 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:

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

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

3. 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 time-scale 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 core-hour to be completed.

## 7 Understanding the bulge dynamics

In this section, we show how the bar pattern speed, dark matter density and in-plane 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 first-guess 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 pre-factor 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 mass-to-clump 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 bar-like 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 side-on 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 mid-plane \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 mass-to-clump 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 mass-to-clump ratio has some degeneracy with both the pattern speed and the central mass: low mass-to-clump 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:

1. 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.

2. 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 mass-to-clump/pattern speed degeneracy and recover the pattern speed of the galactic bar. In Section 8.2, we study the mass-to-clump/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 mass-to-clump 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 mass-to-clump ratios and the corresponding best value of M_{c}. Good fits to the data are found in the range \Omega_{b}\sim 35-42.5\,\rm{km\,s^{-1}\,kpc^{-1}} depending on the mass-to-clump ratio. An increase in mass-to-clump 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 mass-to-clump ratio in the range 900-1100, 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 mass-to-clump 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 mass-to-clump 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.20-0.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 mass-to-clump 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.2-0.3\,\rm{mas\,yr^{-1}} is only slightly larger than the field-to-field 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 in-plane 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 mass-to-clump 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;

• Varying the outer disk scalelength from our fiducial 2.4\,\rm{kpc} to either 2.15\,\rm{kpc} (Bovy & Rix, 2013) or 2.6\,\rm{kpc} (Jurić et al., 2008);

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 non-axisymmetric. 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 late-type 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:

1. 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.

2. Non-axisymmetric long bar from 2D ‘photometry’: An alternative way to define the bar component is to search for non axisymmetries in the face-on 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 non-axisymmetric long bar mass of M_{\rm{non-axi}}=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.

3. Bar-following 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. Bar-supporting 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 Rosetta-like orbits in the bar frame for which f_{r}/f_{x}\neq 2; they build the inner disk. This orbit-based 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 bar-supporting orbit of 1.04\pm 0.06\times 10^{10}\,\,\rm{M}_{\odot}. This estimate misses the non-bar following orbits in the bulge.

W15 determined from a combination of VVV, UKIDSS and 2MASS data the length of the bar and found a half-length of 5.0\pm 0.2\,\rm{kpc}. Since our model is the first non-parametric 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 half-length:

• 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 half-length 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 mass-to-clump 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 20-40\,\rm{kpc} (Bland-Hawthorn & 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 6-8\,\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 power-law 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 \sim-0.6, whose range is shown as the dark grey span in Figure 22. We thus conclude that a central power-law cusp shallower than -\alpha\sim-0.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).

## 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 50-65\,\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 25-40\,\rm{km\,s^{-1}\,kpc^{-1}} from gas dynamics (Weiner & Sellwood, 1999; Rodriguez-Fernandez & 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 gas-dynamical 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 half-length 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 half-length 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 near-exponential vertical profile with scaleheight 45\,\rm{pc} corresponding to a nuclear stellar disk with axial ratio \sim 3-5: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 high-velocity peaks in the line of sight velocity distributions of the APOGEE commissioning data (Nidever et al., 2012). Such a large-scale 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 large-scale 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 pseudo-bulge 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, metal-poor stars with \rm{[Fe/H]}\leq-0.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érez-Villegas 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 25-30\,\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 N-body 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 N-body 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 mass-to-clump ratio of 1014, close to the value of 1000 we measured directly222This corrects Fig. 16 of P15 where all model values of the mass-to-clump 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 N-body 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 (Bland-Hawthorn & 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 pre-existing 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 power-law 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 EROS-II and MOA-II 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 Made-to-Measure method. We first create a set of N-body models of barred disks that broadly matches the bulge, bar and outer disk density by adiabatic adaptation of a initial N-body 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 three-dimensional 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:

1. 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 half-length (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.

2. 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.

3. 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 pseudo-bulge.

4. 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 non-exponential surface density for the disk in the bar region and illustrate the transition between the bar region and the outer disk.

5. 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 stellar-to-dark matter fraction in the bulge and find a preference for a mass-to-clump 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 power-law 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 best-fitting model is the first non-parametric 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
• Bland-Hawthorn & Gerhard (2016) Bland-Hawthorn 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
• Cabrera-Lavers et al. (2008) Cabrera-Lavers A., González-Fernández C., Garzón F., Hammersley P. L., López-Corredoira 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., Faucher-Giguè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 Alma-Ata, 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., Foreman-Mackey 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ópez-Corredoira et al. (2005) López-Corredoira M., Cabrera-Lavers 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
• Martinez-Valpuesta & Gerhard (2011) Martinez-Valpuesta I., Gerhard O., 2011, ApJ, 734, L20
• Martinez-Valpuesta et al. (2006) Martinez-Valpuesta 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., Martinez-Valpuesta 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érez-Villegas et al. (2017) Pérez-Villegas 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., Martinez-Valpuesta 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
• Rodriguez-Fernandez & Combes (2008) Rodriguez-Fernandez 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
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters