<figure><img src="image_1.png" /></figure>

# The Bulge Metallicity Distribution from the APOGEE Survey

Ana E. García Pérez<sup>1,2,3</sup>, Melissa Ness<sup>4</sup>, Annie C. Robin<sup>5</sup>, Inma Martinez-Valpuesta<sup>2,3</sup>, Jennifer Sobeck<sup>6</sup>, Gail Zasowski<sup>7</sup>, Steven R. Majewski<sup>1</sup>, Jo Bovy<sup>8,9</sup>, Carlos Allende Prieto<sup>2,3</sup>, Katia Cunha<sup>10,11</sup>, Léo Girardi<sup>12,13</sup>, Szabolcs Mészáros<sup>14,19</sup>, David Nidever<sup>15</sup>

<sup>1</sup>Department of Astronomy, University of Virginia, Charlottesville, VA 22904-4325, USA

<sup>2</sup>Instituto de Astrofísica de Canarias, E-38205 La Laguna, Tenerife, Spain

<sup>3</sup>Departamento de Astrofísica, Universidad de La Laguna, E-38206 La Laguna, Tenerife, Spain

<sup>4</sup>Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany

<sup>5</sup>Institut Utinam, CNRS UMR 6213, OSU THETA, Université Bourgogne-Franche-Comté, 41bis avenue de l'Observatoire, F-25000 Besançon, France

<sup>6</sup>Department of Astronomy, University of Washington, Seattle, WA 98195-1580, USA

<sup>7</sup>Department of Physics and Astronomy, University of Utah, Salt Lake City, UT 84112, USA

<sup>8</sup>Department of Astronomy and Astrophysics, University of Toronto, 50 St. George Street, Toronto, ON, M5S 3H4, Canada

<sup>9</sup>Dunlap Institute for Astronomy and Astrophysics, University of Toronto, 50 St. George Street, Toronto, ON, M5S 3H4, Canada

<sup>10</sup>Observatório Nacional, São Cristóvão, Rio de Janeiro, Brazil

<sup>11</sup>Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA

<sup>12</sup>Laboratório Interinstitucional de e-Astronomia—Linea, Rua Gal. José Cristino 77, Rio de Janeiro, RJ—20921-400, Brazil

<sup>13</sup>Osservatorio Astronomico di Padova, INAF, Vicolo dell'Osservatorio 5, I-35122 Padova, Italy

<sup>14</sup>ELTE Eötvös Loránd University, Gothard Astrophysical Observatory, Szombathely, Hungary

<sup>15</sup>National Optical Astronomy Observatories, Tucson, AZ 85719, USA

<sup>16</sup>Astrophysics Research Institute, Liverpool John Moores University, 146 Brownlow Hill, Liverpool, L3 5RF, UK

<sup>17</sup>Laboratoire Lagrange, Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Blvd de l'Observatoire, F-06304 Nice, France

<sup>18</sup>University of Texas at Austin, McDonald Observatory, Fort Davis, TX 79734, USA

Received 2017 September 5; revised 2017 November 17; accepted 2017 November 20; published 2018 January 10

## Abstract

The Apache Point Observatory Galactic Evolution Experiment (APOGEE) provides spectroscopic information of regions of the inner Milky Way, which are inaccessible to optical surveys. We present the first large study of the metallicity distribution of the innermost Galactic regions based on high-quality measurements for 7545 red giant stars within 4.5 kpc of the Galactic center, with the goal to shed light on the structure and origin of the Galactic bulge. Stellar metallicities are found, through multiple Gaussian decompositions, to be distributed in several components, which is indicative of the presence of various stellar populations such as the bar or the thin and the thick disks. Super-solar ([Fe/H] = +0.32) and solar ([Fe/H] = +0.00) metallicity components, tentatively associated with the thin disk and the Galactic bar, respectively, seem to be major contributors near the midplane. A solar-metallicity component extends outwards in the midplane but is not observed in the innermost regions. The central regions (within 3 kpc of the Galactic center) reveal, on the other hand, the presence of a significant metal-poor population ([Fe/H] = −0.46), tentatively associated with the thick disk, which becomes the dominant component far from the midplane (|Z| ≥ +0.75 kpc). Varying contributions from these different components produce a transition region at +0.5 kpc ≤ |Z| ≤ +1.0 kpc, characterized by a significant vertical metallicity gradient.

**Key words:** Galaxy: bulge – Galaxy: structure – stars: abundances – stars: atmospheres

## 1. Introduction

In the standard theoretical framework for galaxy formation and evolution, galaxy formation proceeds by hierarchical merging of cold dark matter clumps and their associated baryons. However, the physics that drives the evolution of baryonic matter, critical for realistically modeling the luminous components of galaxies, remains to be understood. Processes such as star formation and feedback work on scales much smaller than the resolution of current galaxy simulations, which limits the generation of robust predictions (e.g., Agertz et al. 2011).

The Milky Way (MW) bulge is an exemplar of a barred bulge (Dwek et al. 1995) with a low Sérsic index (Widrow et al. 2008) and an X-shaped profile (McWilliam & Zoccali 2010; Nataf et al. 2010). N-body simulations of disk galaxies have demonstrated that bar formation and bar instabilities are important for the evolution of central regions

in spiral galaxies (Combes et al. 1990; Raha et al. 1991; Athanassoula 2005). Bars can form in thin disks and then buckle, which explains observations of rotationally supported bars, peanut shapes, and X-shape profiles in the inner regions of galaxies (Bureau & Athanassoula 2005).

Simulations of MW-like galaxies can form bars and reproduce at least some of the observed properties of the MW bulge (e.g., Guedes et al. 2011; Okamoto 2013). However, the direct attribution of MW bulge properties to bar instabilities and buckling has not yet being established. A vertical metallicity gradient, which has been detected in the MW bulge, was originally thought to be unsustainable after bar buckling due to orbital mixing, but, as discussed by Ness et al. (2013), is indeed possible (see also Martinez-Valpuesta & Gerhard 2013).

The bulk of the MW's bulge stellar population is old (10 Gyr, e.g., Ortolani et al. 1995; Clarkson et al. 2008), but observations of microlensed turnoff dwarfs (Bensby et al. 2013), intermediate mass asymptotic giant branch stars (AGB; Uttenthaler et al. 2007), and planetary nebulae (PNe

<footnote n="19"><sup>19</sup>Premium Postdoctoral Fellow of the Hungarian Academy of Sciences.</footnote>
García-Hernández & Górny 2014) provide evidence of a younger (<5 Gyr) population (see also Gesicki et al. 2014).<fnref n="20" /> A wide range of metallicities is observed in the bulge, with the mode around solar. A kinematical investigation of metal-rich M-type giants by the BRAVA survey (Rich et al. 2007; Howard et al. 2009) revealed that the bulge has cylindrical rotation, leaving little room for a hotter kinematical component (a classical bulge; Shen et al. 2010). However, the bulge also has distinct subpopulations that hint at a complex formation history. These multiple populations create metallicity gradients, which are not reproduced by disk galaxy simulations that ignore the mixing of populations.

Zoccali et al. (2008) observed fields along the minor axis (at $b \le -4°$) and measured vertical metallicity variations of $-0.5 \text{ dex kpc}^{-1}$. The outer vertical gradient was later confirmed by e.g., Gaia-ESO (Gilmore et al. 2012; Rojas-Arriagada et al. 2014) and ARGOS surveys (Freeman et al. 2013; Ness et al. 2013), the study of Johnson et al. (2013), etc. Using photometric data from the Vista Variables in the Via Lactea (VVV) program, Gonzalez et al. (2013) created a map, mainly of the southern bulge, showing a smooth metallicity variation with Galactic longitude and a flattening of the vertical gradient in the inner regions ($|b| \lesssim 5°$). This flattening was first found at high spectral resolution by Rich et al. (2012) in a sample of 44 M-type giants. Observing ~430 stars in the red clump with a modest resolving power ($R = 6500$), Babusiaux et al. (2014) confirmed a flattening in the innermost parts, in fields at ($l = -6°, +6°, +10°, b = 0°$), and ($l = 0°, b = +1°$).

