Chemodynamical Modelling of the Galactic Bulge and Bar

TL;DR

This paper introduces a chemodynamical model of the Milky Way's bulge and bar, integrating metallicity data to analyze stellar populations and their dynamics, revealing distinct properties of metal-rich and metal-poor stars.

Contribution

It extends the Made-to-Measure method to include stellar metallicity, enabling detailed chemodynamical modeling of the galactic bulge, bar, and inner disk.

Findings

Metal-rich stars are strongly barred with disk-like dynamics.

Metal-poor stars show diverse kinematic behaviors, indicating multiple stellar populations.

The model reproduces observed vertex deviations in the bulge.

Abstract

We present the first self-consistent chemodynamical model fitted to reproduce data for the galactic bulge, bar and inner disk. We extend the Made-to-Measure method to an augmented phase-space including the metallicity of stars, and show its first application to the bar region of the Milky Way. Using data from the ARGOS and APOGEE (DR12) surveys, we adapt the recent dynamical model from Portail et al. to reproduce the observed spatial and kinematic variations as a function of metallicity, thus allowing the detailed study of the 3D density distributions, kinematics and orbital structure of stars in different metallicity bins. We find that metal-rich stars with [Fe/H] > -0.5 are strongly barred and have dynamical properties that are consistent with a common disk origin. Metal-poor stars with [Fe/H] < -0.5 show strong kinematic variations with metallicity, indicating varying contributions from the underlying stellar populations. Outside the central kpc, metal-poor stars are found to have the density and kinematics of a thick disk while in the inner kpc, evidence for an extra concentration of metal-poor stars is found. Finally, the combined orbit distributions of all metallicities in the model naturally reproduce the observed vertex deviations in the bulge. This paper demonstrates the power of Made-to-Measure chemodynamical models, that when extended to other chemical dimensions will be very powerful tools to maximize the information obtained from large spectroscopic surveys such as APOGEE, GALAH and MOONS.

Figures

Tables

Table 1. Table 1: Parameters of the four Gaussian EMDFs used to describe the particle MDFs in our chemodynamical M2M modelling, based on table 3 of Ness et al. ( 2013a ) .

	$\alpha$	$\beta$	$\gamma$	$\delta$
${\text{M}}_c^0$	0.11	-0.28	-0.69	-1.18
${\sigma }_c$	0.12	0.13	0.15	0.14

Table 2. Table 2: Stellar masses within 5 ​ kpc 5 kpc 5\,\rm{kpc} from the Galactic Centre for all models and metallicity bins, given in units of 10 10 ​ M ⊙ superscript 10 10 subscript M direct-product 10^{10}\,\,\rm{M}_{\odot}

Models	A	B	C	D
Fiducial	1.16	1.11	0.37	0.04
Variation with flat priors	1.20	1.07	0.37	0.04
Variation with BGM priors	1.10	1.16	0.39	0.04
Variation with $\log g\le 1\text{dex}$	1.12	1.14	0.39	0.04
Variation with [Fe/H] offset of $+0.1\text{dex}$	1.19	1.11	0.35	0.04
Variation with [Fe/H] offset of $-0.1\text{dex}$	1.03	1.17	0.42	0.06

Equations

$$y=\int K(\mathbf{x},\mathbf{v})f(\mathbf{x},\mathbf{v})d^3\mathbf{x}{d}^3\mathbf{v}$$

$$y(t)=\sum _{i=1}^N{w}_i(t)\times K({\mathbf{x}}_i(t),{\mathbf{v}}_i(t))$$

$$\tilde{y}(t)={\int }_{\tau }y(t-\tau )\times e^{-\alpha \tau }d\tau$$

$$\frac{d{w}_i}{dt}=\epsilon w_i\frac{\partial F}{\partial w_i}$$

$$y={\int }_{\text{M}}{\int }_{\mathbf{x},\mathbf{v}}K(\mathbf{x},\mathbf{v},\text{M})f(\mathbf{x},\mathbf{v},\text{M})d^3\mathbf{x}{d}^3\mathbf{v}d\text{M}.$$

$$f(\mathbf{x},\mathbf{v},\text{M})=\sum _i{w}_i\times \delta (\mathbf{x}-{\mathbf{x}}_i)\delta (\mathbf{v}-{\mathbf{v}}_i)f_{\text{M},i}(\text{M}).$$

$$f_{\text{M},i}=\sum _{c\in \alpha -\delta }{w}_{i,c}{\mathcal{G}}_{[{\text{M}}_c^0,{\sigma }_c]}$$

$$y_j=\sum _i{w}_i\times \sum _{c\in \alpha -\delta }{K}_{j,c}({\mathbf{x}}_i,{\mathbf{v}}_i)w_{i,c}$$

$$K_{j,c}({\mathbf{x}}_i,{\mathbf{v}}_i)={\int }_{\text{M}}{K}_j({\mathbf{x}}_i,{\mathbf{v}}_i,\text{M})\times {\mathcal{G}}_{[{\text{M}}_c^0,{\sigma }_c]}(\text{M})d\text{M}.$$

$$\frac{d{w}_{i,c}}{dt}={\epsilon }_c\frac{\partial F}{\partial w_{i,c}}$$

$$F=-\frac{1}{2}\sum _j{\left(\frac{y_j-Y_j}{\sigma (Y_j)}\right)}^2$$

$$K_{j,d,m}^{\mathrm{ARGOS},0}({\mathbf{z}}_i)={\delta }_{j,d,m}^{\mathrm{ARGOS}}({\mathbf{z}}_i)$$

$$K_{j,d,m}^{\mathrm{ARGOS},1}({\mathbf{z}}_i)=\frac{{\delta }_{j,d,m}^{\mathrm{ARGOS}}({\mathbf{z}}_i)}{W_{j,d,m}^{\mathrm{ARGOS}}}\times v_i$$

$$K_{j,d,m}^{\mathrm{ARGOS},2}({\mathbf{z}}_i)=\frac{{\delta }_{j,d,m}^{\mathrm{ARGOS}}({\mathbf{z}}_i)}{W_{j,d,m}^{\mathrm{ARGOS}}}\times v_i^2$$

$$W_{j,d,m}^{\mathrm{ARGOS}}=\sum _i{w}_i{\delta }_{j,d,m}^{\mathrm{ARGOS}}({\mathbf{z}}_i)$$

$${\delta }_{j,d,m}^{\mathrm{ARGOS}}=\{\begin{array}{ll}{\int }_{\mu (d)}^{\mu (d+1)}{f}_i(\mu ,\text{M})d\mu & \text{if i∈ field j }\\ & \text{and M in bin m,}\\ 0& \text{otherwise}\end{array}$$

$$\begin{array}{rl}& y_{j,d,m}^{\mathrm{ARGOS},0}=\sum _i{w}_i\times \sum _{c\in \alpha -\delta }{w}_{i,c}\times \\ & {\int }_{\text{M}\in m}{\int }_{\mu (d)}^{\mu (d+1)}{f}_i(\mu ,\text{M})\times {\mathcal{G}}_{[{\text{M}}_c^0,{\sigma }_c]}(\text{M})d\mu d\text{M}\end{array}$$

$$\tan (2{\theta }_{r,l})=\frac{2{\sigma }_{r,l}^2}{|{\sigma }_r^2-{\sigma }_l^2|}$$

Full text

††pubyear: 2017††pagerange: Chemodynamical Modelling of the Galactic Bulge and Bar–References

1 Introduction

In the last two decades, our understanding of the inner Galaxy and in particular its central bulge has dramatically changed, driven by the availability of large-scale surveys. The galactic bulge was initially thought to be a classical bulge, spheroidal remnant of mergers that happened during the hierarchical assembly of the Galaxy. This was supported by its old age (Ortolani et al., 1995) and the presence of a vertical metallicity gradient, which was thought at the time to rule out other formation scenarios (Friedli et al., 1994). The picture changed when non-axisymmetries, found consistently in the HI and CO gas flow (Binney et al., 1991; Englmaier & Gerhard, 1999; Fux, 1999), in the near-infrared light distribution (Blitz & Spergel, 1991; Weiland et al., 1994; Bissantz & Gerhard, 2002) and star counts (Nakada et al., 1991; Stanek et al., 1997; López-Corredoira et al., 2005) revealed the presence of a triaxial and rather short boxy bar structure in the bulge (Dwek et al., 1995; Binney et al., 1997), extending to $2-3\,\rm{kpc}$ from the centre.

Recently, Nataf et al. (2010) and McWilliam & Zoccali (2010) independently discovered the so-called split red clump, caused by the bimodal density of red clump giants (RCGs) in the bulge. This bimodality is best illustrated by the work of Wegg & Gerhard (2013) who were able to reconstruct the 3D density of 8 million RCGs stars identified in the VVV catalogue, by direct deconvolution of their magnitude distributions in about 300 fields throughout the bulge. Their density map reveals very clearly the so-called boxy/peanut (B/P) shape of the galactic bulge. Such B/P structured bulges had already been seen in N-body simulations of barred stellar disks and indicated a different formation scenario, where the galactic bulge formed from disk stars through the vertical instability of a stellar bar (Combes & Sanders, 1981; Raha et al., 1991; Debattista & Sellwood, 2000; Athanassoula, 2005; Martinez-Valpuesta et al., 2006).

Additional pieces of evidence now support this formation scenario for the galactic bulge. First, spectroscopic surveys such as the BRAVA (Kunder et al., 2012) and ARGOS (Freeman et al., 2013) surveys showed that the bulge rotates cylindrically, characteristic of a disk origin. Then, Wegg et al. (2015) showed by combining the VVV, UKIDSS and 2MASS catalogues that the galactic bulge smoothly segues into the flat long bar, that extends to about $5\,\rm{kpc}$ from the Galactic Centre. They find that both the bulge and the long bar appear at a consistent angle with respect to the line of sight towards the Galactic Centre, with no need for a secondary galactic bar component, as proposed earlier (Benjamin et al. 2005; López-Corredoira et al. 2005; Cabrera-Lavers et al. 2008; but see also Martinez-Valpuesta & Gerhard 2011). Henceforth, the galactic bulge is now believed to be the vertically extended part of a buckled long bar component, qualitatively similar to the bar produced in N-body simulations.

The first comprehensive and self-consistent dynamical model of the bulge and bar structure in the Milky Way was constructed recently by Portail et al. (2017, hereafter P17). Using the Made-to-Measure (M2M) method, P17 were able to find a dynamical equilibrium model that reproduces simultaneously the RCG density in the bulge and bar region from a combination of the VVV, UKIDSS and 2MASS surveys together with stellar kinematics from the BRAVA, OGLE and ARGOS surveys.

Based on our new knowledge of the current dynamical state of the inner Milky Way, the next stage is to try to understand its chemodynamical formation history. Although individual orbits of stars are scrambled during the bar formation and subsequent buckling instability, Martinez-Valpuesta & Gerhard (2013) showed that the Jacobi energy of stars is largely conserved during the process. Thus, by advanced modelling of the current state of the Galaxy we can hope to trace back earlier stages (Di Matteo et al., 2014). This is challenging and requires as much data as possible, in particular chemical information on stars. Several studies have tried to infer a formation scenario for the bulge based on measurements of the stellar metallicity and its gradient along the bulge minor axis. Martinez-Valpuesta & Gerhard (2013) showed that the presence of a vertical metallicity gradient in the bulge was in itself not a strong argument for the presence of a classical bulge component as previously thought. They found that an initial radial gradient in a disk can, after bar formation and buckling, result in a vertical metallicity gradient similar to that observed in the galactic bulge (Minniti et al., 1995; Zoccali et al., 2008; Gonzalez et al., 2011b), thus suggesting that the galactic bulge could have a pure disk origin. However, other studies (e.g. Babusiaux et al., 2010; Gonzalez et al., 2011a; Hill et al., 2011; Uttenthaler et al., 2012; Rojas-Arriagada et al., 2014; Zoccali et al., 2016) have found evidence for two stellar populations in the bulge with different kinematics: a metal-rich population centred on $\textrm{[Fe/H]}\sim 0.25$, rapidly rotating and dynamically cold, and a metal-poor population centred on $\textrm{[Fe/H]}\sim-0.3$, dynamically hotter and more slowly rotating. These authors associated the metal-rich component to a B/P bulge with disk origin and suggested that the metal-poor component is an old classical spheroid, which in Baade’s window would account for about half of the stars in the bulge.

This view of a prominent classical bulge in the Galaxy has been challenged by the ARGOS survey and its associated series of papers. Using the full ARGOS sample, Ness et al. (2013a) showed that the metallicity distributions of stars in different fields across the inner Galaxy were well represented by five different components whose relative fractions change from field to field, thus explaining the vertical metallicity gradient of the bulge as a change in the relative fractions of the components. For each of these components, labeled from A to E in the order of decreasing metallicity they speculate a different origin as summarized below:

• •

A $(\textrm{[Fe/H]}\sim 0.1)$: Metal-rich cold component of the B/P bulge, concentrated to the plane and associated to a more recent star formation episode in the early disk than B.

• •

B $(\textrm{[Fe/H]}\sim-0.3)$: Main component of the B/P bulge originating from the thin disk before the bar formation and secular evolution.

• •

C $(\textrm{[Fe/H]}\sim-0.7)$: More metal-poor thicker disk-like population that does not significantly support the B/P shape of the bulge. C is associated to the hot and less dynamically responsive inner thick disk at the time when the bar formed.

• •

D $(\textrm{[Fe/H]}\sim-1.2)$: Metal weak thick disk, similar to the thick disk seen in the solar neighbourhood.

• •

E $(\textrm{[Fe/H]}\sim-1.7)$: Very metal-poor stars, associated to either a small classical bulge component or the inner part of the stellar halo.

The interpretation of the origin of the different ARGOS components has been discussed by Di Matteo et al. (2014) and Di Matteo et al. (2015) who concluded from analysing the formation of B/P bulges in N-body simulations that A and B have both a thin disk origin, with B formed on average at larger radii than A while C is likely to be the old thick disk already in place at the time of bar formation. However, this latter conclusion has been questioned by Debattista et al. (2017) who showed recently that co-spatial disk stellar populations with different radial dispersions can be separated by the bar formation into different orbit distributions.

Stars with $\textrm{[Fe/H]}\leq-1$, such as traced by the RRLyrae stars in the bulge, also have a controversial origin. Bulge RRLyrae have been found by Pietrukowicz et al. (2015) to be bar shaped and centrally concentrated (but see also Dékány et al. 2013 who find instead that they form a spheroid). From spectroscopic measurements, Kunder et al. (2016) found that they form a hot population, with a velocity dispersion of $\sim 130\,\rm{km\,s^{-1}}$ and rotate only slowly or not at all. They conclude that bulge RRLyrae stars are consistent with forming a small classical bulge component. Instead Pérez-Villegas et al. (2017) showed that both their density and kinematics are consistent with what is expected for the inner part of the stellar halo that co-evolved with the galactic bar.

In this context of controversies, the goal of this paper is two-fold. First, we present an extension of the M2M method for constructing chemodynamical models of galaxies by attaching chemical variables to N-body orbits in a self-consistent dynamical model. Such M2M chemodynamical models have the power to extract maximum information from the data and will be a valuable tool in the future for modelling the data provided by ongoing and future large spectroscopic surveys in the Galaxy. We focus here on the metallicity distribution but future applications could extend the modelling to other chemical quantities, such as for example $[\alpha/\rm{Fe}]$. This is described in Section 2. Then we investigate what currently available data can tell us about the present state of the inner Galaxy in the seven-dimensional extended phase-space of position, velocity and metallicity. We use the spatial and kinematic variations of the metallicity of stars seen in the ARGOS and APOGEE DR12 data to construct an M2M chemodynamical model of the galactic bar region. We base the underlying dynamics on the model of P17 which is a very good representation of the inner Galaxy and provides us with a library of N-body orbits on to which we attach and adjust a metallicity distribution. We describe this fitting procedure in Section 3 and show in Sections 4 and 5 the spatial distribution and kinematics of stars in the inner Galaxy for different metallicities. We finally discuss possible formation scenarios for the bulge, bar and inner disk, based on our chemodynamical model and in the light of other works and conclude in Section 6.

2 Made-to-measure chemodynamical modelling

2.1 The Made-to-Measure method

The M2M method is a particle-based modelling technique that allows the creation of constrained equilibrium models (Syer & Tremaine, 1996). Initially developed to tailor N-body initial conditions (Syer & Tremaine, 1996; Dehnen, 2009), the M2M method was adapted by De Lorenzi et al. (2007) to create dynamical models of real galaxies by fitting observational data. The idea behind the M2M method is to slowly adapt the particle weights of a N-body model in order to reproduce a given set of data. The N-body weights, representing mass elements in phase-space in the classical N-body approach are simultaneously seen as weighting factors of the N-body orbits on which the particles are. The M2M method has the advantage of being applicable to complex systems like barred galaxies where other modelling techniques such as moment-based methods (Binney et al., 1990; Cappellari et al., 2009), classical distribution function-based methods (Dejonghe, 1984; Qian et al., 1995), action-based methods (Binney, 2010; Sanders & Binney, 2013) or orbit-based methods (Schwarzschild, 1979; Thomas et al., 2009) would be either not applicable or not currently practical.

Formally, we represent the equilibrium state of the Galaxy by its distribution function $f(\bm{x,v})$ and write any observable $y$ as

[TABLE]

where $\bm{x,v}$ is the phase-space vector and $K$ is called the kernel of the observable. In the M2M method, $f(\bm{x,v})$ is approximated by a discrete sample of $N$ particles with particle weights $w_{i}(t)$. Equation 1 is then evaluated by

[TABLE]

where $(\bm{x}_{i}(t),\bm{v}_{i}(t))$ are the phase-space coordinates of particle $i$ at time $t$. In practice, $y(t)$ is often noisy due to the limited number of particles in the model. A convenient way to significantly reduce the particle noise is to replace the observable $y(t)$ by the temporally smoothed observable $\tilde{y}(t)$, defined as

[TABLE]

where $\alpha$ sets the temporal smoothing time-scale.

The M2M method consists of adjusting the particle weights $w_{i}$ in order to maximize a given profit function $F$. This is done through a simple gradient ascent method where the particle weights evolve with time according to

[TABLE]

with $\varepsilon$ a numerical factor that sets the typical time-scale of the weight evolution. The profit function $F$ usually consists of a $\chi^{2}$ term that drives the model towards some data plus possibly a regularization term for minimizing the broadening of the weight distribution. After a decade of refinement, this modelling technique is now mature and has been used in both the 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; Portail et al. 2015a; P17). Pioneering work in M2M chemodynamical modelling has also been done recently by Long (2016) who modelled absorption line strengths in four elliptical galaxies.