A weaker longitudinal metallicity gradient is present in the inner bulge region, as is clearly seen in the metallicity map of Gonzalez et al. (2013). This behavior with Galactic longitude was confirmed at a higher spectral resolution for $b \sim -3°$ by the GIRAFFE Inner Bulge Survey (GIBS) in Gonzalez et al. (2015).

These metallicity gradients are not a consequence of a single narrow metallicity distribution that shifts in mean metallicity as a function of *b*. Instead, these gradients appear to reflect the varying contribution of different populations (e.g., Hill et al. 2011; Ness et al. 2013; Rojas-Arriagada et al. 2014; Gonzalez et al. 2015). The different scale heights of various metallicity subpopulations are tied to their different kinematics. The most dramatic example of this effect is the X-shaped bulge. Metal-rich stars are preferentially associated with this structure (e.g., Hill et al. 2011; Ness et al. 2013; Rojas-Arriagada et al. 2014; Zasowski et al. 2016), while metal-poor stars are not (e.g., Uttenhaler et al. 2012; Ness et al. 2013). This association may be explained by the way stars are redistributed as a function of their initial birth radius into the bulge (Di Matteo et al. 2014). Some studies assign the metal-poor stars to a spherical component (e.g., Hill et al. 2011; Dékány et al. 2013; Rojas-Arriagada et al. 2014; Zoccali et al. 2017), while others associate them with a disk-like structure (e.g., Ness et al. 2016; Portail et al. 2017).

Many of the observed properties of the bulge described above have been reproduced in recent cosmological simulations of galaxy formation (Inoue & Saitoh 2012; Martig et al. 2012; Okamoto 2013). Simulations are now capable of following the evolution of baryons throughout the history of the universe, and therefore can model long-timescale secular processes, such as the formation of a younger population of

stars in the inner galaxy as the result of gas flows driven by internal dynamical formation processes (Obreja et al. 2013; Ness et al. 2014).

While the recent successes of cosmological models are encouraging, such models rely on simple recipes for handling the sub-grid physics and initial conditions. The improvement of these models can only be accomplished through increasingly detailed observations, which permit the refinement of both the initial conditions and the sub-grid physics.

The MW is an invaluable tool in addressing the complicated problem of correctly simulating spiral galaxies, as it is possible to resolve its constituent stars into separate subpopulations. Quantitative knowledge about vertical and radial metallicity gradients in the MW bulge, particularly at low Galactic latitude, are key. Parameterizing the metallicity gradients in detail across the bulge into the disk, and from the midplane outwards to high latitude, is critical to understanding the bulge's formation history and ultimately being able to produce self-consistent simulations capable of describing the large-scale properties of our Galaxy.

The Apache Point Observatory Galactic Evolution Experiment (APOGEE; Majewski et al. 2017), a program of the Sloan Digital Sky Survey III (SDSS-III; Eisenstein et al. 2011), has produced the most complete chemokinematical database of stars useful for mapping the properties of the inner Galaxy based on high-quality spectra at a resolution of $R = 22,500$. The APOGEE H-band wavelength observations easily penetrate the heavily dust-extincted portions of the MW bulge and disk, and therefore they allow for the study of the metallicity variations not only in the outer bulge but—importantly—in poorly studied low-latitude regions, including the inner bulge and along the Galactic plane. As the APOGEE survey has covered predominantly the northern part of the bulge, it complements Gaia-ESO, BRAVA, ARGOS, and GIBS, which are primarily southern hemisphere surveys. In addition, APOGEE has provided more detailed chemical information, with individual element abundances for around 15 atomic species. This has made the identification of chemically peculiar groups of stars in the Galactic bulge possible (García Pérez et al. 2013; Schiavon et al. 2017).

Kinematical and metallicity 2D map ($l, b$) based on APOGEE DR12 data have been presented in Ness et al. (2016). They show a bulge rotating cylindrically and with small gradients of radial velocities and metallicities in the innermost regions. The present work focuses on the *metallicity distribution function* (MDF) of the inner Galaxy and its 3D variations. We pay careful attention to possible biases in the APOGEE DR12 results, and how they might influence the APOGEE mapping of bulge chemistry. The structure of the paper is as follows: Section 2 provides a description of the observations and summarizes how APOGEE determines stellar metallicities. An assessment of the sample selection effects is given in Section 3, while Section 4 addresses how distances are determined and used to winnow the sample to stars in the central Galaxy. Metallicity maps in Galactic Cartesian coordinates XYZ with origin at the Galactic center and the distribution of individual metallicities are presented in Section 5. The APOGEE results are discussed in terms of bulge structure models in Section 6 and final conclusions are offered in Section 7.

<footnote n="20">²⁰ A young and metal-rich population is also seen in the inner regions of the Andromeda Galaxy (Boyer et al. 2013).</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Distribution of the 7545 stars within 4.5 kpc from the Galactic center included in this study. Fields with $b > 40°$ are not shown for clarity. Each displayed field is 2°–3° in diameter and typically contains ~100 stars. The GIBS fields observed at high spectral resolution (stars) and Babusiaux et al. (2014, tan circles) fields are marked.</figcaption></figure>

## 2. Observations and Metallicities

### 2.1. Observations

The results presented in this paper are for a sample of 7545 giant stars observed in the APOGEE survey between July 2011 and July 2014 that have estimated distances that place them within 4.5 kpc from the Galactic Center (GC; $d_{GC} \le 4.5$ kpc). This choice was motivated by sampling the edges of the long bar (Wegg et al. 2015). The stars are distributed in 83 APOGEE pointings of typically 100 stars each, as illustrated in Figure 1. The individual circular fields, ranging from 1 to 3 degrees in diameter, are identified by their Galactic coordinates $(l, b)$ in degrees, e.g., 00+02. Five fields toward the central Galaxy—the APOGEE field centered on the Galactic Center (000+00), a BRAVA field (000–05), Baade’s Window (001–04), and the Sagittarius fields, SGRC and SGRCM-4—were excluded due to their special target selection criteria. The observations cover both the outer bulge ($|l|$ or $|b| > 4°$) and importantly, the poorly studied low-latitude bulge ($|b| < 4°$). Of the 83 fields, 16 lie in the inner bulge, $|(l, b)| \le 4°$, while 29 outer fields (up to $l = 30°$) have $|b| \le 2°$. Previously such low-latitude regions had very few stars observed with high-resolution spectroscopy: only a few dozen M-type giants with $|b| \le 2°$ and a few hundred G-type and K-type giants. The vast majority of the sample is at distances between 4 and 12 kpc from the Sun, and suffers line-of-sight extinctions between 0.2 and 1 mag.

The APOGEE H-band spectra, acquired with a cryogenically cooled, multi-object spectrograph (Wilson et al. 2012) coupled to the Sloan Foundation 2.5 m telescope (Gunn et al. 2006) at the Apache Point Observatory, in New Mexico and recorded by three HAWAII-2RG detectors, span the wavelength range 1.51–1.69 µm. The instrument has 300 input fibers and, in the standard APOGEE configuration, approximately 230 fibers are assigned to science targets<fnref n="21" /> and 70 are reserved for calibration: 35 targeting hot stars to record the telluric absorption pattern plus 35 sky fibers. To achieve Nyquist sampling at the shortest wavelengths, multiple exposures are taken while dithering the

detector array by half a pixel (see Nidever et al. 2015, for further details).

For the bulge fields, stars were selected from the 2MASS Point Source Catalog (Skrutskie et al. 2006) by color and magnitude, largely adopting $[J - K_s]_0 \ge 0.5$ and $7 \le H \le 11$, although a few fields were considered part of the APOGEE disk sample and slightly different selection criteria, with deeper integrations, were applied for them. The specified color cut was adopted to minimize the contamination by foreground dwarf stars (most prominent at $[J - K_s]_0 < 0.8$), retaining potential low-metallicity giants in the sample. The faint limit of $H < 11$ was set to ensure a minimum signal-to-noise ratio (S/N) of 100 per pixel for any inner Galaxy fields where single, approximately 1 hr visits were used (as opposed to the 3 visits norm for APOGEE fields; for more details, see Zasowski et al. 2013). Only ~687 of the 7545 stars presented here have S/N < 100 and all have S/N ≥ 50. Reddening corrections were estimated by combining near- and mid-IR photometry (2MASS, IRAC, and WISE), using the RJCE method (Majewski et al. 2011) and the Indebetouw et al. (2005) extinction law.

Raw data were processed with APOGEE’s custom data reduction pipeline (Nidever et al. 2015), following a standard procedure: pixel dither combination, spectral extraction, wavelength calibration, sky emission and telluric contamination correction, and (when applicable) visit combination. All the spectra have been publicly released as part of the SDSS Data Release 12 (DR12; Alam et al. 2015).

### 2.2. Metallicity Determination and Sample Selection

The APOGEE Stellar Parameter and Abundance Pipeline (ASPCAP; García Pérez et al. 2016) was employed to determine stellar metallicities ([Fe/H]) simultaneously with the other atmospheric parameters $T_{eff}$, log g, [C/Fe], [N/Fe], and [α/Fe]. ASPCAP relies upon $\chi^2$ minimization to match each star’s entire APOGEE spectrum to a library of precomputed, LSF-convolved (FWHM resolving power $R \equiv \lambda/\delta\lambda \sim 22,500$), and normalized synthetic spectra (Shetrone et al. 2015; Zamora et al. 2015). The microturbulence was tied to the surface gravity value by the relation $\xi_t = 2.478 - 0.325 \log g$, derived from the analysis of a subsample of APOGEE data. The final ASPCAP metallicities were calibrated to well-known values of a sample of globular and open cluster stars (Holtzman et al. 2015). Based on the dispersion around the calibration values, the typical metallicity accuracy is estimated to be about 0.12 dex. However, the precision of the measurements is significantly better, typically about 0.05 dex (Holtzman et al. 2015). This precision is usually enhanced when abundance ratios such as [O/Fe] are considered. In fact, Nidever et al. (2014) found, for red clump stars in the thin disk, a spread in [α/Fe] at any given [Fe/H] between 0.02 and 0.04 dex, for high S/Ns, and Bertrán de Lis et al. (2016) found that stars in open clusters with similar temperatures showed consistent [O/Fe] ratios to within 0.01 dex.

Our sample is dominated by cooler red giants, and consequently, our metallicities are slightly more uncertain than the bulk of the APOGEE sample, around 0.05–0.09 dex. Additional details about the APOGEE DR12 extracted parameters and abundances may be found in Holtzman et al. (2015).

To select stars with reliable ASPCAP parameters, the APOGEE_ASPCAPFLAG bitmask flag was used.<fnref n="22" /> For the

<footnote n="21">²¹ The number of stars per field analyzed here is further reduced by our quality criteria and fiber-to-fiber distance limitations (see Sections 3 and 4).</footnote>

<footnote n="22">²² https://www.sdss3.org/dr12/algorithms/bitmasks.php</footnote>
present study, stars were removed from our sample if any of the following were true: unreliable $T_{\text{eff}}$, log g, or [Fe/H]; large differences between photometric and spectroscopic $T_{\text{eff}}$ estimates; large χ² values; low S/N; and indications of rapid rotation from broadened line profiles. The final sample stars have parameters inside the ranges −0.4 ≲ log g ≲ 4.0, 3600 ≲ $T_{\text{eff}}$ ≲ 5500 K, and −2.7 ≲ [Fe/H] ≲ +0.6, which match those where the DR12 atmospheric parameters are calibrated. Known cluster members, as reflected in the APOGEE_TARGET1 and APOGEE_TARGET2 flags, were removed from the sample. In addition, stars with a radial velocity dispersion (vscatter) larger than 1 km s⁻¹ were excluded, since that is usually an indication of binarity.

### 3. Metallicity Bias

Large stellar samples of giants spanning different regions of the bulge are ideal for exploring this Galactic component. That was the main motivation for including the observations here described in the APOGEE survey. However, such samples can suffer from selection biases associated with target selection and/or limitations in the spectroscopic analysis, which may skew the derived parameters and overall sample statistics. For a typical bulge distance and age (8 kpc and 10 Gyr), the APOGEE database samples only the top of the red and asymptotic giant branches (RGB and AGB, see Figure 2). In the present work, the culled stars comprise ~63% of the 20,707 survey giants (log g ≲ 3.8) in the 83 fields considered. The most common rejection factor after the distance cut, was having poor $T_{\text{eff}}$ estimates, mainly due to the proximity to the cool edge of our model grids ($T_{\text{eff}}$ = 3500 K), which affects 17% of the 20,707 stars. Among these cool stars lost from the sample, we are preferentially missing the most metal-rich stars, which could somewhat distort the high-end of our inferred metallicity distributions.

We make a quantitative estimation of the biases present in the APOGEE bulge sample by using the Chabrier (2001) IMF and integrating it for different Marigo et al. (2008) isochrones to identify which fraction of any given mono-age and mono-metallicity population in the bulge makes the APOGEE cut: $7 \le H \le 11$, ($J - K_s$) ≥ 0.5, 3600 ≲ $T_{\text{eff}}$ ≲ 5500 K, and −0.4 ≲ log g ≲ 4.0. The integral is computed for multiple isochrones with a relevant range in age (5 and 10 Gyr), and metallicity (approximately between −2.0 and +0.5). Variations in extinction ($A_K$ = 0 and 1) and distance (5 and 8 kpc) are also considered. The fraction of stars observed for each case is computed as the fraction between the integral over the part of the isochrone that satisfies the APOGEE cut ($\int_{c} \xi dM$) and the same integral over the entire range of stellar masses for RGB and AGB stars ($\int_{g} \xi dM$).

A few examples of our APOGEE targeting efficiency estimates are shown in Figure 3. The top row of panels shows the integral of the IMF over the window defined by the APOGEE cut. The bottom row of panels in the figure shows the relative fraction of stars that make the cut. More complete coverage of the entire RGB–AGB is achieved at shorter distances, and the same is true for lower values of the interstellar extinction.

As discussed above, metal-poor stars are expected to have a better sampling than metal-rich stars because of the $T_{\text{eff}}$ cut: at the highest metallicities ([Fe/H] ≳ +0.2), the brightest parts of

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Selected sample (black) is compared with theoretical isochrones (Marigo et al. 2008) for an age of 10 Gyr and five different metallicities. The gray regions indicate the selected parameter space, defined by brightness (for a typical bulge distance of 8.0 kpc from the Sun) and color limits, as well as calibration restrictions.</figcaption></figure>