2.2 Modelling the metallicity of stars

In regular M2M modelling, $f$ is simply a function of phase-space. The method can be extended to chemodynamical models by extending the phase-space to chemical dimensions. Introducing the metallicity M into the phase-space, Equation 2 becomes

[TABLE]

In all the following, we use the term ‘metallicity’ to mean the iron abundance of a star, denoted alternatively by [Fe/H] or shorter M for readability of the equations. Several choices can be made in constructing a particle-based representation of $f(\bm{x,v},\textit{M})$. One possibility is to assign a single metallicity value $\textit{M}_{i}$ to each particle, and approximate $f$ as the sum of delta functions $f(\bm{x,v},\textit{M})=\sum_{i}w_{i}\times\delta(\bm{x}-\bm{x}_{i})\,\delta(\bm{v}-\bm{v}_{i})\,\delta(\textit{M}-\textit{M}_{i})$. However, in the overall potential, the particles represent orbits that can in reality be populated by many stars of different metallicities. In such a description, modelling a wide distribution of metallicity populating near-identical orbits requires many particles, with equivalent phase-space coordinates but different metallicities. Instead, it is more efficient to adopt for each particle a metallicity distribution function (MDF), denoted $f_{\textit{M},i}$, and write the distribution function as

[TABLE]

The space of all possible MDFs is of infinite dimension and has to be approximated by a finite dimensional space. We model $f_{\textit{M},i}(\textit{M})$ as an expansion on some mathematical set of elementary metallicity distribution functions (EMDFs), common for all particles, and thus represent $f_{\textit{M},i}(\textit{M})$ by a discrete set of numbers. The choice of the number and form of the EMDFs is constrained by both physical and technical considerations. First, as discussed later in Section 3.2, stars with different metallicities have different luminosity functions and thus usually have a different probability to be detected in a given survey depending on its selection criteria. The set of EMDFs should allow a faithful representation of the full MDF of stars in the Galaxy, and not for example only in a limited metallicity range. Secondly, the number of parameters used to represent the particle MDFs should be large enough to capture all relevant variations with metallicity shown by the fitted data, but also be as small as possible to avoid degeneracies and computational overhead.

From the ARGOS data, Ness et al. (2013a) showed that the MDF and its field-to-field variations across the bar region could be well represented by only five numbers, representing the varying amplitudes of five Gaussian MDFs with approximately fixed means and dispersions. Hence those Gaussian MDFs, denoted $\alpha$, $\beta$, $\gamma$, $\delta$ and $\epsilon$ in the order of decreasing metallicity form a suitable set of EMDFs for modelling the full MDF of any particle. Note that we do not need to assume that the Gaussian EMDFs of Ness et al. (2013a) are the MDFs of actual physical components of the Galaxy with different origin and history. Their use is purely restricted to a mathematical set of EMDFs that allow modelling reasonable particle MDFs with only a few parameters per particle, as demonstrated later in Figure 3.

For efficiency we neglect the Gaussian $\epsilon$ since it represents only about half a percent of all stars in the ARGOS survey and is thus very poorly constrained by the data. Note that the Gaussian EMDFs $\alpha$-$\delta$ were originally named A-D in Ness et al. (2013a). We chose here a different naming convention to highlight the difference between the Gaussian EMDFs $\alpha$-$\delta$, and the metallicity bins A-D introduced by Ness et al. (2013b) and used later in Section 3.

The MDF of each particle is then modelled as

[TABLE]

where $w_{i,c}$ are the chemical weights, i.e. weights of the Gaussian EMDFs $c$ in the MDF of particle $i$, $\textit{M}^{0}_{c}$ and $\sigma_{c}$ are the mean and standard deviation of the EMDFs and $\mathcal{G}_{[\mu,\sigma]}$ is the Gaussian distribution of mean $\mu$ and standard deviation $\sigma$. Adopted values for $\textit{M}^{0}_{c}$ and $\sigma_{c}$ are quoted in Table 1 and correspond to the means of the values determined by Ness et al. (2013a) (see their table 3).

From the MDF of every particle, we can evaluate any metallicity dependent observable $y_{j}$ on the model by

[TABLE]

where

[TABLE]

The purpose of the new kernel function $K_{j}(\bm{x}_{i},\bm{v}_{i},\textit{M})$ is to limit the observable $y_{j}$ to a given volume in phase-space and apply on the model any observational bias that can exist on the real data which the model is compared to. We describe in detail how we model these observational biases for the ARGOS and APOGEE surveys in Section 3.

2.3 M2M fit of the MDF

The M2M adaptation of the particle chemical weights $w_{i,c}$ is done in a similar fashion as the M2M evolution of the physical weights of the particles. The chemical weights $w_{i,c}$ are evolved in the direction of the gradient of the profit function $F$ according to

[TABLE]

where $\varepsilon_{c}$ set the strength of the MDF fitting.

As a profit function we adopt a simple chi-square given by

[TABLE]

where $y_{j}$, $Y_{j}$ and $\sigma(Y_{j})$, defined in more detail in Section 3, are respectively the model observables, the data and its associated errors. The $w_{i,c}$ are additionally normalized to $\sum_{c}w_{i,c}=1$ at every step for every particle. Note that the weight adaptation described in Equation 10 requires some initial values for the chemical weights $w_{i,c}$ to start with, as discussed in more detail in Section 3.1. For simplicity, we keep the particle masses $w_{i}$ constant throughout the M2M fit and only adapt the chemical weights $w_{i,c}$. This is justified since our initial model is already a good dynamical model of the entire galactic bar region, as described in the next section.

The M2M fit is performed in three successive phases. First, the model observables are initialized by evolving the N-body model for a time $\rm{T_{smooth}}$ and computing the temporally smoothed model observables. Then the model is evolved while adapting the chemical weights $w_{i,c}$ for a time $\rm{T_{M2M}}$. Finally the model is relaxed for a time $\rm{T_{relax}}$, to test the stability of the M2M fit. By using an internal time unit (iu) in which a circular orbit at $4\,\rm{kpc}$ has an orbital period of unity, we adopt for $\rm{T_{smooth}}$, $\rm{T_{M2M}}$ and $\rm{T_{relax}}$ the respective values of 4, 40 and 16 iu. The temporal smoothing time-scale is set to $\alpha=1\,\rm{iu}^{-1}$ and the strength of the MDF fitting to $\varepsilon_{c}=10\,\rm{iu}^{-1}$.

The weight evolution described by Equation 10 is performed iteratively by a first-order integration with small time step $\delta t=10^{-3}\,\rm{iu}$. Between two iterations of the weight adaptation procedure, the N-body model is integrated forward using an adaptive drift-kick-drift leap-frog algorithm. The gravitational potential is computed directly from the particle distribution (stellar and dark matter particles) using the hybrid grid method of Sellwood & Valluri (1997), modified in P17 for allowing better vertical resolution in the bar region. We assume that the gravitational potential rotates at a constant pattern speed of $40\,\rm{km\,s^{-1}\,kpc^{-1}}$ corresponding to the pattern speed of the bar in the model of P17 that matches both the BRAVA data in the bulge and the ARGOS data in the bar region.

3 Dynamical model and data constraints on the MDF

In this section, we describe in detail the initial dynamical model used to provide the N-body orbits on to which we fit the particle MDFs in Section 3.1, together with the two spectroscopic data sets, the ARGOS survey in Section 3.2 and APOGEE surveys in Section 3.3. We then describe and compare our fiducial model to the data, and present variational models used to quantify systematic uncertainties.

3.1 Dynamical model of the galactic bar region

In P17, we presented the first non-parametric dynamical model of the entire bar region. This model reproduces the 3D density of RCGs in the bulge from Wegg & Gerhard (2013), the magnitude distributions of RCGs across the bar region from combination of the VVV, UKIDSS and 2MASS surveys (Wegg et al., 2015), together with stellar kinematics from the BRAVA, OGLE and ARGOS surveys. Including the kinematics as a function of distance from the ARGOS survey we were able to recover the bar pattern speed and the mass distribution of both stars and dark matter in the inner $5\,\rm{kpc}$ of the Galaxy, hence constraining the full effective potential in the galactic bar region. We base our work on this model that provides a gravitational potential and a self-consistent library of N-body orbits. Note that throughout the chemodynamical M2M fit the particles masses and the gravitational potential in the bar frame are both kept fixed.

Figure 1 shows the face-on surface density of this model. The bar rotates at $\Omega=40\,\rm{km\,s^{-1}\,kpc^{-1}}$, as required to fit both the BRAVA kinematics in the bulge and the ARGOS radial velocity field in the bar region, and has in the model a half-length of $5.30\pm 0.36\,\rm{kpc}$, in good agreement with the determination of bar the half-length in the Milky Way of $5.0\pm 0.2\,\rm{kpc}$ by Wegg et al. (2015).

We introduce the metallicity of stars in the dynamical model by specifying for each particle an initial value for the four chemical weights $w_{i,c}$. We call these initial weights the “prior” MDF although they do not correspond to a formal Bayesian prior. It is a priori unclear how much the final result of a chemodynamical M2M fitting depends on the assumed prior MDF, as the chemical weights $w_{i,c}$ are only allowed to change when travelling through a portion of space constrained by the data. Since we focus on the bar region that is well covered by the ARGOS survey, our fiducial choice for the prior MDF is the fractions that correspond to the overall MDF of all the ARGOS stars: $w_{i,\alpha}=0.23$, $w_{i,\beta}=0.43$, $w_{i,\gamma}=0.29$ and $w_{i,\delta}=0.05$. We consider alternative priors and quantify their effects on the final results in Section 3.4.

The bar is oriented at an angle of $28\degree$ to the line of sight towards the Galactic Centre (Wegg et al., 2015). Following the recommendations of Bland-Hawthorn & Gerhard (2016), we place the Sun at $R_{0}=8.2\,\rm{kpc}$ and assume that the Local Standard at Rest (LSR) is on a circular orbit at $V(R_{0})=238\,\rm{km\,s^{-1}}$. We adopt for the peculiar motion of the Sun in the LSR the values of $(U,V,W)=(11.1,12.24,7.25)\,\rm{km\,s^{-1}}$ as measured by Schönrich et al. (2010). Altogether these assumptions lead to a proper motion of the Galactic Centre of $30.5\,\rm{km\,s^{-1}\,kpc^{-1}}$ consistent with the total angular velocity on the Sun of $30.57\pm 0.43\,\rm{km\,s^{-1}\,kpc^{-1}}$ measured from proper motion of masers by Reid et al. (2014). All radial velocities quoted in this work are expressed in the Galactocentric inertial frame. We re-transformed the radial velocity measurements of the ARGOS and APOGEE surveys using our set of assumptions since they differ from the traditional transformations used by many spectroscopic surveys. In order to predict the selection function of the surveys and apply them when observing the model, we use the PARSEC isochrones (Bressan et al., 2012; Chen et al., 2014; Tang et al., 2014) and assume a $10\,\rm{Gyr}$ old population for a Kroupa IMF.

3.2 ARGOS as a function of distance and metallicity

The Abundance and Radial velocity Galactic Origin Survey (ARGOS) is a large spectroscopic survey of about $28000$ stars towards the galactic bulge and inner disk, designed to sample RCGs all the way from the nearby disk ($\sim 4.5\,\rm{kpc}$ from the Sun) to the far side of the Galactic Centre ($\sim 13\,\rm{kpc}$ from the Sun). The survey covers lines of sight at longitudes between $l=-20\degree$ and $l=20\degree$ and latitudes between $b=-5\degree$ and $-10\degree$, as well as a few northern fields, as shown in Figure 2. The ARGOS stars are selected from the 2MASS catalogue according to a procedure fully described in Freeman et al. (2013) for which we modelled the selection function in P17. In brief, the ARGOS selection function consists of (i) a weight for each star and (ii) a function of the extinction corrected magnitude $C({K_{s}}_{0})$ for each field. The set of weights comes from the fact that for each field an equal number of stars were selected in three magnitude bins in order to sample the full range of distance from the nearby disk to the far side of the bulge. The selection function as a function of magnitude represents the probability for a star to pass the various selection criteria of the survey, thus taking into account the magnitudes and colour cuts, high-quality imaging criteria and incompleteness of the input catalogue that altogether bias the survey towards nearby stars. Details about the modelling of the ARGOS selection function can be found in P17.

From the medium resolution spectra of the ARGOS stars, Ness et al. (2013a) estimated various stellar parameters including the radial velocity and the metallicity with a typical error of $1\,\rm{km\,s^{-1}}$ and $0.1\,\rm{dex}$ respectively. Since the ARGOS stars are primarily RCGs which are approximate standard candles, we can infer accurate stellar distances from their photometry. We adopt an absolute magnitude of $M_{K_{s},\rm{RCG}}=-1.72$ as used by Wegg & Gerhard (2013) for the bulge, and bin stars as a function of distance. In P17, we constructed the ARGOS radial velocity field by computing for each survey field the mean velocity and velocity dispersion of the stars in several distance modulus bins of $0.25\,\rm{mag}$ along the line of sight. We adopt here the same distance bins as P17.

Since the number of stars per distance bin is typically $\sim 80-100$ in each field, we can split the bins further in the metallicity space. Following Ness et al. (2013b), we define four metallicity bins of $\Delta\textrm{[Fe/H]}=0.5\,\textrm{dex}$ called A-D by $0.5\geq\textrm{[Fe/H]}>0$ (A), $0\geq\textrm{[Fe/H]}>-0.5$ (B), $-0.5\geq\textrm{[Fe/H]}>-1$ (C) and $-1\geq\textrm{[Fe/H]}>-1.5$ (D). Binning the stars avoids correlations between data points and is therefore well suited to the $\chi^{2}$ minimization of the M2M method.

The difference between the Gaussian EMDFs $\alpha-\delta$ introduced in Section 2.2 and the metallicity bins A-D introduced here is best described by Figure 3. In this figure, the histogram shows a typical MDF observed in a field towards the bulge. By binning the stars in the four bins A-D, we obtain a low resolution version of the MDF. However, this low resolution MDF is not a faithful representation of the real underlying MDF, which for example in bin A has many more stars close to the lower end of the bin than stars close to the higher end of the bin. A more faithful representation of the MDF is obtained using the set of EMDFs $\alpha-\delta$. Indeed, from the counts in the four bins A-D, we can determine four unique numbers $w_{c}$, corresponding to the weights of the Gaussian EMDFs $\alpha$-$\delta$, for which the integral of the total MDF has the same counts as the stars in each of the metallicity bins. Summing up the Gaussian EMDFs, we obtain a full MDF, using only these four numbers, that is a very good representation of the underlying histogram as shown in Figure 3. The relation between the counts in bins A-D and the chemical weights of Section 2.2 is given by a $4\times 4$ matrix, whose elements correspond to the contribution of each Gaussian EMDFs to each bin. Given our choices of bins and Gaussian EMDFs, this matrix is mostly diagonal.

For a field indexed by $j$, distance bin indexed by $d$ and metallicity bin indexed by $m$ we compute the following three quantities:

The fraction $f_{j,d,m}$, of the observed stars in distance bin $d$ of field $j$ that have metallicities within bin $m$ 2. 2.

Mean velocity $v_{j,d,m}$ of the stars considered above 3. 3.

Velocity dispersion $\sigma_{j,d,m}$ of the stars considered above

Errors in $f_{j,d,m}$, $v_{j,d,m}$ and $\sigma_{j,d,m}$ are estimated respectively as $1/\sqrt{n}$, $\sigma_{j,d,m}/\sqrt{n}$ and $\sigma_{j,d,m}/\sqrt{2\,n}$ for $n$ the number of stars in bin $(j,d,m)$. Since $\sigma_{j,d,m}$ is biased towards lower values for small number statistics, we exclude from further considerations the kinematics of bins with less than five stars. In some fields and metallicity bins, there are not enough stars in the ARGOS sample to also bin in distance. Hence for each field $j$ and metallicity bin $m$, we compute in addition the three similar quantities $f_{j,m}$, $v_{j,m}$ and $\sigma_{j,m}$ corresponding to the fraction, mean velocity and velocity dispersion of stars integrated along the line of sight. $f_{j,d,m}$, $v_{j,d,m}$ and $\sigma_{j,d,m}$ are shown together with our fiducial chemodynamical model described in Section 3.4.

In order to compare the model to the data, we need to form model observables that are the model equivalent quantities of $f_{j,d,m}$, $v_{j,d,m}$ and $\sigma_{j,d,m}$. To do so, we first turn the model particles into mock stars, then apply the ARGOS selection function and finally estimate the distance of the remaining mock stars assuming they are RCGs, as done for the observed stars. This is described in detail in P17 and will only be outlined here. For a given particle, we simulate a mock $10\,\rm{Gyr}$ old stellar population for a Kroupa IMF and use the PARSEC isochrones in combination with the selection function to predict the magnitude distribution of ARGOS mock stars, i.e. RCGs on top of a background of giant stars. We then apply the distance estimation based on the RCG magnitude as we do for the ARGOS stars and compute the distance distribution of observable mock stars drawn from the particle. This distribution, called $f_{i}(\mu)$ in P17, encloses all selection and stellar population effects but does not depend on metallicity. Here, we extend $f_{i}(\mu)$ to $f_{i}(\mu,\textit{M})$ taking into account the fact that the luminosity function of different metallicities is slightly different: more metal-rich RCGs are brighter in the near-infrared than more metal-poor RCGs. The effect of metallicity is illustrated in Figure 4 where we show the distance distributions of observable mock stars corresponding to a particle located at $\mu=14.0$ for different metallicity distributions corresponding to the four Gaussian EMDFs $\alpha-\delta$. More metal-rich stars appear at closer distances than more metal-poor stars.

Motivated by the general form of an observable given in Equation 8, we finally define our NMAGIC observables as the number count and first and second mass-weighted velocity moments. The NMAGIC kernel of the observables to be used in Equation 9 are then given by the following equations for a field $j$, a distance bin $d$ and a metallicity bin $m$:

[TABLE]

[TABLE]

and

[TABLE]

where $\bm{z}_{i}=(\bm{x}_{i},\bm{v}_{i},\bm{M}_{i})$ is the extended phase-space vector of particle $i$ and ${v}_{i}$ its radial velocity. $W_{j,d,m}^{\rm{ARGOS}}$ is given by

[TABLE]

and $\delta_{j,d,m}^{\rm{ARGOS}}(\bm{z}_{i})$ by

[TABLE]

with $\mu(d)$ and $\mu(d+1)$ the boundaries of the distance modulus bin $d$. NMAGIC kernels for data integrated along the line of sight are defined in a similar way.

These kernels are then inserted into equations (8) and (9) to give the full expression for a model observable. In the case of the model observable ‘number count in field $j$, distance bin $d$ and a metallicity bin $m$’, denoted $y_{j,d,m}^{\rm{ARGOS},0}$, this insertion sums up to

[TABLE]

and similarly for the first and second mass weighted velocity moments.

3.3 APOGEE as a function of metallicity

The Apache Point Observatory Galactic Evolution Experiment (APOGEE) is a high resolution spectroscopic survey designed to sample a large number of stars ($\sim 1.5\times 10^{5}$) in all possible galactic environments, from the stellar halo to the disk and the bulge (Eisenstein et al., 2011; Majewski, 2012; Nidever et al., 2012). APOGEE operates in the near-infrared and is therefore able to pierce through highly dust obscured regions, making it particularly useful close to the galactic plane. Among the $\sim 1.5\times 10^{5}$ stars about $10^{4}$ are located in fields towards the bulge and long bar, some of which are directly in the galactic plane. We thus use APOGEE to complement ARGOS in the galactic plane, in the fields shown in Figure 2. We use the stellar parameters and spectroscopic distances of stars in the APOGEE bulge and bar fields determined by Ness et al. (2016) from the publicly available spectra of DR12 (García Pérez et al., 2013; Holtzman et al., 2015; Shetrone et al., 2015) using the data-driven method The Cannon (Ness et al., 2015). The APOGEE bar and bulge fields do not contain enough stars for binning both in distance and metallicity. We thus consider for each field only a single distance bin ranging from $4\,\rm{kpc}$ to $12\,\rm{kpc}$ from the Sun and bin the stars in the four metallicity bins A-D. As in Section 3.2, we construct for each APOGEE field the quantities $f_{j,m}$, $v_{j,m}$ and $\sigma_{j,m}$ corresponding respectively to the fraction of observed stars of field $j$ that have metallicities within bin $m$, their mean velocity and their velocity dispersion.

The APOGEE survey is targeted at bright red giants selected from the 2MASS catalogue according to a complex procedure fully described in Zasowski et al. (2013). To fairly compare the model to the APOGEE data, we need to model the distance and metallicity biases introduced by the APOGEE selection criteria. Luckily Bovy (2016) and Bovy et al. (2016a) computed and made available the APOGEE selection function for any kind of tracer stellar population. Using the PARSEC isoschrones to simulate a mock $10\,\rm{Gyr}$ old stellar population with a Kroupa IMF at a given metallicity M, the APOGEE selection function provides us with the fraction of stars observable by the APOGEE survey as a function of distance $\mu$. This fraction corresponds to the quantity $f_{i}(\mu,\textit{M})$ already used to define the ARGOS kernels in Equation 16. Hence, our APOGEE kernels are basically the same as Equation 16 with the exception that we have only one large bin in distance covering all the way from $4\,\rm{kpc}$ to $12\,\rm{kpc}$ along the line of sight. $f_{j,m}$, $v_{j,m}$ and $\sigma_{j,m}$ are shown together with our fiducial M2M chemodynamical model in the next subsection.

3.4 Fiducial and variation models

Our fiducial model is obtained by fitting the best-fitting dynamical model of P17 to the ARGOS and APOGEE (DR12) data using the M2M chemodynamical method described in Section 2. ARGOS and APOGEE have a different spatial coverage and complement each other in constraining the entire galactic bar region, as shown by Figure 2. The evolution of $\chi^{2}$ during the M2M fit is shown in Figure 5, for observables in the four metallicity bins separately. It shows the typical evolution of $\chi^{2}$ in M2M modelling, staying constant during $\rm{T_{smooth}}$, decreasing during $\rm{T_{M2M}}$ while fitting, to finally slightly increase again and stabilize during $\rm{T_{relax}}$.

The three panels of Figure 6 show respectively the comparison between our fiducial model and the ARGOS data for the fraction of stars in each metallicity and distance bin $f_{j,d,m}$, their mean radial velocity $v_{j,d,m}$ and radial velocity dispersion $\sigma_{j,d,m}$. The distance distribution of the metallicity bins B and C are remarkably well fitted. Our model slightly overpredicts A in the central field along the minor axis at $b=-5\degree$, which is also the field most affected by selection effects due to the important extinction and crowding. The kinematics is also well fitted and reproduces the main trends as a function of metallicity already described by Ness et al. (2013b): B is a hotter and slightly faster rotating replica of A while C and D are hotter and slower than both B and A at $|b|=5\degree$. We note however that the model systematically overpredicts the dispersion of the very cold metallicity bin A by about $10\,\rm{km\,s^{-1}}$.

Similarly, the three panels of Figure 8 show respectively the comparison of our model with the fraction $f_{j,m}$, mean radial velocity $v_{j,m}$ and radial velocity dispersion $\sigma_{j,m}$ as computed from the APOGEE data. In all metallicity bins, the model performs well at reproducing the data. We note however in metallicity bin A that the model tends to systematically underpredict the fraction of stars, an opposite trend to that found in the ARGOS data. This disagreement is probably due to a slight inconsistency between the two data sets, possibly originating from systematic effects in the selection function or from a difference in the upper end of the metallicity scales between ARGOS and APOGEE.

Our fiducial model is based on modelling assumptions whose effects need to be quantified before analysing the spatial and kinematic distribution of stars in different metallicity bins. To this end we repeated the modelling for five alternative models with different assumptions as described below:

(i) and (ii):

As already mentioned in Section 3.1, part of the Galaxy is not constrained by our data sets and therefore the model partly depends on the prior MDF distribution assumed. We thus consider two variations of the prior distribution, the first corresponding to “flat-priors” where we initialize all $w_{i,c}$ to $0.25$ and a second more complex corresponding for each particle to the values determined from the MDF at the position of the particle in the Besançon Galaxy Model (Robin et al., 2003).

(iii):

The APOGEE sample is subject to a severe metallicity bias in the inner Galaxy, as discussed in the appendix of Hayden et al. (2015) and shown very clearly in fig. 2 of Ness et al. (2016). Only a few stars in the APOGEE fields towards the bulge are actually in the bulge, and those are preferentially metal-poor stars. This bias against metal-rich stars is already modelled by the APOGEE selection function but is also somewhat uncertain since it relies on stellar population models and assumptions about the IMF and age of the stellar population in the inner part of the Galaxy. Hayden et al. (2015) showed that this bias could be reduced by removing from the sample all stars with $\log g\leq 1\,\textrm{dex}$. We thus repeat the fitting using this cut in surface gravity and modify the selection function accordingly.

(iv) and (v):

When using APOGEE in combination with ARGOS we might expect a systematic difference in the metallicity scales of the two surveys. Such differences in scale are always present between different spectroscopic surveys and Smiljanic et al. (2014) showed that even when analysing the exact same spectra, different working groups would measure [Fe/H] with systematic variations of $\sim 0.1\,\textrm{dex}$. García Perez et al. (in preparation) find a possible indication for a difference in metallicity scale between ARGOS and APOGEE, based on about 100 stars in common between the two surveys. However, this comparison is limited to the brightest ARGOS stars in only one field and thus does not allow the precise characterization of the overall relation between the metallicity scales of the two surveys. Thus, we limit ourselves here to considering two alternative APOGEE data sets where values for [Fe/H] are offsets by $\pm 0.1\,\textrm{dex}$.

The difference between these variation models is found to be relatively small in the region covered by our data, as shown by the coloured areas in Figure 6. In all the following, we will use the standard deviation of the results provided by the fiducial model plus its five variations as a measure of the systematic uncertainty affecting our results. Table 2 gives an overview of all models and provides the stellar mass within $5\,\rm{kpc}$ from the Galactic Centre in all metallicity bins for the fiducial and variations models.

4 Spatial distribution of the metallicity in the inner Galaxy

In this section, we show how the spatial distributions of stars in the Galaxy vary as a function of metallicity. Since all particles carry a full continuous MDF, we can in principle observe the model at any metallicity. To ease the comparison with other work, we present our model in the four metallicity bins A - D. To recall, the bins A, B, C and D correspond respectively to $0.5\geq\textrm{[Fe/H]}>0$ (A), $0\geq\textrm{[Fe/H]}>-0.5$ (B), $-0.5\geq\textrm{[Fe/H]}>-1$ (C) and $-1\geq\textrm{[Fe/H]}>-1.5$ (D).

4.1 The 3D distribution of metallicity