the RGB and AGB become cooler than the specified $T_{\text{eff}}$ limit. For a 10 Gyr old bulge population, the APOGEE sample would consist of 0.9–1 $M_{\odot}$ stars in the RGB and AGB phases. For a 5 Gyr old bulge population, the APOGEE observations would include instead 1.0–1.3 $M_{\odot}$ stars. The APOGEE bulge sample cuts favor low ([Fe/H] ~ −1.3) over high metallicities ([Fe/H] > 0).

Stars in the red giant branch are statistically better represented at younger ages, closer distances, and/or lower extinctions. Overall, the highest metallicities will be under-represented by some fraction, around 30%, relative to the stars at [Fe/H] ~ −1 (up to 90% for the most distant regions with high extinction). Our APOGEE-based RGB–AGB MDFs (Section 5.3) will be affected by these issues, and the derived metallicity distributions will be distorted, but since departures from the truth distributions should be similar across regions at similar distances, relative variations across the bulge are much more robust, and we focus on those in the present analysis.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Integrated IMF over the isochrones (top) in the APOGEE brightness, color, $T_{eff}$, and log g ranges, and its fraction (bottom) in the entire RGB-AGB range. Values are shown as a function of metallicity for different combinations of distance, extinction, and age. The results of the integrations are normalized to a stellar population of one solar mass.</figcaption></figure>

4. Distances

For the definition of the bulge sample, we adopted a solar Galactocentric distance of 8 kpc, and the distance estimates from Hayden et al. (2015) with a limit of 4.5 kpc in Galactocentric distance ($d_{GC}$). This limit is intended to restrict the sample to the inner Galaxy. Hayden et al. (2015) used a Bayesian method that assumed three stellar density priors (bulge, disk, and halo) and relied upon reddening estimates (Zasowski et al. 2013) as well as Padova isochrones (Bressan et al. 2012) to generate stellar distances of ~20% accuracy. Specifically, Hayden et al. computed probability density functions (PDFs) for various combinations of metallicity, mass, and age based on the probability of belonging to the triaxial bulge, disk, and halo, and taking into account the initial mass function. With these distance determinations, ~50% of the stars observed in the 83 fields have $d_{GC} \le 4.5$ kpc and ~7000 are within 1.25 kpc of the plane. Note that some foreground contamination is present in each of the fields.

Figure 4 shows the spatial distribution in Galactic Cartesian coordinates of the stellar sample, while Figure 5 displays the distribution of their distances for different Galactic longitudes and derived heights from the plane ($Z = d \sin b$). Only regions with more than 20 stars are considered in the latter figure. We find that APOGEE-1 samples mostly the near side of the bulge ($d \lesssim 8$ kpc). To study variations with metallicity, stars were placed into three groups: high metallicity ($[\text{Fe/H}] \ge 0.0$), intermediate metallicity ($-0.5 \le [\text{Fe/H}] \le +0.0$), and low metallicity ($[\text{Fe/H}] \le -0.5$). This grouping scheme was informed by the MDFs described in Section 5.3. As shown in Figure 5, in the midplane, the most metal-poor population tends to be more distant than the metal-rich population and covers a wider range of distances at low Galactic longitudes. The intermediate metallicity population has a median in distance distribution between these two. As Galactic longitude increases or proceeds higher in the bulge, the separation in distance between the high and low-metallicity groups seems to grow smaller. Some of this departure at low $l$ may be attributed to the bias toward low-metallicity stars at large distances and higher

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Spatial distribution of the stellar sample in Galactic Cartesian coordinates XYZ.</figcaption></figure>

extinction, as described above. Nevertheless, some of this difference may also be intrinsic to the structure of the inner Galaxy, e.g., metallicity populations of different scale heights (Robin et al. 2012, 2014).

5. Metallicity Maps

The existing APOGEE observations provide larger bulge coverage at high spectral resolution than any other existing data set. Variations in the metallicity distribution across the bulge inform about possible formation scenarios for the bulge and shed light on the link between these regions with other Galactic components/stellar populations (e.g., disk, bar, and halo). We have used the APOGEE metallicities to calculate the median value at various positions along the inner Galaxy. This is done
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Distance distribution (from the Sun) of the bulge sample ($d_{GC} \le 4.5$ kpc) in bins of 0.5 kpc and separated by Galactic longitude (10° bins), heights from the midplane (1 kpc bins), and metallicities. Histograms are normalized to peak values and the different colors represent the three different metallicity groupings.</figcaption></figure>

after binning the sample in Galactocentric (X, Y, Z) with sizes of 0.5 kpc, 0.5 kpc, and 0.25 kpc, respectively. Previous results for the bulge have suggested the presence of multiple metallicity components of different relative contributions (e.g., Ness et al. 2013). Therefore, characterizing these distributions by medians rather than means was adopted because the former have lower sensitivity to outliers and thus give a more robust representation of the contributing metallicities.

Figure 6 shows (X, Y) maps of the different median metallicity as a function of height (Z) from the Galactic plane. Typical errors for the median metallicities were found to be ~0.05–0.10 dex and were estimated from bootstrapping simulations. In Figure 6, the number of stars does vary with position (from a typical value of five beyond the GC or far from the midplane, to a few tens at our side of the bulge).

## 5.1. A Metal-rich Bulge at Low Heights

Our results, illustrated in Figure 6, present a metal-rich bulge at low heights ($|Z| \le 0.50$ kpc). This has been suggested by previous studies, but those were based on 2D maps (l, b), without spatial resolution along the line of sight. Our 3D maps show significant variations in metallicity with position within the bulge. The side of the bulge closest to the Sun (Y < 0) appears to be metal-rich ($[\text{Fe/H}] \sim +0.2$), while the more distant parts (Y > 0) seem to be more metal-poor ($[\text{Fe/H}] \lesssim -0.2$).

Figure 7 collapses the information on that axis to offer a different perspective of the median metallicity as a function of Galactocentric distance ($R_{GC} = \sqrt{X^2 + Y^2}$, with $R_{GC}$ set to negative values at Y < 0). Changing the sign of $R_{GC}$ depending on whether a location is before or beyond the GC. The Galactic center is very useful to separately consider the more distant

regions, which are more prone to systematic effects (see Section 3). The part of the bulge that is closer to the Sun seems to be more homogeneous in metallicity, although, low-metallicity regions are also observed. The far/distant side of the bulge exhibits, in general, lower metallicities, but those low-metallicity regions can be followed, at intermediate Galactic longitudes ($l \sim 15°$), by regions of higher metallicities.

The APOGEE survey was conducted from the northern hemisphere, but did manage to observe some lower latitude regions of the southern Galactic hemisphere, albeit with overall poorer statistics, a situation now being remedied by data acquisition with the southern-hemisphere-based spectrograph of APOGEE-2. In general, data from the northern Galactic latitudes appear fairly similar to those obtained from regions located south of the Galactic plane.

As discussed in previous sections, the observed variation in metallicity with heliocentric distance is suggestive of biases in the stellar sample due to the cool limits of the model atmospheres used in the spectral analysis. The Galactic bar can contribute to the observed variations; however, its effect should have a marked dependence on Galactic longitude, which is not observed, and a symmetry in metallicity with respect to the bar location would be expected, which is not apparent in the data, probably because of our sample selection (See Section 6). The parts of the bulge closer to the Sun are the ones less affected by sample biases, and we will focus on those for the reminder of the paper.

## 5.2. Vertical Gradient

Far from the midplane ($|Z| \ge +0.75$ kpc), Figure 7 shows a more homogeneous bulge dominated by stars with relatively low metallicity ($[\text{Fe/H}] \lesssim -0.5$ on average). Interestingly,
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Metallicity maps for the APOGEE bulge sample at different heights from the midplane and in X- and Y-bins of 0.5 kpc. The Sun is at (X, Y) = (0, −8). The color denotes metallicity and the dotted circles represent R<sub>GC</sub> of 2 kpc and 4.5 kpc. A minimum metallicity of −0.7 is assigned for distinction.</figcaption></figure>