The surface density in face-on and side-on projections of our fiducial model in the four metallicity bins A-D are shown in Figure 9. White regions indicate parts of the Galaxy where the systematic uncertainties are larger than $30\%$, or outside the central $5.5\,\rm{kpc}$. Bin A stars are concentrated to the plane and are strongly bar shaped. The hammer-like shape of the bar ends is similar to the shape of the superthin bar component, discovered in the VVV + UKIDSS + 2MASS data by Wegg et al. (2015). This hammer-like shape is confined to the plane and although being very clear in the mid-plane density, it is usually absent in the total surface density because of the mass distribution above and below the plane (see Figure 1 for example). We thus speculate that the superthin bar component is significantly metal rich. Bin B stars are also strongly bar shaped but are less concentrated to the plane and do not significantly contribute to the superthin bar.

For metallicity bins C and D, the stellar densities are much more extended vertically, are both less bar shaped but do not show any significant B/P shape in the bulge. Outside the central $\,\rm{kpc}$, from its side-on projection in Figure 9, the stellar density in bin C is found to have a vertical flattening of $\sim 0.3$. Hence, at these radii, it has the shape of a thick disk. However in the central $\,\rm{kpc}$, stars in bins C and D appear to be consistently more centrally concentrated than the inward extrapolation of their outer exponential density profiles, as discussed in more detail in Section 4.3. Note also that only a small fraction of stars from our data sets belong to bin D (<5% of the stars in ARGOS) and therefore our results for this metallicity bin are more uncertain.

Another perspective on the spatial distribution of the different metallicities in the Galaxy is given in the galactic coordinates projections of Figure 10. Stars in A and B are found to dominate in the bulge and close to the plane and the distribution of stars in B is found to resemble a thicker replica of the distribution of stars in A. This figure clearly shows the vertical thickening of the stellar density with decreasing metallicity, resulting in a vertical metallicity gradient as discussed in the next subsection.

4.2 Metallicity gradients in the bar region

The presence of a vertical metallicity gradient along the minor axis of the galactic bulge was first reported by Minniti et al. (1995). They derived the metallicity of many K giants at $1.5-1.7\,\rm{kpc}$ from the galactic plane and, by comparing to previous metallicity measurements in Baade’s window, concluded the presence of a vertical gradient estimated to $-0.35\,\textrm{dex}.\,\rm{kpc}^{-1}$. Later studies confirmed the presence of a gradient although its magnitude has a large scatter in the literature. Zoccali et al. (2008) observed a gradient of $-0.6\,\textrm{dex}.\,\rm{kpc}^{-1}$ from bulge K giants between $b=-4\degree$ and $b=-6\degree$, while Ness et al. (2013a) finds a gradient of $-0.45\,\textrm{dex}.\,\rm{kpc}^{-1}$ from 2000 ARGOS stars between $b=-5\degree$ and $b=-10\degree$. In addition, Rich et al. (2007) found evidence for a flattening of the metallicity gradient between $b=-1\degree$ and $b=-4\degree$. This flattening is also well recovered by the first continuous metallicity maps of the bulge of Gonzalez et al. (2011b, 2013) who derived photometric metallicities of red giants across the entire bulge from 2MASS and VVV, finding an overall gradient of $-0.28\,\textrm{dex}.\,\rm{kpc}^{-1}$ for latitudes between $b=-2\degree$ and $b=-10\degree$. This flattening of the metallicity gradient for $|b|\leq 4\degree$ is probably at the origin of the large scatter of the measurements found in the literature, which all refer to different latitude ranges.

In Figure 11 we show maps of the mean metallicity and metallicity dispersion of our fiducial model in both galactic and side-on projections. In the plane for $|b|<5\degree$, the mean metallicity of the model is quite uncertain since in addition to having a broad MDF with a dispersion of $\sim 0.35\,\textrm{dex}$, it is only constrained by the APOGEE data which struggle to reach the bulge. Outside $|b|=5\degree$ the presence of a vertical gradient is very clear. Along the minor axis of the bulge, including only stars between $6\,\rm{kpc}$ and $10\,\rm{kpc}$ along the line of sight, we find a vertical gradient of the median metallicity in the bulge between $|b|=5\degree$ and $|b|=10\degree$ of $-0.45\,\textrm{dex}.\,\rm{kpc}^{-1}$. This is in good agreement with Ness et al. (2013a), which is expected since the model fits the ARGOS data well. This gradient is however larger than the photometric gradient measured by Gonzalez et al. (2013) although the latitude range in the latter study includes part of the strip at $|b|\leq 4\degree$ where the metallicity gradient flattens. Note that the gradient of the mean metallicity differs from that of the median; for our fiducial model, the former is only $-0.33\,\textrm{dex}.\,\rm{kpc}^{-1}$ in the same region of the bulge.

Interestingly, the vertical gradient of the median metallicity in the bulge between $|b|=5\degree$ and $|b|=10\degree$ of $-0.45\,\textrm{dex}.\,\rm{kpc}^{-1}$ drops only to $-0.37\,\textrm{dex}.\,\rm{kpc}^{-1}$ when limited to stars with $\textrm{[Fe/H]}\geq-0.5$. Hence, the vertical metallicity gradient in the bulge is mostly produced by a gradual change in the metallicity of metal-rich stars (i.e stars in A and B), and less by a gradual change in the relative amplitude of a unique metal-rich B/P bulge (stars in A and B together) to a more metal-poor surrounding (stars in C and D together) such as a thick disk (Bekki & Tsujimoto, 2011) or a classical bulge (Hill et al., 2011). This is also illustrated by Figure 11 in which the shape of the B/P bulge, supported by metal-rich stars, remains visible at all heights.

Several interpretations for the origin of this gradient can still be given, from the presence of physically distinct stellar populations (Ness et al., 2013a; Fragkoudi et al., 2017), from an initial radial gradient in the early disk (Martinez-Valpuesta & Gerhard, 2013; Di Matteo et al., 2014) or from a continuum of disk populations separated by the bar formation (Debattista et al., 2017).

Outside the bulge, the mean metallicity only slowly changes with radius, with slight signs of flaring for galactocentric distances larger than $3-4\,\rm{kpc}$, as also found by Bovy (2016) who showed from APOGEE RCGs that the scaleheight of more metal-poor stars increases with increasing radius.

4.3 A central concentration of metal-poor stars: classical bulge or inner stellar halo

Figure 9 shows that the stellar densities in both bins C and D are more centrally concentrated in the inner $\,\rm{kpc}$ of the Galaxy than they are on larger scales. This is shown in more detail in Figure 12 where we plot the vertical density at different points along the major axis of the bar for all four metallicity bins. In this figure, stars in bins C and D have similar vertical density profiles for all cuts. For cuts located at major axis distance $x$ larger than a $\,\rm{kpc}$, they share the vertical exponential profile of a thick component with scaleheight $\sim 500\,\rm{pc}$. However, in the inner $\,\rm{kpc}$ of the Galaxy, an extra component with short scaleheight rises towards the Galactic Centre and is seen consistently only in the metal-poor bins C and D. For the central cut, the scaleheight of this inner metal-poor component is only $\sim 70\,\rm{pc}$, and it vanishes below the large-scale thick disk density at $|z|\geq 600\,\rm{pc}$.

The presence of a central concentration of metal-poor stars in the inner $\,\rm{kpc}$ of the Galaxy has already been reported by several authors using different data sets. Both Pietrukowicz et al. (2015) and Dékány et al. (2013) found, respectively from OGLE and VVV data, that RRLyrae stars in the bulge were centrally concentrated, although they disagree about the shape of the RRLyrae density distribution. In addition, Schultheis et al. (2015) reported the presence of a peak of metal-poor stars in the APOGEE data within the central degree of the Galaxy. Finally, Zoccali et al. (2016) combined spectroscopic data from the recent GIBS survey with the bulge stellar density from Valenti et al. (2016) and find a central concentration of metal-poor stars at $b=-2\degree$ that vanishes at $b=-4\degree$.

Interpreting the 3D distribution and central concentration of stars in bins C and D in terms of physical components requires more data and is beyond the scope of this paper. Theoretically at least three metal-poor components are expected or suspected to contribute to the metallicity bins C and D in the inner Galaxy. The thick disk, locally metal-poor with a flat radial metallicity gradient (Bensby et al., 2007; Cheng et al., 2012; Bovy et al., 2016b) is expected to contribute to both C and D in the inner Galaxy, and could possibly be centrally compressed because of the deep gravitational potential of the nuclear bulge (Launhardt et al. 2002; P17). In addition, the stellar halo with its steep radial profile (Wetterer & McGraw, 1996; Jurić et al., 2008) is also a possible candidate (Pérez-Villegas et al., 2017), although it is not clear that its mass in the bulge would be sufficient to provide a significant contribution to C and D. Finally, the central concentration of metal-poor stars could also indicate the presence of a classical bulge component, a fossil of early mergers but whose existence has not yet been established.

To conclude, from the 3D mass distributions of stars with metallicity in the bins C and D obtained by our modelling, it seems that on large-scale bin C stars are mainly distributed as a thick disk. This is supported by its flattening of only $\sim 0.3$, quite extreme for a significant contribution from a classical bulge or from the inner stellar halo, according to both observations of external galaxies (Kormendy & Kennicutt, 2004) and simulations (flattening $\sim 0.5$ for low-mass classical bulge in Saha et al. 2012; or $\sim 0.6$ for the inner stellar halo in Pérez-Villegas et al. 2017). Inside the central $\,\rm{kpc}$, an extra metal-poor component dominates but a number of interpretations are possible. More data would be required to investigate this further.