some locations show super-solar metallicity, and a couple quite
low metallicities ([Fe/H] < -1). APOGEE's increased cover-
age of the low-latitude bulge allows us to firmly establish the
presence of a vertical metallicity gradient, consistent with the
findings of Zoccali et al. (2008), Gonzalez et al. (2013), and
Ness et al. (2013). The gradient is no longer evident on the
distant part of the bulge, but, as discussed in the previous
section, our sample lacks metal-rich stars in those regions,
especially closer to the plane where the extinction is stronger.

Figure 8 shows the median metallicity in the region located
in the part of the bulge closer to the Sun (−5 ≤ R<sub>GC</sub> ≤ 0) as a
function of distance from the plane. The slope of the vertical
gradient is not constant, but it appears to be the steepest at

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Median metallicities as a function of Galactocentric distance projected in the midplane (with a typical error bar shown). Negative values of distances were adopted for the side of the Galactic center closest to the Sun to distinguish between the near- and far-sides. The color bar indicates the Galactic longitude.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Median values of the metallicities in Figure 7 at different heights. Only values for regions of $-5 \le R_{GC}(\text{kpc}) \le 0$ are considered.</figcaption></figure>
intermediate heights from the Galactic plane, at $+0.50 \le |Z| \le +1.00$ kpc, with changes in metallicity of approximately $-0.2$ dex in 0.25 kpc. These regions have a sparser coverage in the southern Galactic hemisphere than in the northern one; however, the results for both are fully consistent.

An inner flattening of the metallicity gradient was suggested in earlier studies (e.g., Ramírez et al. 2000; Rich et al. 2012; Babusiaux et al. 2014), which had data for only a few (l, b) locations. We do not only confirm the flattening, but also show the presence of a transition region at intermediate heights, and a flattening beyond $|Z| > 1$ kpc.

### 5.3. Metallicity Distribution Functions

There are only a few studies of the bulge MDF that include low-latitude regions, and these have been restricted to a narrow range in (l, b) (e.g., Rich et al. 2007, 2012; Gonzalez et al. 2015; Zoccali et al. 2017). Babusiaux et al. (2014) had observations in fields at $b = 0°$, but their results were based on optical spectra acquired at a lower spectral resolution ($R = 6500$). The APOGEE database now permits the most complete and accurate study of the distribution of individual metallicities for stars in the Galactic midplane and inner bulge.

The metallicity distributions at different projected Galactocentric radii and heights, in spatial bins of 2 kpc and 0.25 kpc respectively, are shown in Figure 9. Only spatial bins including more than 30 stars are presented. Normalized MDFs are displayed in 0.15 dex metallicity bins, which are twice as large as the typical metallicity uncertainty for stars in the sample.

In the midplane, corresponding to the bottom row of panels in Figure 9, stars from low to super-solar metallicities are observed ($[\text{Fe/H}] \sim -1$ to $+0.5$), but the metal-rich stars are the dominant components. This is particularly true in the regions on the bulge quadrants closer to the Sun, at $R_{GC} < 0$, which we trust to be the ones less affected by sample biases. This metallicity range is very similar to that reported by Gonzalez et al. (2015), but it does reach significantly lower metallicities than in the studies by Rich et al. (2007, 2012), most likely due to a smaller sample size that makes them miss the rare very low-metallicity stars. On the other hand, stars of lower metallicity are major contributors far from the plane ($|Z| > 0.75$). Note the presence of a significant metal-poor contribution in the central regions, as previously seen in APOGEE data for the GC by Schultheis et al. (2015), and reported earlier for a few inner fields by Babusiaux et al. (2014).

### 5.4. MDF Decomposition

The detection of different metallicity distributions in the inner Galaxy can be interpreted in terms of density variations in multiple overlapping metallicity components (e.g., disk, bar, classical bulge, and inner halo), as suggested by Ness et al. (2013), rather than a bulk change in the overall population metallicity. For each of the selected regions in Figure 9, the distribution of the metallicities contains information about such components.

Ness et al. (2013) concluded that a minimum sample of ~500 ARGOS survey stars were required to detect multiple metallicity components. The minimum number of stars per bin must be lower for the APOGEE sample, due to its greater metallicity precision: 0.05–0.09 dex versus 0.13 dex for ARGOS. Lindegren & Feltzing (2013) found that the minimum sample size required for resolving two different chemical

distributions separated by $r$ times the standard deviation (i.e., the measurement uncertainty) could be approximated by the expression $N_{min} \approx \exp(0.6 + 13 r^{-0.8})$. This means that two populations whose metallicities differ by 0.32 dex could be resolved only with a sample of ~500 stars measured with the precision of the ARGOS observations, while only ~60 would be required at the typical precision of the APOGEE bulge sample. In our analysis, we typically have more than 100–200 stars per spatial bin, and, for the most part, we find a smooth variation of the distributions across neighboring regions.

The metallicity distributions often exhibit multiple peaks, and vary with position. A three Gaussian (3G) decomposition of the MDFs based on a maximum likelihood estimator and an analysis of jackknife samples (Bovy et al. 2011) returned components at four different metallicities, marked with vertical lines in Figure 9: $+0.32$ (metal-rich), $+0.00$ (intermediate metallicity), and $-0.46$ and $-0.83$ (metal-poor). It should be noted that the separation of the multiple components is larger than the uncertainties in their fitted mean metallicities ($\sigma \le 0.15$ dex). The four different values are related to the four distinct centers found from the various 3G decompositions, with generally three of these four discriminated at each location.

There is an indication of a metal-poor component at $[\text{Fe/H}] = -1.22$ far from the plane ($|Z| > +1.00$ kpc), which may be connected to the stellar halo (Allende Prieto et al. 2014). The fraction of very metal-poor stars detected in this study is larger than that found in the ARGOS survey: 1.47% versus 0.07%, respectively, for a metallicity of $[\text{Fe/H}] \le -1.5$. Most of the ARGOS metal-poor stars are seen at high Galactic latitudes ($b \ge 6°$). Possible explanations for the greater numbers of the current work include a larger presence of this population at low heights and/or a metallicity bias.

The components peaking at $[\text{Fe/H}] \sim -0.83$ and $-0.46$, more prominent far from the plane, may be related to the thick disk (see, e.g., Lee et al. 2011). For a more appropriate comparison in terms of homogeneity and proximity, compare with disk values in Hayden et al. (2014). The contribution of the component centered at $[\text{Fe/H}] = -0.46$ is significant in the central regions of the low bulge ($0 \le |Z| \le +0.25$ kpc). In fact, the fraction of metal-poor stars ($[\text{Fe/H}] \le -0.3$) along the midplane ranges from 7% at the nearest location to 41% in the GC. Such differences are large in comparison to the noise.

The components at metallicities ~0.2–0.3 resemble the distributions reported for the central parts of the thin disk (see, e.g., *R* ~ 5 kpc in Figure 7 of Hayden et al. 2014). This may be indicative of a bulge with a disk origin. The solar and super-solar metallicity components have a larger contribution at low $|Z|$, becoming the major contributors (especially the most metal-rich) on the side closer to the Sun. Interestingly, the solar-metallicity component extends to 3 kpc in radius (cylindrical coordinates). This component is present at low heights, independent of the set of heliocentric distance estimates employed (estimates other than those of the current work were also investigated); however, it is not visible in the most central regions. More uncertain is its vertical extent, whose contribution extends significantly beyond the intermediate heights depending on the adopted set.

The values of the metallicity components found in this study ($+0.32$, $+0.00$, $-0.46$, and $-0.83$) are consistent with those previously reported in the literature. The first and third most metal-rich components are in good agreement with the
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Metallicity distribution functions (tan histograms) in bins of 0.15 dex arranged by projected Galactocentric distance and distance from the Galactic midplane (including only those regions having samples of *N* > 30). Histograms are normalized to their area. Southern data are shown in brown and their associated number of stars is given after the slash. The three Gaussian decomposition (blue curves) is displayed only for northern regions. The vertical lines show the mean metallicity of the individual Gaussians and the median metallicity of the MDF (dashed).</figcaption></figure>

high-spectral-resolution results from GIBS (+0.26 and -0.31; Gonzalez et al. 2015) and the Gaia-ESO survey (+0.18 and -0.50; Rojas-Arriagada et al. 2014). The agreement with ARGOS (+0.10, -0.28, -0.68, and -1.18; Ness et al. 2013) is slightly worse, possibly due to differences with their metallicity scale. Based on a common stellar sample, their estimates in the super-solar metallicity regime are lower than the APOGEE values, and that explains some of the different components identified. However, in the low-metallicity regime, there are discrepancies between the identified populations that cannot be explained by a metallicity offset. Some differences with the study of microlensed dwarfs (Bensby et al. 2017) are also observed. We note that both Ness et al. (2013) and Bensby et al. (2017) cover different parts of the bulge than this study. Both studies find more peaks than those obtained in our 3G decompositions, despite the smaller stellar sample in Bensby et al. (2017). A metal-rich component has also been detected in the midplane by Babusiaux et al. (2014; +0.20). The good agreement demonstrated with various literature studies as well as the small derived uncertainties in the metallicity

decomposition offer further support for the distributions we identify in the APOGEE data.

## 6. Model Comparisons

The metallicity results of this study are compared with two different models: the *N*-body dynamical simulation of Martinez-Valpuesta & Gerhard (2013; MVG hereafter) and the population synthesis model from the Besançon Galaxy Model (Robin et al. 2012, 2014; BGM hereafter). The latter model relies upon more assumptions regarding the Galactic gravitational potential and directly aims to reproduce the observed properties of the stellar populations.

### 6.1. MVG Simulations

The MVG simulation consists of a boxy bulge that evolved from an exponential disk ($Q = 1.5$, scale-length of 1.29 kpc, and scale-height of 0.225 kpc) embedded in a live dark matter halo, and that suffered from instabilities and bar buckling (see Martinez-Valpuesta & Gerhard 2011). The resulting bar has a half-length of 4.5 kpc and an orientation of 25° between the bar
<figure><img src="image_11.png" /><figcaption>**Figure 10.** MVG mean metallicity distribution at $1.50 \le |z| \le 1.00$ (top), $0.50 \le |z| \le 1.00$ (middle), and $0.00 \le |z| \le 0.50$ (bottom) for $4 \le d \le 12.0$ kpc. Metallicities are color-coded and the dashed lines represent $l = 0°, \pm 5°, \pm 10°$, and $\pm 15°$, and projected radii of 3.5 kpc. The orientation of the bar is indicated by a solid line.</figcaption></figure>

major axis and the Sun-GC axis. Metallicity was added to the simulation by assigning a radial metallicity gradient ([Fe/H] = 0.6 - 0.4 $R_{GC}$); $R_{GC}$ in units of kiloparsecs) to the initial disk, which was chosen to reproduce the vertical metallicity variations observed by Gonzalez et al. (2013). According to the model, low-metallicity stars from the outer disk are mapped inwards to high latitudes producing the vertical gradients. The simulation snapshot we adopted was taken after the system was relaxed, $t \sim 1.9$ Gyr.

The simulation shows a metal-rich inner bulge elongated along the bar and surrounded by a metal-poor disk (see Figure 10) due to the initial setup. Low-metallicity stars come from the outer disk (no thick disk is included). The model does not show the asymmetry in the metallicity distributions we observe between the quadrants closer to the Sun and those beyond the Galactic center, in line with our conclusion that those are the result of a bias in our sample. On the other hand, we would have expected to observe the symmetry around the bar position shown by the simulations, but we do not.

The simulation cannot reproduce the high metallicities observed in the solar neighborhood nor of the inner disk. Furthermore, Hayden et al. (2014) show a quite flat radial gradient for the thin disk near the bulge at $|Z| < 0.25$ kpc. This is in contrast with the larger gradient adopted in MVG. Milder metallicity gradients, such as those observed near the Sun, may reproduce better the high metallicities we observe in the midplane.

## 6.2. BGM Model

This model consists of a mixture of multiple stellar populations: bar, thin and thick disks, and halo. Specific properties are assigned as follows:

1. A thin disk with ages from 0 to 10 Gyr, with an age-metallicity relation from Haywood (2008) in the solar neighborhood, and a radial metallicity gradient of $-0.07 \text{ dex kpc}^{-1}$. Its scale-length has been constrained from a study of 2MASS star counts presented in Robin et al. (2012).

2. A bar with an age of 8 Gyr, an average solar metallicity, and no gradients. The shape of the bar has been determined from 2MASS color-magnitude diagrams (Robin et al. 2012).

3. A thick disk having two epochs of star formation at ages of 10 and 12 Gyr. Its characteristics have been determined in Robin et al. (2014). The mean metallicities are $-0.5$ and $-0.8$, respectively, and no metallicity gradients are assumed.

4. A stellar halo, with an age of 14 Gyr, a mean metallicity of $-1.5$, and no metallicity gradient.

The kinematics for each population are computed mainly as described in Robin et al. (2003) for the thin and thick disks, and for the stellar halo, as given after the updates on the age velocity dispersion relation coming from the fit to RAVE and Gaia TGAS data (Robin et al. 2017). For the bar, the full 3D velocity field is computed using an N-body model from Debattista et al. (2006), scaled to fit BRAVA's data (Gardner et al. 2014; Robin et al. 2014).

This model has been observationally constrained, but in that exercise no APOGEE data were used. In the comparison below, APOGEE data are simulated by applying the selection criteria introduced during the survey targeting process. The number of targets in each field are selected exactly as was done for actual APOGEE observations. Further cuts are applied to remove regions of the log g-versus-$T_{\text{eff}}$ plane compromised by ASPCAP's limitations.

The sample extracted from this model is therefore restricted to $4000 \le T_{\text{eff}} (\text{K}) \le 4500$. Cuts in distances are not applied to avoid introducing uncertainties associated with the observed distance estimates. However, the high $T_{\text{eff}}$ cut provides a natural culling of most of the foreground giants.

Observed and simulated MDFs are compared in three latitude bins in Figure 11. The APOGEE observations are overall well fitted by the simulations. Nonetheless, there are some differences, e.g., the super-solar metallicity contribution is overestimated in the model. The variation in metallicity as a function of Galactic latitude is produced by the different proportions of the populations included in the model, distorted by the APOGEE selection function. In these simulations, the dominant populations are the thin disk and the bar at low latitudes, the thick disk at high latitudes, and a combination of
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Observed (black) and Besançon simulated (red) metallicity distribution functions in bins of 0.25 dex. Only bulge fields with $H \le 11$ and only stars with $4000 \le T_{\text{eff}} (\text{K}) \le 4500$ are considered. The MDFs are arranged by Galactic latitude. Top: $10^\circ \le b \le 17^\circ$, middle: $4^\circ \le b \le 10^\circ$, bottom: $b \le 4^\circ$.</figcaption></figure>

both components in between. There is no need to include a specific bulge component to reproduce the observed distributions.