4.4 Support to the B/P bulge and bar

In the ARGOS bulge sample, only stars in A and B contribute to the B/P shape. This is demonstrated by Figure 13 where we show the stellar density times the cube of the line-of-sight distance ($\rho\times r^{3}$), as a function of distance modulus for the four metallicity bins in a minor axis field at $b=-5\degree$. The quantity plotted in this figure $\rho\times r^{3}$ is proportional to the number count of stars as a function of distance modulus for a standard candle stellar population. This figure can thus be compared to the similar diagram of Ness et al. (2012) where they show from the raw ARGOS data directly that only stars in A and B have a bimodal magnitude distribution, the so-called ‘split red clump’, and that therefore only stars in A and B support the B/P shape of the bulge. Our model confirms these results, taking also into account the ARGOS selection function.

From the face-on views of Figure 9 it is clear that stars of all metallicities support the bar, but in terms of mass, most of the support to the bar is provided by metal-rich stars. This is shown by the upper panel of Figure 14 where we plot the radial profiles of the $m=2$ Fourier coefficients in the four metallicity bins relative to the overall $m=0$ Fourier coefficient. Metal-poor stars in bins C and D provide only a weak mass support to the bar, mostly because they represent only a small fraction of the mass overall (see Table 2). However, there is still a significant fraction of the metal-poor stars that does support the bar, as shown by the lower panel of Figure 14 where we plot the radial profiles of the $m=2$ Fourier coefficients in the four metallicity bins relative to the $m=0$ Fourier coefficients of the respective metallicity bin.

Stars that do not strongly support the bar are instead building the inner disk. To isolate the bar from the inner disk we classify the N-body orbits of the model of P17 into two classes, the bar-supporting orbits and the not bar-supporting orbits. Following the method of Portail et al. (2015b), we identify the bar-supporting orbits as the orbits for which $f_{r}/f_{x}=2\pm 0.1$, where $f_{r}$ and $f_{x}$ are respectively the dominant frequency of the oscillations in cylindrical radius and in the bar major axis coordinate along the orbits.

In Figure 15, we show for bins A and B the fraction of the surface density in galactic coordinates that originates from bar-supporting orbits. Both metallicity bins exhibit an increase of bar-supporting orbits towards the corners of the bulge, as expected for orbits supporting the B/P shape of the bulge. Note that for bin A about half the stars in this “orbital X-shape” support the bar.

How are the bar and inner disk orbits divided between the different metallicity bins? Using our orbit classification, we find that the galactic bar is composed of $52\%\pm 3\%$ of stars with metallicities within bin A, $34\%\pm 1\%$ of stars with metallicities within bin B, and only $12\%\pm 1\%$ and $2\%\pm 1\%$ of stars with metallicities within bin C and D, respectively. For the inner disk within corotation, the fractions of stars in A, B, C and D are respectively $38\%\pm 3\%$, $47\%\pm 2\%$, $14\%\pm 1\%$ and $1\%\pm 1\%$. Note that although B is found in Ness et al. (2013b) to be the most populated metallicity bin in the ARGOS data, we find here that bin A stars give the main mass support to the bar and have a similar role as bin B stars in building the inner disk.

5 Kinematic distribution of the metallicity in the Galaxy

In this section, we show how stars in the different metallicity bins differ by their kinematics and orbital properties. As already pointed by Babusiaux et al. (2010) and Ness et al. (2013b), metal-rich stars with $\textrm{[Fe/H]}\geq-0.5$ (i.e, in bins A and B) have very different kinematics than more metal-poor stars with $\textrm{[Fe/H]}\leq-0.5$ (i.e in bins C and D). We now analyse in detail these two main groups separately using our fiducial chemodynamical model of Section 3.4.

5.1 Bar-like kinematics of metal-rich stars

The kinematics of stars in bin A and B are found in the ARGOS data to be similar to each other, as shown by Ness et al. (2013b). They are both rapidly and cylindrically rotating, i.e with only weak variations as a function of latitude. They have similar dispersion profiles, both steeply rising towards the centre. In the ARGOS data, bin B stars are found to rotate slightly faster than bin A stars, with a velocity dispersion about $20\,\rm{km\,s^{-1}}$ higher at the centre.

The similarities between the kinematics of stars in these two metallicity bins have been interpreted as an indication of a common thin disk origin by Di Matteo et al. (2014). In their N-body simulations of bar forming disks, stars that formed on average at larger radii are mapped after bar formation and buckling into a population that exhibits faster rotation and higher dispersion. Thus, with an initially decreasing metallicity gradient in the disk, birth radius and metallicity become correlated and naturally reproduce the kinematic behaviour of bin A and B stars, thus arguing for a common disk origin for these stars.

This simple picture however does not entirely capture the trend we see in our fiducial chemodynamical model.

Figure 16 shows maps in galactic coordinates of the mean radial velocity and velocity dispersion for stars in bins A and B, integrating along the line of sight between $3$ and $12\,\rm{kpc}$ from the Sun. We see that B is indeed faster than A for $b\geq 5\degree$ where the ARGOS fields are located, but not in the plane for $b\leq 5\degree$ where we find the inverse trend, with A being about $10\,\rm{km\,s^{-1}}$ faster than B at $l=20\degree$. This result likely originates from the complex orbital structure of the in-plane bar, mostly populated by bin A stars. Although the study of the orbital structure of the bulge and bar in the dynamical model of P17 is beyond the scope of this paper, we can investigate the origin of this separate relative behaviour of A and B in the plane by experimenting with various cuts in the model. Doing so we find that the particles that are responsible for A being faster than B in the plane are located less than $2\,\rm{kpc}$ away from the bar major axis, and predominantly along the rear-side of the bar. These particles produce a high-velocity tail in the line-of-sight distributions, thus shifting the mean velocity in bin A to higher values. We thus speculate that the faster rotation of in-plane stars in bin A is likely to be caused by the complex orbital structure of the bar, in possible relation with the orbits suggested by Aumer & Schönrich (2015) to be at the origin of the high-velocity peak observed in the APOGEE data for $l=4\degree-14\degree$ and $|b|\leq 2\degree$ (Nidever et al., 2012). Note however that the model only shows a high-velocity tail, with no sign for a separate high-velocity peak. It appears only for $l\geq 8\degree$ and is thus unrelated to the possible presence of a $\sim\,\rm{kpc}$ scale nuclear disk as suggested by Debattista et al. (2015). A detailed orbit analysis would be required to investigate the cause of the fast rotation of bin A stars further.

5.2 Kinematics of the metal-poor stars

The metal-poor stars, often defined by $\textrm{[Fe/H]}\leq-0.5$, have been shown to be kinematically distinct from their metal-rich counterparts. They generally rotate more slowly and have a higher dispersion, showing that metal-poor stars cannot be seen as just the metal-poor tail of a single population with a broad MDF.

This is illustrated in more detail in Figure 17 where we show maps of the mean radial velocity and velocity dispersion of stars in bins C and D from our fiducial model. By comparing with Figure 16, both C and D appear to have a higher dispersion than their surroundings, except in the central few degrees where the metal-rich stars have larger dispersions, caused by their bar-supporting orbits. Debattista et al. (2017) showed recently that hotter stellar populations are less responsive to the bar formation and buckling. Similar conclusion is obtained by Di Matteo (2016) for vertically thicker stellar populations. Thus if C and D were initially populated by a hot and/or thick stellar population, this would explain their relative weak support to the bar and their absence of bar-like kinematic features. We note however that metal-poor stars with $\textrm{[Fe/H]}\leq-0.5$ do not share a unique kinematic behaviour but rather show smoothly evolving trends with metallicity. In the ARGOS data, bin C stars at $b=-5\degree$ reach a mean velocity of $\sim 120\,\rm{km\,s^{-1}}$ at $l=20\degree$ while bin D stars reach only $\sim 70\,\rm{km\,s^{-1}}$.

5.3 Vertex deviation in the bulge

Because the physical components present in the inner Galaxy (at least thick disk, stellar halo and possibly classical bulge) are likely to overlap in metallicity space, we expect smooth kinematic transitions as a function of metallicity. This is illustrated by Figure 18 where we show the vertex deviation in Baade’s window $\theta_{r,l}$ of our model compared to the data from Soto et al. (2007) and Babusiaux et al. (2010). The vertex deviation $\theta_{r,l}$ is defined by

[TABLE]

where $\sigma^{2}_{r,l}$, $\sigma_{r}^{2}$ and $\sigma_{l}^{2}$ are respectively the covariances of the radial velocity and proper motion dispersion in the longitudinal direction, and their two associated variances. The vertex deviation represents the angle of the principal axis of the velocity ellipsoid in non-isotropic cases. We selected only stellar particles between $6\,\rm{kpc}$ and $10\,\rm{kpc}$ away along the line of sight within a cone of half a degree centred on $l=1\degree$, $b=-3.9\degree$ and evaluated the vertex deviation of the model as a function of metallicity on a grid with a resolution of $0.1\textrm{dex}$.