Our observed metal-rich components and the component at $-0.46$ would be associated with the bar+thin disk and the thick disk, respectively, in the model. The inner Galaxy shows a vertical transition from metal-rich to metal-poor brought by a changeover from a region dominated by a bar+thin disk to a thick-disk one, in line with our data.

Still, the thick disk would have a significant concentration in the central regions and would be the main sampled population at the far side of the bulge. This is caused by our target and field selection, and the shortest scale-length of the thick disk in comparison with that of the thin disk (Bensby et al. 2011; Bovy et al. 2012).

Note that the chemodynamical model of Portail et al. (2017) suggested that stars with metallicities as low as [Fe/H] $\sim -0.5$ are strongly barred. A boxy/peanut-like structure was also assigned to stars of such low metallicities in Ness et al. (2016). Though, both studies are based on a metallicity grouping based on the ARGOS components. The ARGOS and APOGEE surveys are not necessarily on the same metallicity scale (Schultheis et al. 2017); therefore, the question is whether the findings of Ness et al. (2016) and Portail et al. (2017) are robust

enough for the analysis assumptions; e.g., APOGEE versus their assumed ARGOS metallicity grouping.

The interpretation of the solar-metallicity component in terms of the BGM model is more challenging. In the simulation, the bar stops at 3.5 kpc from the GC (the thin and thick disks extend beyond), while the observed solar-metallicity component extends farther. An association of this component with the bar is not straightforward, because the bar and thin disk may not be chemically distinct. However, should the association be confirmed (e.g., using kinematics), our observations would give further support to the existence of a long bar ($\sim 4.5$ kpc), for which additional recent support has been offered (Wegg et al. 2015).

The BGM model has indications of a very metal-poor component (old thick disk and halo) everywhere, but with a limited contribution at low heights and large Galactocentric radii, consistent with our non-detections.

## 7. Conclusions

Spectroscopic observations in the IR of the central Galactic regions contain precious information relevant to the physical processes that participated in the formation of the bulge. That the metal-rich stars there are associated with a pseudo bulge is largely based on observations at intermediate and high Galactic latitudes, rather than at the latitudes typical of the bar. The nature of the metal-poor stars is somewhat more uncertain, with several proposed scenarios (e.g., a classical bulge or a thick disk).

Our study, based on high-quality APOGEE data, is unique in spatial coverage, allowing us to carry out a thorough, in situ, investigation of the connection between the bulge and the bar. Ours is the first large-scale 3D map that combines mean metallicities and MDFs based on spectra delivering $\sigma_{[Fe/H]} \lesssim 0.05-0.09$ dex uncertainties for stars across the inner bulge. The study comprises $\sim 7545$ stars in 83 fields, largely with $|b| \le 4^\circ$, and over longitudes from $l = -5^\circ$ to $l = 32^\circ$.

Stars from low to super-solar metallicities are observed in all regions. At low heights ($<0.5$ kpc) the APOGEE data show an overall super-solar metallicity bulge ($\sim +0.2$), and a metal-poor ($\sim -0.4$) population far from the plane ($|Z| > 1.00$ kpc) with a smooth transition in between. The largest vertical metallicity gradients are observed at intermediate heights, with shallower slopes on both ends. The far side of the bulge appears metal-poor through almost all heights, but after detailed evaluation we conclude that this effect is merely an artifact of the selection and analysis biases.

We make decompositions of the MDFs at different locations within the bulge into multiple Gaussian components, supported by maximum likelihood and jackknife techniques. This analysis suggests the presence of four metallicity components at $+0.32$ (super-solar), $+0.00$ (solar), and $-0.46$ and $-0.83$ (metal-poor), which are of different strength across the bulge. The two metal-rich components are observed at low and intermediate heights, but only one of them (super-solar) is observed in the most central regions. The solar component extends more than 3 kpc in the direction of the Sun and beyond the region where we find the metal-poor components. The metal-poor component at $-0.46$, which is also centrally present at low heights, dominates at greater heights.

A possible interpretation of these components, based on their metallicity and model predictions, is their association with the bar, the thin and thick disks. A comparison with the Besançon
model indicates that the bar+thin disk, and the thick disk, contribute mostly at low- and at high-Z distances, respectively, with a smooth transition in between. Changing contributions of the different populations provide a simple explanation for the flattening of the vertical metallicity gradient in the inner regions. Another possible interpretation (motivated by the MVG model) is that the bar changes the stellar orbits of the low-metallicity stars in height (Z) and radius, introducing chemical gradients far from the midplane. Our main discrepancy with this model is our lack of observed metal-poor regions in the midplane on the near side of the bulge, which may be indicative of an inappropriate model. Models with star formation in situ, which may remedy the problem, are under construction.

The combination of chemistry and kinematics brings an improved characterization of the MW central regions. Further progress will be possible in the near future with an expanded stellar sample from the ongoing APOGEE-2 survey, including observations from the southern hemisphere, which offers a much better view of the central parts of the Galaxy. New data and the associated improved statistics and coverage will be invaluable for disentangling the nature of the complex metallicity variations discussed in this work.

Support for A.E.G.P. was provided by SDSS-III/APOGEE. C.A.P. is grateful for support from MINECO for this research through grant AYA2014-56359-P. Sz.M. has been supported by the Premium Postdoctoral Research Program of the Hungarian Academy of Sciences, and by the Hungarian NKFI Grants K-119517 of the Hungarian National Research, Development and Innovation Office. BGM simulations were executed on computers from the Utinam Institute of the Université de Franche-Comté supported by the Région de Franche-Comté and Institut des Sciences de l'Univers (INSU).

Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. The SDSS-III website is http://www.sdss3.org/. SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, University of Cambridge, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofísica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University.

## ORCID iDs

Ana E. García Pérez https://orcid.org/0000-0003-2184-6198

Annie C. Robin https://orcid.org/0000-0001-8654-9499

Inma Martinez-Valpuesta https://orcid.org/0000-0003-2143-2446

Steven R. Majewski https://orcid.org/0000-0003-2025-3147

Jo Bovy https://orcid.org/0000-0001-6855-442X

Carlos Allende Prieto https://orcid.org/0000-0002-0084-572X

Léo Girardi https://orcid.org/0000-0002-6301-3269

David Nidever https://orcid.org/0000-0002-1793-3689

Mathias Schultheis https://orcid.org/0000-0002-6590-1657

Matthew Shetrone https://orcid.org/0000-0003-0509-2656

Verne V. Smith https://orcid.org/0000-0002-0134-2024

## References

Agertz, O., Teyssier, R., & Moore, B. 2011, MNRAS, 410, 1391

Alam, S., Albareti, F. D., Allende Prieto, C., et al. 2015, ApJS, 219, 12

Allende Prieto, C., Fernández-Alvar, E., Schlesinger, K. J., et al. 2014, A&A, 568, A7

Athanassoula, E. 2005, MNRAS, 358, 1477

Babusiaux, C., Katz, D., Hill, V., et al. 2014, A&A, 563, A15

Bensby, T., Alves-Brito, A., Oey, M. S., Yong, D., & Meléndez, J. 2011, ApJL, 735, L46

Bensby, T., Feltzing, S., Gould, A., et al. 2017, A&A, 605, A89

Bensby, T., Yee, J. C., Feltzing, S., et al. 2013, A&A, 549, A147

Bertran de Lis, S., Allende Prieto, C., Majewski, S. R., et al. 2016, A&A, 590, A74

Bovy, J., Hogg, D. W., & Roweis, S. T. 2011, AnApS, 5, 1657

Bovy, J., Rix, H.-W., Liu, C., et al. 2012, ApJ, 753, 148

Boyer, M. L., Girardi, L., Marigo, P., et al. 2013, ApJ, 774, 83

Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127

Bureau, M., & Athanassoula, E. 2005, ApJ, 626, 159

Chabrier, G. 2001, ApJ, 554, 1274

Clarkson, W., Sahu, K., Anderson, J., et al. 2008, ApJ, 684, 1110

Combes, F., Debbasch, F., Friedli, D., & Pfenniger, D. 1990, A&A, 233, 82

Debattista, V. P., Mayer, L., Carollo, C. M., et al. 2006, ApJ, 645, 209

Dékány, I., Minniti, D., Catelan, M., et al. 2013, ApJL, 776, L19

Di Matteo, P., Haywood, M., Gómez, A., et al. 2014, A&A, 567, A122

Dwek, E., Arendt, R. G., Hauser, M. G., et al. 1995, ApJ, 445, 716

Eisenstein, D. J., Weinberg, D. H., Agol, E., et al. 2011, AJ, 142, 72

Freeman, K., Ness, M., Wylie-de-Boer, E., et al. 2013, MNRAS, 428, 3660

García Pérez, A. E., Allende Prieto, C., Holtzman, J. A., et al. 2016, AJ, 151, 144

García Pérez, A. E., Cunha, K., Shetrone, M., et al. 2013, ApJL, 767, L9

García-Hernández, D. A., & Górny, S. K. 2014, A&A, 567, AA12

Gardner, E., Debattista, V. P., Robin, A. C., Vásquez, S., & Zoccali, M. 2014, MNRAS, 438, 3275

Gesicki, K., Zijlstra, A. A., Hajduk, M., & Szyszka, C. 2014, A&A, 566, A48

Gilmore, G., Randich, S., Asplund, M., et al. 2012, Msngr, 147, 25

Gonzalez, O. A., Rejkuba, M., Zoccali, M., et al. 2009, ApJL, 702, L153

Gonzalez, O. A., Zoccali, M., Vasquez, S., et al. 2015, A&A, 584, A46

Guedes, J., Callegari, S., Madau, P., & Mayer, L. 2011, ApJ, 742, 76

Gunn, J. E., Siegmund, W. A., Mannery, E. J., et al. 2006, AJ, 131, 2332

Hayden, M. R., Bovy, J., Holtzman, J. A., et al. 2015, ApJ, 808, 132

Hayden, M. R., Holtzman, J. A., Bovy, J., et al. 2014, AJ, 147, 116

Haywood, M. 2008, MNRAS, 388, 1175

Hill, V., Lecureur, A., Gómez, A., et al. 2011, A&A, 534, A80

Holtzman, J. A., Shetrone, M., Johnson, J. A., et al. 2015, AJ, 150, 148

Howard, C. D., Rich, R. M., Clarkson, W., et al. 2009, ApJL, 702, L153

Indebetouw, R., Mathis, J. S., Babler, B. L., et al. 2005, ApJ, 619, 931

Inoue, S., & Saitoh, T. R. 2012, MNRAS, 422, 1902

Johnson, C. I., Rich, R. M., Kobayashi, C., et al. 2013, ApJ, 765, 157

Lee, Y. S., Beers, T. C., An, D., et al. 2011, ApJ, 738, 187

Lindegren, L., & Feltzing, S. 2013, A&A, 553, A94

Majewski, S. R., Schiavon, R. P., Frinchaboy, P. M., et al. 2017, AJ, 154, 94

Majewski, S. R., Zasowski, G., & Nidever, D. L. 2011, ApJ, 739, 25

Marigo, P., Girardi, L., Bressan, A., et al. 2008, A&A, 482, 883

Martig, M., Bournaud, F., Croton, D. J., Dekel, A., & Teyssier, R. 2012, ApJ, 756, 26

Martinez-Valpuesta, I., & Gerhard, O. 2011, ApJL, 734, L20

Martinez-Valpuesta, I., & Gerhard, O. 2013, ApJL, 766, L3

McWilliam, A., & Zoccali, M. 2010, ApJ, 724, 1491

Nataf, D. M., Udalski, A., Gould, A., Fouqué, P., & Stanek, K. Z. 2010, ApJL, 721, L28

Ness, M., Debattista, V. P., Bensby, T., et al. 2014, ApJL, 787, L19
Ness, M., Freeman, K., Athanassoula, E., et al. 2013, MNRAS, 430, 836

Ness, M., Zasowski, G., Johnson, J. A., et al. 2016, ApJ, 819, 2

Nidever, D. L., Bovy, J., Bird, J. C., et al. 2014, ApJ, 796, 38

Nidever, D. L., Holtzman, J. A., Allende Prieto, C., et al. 2015, AJ, 150, 173

Obreja, A., Domínguez-Tenreiro, R., Brook, C., et al. 2013, ApJ, 763, 26

Okamoto, T. 2013, MNRAS, 428, 718

Ortolani, S., Renzini, A., Gilmozzi, R., et al. 1995, Natur, 377, 701

Portail, M., Wegg, C., Gerhard, O., & Ness, M. 2017, MNRAS, 470, 1233

Raha, N., Sellwood, J. A., James, R. A., & Kahn, F. D. 1991, Natur, 352, 411

Ramírez, S. V., Stephens, A. W., Frogel, J. A., & DePoy, D. L. 2000, AJ, 120, 833

Rich, R. M., Origlia, L., & Valenti, E. 2012, ApJ, 746, 59

Rich, R. M., Reitzel, D. B., Howard, C. D., & Zhao, H. 2007, ApJL, 658, L29

Robin, A. C., Bienaymé, O., Fernández-Trincado, J. G., & Reylé, C. 2017, A&A, 605, A1

Robin, A. C., Marshall, D. J., Schultheis, M., & Reylé, C. 2012, A&A, 538, A106

Robin, A. C., Reylé, C., Derrière, S., & Picaud, S. 2003, A&A, 409, 523

Robin, A. C., Reylé, C., Fliri, J., et al. 2014, A&A, 569, A13

Rojas-Arriagada, A., Recio-Blanco, A., Hill, V., et al. 2014, A&A, 569, A103

Schiavon, R. P., Zamora, O., Carrera, R., et al. 2017, MNRAS, 465, 501

Schultheis, M., Cunha, K., Zasowski, G., et al. 2015, A&A, 584, A45

Schultheis, M., Rojas-Arriagada, A., García Pérez, A. E., et al. 2017, A&A, 600, A14

Shen, J., Rich, R. M., Kormendy, J., et al. 2010, ApJL, 720, L72

Shetrone, M., Bizyaev, D., Lawler, J. E., et al. 2015, ApJS, 221, 24

Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163

Utenthaler, S., Lebzelter, T., Palmerini, S., et al. 2007, A&A, 471, L41

Utenthaler, S., Schultheis, M., Nataf, D. M., et al. 2012, A&A, 546, A57

Wegg, C., Gerhard, O., & Portail, M. 2015, MNRAS, 450, 4050

Widrow, L. M., Pym, B., & Dubinski, J. 2008, ApJ, 679, 1239

Wilson, J. C., Hearty, F., Skrutskie, M. F., et al. 2012, Proc. SPIE, 8446, 84460

Zamora, O., García-Hernández, D. A., Allende Prieto, C., et al. 2015, AJ, 149, 181

Zasowski, G., Johnson, J. A., Frinchaboy, P. M., et al. 2013, AJ, 146, 81

Zasowski, G., Ness, M. K., García Pérez, A. E., et al. 2016, ApJ, 832, 132

Zoccali, M., Hill, V., Lecureur, A., et al. 2008, A&A, 486, 177

Zoccali, M., Vasquez, S., Gonzalez, O. A., et al. 2017, A&A, 599, A12