It is clear in Figure 18 that $\theta_{r,l}$ is significantly non-zero for stars with $\textrm{[Fe/H]}\geq-0.5$ and gradually transitions to zero for more metal-poor stars. Our fiducial model reproduces the observed evolution of $\theta_{r,l}$ with metallicity well, providing a good fit to the data although not directly fitted to it.

The absence of a significant vertex deviation for $-1\leq\textrm{[Fe/H]}\leq-0.5$ (C) has been interpreted by Soto et al. (2007) and Babusiaux et al. (2010) as a signature of the presence of an old metal-poor spheroid. Our model shows instead that bin C stars have no significant vertex deviation $\theta_{r,l}$ in Baade’s window, which is shown in Figure 12 to be sufficiently far away from the Galactic Centre to be dominated there by the large-scale disk-like density of bin C stars found in Section 4.1. In addition, C is shown in Figure 17 to be cylindrically rotating with a rotation only slightly slower than that of the metal-rich stars.

We thus conclude that our interpretation for stars in the metallicity bin C forming a thick disk outside the central $\,\rm{kpc}$ is consistent with its hot kinematics and relatively fast rotation. The lack of significant vertex deviation for metal-poor stars in Baade’s window is caused by the large-scale thick disk distribution of bin C stars, with no need for a large classical bulge component.

Little can be reliably said from our model about the metallicity bin D since it is only poorly constrained by the data. D is found in the model to be cylindrically rotating, but at a slower rate than C and has a dispersion that varies only slightly in the inner $10\degree$ of the Galaxy. Differences between C and D are likely to highlight different contributions of the underlying metal-poor stellar populations described above. For example, the lower rotation of D with respect to C could be explained by a higher contribution in D of a non-rotating stellar halo.

6 Discussion and conclusion

We presented an extension of the Made-to-Measure (M2M) method for modelling self-consistent equilibrium chemodynamical models of galaxies. Its application in the inner Milky Way results in the first chemodynamical particle model fitted to chemo-kinematic data in the bulge and bar. Such chemodynamical models will be invaluable in the near-future for interpreting chemical abundances and stellar parameters provided by the many on-going and future large-scale spectroscopic surveys of the Galaxy. Indeed, M2M chemodynamical models allow to maximize the information extracted from spectroscopic data by constraining at once the spatial distribution and kinematics of stars with the orbital distribution of a self-consistent dynamical model. Since it is built upon the standard M2M method, it is flexible enough to be applicable to complex systems such as the inner barred galactic disk, where other modelling technique would be either not applicable or too complex. We applied the method based on the dynamical model of the bar region recently made by Portail et al. (2017), and successfully reproduced the spatial and kinematic variations of the metallicities of stars seen in the ARGOS and APOGEE (DR12) surveys. Our chemodynamical model allows us to study the dynamics of stars with different metallicities in great detail, unveiling their distinct 3D density distributions, kinematics and orbital structures. Summarizing the properties of stars in the different metallicity bins in the inner $5\,\rm{kpc}$ of the Galaxy, as found in our model, but also in the light of the work of Ness et al. (2013a, b) and Di Matteo et al. (2014, 2015), we conclude the following:

** A ($0.5\geq\textrm{[Fe/H]}>0$):**

We find that A is the metallicity bin that contributes most to the galactic bar and bulge structure in its current state. A is strongly barred and hosts $52\%\pm 3\%$ of the bar-supporting orbits. It is generally thin, dominates over the other metallicity bins in the plane and provides a very strong support to the B/P shape of the bulge. Its hammer-like shape in the outer parts of the bar indicates that a significant fraction of A belongs to the superthin bar. A has a very low dispersion of only about $\sigma\approx 60-80\,\rm{km\,s^{-1}}$ that steeply rises in the central $5\degree$, caused by the bar-supporting orbits. Stars in A are found to be more rapidly rotating than in the other metallicity bins for longitudes between $10$ and $30\degree$ at $|b|\leq 5\degree$ and show a high-velocity tail in the line-of-sight velocity distribution there. We suggest that this tail, shifting the mean velocity to a higher value, is caused by the complex orbital structure of the bar in the plane. A also contributes in our model $38\%\pm 3\%$ of the inner disk within corotation. In good agreement with Ness et al. (2013b) and Di Matteo et al. (2015), bin A is thus consistent with being partly populated by stars with an early thin disk origin. However, since the superthin bar is probably associated with a more recent star formation episode (Wegg et al., 2015), bin A stars are probably a composite stellar population, with old and metal-rich stars formed in the early disk and later mapped into the bar, together with stars formed more recently in the superthin bar or trapped by the bar.

** B ($0\geq\textrm{[Fe/H]}>-0.5$):**

We find that B is the second main metallicity bin of the galactic bulge and bar structure, which dominates over A at $|b|\geq 5\degree$. B is also strongly barred, contains $34\%\pm 1\%$ of the bar-supporting orbits and also contributes to the B/P shape of the bulge. B is found to be more extended vertically than A, in good agreement with Ness et al. (2013a). It is generally hotter than A and also shows a steep rise in the velocity dispersion towards the centre, a signature of its barred nature. B rotates generally faster than A at $|b|\geq 5\degree$, except in the plane where it does not share the complex orbital structure of A. It accounts for $47\%\pm 2\%$ of the inner disk, a similar share to that of bin A. B thus appears to be consistent with being populated by stars with a disk origin, formed from stars initially located at larger radii than bin A stars, as found by Di Matteo et al. (2014) to explain the faster rotation and larger dispersion.

** C ($-0.5\geq\textrm{[Fe/H]}>-1$):**

C shows significant differences from A and B. It is thicker, and has a slightly slower rotation and a slightly higher dispersion than the more meal-rich stars. C is found to contribute weaker support to the bar, hosting only $12\%\pm 1\%$ of the bar-supporting orbits. Outside the central $\,\rm{kpc}$, C is strongly vertically flattened ($\sim 0.3$) and has a thick disk like density with a fast and cylindrical rotation, while in the inner $\,\rm{kpc}$, it becomes significantly more centrally concentrated. Several metal-poor components are expected to be present in the inner Galaxy, thick disk, stellar halo and possibly a classical bulge which would all be expected to contribute to the metallicity bin C at some level. From both its spatial distribution and kinematics we conclude that outside the central $\,\rm{kpc}$ C is mainly populated by thick disk stars. Whether or not bin C stars were already in a thick disk before bar formation, as suggested by Di Matteo et al. (2015), is however still unclear. Debattista et al. (2017) showed instead that a thin but initially radially hot disk population could also be mapped into a thick distribution after the bar formation. In the central $\,\rm{kpc}$ of the Galaxy, the concentration of metal-poor stars can be interpreted in various ways and more data would be required to investigate this further.

** D ($-1\geq\textrm{[Fe/H]}>-1.5$):**

D is only weakly constrained by the ARGOS data, and even less by the APOGEE data since the fraction of metal-poor stars increases with height above the plane. In our models, D is found to host a very thick, high dispersion and slowly rotating stellar population. It is also centrally concentrated in the inner $\,\rm{kpc}$, similarly to C. However, D is different from a metal-poor tail of C, with a rotation that is not cylindrical and a nearly constant dispersion in the inner $10\degree$ of the Galaxy, possibly indicating the contribution of an extended non-rotating and metal-poor component which we speculate could be the extended stellar halo.

To conclude, most of the metal-rich stars (B and most of A) seem to be consistent with a common and relatively simple origin, the early thin disk, which after bar formation and buckling became mapped into different spatial distributions and a rather complex orbital structure. Part of A could also be associated to a more recent star formation episode, as indicated by the similarities between the superthin bar and the shape of the bar end seen in this metallicity bin. On the contrary, metal-poor stars (C and D stars) appear to presently have a relatively simple structure, only weakly bar-following but they are probably composed of a complex superposition of components: thick disk (either originally thick or thickened by the bar formation), inner halo and possibly a classical bulge. Note however that in our model, the lack of significant vertex deviation for metal-poor stars in Baade’s window is naturally reproduced by the thick disk distribution of bin C stars and therefore does not argue for the presence of a classical bulge in the Milky Way.

Our M2M chemodynamical model of the bar region, although the first of its kind, represents a significant step towards more complete models, that by including additional chemical dimensions will be able to better constrain the formation history of the Galaxy. On-going and future spectroscopic surveys such as APOGEE, GAIA-ESO, GALAH and MOONS will provide stellar abundances of 10-20 elements for more than $10^{5}$ stars each. This profusion of data allows extending the dimensionality of the chemical phase-space in M2M chemodynamical models to other elements that better constrain the origin of stars. The abundance in $\alpha$-elements is of primary importance, since it is a natural chemical “clock” giving the time-scale of stars formation of a given stellar population. We plan in the near future to first extend the modelling of the bar region to a two dimensional chemical space, including [Fe/H] and $\rm{[\alpha/Fe]}$, and later to other elements and stellar ages.

Acknowledgments

We acknowledge the use of the p3lib potential solver, made available to us by Jerry Sellwood and Monica Valluri. We are grateful to Jo Bovy for publicly releasing his tools dealing with APOGEE data, and particularly its selection function. We thank the organizers of the conferences “The Milky Way and its environment” in Paris and GASP 2016 in Canberra for the opportunity to present this work before publication. We thank the anonymous referee for precise reading and helpful comments.