The Gaia-ESO Survey: Exploring the complex nature and origins of the Galactic bulge populations

TL;DR

This study uses Gaia-ESO survey data to analyze the chemical, kinematic, and spatial properties of the Galactic bulge, revealing its complex bimodal metallicity distribution and linking different populations to formation scenarios.

Contribution

It provides the first large-area, homogeneous spectroscopic analysis of bulge MDF and alpha-abundances, offering new insights into bulge formation and evolution.

Findings

Bimodal metallicity distribution confirmed across the bulge.

Metal-rich stars show bar-like kinematics and X-shape association.

Metal-poor stars have isotropic kinematics and do not participate in the X-shape.

Abstract

Abridged: We used the fourth internal data release of the Gaia-ESO survey to characterize the bulge chemistry, spatial distribution, kinematics, and to compare it chemically with the thin and thick disks. The sample consist on ~2500 red clump stars in 11 bulge fields ($-10^\circ\leq l\leq+8^\circ$ and $-10^\circ\leq b\leq-4^\circ$), and a set of ~6300 disk stars selected for comparison. The bulge MDF is confirmed to be bimodal across the whole sampled area, with metal-poor stars dominating at high latitudes. The metal-rich stars exhibit bar-like kinematics and display a bimodality in their magnitude distribution, a feature which is tightly associated with the X-shape bulge. They overlap with the metal-rich end of the thin disk sequence in the [Mg/Fe] vs. [Fe/H] plane. Metal-poor bulge stars have a more isotropic hot kinematics and do not participate in the X-shape bulge. With similar Mg-enhancement levels, the position of the metal-poor bulge sequence "knee" is observed at [Fe/H]$_{knee}=-0.37\pm0.09$, being 0.06 dex higher than that of the thick disk. It suggests a higher SFR for the bulge than for the thick disk. Finally, we present a chemical evolution model that suitably fits the whole bulge sequence by assuming a fast ($<1$ Gyr) intense burst of stellar formation at early epochs. We associate metal-rich stars with the B/P bulge formed from the secular evolution of the early thin disk. On the other hand, the metal-poor subpopulation might be the product of an early prompt dissipative collapse dominated by massive stars. Nevertheless, our results do not allow us to firmly rule out the possibility that these stars come from the secular evolution of the early thick disk. This is the first time that an analysis of the bulge MDF and $\alpha$-abundances has been performed in a large area on the basis of a homogeneous, fully spectroscopic analysis of high-resolution, high S/N data.

Figures

Tables

Table 1. Table 1: Characterization of the observed fields. Signal-to-noise ratios are the field average. E(J-K) G11 corresponds to the reddening as computed from the extinction maps of Gonzalez et al. ( 2011 ) in a box of 30 arcmin per side centered in the respective ( l 𝑙 l , b 𝑏 b ) coordinates. E(J-K) fit values are those estimated in Sect.  3 from isochrone fitting. Finally, N G G {}_{\textmd{G}} /N fld fld {}_{\textmd{fld}} provides the ratio between giant ( log ⁡ ( g ) < 3.5 𝑔 3.5 \log(g)<3.5 ) and the total number of stars per field.

Field	$l$	$b$	$\overline{\text{SN}}$	E(J-K)G11	E(J-K)fit	N${}_{\text{G}}$/N${}_{\text{fld}}$
name
p1m4	1.00	-3.97	364	0.26	0.20	359/369
p0m6	0.18	-6.03	254	0.14	0.17	180/204
m1m10	-0.74	-9.45	340	0.03	0.06	131/187
p7m9	6.85	-8.87	244	0.10	0.11	200/221
m10m8	-9.78	-8.09	347	0.03	0.08	234/310
m4m5	-3.72	-5.18	192	0.19	0.18	89/94
m6m6	-6.57	-6.18	168	0.13	0.13	189/206
p0m8	0.03	-8.06	310	0.06	0.07	81/98
p2m9	1.71	-9.22	362	0.08	0.08	81/105
p8m6	7.63	-5.86	250	0.22	0.20	284/302
p6m9	6.01	-9.62	279	0.09	0.08	159/196

Table 2. Table 2: Line-of-sight Galactocentric radial velocities of stars located in the bright and faint peaks of the metal-rich magnitude distribution. Units are in km s − 1 1 {}^{-1}\leavevmode\nobreak\ .

		V${}_{\text{GC}}$ bright	V${}_{\text{GC}}$ faint
p0m6	Mean	$-37.7\pm 17.9$	$12.0\pm 19.6$
	$\sigma$	$91.1\pm 12.6$	$94.0\pm 13.9$
	Number	26	33
p0m8	Mean	$6.0\pm 21.6$	$18.5\pm 18.1$
	$\sigma$	$52.9\pm 15.3$	$51.1\pm 12.8$
	Number	6	8
p2m9	Mean	$-35.7\pm 27.5$	$-18.0\pm 18.7$
	$\sigma$	$61.4\pm 19.4$	$52.8\pm 13.2$
	Number	5	8
m1m10	Mean	$-10.2\pm 14.7$	$-33.6\pm 13.1$
	$\sigma$	$54.5\pm 10.3$	$49.1\pm 9.3$
	Number	14	14

Table 3. Table 3: Parameters of the best GMM model for each field where two components were found. The subscripts “MP” and “MR” refer to the metal-poor and metal-rich Gaussian components, respectively.

Field	${\overline{\text{[Fe/H]}}}_{MP}$	${\sigma }_{MP}$	WeightMP	${\overline{\text{[Fe/H]}}}_{MR}$	${\sigma }_{MR}$	WeightMR
name
m4m5	-0.09	0.36	0.65	0.51	0.13	0.35
p1m4	-0.33	0.31	0.64	0.43	0.15	0.36
m6m6	-0.30	0.27	0.60	0.40	0.20	0.40
p0m6	-0.27	0.31	0.64	0.49	0.15	0.36
p8m6	-0.25	0.35	0.63	0.44	0.15	0.37
m10m8	-0.40	0.29	0.78	0.33	0.17	0.22
p7m9	-0.44	0.25	0.49	0.30	0.23	0.51
m1m10	-0.44	0.26	0.68	0.39	0.19	0.32
p6m10	-0.43	0.26	0.70	0.38	0.20	0.30

Equations

$$d(a.m)$$

$$+\frac{{\left[{\text{[Fe/H]}}_{\ast }-\text{[Fe/H]}(a,m)\right]}^2}{{\sigma }_{\text{[Fe/H]}\ast }^2},$$

$$W(a,m)=P_m{P}_{IMF}\left[e^{-d(a,m)}\right].$$

$$V_{GC}$$

$$+16.5\left[\sin (b)\sin (25+\cos (b)\cos (25)\cos (l-53)\right]).$$

$${\left(\frac{d{\sigma }_{\text{gas}}}{dt}\right)}_{\text{infall}}=A(r)X_{\text{i}}{e}^{-t/T_{\text{inf}}},$$

$$\psi (t)=\nu {\sigma }_{\text{gas}}^k$$

Full text

*Aims. *We used the fourth internal data release of the Gaia-ESO survey to characterize the bulge metallicity distribution function (MDF), magnesium abundance, spatial distribution, and correlation of these properties with kinematics. Moreover, the homogeneous sampling of the different Galactic populations provided by the Gaia-ESO survey allowed us to perform a comparison between the bulge, thin disk, and thick disk sequences in the [Mg/Fe] vs. [Fe/H] plane in order to constrain the extent of their eventual chemical similarities.

*Methods. *We obtained spectroscopic data for $\sim 2500$ red clump stars in 11 bulge fields, sampling the area $-10^{\circ}\leq l\leq+8^{\circ}$ and $-10^{\circ}\leq b\leq-4^{\circ}$ from the fourth internal data release of the Gaia-ESO survey. A sample of $\sim 6300$ disk stars was also selected for comparison. Spectrophotometric distances computed via isochrone fitting allowed us to define a sample of stars likely located in the bulge region.

*Results. *From a Gaussian mixture models (GMM) analysis, the bulge MDF is confirmed to be bimodal across the whole sampled area. The relative ratio between the two modes of the MDF changes as a function of $b$, with metal-poor stars dominating at high latitudes. The metal-rich stars exhibit bar-like kinematics and display a bimodality in their magnitude distribution, a feature which is tightly associated with the X-shape bulge. They overlap with the metal-rich end of the thin disk sequence in the [Mg/Fe] vs. [Fe/H] plane. On the other hand, metal-poor bulge stars have a more isotropic hot kinematics and do not participate in the X-shape bulge. Their Mg enhancement level and general shape in the [Mg/Fe] vs. [Fe/H] plane is comparable to that of the thick disk sequence. The position at which [Mg/Fe] starts to decrease with [Fe/H], called the “knee”, is observed in the metal-poor bulge at $\textmd{[Fe/H]}_{\textmd{knee}}=-0.37\pm 0.09$, being 0.06 dex higher than that of the thick disk. Although this difference is inside the error bars, it suggest a higher star formation rate (SFR) for the bulge than for the thick disk. We estimate an upper limit for this difference of $\Delta\textmd{[Fe/H]}_{knee}=0.24$ dex. Finally, we present a chemical evolution model that suitably fits the whole bulge sequence by assuming a fast ($<1$ Gyr) intense burst of stellar formation that takes place at early epochs.

*Conclusions. *We associate metal-rich stars with the bar boxy/peanut bulge formed as the product of secular evolution of the early thin disk. On the other hand, the metal-poor subpopulation might be the product of an early prompt dissipative collapse dominated by massive stars. Nevertheless, our results do not allow us to firmly rule out the possibility that these stars come from the secular evolution of the early thick disk. This is the first time that an analysis of the bulge MDF and $\alpha$-abundances has been performed in a large area on the basis of a homogeneous, fully spectroscopic analysis of high-resolution, high S/N data.

1 Introduction

The Galactic bulge is the Rosetta stone for our understanding of galaxy formation and evolution. Being a major Galactic component, comprising around a quarter of the Milky Way stellar mass ($M_{\textmd{bulge}}=2.0\pm 0.3\times 10^{10}M_{\odot}$, Valenti et al. 2016), and covering around 500-600 square degrees in the sky, the Galactic bulge provides us with the closest example of this kind of frequent galactic structure. As a predominantly old stellar population (Zoccali et al. 2003; Clarkson et al. 2008)111Clarkson et al. (2011) studied the CMD of proper motion selected bulge stars with HST photometry, and found that only $\leq 3.4$% of the bulge population can be younger than 5 Gyr. However, from the spectroscopic analysis of a sample of lensed bulge dwarfs, Bensby et al. (2013) found that nearly 22% are younger than 5 Gyr., it witnessed the very early formation history of the Milky Way, and its stars contain a detailed record of the past chemodynamical events that shaped its current observable properties. This valuable information can be read from photometric and/or spectroscopic observations of its resolved stellar populations. The star-by-star study of its stellar content, together with the great degree of detail that is possible to achieve with the current large aperture telescopes and multiobject spectroscopy, have turned the bulge into an opportunity to perform near field cosmology in order to test any envisaged scenario of galaxy formation.

Currently, there are two broad scenarios of bulge formation. The first assumes an early prompt formation, whether through a dissipative collapse of a primordial cloud contracting in a free-fall time (Eggen et al. 1962) or through the accretion of substructures, disk clumps, or external building blocks in a $\Lambda\textmd{CDM}$ context (Scannapieco & Tissera 2003; Immeli et al. 2004). The predicted outcome of this process is a classical bulge, a centrally concentrated spheroidal structure, predominantly made up of old stars and dynamically sustained by isotropic random orbital motions. The second scenario conceives the bulge formation as the product of secular internal evolution of the early disk over longer timescales. In this case, dynamical instabilities of the early inner disk lead to the formation of a bar, a structure which subsequently undergoes vertical instabilities, buckling, and redistributing disk angular momentum in the vertical direction. The resulting structure – which has a characteristic boxy peanut (B/P) or, in extreme cases, an X-shaped morphology – is commonly called a pseudobulge.

In the last decade, the study of the Milky Way bulge has experienced a revolution, mainly driven by technical improvements in instrumentation and telescope aperture, allowing the execution of several mid- and large-scale spectroscopic and photometric surveys of the central Galactic region. The complex picture that has emerged from this very active research makes it evident that the Galactic bulge can no longer be considered a simple homogeneous structure.

The Galactic bulge hosts a bar (e.g., de Vaucouleurs 1964; Liszt & Burton 1980; Weiland et al. 1994), currently characterized as a triaxial structure of $\sim 3.5$ kpc in length flaring up into an X-shape structure (Wegg & Gerhard 2013; Ness & Lang 2016). This configuration is predicted as an outcome of secular disk evolution.

On the other hand, the metallicity distribution function (MDF) study by Zoccali et al. (2008) demonstrated the existence of a vertical metallicity gradient along the bulge minor axis in the range $b=[-4:-12]^{\circ}$. This gradient, already suggested by Minniti et al. (1995), was interpreted as the signature of classical bulge formation. Using the same sample, Babusiaux et al. (2010) showed that metal-rich stars present a vertex deviation compatible with bar driven kinematics. Instead, the metal-poor component exhibits isotropic kinematics, as expected for a classical spheroid. The work of Hill et al. (2011) on Baade’s window, revealed that these kinematical signatures can be correlated with a bimodal nature of the MDF, which is also found in other fields (Uttenthaler et al. 2012; Rojas-Arriagada et al. 2014; Gonzalez et al. 2015; Zoccali et al. 2017). The work of Ness et al. (2013a) challenged this picture from their analysis of $\sim 10200$ likely bulge stars from the ARGOS survey. In fact, their MDFs from $l=\pm 15^{\circ}$ strips at $b=-5^{\circ},\,-7.5^{\circ},\,-10^{\circ}$ are trimodal. They related the double red clump feature, a signature of the B/P bulge, only with $\textmd{[Fe/H]}\geq-0.5$ stars, and the vertical metallicity gradient with a change in the relative size of the metallicity components. The determination of the intrinsic shape of the bulge metallicity distribution function is fundamental because its exact multimodal shape can be related with a number of different bulge formation channels.

In this general context, attempts to conciliate morphological, chemical, and kinematical evidence argue for a composite nature of the bulge. Recent research seems to agree on the bar-driven secular origin of the metal-rich bulge. Secular evolution through disk instability is able to reproduce the chemical, morphological, and kinematic properties displayed by bulge stars in this metallicity range. Instead, there is less consensus on the origin of the metal-poor bulge. Its spatial distribution seems to be uncorrelated with the bar position, appearing as an extended, centrally concentrated and possibly spheroidal component. This is supported by the distribution found for other tracers of metal-poor old populations such as RR Lyrae stars (Pietrukowicz et al. 2012, Dékány et al. 2013, Kunder et al. 2016, Gran et al. 2016, but see also Pietrukowicz et al. 2015). On the chemical abundance side, $\alpha$-abundance ratios with respect to iron are systematically enhanced over its whole metallicity range.

Detailed comparisons between bulge and thick disk samples in the [$\alpha$/Fe] vs. [Fe/H] plane provide a direct way to try to understand the origin of the metal-poor bulge. Early attempts in this direction (Zoccali et al. 2006; Lecureur et al. 2007; Fulbright et al. 2007) claimed that the bulge presents higher $\alpha$-enhancements relative to the thick disk. Meléndez et al. (2008) and Alves-Brito et al. (2010) attributed this result to systematic effects arising from the comparison of giant and dwarf samples given their different temperature and gravity regimes. Their homogeneous sample of bulge and local thick disk giants display chemical similarities, with similar trends in the [$\alpha$/Fe] vs. [Fe/H] plane, and presumably a comparable location of the so-called “knee” in the sequences of both populations. Similarities between the bulge and the thick disk have also been suggested using dwarf stars (Bensby et al. 2013, 2014). The study of the detailed chemical abundance patterns from statistically significant homogeneously analyzed samples can shed light on the initial conditions, physical processes, and relative timescales characterizing formation and evolution of the bulge and thick disk populations.

All in all, the puzzle of bulge formation has many pieces, and not all of them are currently in their definitive place. In this paper, we provide new evidence on some of the issues discussed above. To this end, we made use of data coming from the fourth internal data release of the Gaia-ESO survey (iDR4). The Gaia-ESO survey is a large ongoing public spectroscopic survey (300 nights from the end of 2011 to the end of 2016) targeting $\sim 10^{5}$ stars distributed in all the main components of the Milky Way: the halo, bulge, and the disk system (Gilmore et al. 2012). The present study is an extension of our previous work (Rojas-Arriagada et al. 2014), which was based on a subset of the fields studied here and not including the analysis of individual abundances. The structure of the paper is as follows. In Sect. 2 the data are presented, the selection function of the Gaia ESO survey described, and the data processing outlined. In Sect. 3 we present the method and the results obtained for stellar distances and reddening determinations from an isochrone fitting procedure. The bulge metallicity distribution function is presented in Sect. 4, while the trends in the [Mg/Fe] vs. [Fe/H] and correlations with kinematics in Sect. 5. A search for chemical similarities between the bulge and the thick disk is presented in Sect. 6. A comparison with a chemical evolution model is presented in Sect. 7. Finally, the discussion and our conclusions are drawn in Sect. 8.

2 Data

In the present study, we made use of data coming from the fourth internal data release of the Gaia-ESO survey. The Gaia-ESO survey consortium is based on working groups in charge of the different tasks, from target selection and observation to the derivation of the different fundamental parameters and abundances required to achieve the scientific goals of the survey. A general description of the survey can be found in Gilmore et al. (2012), while a description of the data processing flow is briefly outlined below.

We work with a sample of 2320 red clump stars from observations collected up to the iDR4 of the Gaia-ESO survey. They are distributed in 11 pointings toward the bulge region. The positions of the observed fields222We refer to the fields throughout the paper by a name convention using their Galactic longitude $l$ and latitude $b$, and the p/m letter coding the $\pm$ sign, to assemble their names. For example, the field at $l=7$ $b=-9$ is named p7m9. are illustrated in Fig. 1 overplotted on top of an extinction map of the bulge region (using data from the extinction maps of Schlegel et al. 1998). Five of the fields were already observed during the first nine months of the Gaia-ESO survey project, and released in the iDR1. They were analyzed in a previous Gaia-ESO survey publication (Rojas-Arriagada et al. 2014), although the $\alpha$-abundances were not included in iDR1. For comparison purposes, a sample of 228 red giant branch (RGB) and red clump stars in Baade’s window was adopted from Zoccali et al. (2008) and Hill et al. (2011). These comparison stars were reobserved and analyzed in the same way as the rest of the Gaia-ESO survey bulge targets, and added to the main sample making a total of 2548 stars. In addition, a set of spectroscopic benchmark stars were observed to calibrate the computed spectroscopic parameters. Spectra were obtained with the ESO/VLT/FLAMES facility (Pasquini et al. 2000) in the MEDUSA mode of the GIRAFFE multi-object spectrograph. Only the HR21 setup was employed (except for half of the stars in the comparison sample, which were also observed with the HR10 setup), providing a spectral coverage spanning from 8484 to 9001 Å with a resolving power of $R\sim 16200$. The general quality of the obtained spectra is quite good: the average signal-to-noise ratio (S/N) is 290 and no spectrum has less than 80 per resolution element.

2.1 Target selection

The targets were selected with a photometric selection function specifically designed for the bulge portion of the Gaia-ESO survey. It made use of $J$ and $K_{s}$ photometry available from the Vista Variables in the Via Lactea project (VVV; Minniti et al. 2010). This selection is illustrated in Fig. 2. A generic color cut selects stars with $(\textmd{J-K}_{s})_{0}>0.38$ mag, which is imposed on the dereddened photometry in each field according to the values estimated from the reddening map of Gonzalez et al. (2011)333These maps, derived from VVV and 2MASS photometric data, are accessible at http://mill.astro.puc.cl/BEAM/calculator.php.. This cut, defining the left border of the selection box in Fig. 2, is blue enough to allow metal-poor bulge stars to be included in the sample, but has the drawback of including a number of foreground dwarf main-sequence stars. This sample contamination enters in a variable proportion according to the field extinction. The latter because the dwarf disk stars are distributed in the CMD mostly in a vertical band at the blue side of the bulge RC. This blue plume is on average less affected by the reddening than the bulge RC, so that the difference in color between the two features depends on the specific field extinction. In Fig. 2, the dwarf thin disk plume is visible at $J-K_{0}\sim 0.35$ mag, while that corresponding to the disk RC at $J-K_{0}\sim 0.65$ mag. Since RC stars are good standard candles, the RC sequence clumps in magnitude whenever these stars clump spatially. This happens at $J_{0}=13.5$ mag which, in fact, corresponds to the mean apparent magnitude of a RC star located in the Galactic bulge. On the other hand, a generic magnitude cut selects stars with $(12.9<\textmd{J}_{0}<14.1$ mag. This 1.2 mag interval is in general large enough to select stars located in the bulge RC peak of the field luminosity function, accounting for the spatial distance spread of the bar and the change in mean magnitude with longitude because of the bar position angle. In a number of fields where a double RC is observed in the luminosity function, the magnitude cut would not fit the entire magnitude extension of the bar. In these cases, an extension of the magnitude limit was allowed to include up to 30% of the targets in an extra 0.3 mag below the nominal cut.

The above selection function draws the main sample of 2320 RC stars. Instead, the sample of 228 comparison stars have selection functions described in Zoccali et al. (2008) and Hill et al. (2011).

2.2 Radial velocities, stellar parameters, and individual abundances

Radial velocities are measured by the Gaia-ESO survey with a dedicated pipeline by cross-correlation against real and synthetic spectra (Koposov et al., in prep). In our sample, velocity uncertainties are lower than $0.4$ km s{}^{-1}\leavevmode\nobreak\. The determination and compilation of a recommended set of atmospheric parameters and elemental abundances is performed by the Gaia-ESO survey working group 10 (WG10) for all the F-, G-, and K-type stars observed with GIRAFFE. A detailed description of the process will be published in Recio-Blanco et al. (in prep). In short, the individual spectra are analyzed using three independent approaches: Spectroscopy Made Easy (SME; Valenti & Piskunov 1996), FERRE (Allende Prieto et al. 2006), and MATISSE (Recio-Blanco et al. 2006). This is performed in a model-driven way by comparing the observed spectra against synthetic templates, whether interpolated from a dense grid or computed on the fly. In this way T_{\rm eff}\leavevmode\nobreak\, \log{g}\leavevmode\nobreak\, [M/H], and [$\alpha$/Fe] are determined by the three nodes. A set of spectroscopic benchmark stars (Jofré et al. 2015) is analyzed in the same way. For each node, the differences between the calculated and the nominal fundamental parameters are estimated for the set of benchmark stars. Using these values, the node results for a given program star are bias corrected into the astrophysical scale given by the benchmark stars, and then combined in average to produce a unique set of atmospheric parameters while reducing the random errors of individual determinations. The corresponding errors are computed as the node-to-node dispersion in order to properly account for large node-to-node discrepancy in low-quality parametrization. They constitute the recommended set of model-driven, multi-method fundamental parameters by the Gaia-ESO survey consortium.

This set of parameters, the Gaia-ESO survey linelist used to compute the synthetic spectra in the previous step (Heiter et al. 2015), and the MARCS model atmospheres (Gustafsson et al. 2008), are adopted to determine the elemental abundances of $\alpha$- and iron-peak elements (including the iron and magnesium abundances used in this work) using SME and an automated spectral synthesis method (Mikolaitis et al. 2014). The results from the two methods compare well, and only small bias corrections are needed. The final abundances for each element are calculated as the average of the two individual determinations, while errors are taken proportional to the absolute difference between them. Finally, abundances relative to the Sun are derived by adopting the solar composition of Grevesse et al. (2007). They constitute the recommended set of abundances by the Gaia-ESO survey consortium.

It is worth highlighting here that, contrary to the Gaia-ESO survey iDR1 (used in our previous bulge study Rojas-Arriagada et al. 2014), the procedure described above includes three improvements: (1) the use of three codes instead of one to compute the fundamental parameters, thus providing final results with smaller statistical, and hopefully, systematic errors; (2) the availability of elemental abundances which enable us to perform a more detailed analysis than that presented in Rojas-Arriagada et al. (2014); and (3) a more robust calibration of both the stellar parameters and abundances thanks to a larger sample of observed benchmarks.

Although the present sample contains a number of fields already studied from the iDR1, the fundamental parameters and abundances adopted here come from the iDR4, as is true for the rest of the sample.

In Fig. 3, we display the Hertzsprung-Russell (HR) diagram, using the fundamental parameters of bulge stars for which the iron determinations from FeI lines are available (almost all of which also have Mg measurements). We can verify the general good quality of the stellar parametrization because the main HR features, main sequence, turn-off and red clump are clearly distinguishable. It is also apparent that the nature of the Gaia-ESO survey selection function leads, as anticipated, to a sample with some contamination from dwarf main-sequence stars. We select stars with $\log(g)<3.5$ (and $\log(g)>1.5$ to avoid giants for which stellar parametrization could suffer from modeling uncertainties) as our RC bulge sample. The dichotomy between the RC and dwarf stars is explicitly shown in the CMD diagram in Fig. 4. The figure clearly shows that the dwarf contaminants are preferentially located on the blue side of the CMD toward the locus where the blue plume of disk dwarf stars is visible in a general bulge field CMD (cf. Fig. 2). A small number of dwarf stars are visible at (J-Ks)${}_{0}\gtrsim 0.65$. They correspond to a fraction of the stars with $\log(g)$ values that are slightly higher than the cut at $\log(g)=3.5$ dex. Although they could correspond to RC members according to their colors, we adopt their spectroscopic classification. This good general correspondence between stars in the HR and CMD diagrams constitutes a sanity check on the internal consistency of the stellar parametrization.

In the following we made use of the giant/RC sample defined above. It contains mostly RC with a contribution of RGB stars, and is composed of 1987 stars (including stars from the comparison sample).

3 Distance and reddening estimations

We calculated individual line-of-sight spectrophotometric distances and reddenings for the whole sample with available FeI measurements (2273 stars). The adopted procedure made use of the fundamental parameters T_{\rm eff}\leavevmode\nobreak\; \log{g}\leavevmode\nobreak\; [Fe/H] (from FeI lines); and VISTA $J$, $H$, and $K_{s}$ photometry and associated errors to compute simultaneously the most likely line-of-sight distance and reddening by isochrone fitting with a set of PARSEC isochrones444Available at http://stev.oapd.inaf.it/cgi-bin/cmd. The general approach, rather similar to other methods in the literature (e.g., Zwitter et al. 2010; Ruchti et al. 2011; Kordopatis et al. 2011), is outlined below.

We consider a set of isochrones spanning ages from 1 to 13 Gyr in steps of 1 Gyr and metallicities from $-2.2$ to $+0.5$ dex in steps of 0.1 dex. In practice, for a given age and metallicity, each isochrone consists of a sequence of model stars with increasing mass located along a track in the T_{\rm eff}\leavevmode\nobreak\vs. \log{g}\leavevmode\nobreak\plane from the main-sequence to the AGB. Each model star is characterized by theoretical values of the absolute magnitudes $M_{J}$, $M_{H}$, and $M_{K_{s}}$. On the other hand, an observed star is characterized by a vector containing a set of fundamental parameters and observed passband magnitudes $\{T_{\textmd{eff}},\log(g),\textmd{[Fe/H]},J,H,K_{s}\}$, together with their associated errors. Given the three fundamental parameters T_{\rm eff}\leavevmode\nobreak\, \log{g}\leavevmode\nobreak\, and [Fe/H], a star can be placed in the isochrone T_{\rm eff}\leavevmode\nobreak\-\log{g}\leavevmode\nobreak\-[Fe/H] space.

• 2.

We compute the distance from this observed star to the whole set of model stars considering all the isochrones. To this end, we adopt the metric

[TABLE]

where $T_{\textmd{eff}}(a,m)$, $\log(g)(a,m)$, and $\textmd{[Fe/H]}(a,m)$ are the fundamental parameters, depending on the age $a$ and mass $m$, characterizing the isochrone model stars. The quantities with a star subscript stand for the fundamental parameters and errors ($\sigma_{T_{\textmd{eff}}\ *},\ \sigma_{\log(g)\ *}$) of the observed star.

• 3.

Using this metric, we compute weights associated with the match of the observed star with each point of the isochrone collection

[TABLE]

This weight is composed of three factors:

• a.

$P_{m}$ accounts for the evolutionary speed of the model stars along the isochrone. The isochrones are constructed in order to roughly distribute their model stars uniformly along them. This means that a simple unweighted statistic using all the model stars will lead to overweight short evolutionary stages and not long-lived ones. A way to correct for this effect is to include a weight $P_{m}$ proportional to the $\Delta m$ between contiguous model stars in order to assign more weight to the long-lived evolutionary stages where a randomly selected star is more likely to be;

• b.

$P_{\tiny{IMF}}$ accounts for the fact that, given a stellar population, the number of stars per mass interval $dN/dm$ is not uniform. In fact, this distribution is given by the initial mass function (IMF)555In practice, we made use of the PARSEC isochrone quantity int_IMF, which is the cumulative integral of the IMF along the isochrone. In fact, following Girardi et al. (2000), we assume that “the difference between any two values is proportional to the number of stars located in the corresponding mass interval”.;

• c.

The third factor is an exponential weight associated with the distance of the observed star with respect to each model star, given the adopted metric. We can use the weights $W(a,m)$ to compute any kind of weighted statistics.

• 4.

We calculate for a given observed star the likely values of its absolute magnitudes $M_{J}$, $M_{H}$, and $M_{K}$ from the set of isochrones.

• 5.

We compute the line-of-sight reddening by comparing the theoretical color with the observed color $E(J-K)=(J_{obs}-K_{obs})-(M_{J}-M_{K})$.

• 6.

Finally, from these values, by considering the observed photometry $J$, $H$, $K_{s}$ and the estimated reddening, we compute the distance modulus and then line-of-sight distances.

We computed distances and reddening values (field averages are quoted in Col. 6 in Table 1) for the whole bulge sample with available [Fe/H] values. Typical internal errors in distance are about 25-30%. Using the ($l$,$b$) star positions, we also compute the Galactocentric Cartesian coordinates $X_{GC}$, $Y_{GC}$, and $Z_{GC}$, and the cylindrical Galactocentric radial distance $R_{GC}=\sqrt{X_{GC}^{2}+Y_{GC}^{2}}$. The distribution of the latter is shown in Fig. 5, separately for the giant and dwarf portions of our sample. We can see how the stars we found to be foreground dwarf contaminants based on their \log{g}\leavevmode\nobreak\values are in fact located mainly at 7-8 kpc, in the solar neighborhood. On the other hand, the presumed RC bulge stars are found in a narrow distribution with a peak at $\sim 1.5$ kpc. We did not expect this maximum to be at $R_{GC}=0$ kpc given that most of our fields are several degrees apart from the Galactic plane. The shape of the $R_{GC}$ distribution led us to introduce a radial distance cut, defining a working sample of likely bulge stars. To this end, we adopted the criterion $R_{GC}\leq 3.5$ pc. This restriction was applied to the giant sample already defined from their \log{g}\leavevmode\nobreak\values. The resulting working sample is composed of 1583 stars.

4 Metallicity distribution function

We studied the shape of the MDF from our working sample of likely bulge stars, excluding the comparison stars because their different selection function might bias the MDF toward high metallicity. As a first glimpse of the bulge MDF, we split the sample into two groups of fields which are close to or far from the plane. They are a combination of fields located at $b>-7^{\circ}$ and $b<-7^{\circ}$, respectively (the horizontal dashed gray line in Fig. 1). In this way, each half contains a similar number of fields. While it is true that this exercise can blur specific MDF field-to-field variations, it allowed us to increase the number statistics to investigate the general characteristics of the bulge MDF. The two subsamples are displayed in the upper and middle panels of Fig. 6. Two things are immediately apparent. First, the MDFs present a clear bimodal distribution with a narrow metal-rich component peaking at super-solar metallicities and another broader and metal-poor component peaking at $\textmd{[Fe/H]}\approx-0.4/-0.5$ dex (in agreement with Hill et al. 2011 and Gonzalez et al. 2015, but in contrast with the trimodal MDF of Ness et al. 2013a). Second, the relative proportion of stars comprising the two peaks changes with Galactic latitude. In fact, the size of the metal-rich component decreases with respect to the metal-poor one while going far from the Galactic plane. Broadly speaking, our metal-poor and metal-rich MDF components encompass the metallicity ranges ${\rm-1.0\leq[Fe/H]\leq 0.0}$ dex and ${\rm 0.0\leq[Fe/H]\leq 1.0}$ dex. The incidence of stars with ${\rm[Fe/H]<-1.0}$ dex is low ($1.7\%$ of our sample), and given its small number we do not attempt here a detailed analysis of its properties. Accounting for our distance cut to select likely bulge members, these stars might be a combination of halo passing-by stars and the metal-poor tail of the endemic bulge population.

To quantify these facts, as we did in Rojas-Arriagada et al. (2014), we performed a Gaussian mixture models (GMM) decomposition666See Ness et al. (2013a) and Rojas-Arriagada et al. (2016) for a mathematical description of the procedure and its application to the analysis of chemical distributions of stellar populations. on the two MDFs. In both cases, the Akaike information criterion, used for model selection, gave preference to a two-component solution with a high relative probability. Close to the plane, the narrow metal-rich component ($\sigma=0.16$ dex) encompasses 36% of the probability density of the model, while the broader metal-poor component ($\sigma=0.33$ dex) the remaining 64%. On the other hand, far from the plane, the metal-rich ($\sigma=0.35$ dex) and metal-poor ($\sigma=0.29$ dex) components account for 30% and 70% of the relative weights, respectively.

As a qualitative comparison, in the lower panel of Fig. 6 we display the MDF of the sources classified as dwarfs according to their \log{g}\leavevmode\nobreak\values, which are mostly solar neighborhood members (Fig. 5). Their distribution resembles what it is observed in the solar neighborhood, for example by the Geneva-Copenhagen survey (e.g., Casagrande et al. 2011). It is clear that these stars have a MDF with a significantly different shape with respect to the bulge sample. Their MDF has a long tail toward low metallicity (partially due to the contribution of the local thick disk) and a sharp decline toward $\textmd{[Fe/H]}=0.4$ dex. The distribution presents a strong peak at solar metallicity, precisely at the locus where the dip in the bimodality of the bulge MDF is located.

The individual MDFs of the 11 bulge fields analyzed in this work are shown in Fig. 7. Individual GMM decompositions were attempted in each field (parameters of the best GMM fits in Table 3). In agreement with Fig. 6, the preferred GMM model has two components, except in two fields (p0m8 and p2m9) where the lower number of stars prevents the GMM from giving strong statistic assessments. When comparing MDF decompositions between strips that are at a similar latitude (rows in Fig. 7), a decline in the number of metal-rich stars in favor of metal-poor stars with increasing distance from the Galactic plane is visible (as seen also in Zoccali et al. 2008; Ness et al. 2013a). On the other hand, while comparing fields at similar latitude, those located at positive longitudes tend to have a more enhanced metal-rich component. This asymmetry with respect to the minor axis was already characterized in the photometric metallicity map of Gonzalez et al. (2013). As described there, it is just a perspective effect due to the bar position angle; at positive longitudes the line of sight intersects the bar at shorter distance from the plane than at negative longitudes. This means that at positive longitude our lines of sight sample regions with higher dominance of metal-rich stars than at the symmetric fields at negative longitude; consequently, the relative size of the metal-rich peak is higher, as observed in Fig. 7.

4.1 Quantification of metallicity gradients

From the GMM profiles, we first determined the metallicity at which the peaks of the two populations are located in each field (with the exception of p0m8 and p2m9). Then we computed metallicity gradients with $l$ and $b$ independently for the two populations. We also computed mean field metallicity gradients with $l$ and $b$. For both metal-rich and metal-poor populations, we found negligible gradients with $l$ but noticeable variations with $b$ (gradients of $-0.18\ \textmd{dex/kpc}$ and $-0.31\ \textmd{dex/kpc}$, respectively). A gradient of $-0.24\ \textmd{dex/kpc}$ was found for the variation of the mean field metallicity with $b$. These values were computed by assuming all the fields centers projected on a plane at 8 kpc (to be consistent with other studies and to allow comparison). Our results are compatible with the presence of internal vertical gradients in both metallicity populations, with the gradient of the metal-poor fraction being $\sim 60$ percent higher than that displayed by the metal-rich stars. In this sense, the global metallicity gradients, traditionally measured from the mean field metallicity variations with $b$, can be interpreted as the interplay of two effects: the variation of the relative proportion in which both populations contribute to the global field MDF plus the presence of internal gradients in both components. As a reference, if we compute the vertical gradient in similar fashion, but using the results for fields at $b=-4^{\circ},-6^{\circ},-12^{\circ}$ from Zoccali et al. (2008), we find a gradient of $-0.24\ \textmd{dex/kpc}$, in excellent agreement with the value we derived from our fields. Also, the photometric metallicity map of Gonzalez et al. (2013) indicates a vertical gradient of $-0.28\ \textmd{dex/kpc}$, again in agreement with the global gradient reported here.

4.2 Spatial distribution of the subcomponents

In Fig. 8 we display the generalized histograms of the VVV $K_{s}$ (reddening corrected) magnitude distributions of fields where the double RC feature is present according to the density maps of Wegg & Gerhard (2013). The upper and lower panels show the magnitude distributions of metal-rich and metal-poor stars in each field. From the comparison of the two sets of profiles, it is clear that an enhanced bimodality is drawn by the metal-rich stars. The difference in magnitude between the two peaks changes from field to field, being smaller closer to the plane, thus tracing the distance between the near and far arms of the X-shape bulge. On the other hand, metal-poor stars present nearly flat magnitude distributions, with some tendency, especially in the outermost fields, to have a peak at faint magnitudes. This occurs because the volume observed is bigger at greater distances, due to the cone effect.

It has been suggested that an enhanced bimodality for metal-rich stars can arise or be inflated by stellar evolutionary effects (Nataf et al. 2014). The RGB is redder than the RC, but both become bluer with decreasing metallicity. This implies that the relative contamination of the RC sample with RGB members can increase as a function of metallicity given a color cut in the survey selection function. From a PARSEC isochrone of 10 Gyr and $\textmd{[Fe/H]}=-1.5$ dex (so at the metal-poor end of the bulge MDF), the RC lies at $J-K=0.40$ mag, redder than the GES color cut at $J-K=0.38$ mag. Consequently, our sample should be free of this potential bias. On the other hand, the ratio of RC relative to RGB stars is an increasing function of metallicity, meaning that for example a sample with $\textmd{[M/H]}\sim-1.3$ dex should be 1.75 times larger than one at $\textmd{[M/H]}\sim 0.4$ dex to display features with the same statistical significance. In the combined set of stars from the p0m6, p0m8, p2m9, and m1m10 fields, the ratio between stars with metallicity lower and higher than solar is 1.65, which ensures that this bias source might not be relevant in our case. A third potential bias comes from a metallicity dependence of the magnitude and the strength of the red giant branch bump. These factors can conspire to increase the signal of the faint magnitude peak at high metallicity. While it is true than the exact modeling of the impact of this effect is complicated, it should just increase the difference between the peaks, and does not necessarily invalidate the qualitative presence of two peaks in the magnitude distribution.

In line with previous studies in the literature (De Propris et al. 2011; Uttenthaler et al. 2012; Vásquez et al. 2013), we attempt to characterize the stream motions in the X-shape bulge by comparing the line-of-sight radial velocities of stars around the peaks of the metal-rich magnitude distribution. Given the size of our sample, this exercise may suffer from low number statistics, as evidenced by the relative size of the error bars. The results for the four studied fields are in Table 2. With the exception of p0m6, there are no statistically significant differences in velocity for the bright and faint groups of metal-rich stars. These results are in agreement with previous works for p0m8 (De Propris et al. 2011) and m1m10 (Uttenthaler et al. 2012). The structure of the X-shape bulge is complex; it is composed of the superposition of several stable family orbits. Radial velocity measurements on a larger number of fields might help us to unravel the nature and spatial distribution of these orbit streams.

The above analysis reinforces the bimodal nature of the MDF throughout the bulge area sampled by our fields. In the following, we aim to further characterize the MDF metallicity groups by including $\alpha$-abundances and kinematics into the analysis.

5 Bulge trends in the [Mg/Fe] vs. [Fe/H] plane

Beyond the study of the MDF, the availability of elemental abundances from high-resolution spectroscopy provides us with an important tool to understand the bulge nature. In fact, the trends displayed by stars of any stellar population in the [$\alpha$/Fe] vs. [Fe/H] plane encode important information regarding its IMF and the star formation history. This is particularly critical in Galactic bulge studies as it has been used in attempts to associate the bulge with other Galactic components, in particular with the thick disk (e.g., Zoccali et al. 2006; Fulbright et al. 2007; Alves-Brito et al. 2010; Bensby et al. 2013).

The Gaia-ESO survey iDR4 provides abundances for several species. We focus here on the distribution in the [Mg/Fe] vs. [Fe/H] plane. We adopted magnesium because its abundance determination seems to be less affected by errors in stellar parameters and because its spectral lines are more clearly defined in the GIRAFFE HR21 setup domain than those of the other available $\alpha$-elements (Mikolaitis et al. 2014). Moreover, like oxygen, magnesium is expected to be produced exclusively by SN II explosions, while other alphas have more than one nucleosynthesis channel.

In Fig. 9, we display the [Mg/Fe] vs. [Fe/H] distribution of our working sample in the different fields. Here we include in p1m4 the comparison RGB and RC stars discarded while studying the MDF (since here we are interested in the trends and not in the density distribution). The upper right panel of Fig. 9 shows the whole sample, together with a median profile and $1\sigma$ dispersion band calculated over several small bins in metallicity. This fiducial trend, is then overplotted on the individual field distributions in the remaining panels. The different field samples compare well with the fiducial trend; there are no strong deviations throughout the bulge region.

Moreover, Fig. 9, shows that in every field the curve tends to flatten at metallicity lower than $\sim-0.4$ dex. This is an expected feature from the time-delay model, according to which the $\alpha$-enhancement levels start to strongly decline with [Fe/H] after the maximum of the rate of supernovae Ia explosions is reached. This produces a knee in the [$\alpha$/Fe] vs. [Fe/H] trend whose location provides constraints on the formation timescale estimate of the stellar system. In Fig. 10, we present the whole bulge working sample, together with a best fit bilinear model. The model used to fit the data consists of two linear trends sharing a common point, i.e., the knee, and leaves the other parameters free. The fit is performed in the range $-1.5\leq\textmd{[Fe/H]}\leq+0.1$, covering the metallicity range of the metal-poor bulge component. The fit is performed by means of a $\chi^{2}$ minimization, and errors are taken into account by performing 1000 Monte Carlo samplings from the individual errors in [Mg/Fe]. We can see in Fig. 10 that given the size of the sample and the data dispersion in [Mg/Fe], we cannot constrain the knee position better than $\sim 0.1$ dex, with the resulting value being $\textmd{[Fe/H]}_{\textmd{knee}}=-0.37\pm 0.09$ dex.

The median trend of our bulge stars in the [Mg/Fe] vs. [Fe/H] plane compares well with the distribution of inner disk stars in the [$\alpha$/Fe] vs. [Fe/H] plane presented in Hayden et al. (2015) (their Fig. 4, leftmost panels for $3<R_{GC}<5$ kpc). In both cases, the sequence starts from the locus of high-$\alpha$ metal-poor stars and ends in that of low-$\alpha$ metal-rich ones. In the case of disk stars, as seen by APOGEE, a vertical step in the sequence is visible at ${\rm[Fe/H]\sim-0.1}$ dex, which is not evident in our bulge sample. Except for this, a general similarity between the stellar distribution of bulge and disk(s) samples in the $\alpha$-abundance vs. metallicity plane can be suggested. Nevertheless, a more detailed quantitative comparison is not possible here since there is no guarantee that both surveys are in the same abundance scale. A set of common stars to cross-calibrate them is needed and awaited.

Based on the conclusions drawn in Sect. 4 regarding the bimodal nature of the bulge MDF, we split the sample into metal-rich and metal-poor stars. To this end, we adopted the limits $\textmd{[Fe/H]}=+0.15$ and $+0.10$ dex for the fields close to ($b>7^{\circ}$) and far from ($b<7^{\circ}$) the plane, respectively. In Fig. 11, we display the Galactocentric velocity dispersion777Galactocentric velocity conceptually corresponds to the line-of-sight radial velocity that would be observed by an stationary observer at the Sun’s position. It is calculated as

$\displaystyle V_{GC}$ $\displaystyle=V_{HC}+220\sin(l)\cos(b)$

$\displaystyle\hskip 10.0pt+16.5\left[\sin(b)\sin(25+\cos(b)\cos(25)\cos(l-53)\right]).$

trends of the fields color-coded according to their Galactic latitude (for a comparison of line of sight distance distributions with simulations, see Williams et al. 2016). Given that the individual radial velocity uncertainties are small compared with the field dispersions, the error in the velocity dispersion can be taken as $\sigma/\sqrt{2N}$. The bulge metal-poor components appear to be kinematically hot throughout the whole sampled area, with values around $\sigma V_{GC}=100\textmd{km}\ \textmd{s}^{-1}$. Instead, the metal-rich components present velocity dispersions higher close to the plane, and decrease systematically with $b$.

To see these results more in perspective, we display in Fig. 12 the [Mg/Fe] vs. [Fe/H] distributions of fields close to and far from the plane. Each subsample is split into small areas, color-coded according to their velocity dispersion. On average, metal-rich and metal-poor parcels are kinematically homogeneous in inner fields, while for the outer ones the metal-rich end is kinematically colder. It is worth noting that, according to this figure, there is no evidence of kinematic variations with [Mg/Fe] at fixed metallicity.

Figures 11 and 12 show that the metal-poor bulge component seems to be more kinematically homogeneous than the metal-rich one in the surveyed area. We attempt to test in detail the kinematics of the metal-poor stars by using all of them to construct the velocity dispersion profile displayed in Fig. 13. A $1\sigma$ error band is displayed as a shaded area. An interesting trend is clearly visible: the velocity dispersion increases and then decreases symmetrically around $\textmd{[Fe/H]}\sim-0.4$ dex, which – curiously – is roughly the metallicity where the [Mg/Fe] vs. [Fe/H] knee is located. This is illustrated by the dashed gray line and shaded area depicting the knee’s metallicity position and error. This behavior is different from that displayed, over the same metallicity range, by ARGOS data (cf. Ness et al. 2013b, their Fig. 7; ${\rm-0.8\leq[Fe/H]\leq 0.0}$ dex). In this sense, it is not fully clear whether the velocity dispersion of metal-poor bulge stars increases steadily as a function of decreasing metallicity or presents a more complex behavior, such as that suggested by Fig. 13. We expect to be able to tackle this issue with the next internal data release of the Gaia-ESO survey as its larger spatial sampling will allow us to compare trends with enough statistics at different small regions in ($l,b$). It is important to fully characterize this behavior since it might provide an important observational constraint, and an interesting fact to be explained by chemodynamical numerical models of Milky Way formation.

6 Chemical similarities between the thick disk and the bulge

A major part of the Gaia-ESO survey pointings are devoted to characterizing the disk populations. We take advantage of this sample to chemically compare the disk(s) and the bulge on the basis of a large homogeneous sample.

The disk samples of Gaia-ESO survey are observed with both the HR10 and HR21 GIRAFFE setups. A careful fundamental parameter homogenization, based on benchmark stars, ensures compatibility between the parameters and elemental abundances derived from the HR10+HR21 setup combination (disk) and the HR21 alone (bulge). In Fig. 14 we display the comparison of HR10+HR21 and HR21 iron and magnesium abundances derived for a sample of 144 bulge stars (half of the comparison sample presented in Sect. 2). A very good agreement between the two sets of measurements is visible.

From the whole disk sample, we selected stars satisfying $\textmd{S/N}\geq 45$, ${\rm\Delta T_{teff}\leq 150\,K}$, $\Delta\log(g)\leq 0.23$ dex, $\Delta\textmd{[M/H]}\leq 0.20$ dex, $\Delta\textmd{[Fe/H]}\leq 0.1$ dex, and $\Delta\textmd{[Mg/H]}\leq 0.08$ dex. In this way, we defined a clean disk sample composed of 6313 stars. A separation of thin and thick disk stars in the [Mg/Fe] vs. [Fe/H] plane was performed by following the dip in [Mg/Fe] distribution in several narrow metallicity bins. The separated subsamples are shown in Fig. 15.

As a first qualitative comparison between the disks and the bulge in the [Mg/Fe] vs. [Fe/H] plane, we constructed median curves and dispersion bands for the thin and thick disk sequences. We overplotted the resulting profiles on top of the bulge sample distribution in Fig. 16. We can see that bulge and thick disk stars have comparable [Mg/Fe] enhancement levels over the whole metallicity range spanned in common. Nevertheless, a larger dispersion in [Mg/Fe] of bulge stars relative to the thick disk is apparent along the whole metallicity range. Although this can be a real feature that reveals differences in chemical evolution between the two populations, we cannot rule out the possibility that this effect is the result of the lack of spectral information available from the HR21 setup for the bulge compared to the HR10+HR21 available for the thick disk sample. On the other hand, the thin disk sequence runs under the bulge one and matches it at $\textmd{[Fe/H]}>0.1$ dex. In this way, a chemical similarity between the metal-poor bulge and the thick disk, and between the metal-rich bulge and the thin disk are apparent. This has the important implication that if we want to explain the bulge as the product of secular evolution, we have to include both the thin disk and the thick disk to properly account for the chemical properties of the bulge sequence (in line with the recent claim of Di Matteo et al. 2015). Current suggestions (Shen et al. 2010; Martinez-Valpuesta & Gerhard 2013) include just the thin disk, which is not consistent with the chemical evidence presented here.

We attempt to make more detailed assessments of the chemical similarity between the bulge and the thick disk by comparing the metallicity location of the knee in the two sequences. Unlike previous attempts in this direction, our thick disk sample spans a broader extent in Galactocentric radii, with a significant number of stars observed down to 4 kpc. We selected stars with $|Z_{GC}|\leq 3$ kpc to ensure a nearly homogeneous distribution of $Z_{GC}$ along the sampled radial range.

We split the thick disk sample in five radial portions of approximately the same number of stars in order to probe potential radial variations of the knee position. As we did for the bulge sample, we fit a bilinear model to the thick disk sequence in each radial bin. We use stars in the range $-1.0\leq\textmd{[Fe/H]}\leq+0.1$ dex to avoid the undersampled metal-poor end and the region where the thin and thick disk sequence separation is more uncertain (i.e., around solar metallicity). The results for the five radial bins are displayed in panels a-e of Fig. 17. The metallicity at which the knee is located, and the respective error from 1000 Monte Carlo samplings on the individual [Mg/Fe] errors, are quoted in each panel. We can see that, accounting for the error bars, the position of the thick disk knee does not change through the sampled radial region. This is explicitly shown in panel f, where – except for the last distance bin (with lower number statistics) – the different ${\rm[Fe/H]_{knee}}$ measurements are consistent with being flat with respect to $R_{GC}$. A radial decrease in the knee metallicity position with $R_{GC}$ would imply an inside-out formation for the thick disk, which would conflict with the observed absence of a radial metallicity gradient (Mikolaitis et al. 2014). Instead, the constant ${\rm[Fe/H]_{knee}}$ we found here might imply a formation given by a single star burst in an initially well-mixed media.

A similar shape of the thick disk trend in the [$\alpha$/Fe] vs. [Fe/H] plane for all $R_{GC}$ has been also qualitatively suggested by the APOGEE data (Nidever et al. 2014; Hayden et al. 2015). As already mentioned, a quantitative detailed comparison between the GES and APOGEE results is not possible because of the unavailability of a set of common stars for cross-calibrating their abundance scales. The trends of low- and high-$\alpha$ stars displayed in Fig. 17 and those of Hayden et al. (2015) (the middle and lower rows of their Fig. 4) are comparable: the two disk sequences intersect each other at solar metallicity. In the inner distance bins, rather than a single sequence of disk stars, as suggested in Hayden et al. (2015), both sequences are visible in GES data but with a lack of metal-poor thin disk stars. This is expected if those stars constitute a different outer disk population, as has recently been suggested (Haywood et al. 2013; Rojas-Arriagada et al. 2016).

Given the radial constancy of ${\rm[Fe/H]_{knee}}$ of thick disk stars, we attempt to increase the accuracy of its determination by performing a bilinear fit on the whole thick disk sample (mean Galactocentric radius $R_{GC}=7.1$ kpc). We obtained a value of $\textmd{[Fe/H]}_{knee}=-0.43\pm 0.02$ dex, which we can consider as representative of the whole thick disk (panel g). In panel h, we display a fit performed just considering RC thick disk stars. The resulting $\textmd{[Fe/H]}_{knee}=-0.44\pm 0.04$ dex is in agreement with the figure derived from the whole sample. This demonstrates that no systematics are likely to be introduced in our analysis by using results coming from the combination of dwarf and giant stars.

Finally, we compare the metallicity knee position of the thick disk and bulge sequences. A difference of $\Delta\textmd{[Fe/H]}=0.06$ dex is found. This difference is relatively small with respect to the size of the error bars of both determinations (0.02 and 0.09 dex, respectively). Unfortunately, the uncertainty levels of our abundance measurements prevents us from making a strong assessment on the statistical significance of a null difference. However, assuming the plausible scenario of a nonzero difference, it would have an upper limit of $\Delta\textmd{[}Fe/H]_{knee}=0.24\ \textmd{dex}$, considering its 95% confidence interval.

In summary, we found evidence of a constant SFR with Galactocentric distance for the thick disk formation. In addition, a chemical similarity between the bulge and the thick disk is suggested by the data. A fine-tuned compatibility between the detailed properties of the two sequences is beyond the statistical resolution of the present sample. Nevertheless, some caution should be taken when considering the facts exposed here; although similar enhancement levels are found for the two populations, indicating a similar IMF, the bulge exhibits a larger dispersion in [Mg/Fe] around the mean, a result that needs to be confirmed with a more homogeneous data set. And similarly, although the knee metallicity positions of the two sequences are comparable within the errors, a plausible difference as large as 0.24 dex suggests a difference in the characteristic SFR of the two populations, i.e., the bulge formed on a shorter timescale than the thick disk.

7 Comparison with a chemical evolution model

The modeling of observational data by means of chemical evolution models provides an interesting opportunity to put constraints on the formation timescale of a stellar system. We attempt here to constrain the bulge formation timescale by adopting a model for a bulge formed at early epochs from the dissipative collapse of a cloud accompanied by a strong burst of star formation. To this end, we adopted the model of Grieco et al. (2012). In this work, two bursts of star formation are invoked to model the metal-rich and metal-poor modes of the bulge MDF. We adopt here the model corresponding to the metal-poor bulge. The model assumes a gas infall law given by

[TABLE]

where $X_{\textmd{i}}$ is the abundance of a generic chemical element $i$ in the infall gas, whose chemical composition is assumed to be primordial or slightly enhanced from the halo formation; $T_{\textmd{inf}}$ is the infall timescale, fixed by reproducing present day abundances (MDF), SFR, and stellar mass; and $A(r)$ is a parameter fixed by reproducing the current average total bulge surface mass density. The parametrization of the star formation rate is adopted as a Schmidt-Kennicutt law:

[TABLE]

with $k$ the law index and $\nu$ the star formation efficiency (i.e., the star formation rate per unit mass of gas). The model includes a Salpeter IMF, constant in space and time, which allows the MDF of the metal-poor bulge population to be correctly reproduced. The set of yields are adopted from Romano et al. (2010).

We ran several models, adjusting the parameters to better reproduce the data. Our best model is displayed in Fig. 18, where it is compared to the whole bulge working sample. This model assumes a short timescale for the gas infall $T_{\textmd{inf}}=0.1\ \textmd{Gyr}$ and a very efficient star formation, with $k=1$ and $\nu=25\ \textmd{Gyr}^{-1}$.

The main characteristics of the bulge sequence (enhancement levels, qualitative location of the knee) are well reproduced by this model. We can see that, for $\textmd{[Fe/H]}\geq-1.5$ dex, the predicted [Mg/Fe] abundance ratio steadily decreases with metallicity. This behavior increases from $\textmd{[Fe/H]}\gtrsim-0.4$ dex, which is comparable with the observed locus of the knee, as determined in Sect. 5. The overall formation timescale, as read from the upper axis of Fig. 18, indicates a rather rapid chemical enrichment of the bulge taking place on a timescale of 0.5-0.7 Gyr. Such a short timescale is compatible with a monolithic assembly of the metal-poor bulge. It is worth noting that the model formally shows a continuous enrichment of the gas up to super-solar metallicity. However, the number fraction of such super-solar stars is low compared to populations with lower metallicities. Thus, from the point of view of the chemical evolution model, it is not excluded that the bulge contains a fraction of metal-rich old stars born in situ. In this sense, metal-rich bulge stars might constitute a composite population of stars formed in situ plus a larger fraction of stars with disk origin currently located in this region as the product of secular bar-driven dynamics. This might be compatible with the relatively flat or even bimodal age distribution of metal-rich bulge stars, as suggested from the analysis of microlensed dwarfs (Bensby et al. 2013) and APOGEE giants (Schultheis et al. 2017), respectively.

8 Discussion

In this work, we have made use of the fourth internal data release of the Gaia-ESO Survey to perform a fully spectroscopic analysis of the bulge in the perspective of other Galactic components. The evidence presented here leads us to consider the Galactic bulge as a composite structure, due to the coexistence of two main stellar populations.

From a Gaussian mixture models analysis, the bulge MDF appears as a bimodal distribution comprising a narrow super-solar metal-rich component and a broad metal-poor component. This bimodal nature is verified in all the individual fields, except in those limited by small number statistics. The relative proportion of stars belonging to each of the two populations changes, with metal-poor stars dominating far from the Galactic plane. The line-of-sight Galactocentric velocity dispersion correlates with metallicity, further stressing the likely different nature of the two populations. Metal-poor stars display a high-velocity dispersion around 100 km s{}^{-1}\leavevmode\nobreak\and nearly independent of ($l$, $b$). Instead, metal-rich stars present a more complex behavior; the stars close to the plane are as kinematically hot as the bulge and decrease systematically with $b$ toward disk values. An additional correlation with metallicity appears when considering the bimodal nature of the bulge RC magnitude distribution. In the fields where this feature is visible, the distinction between the magnitude peaks is enhanced if just metal-rich stars are considered, while metal-poor ones display flatter distributions.

This bimodality contrasts with the trimodal ARGOS MDFs found by Ness et al. (2013a). Their distributions are derived from a larger sample of stars ($\sim 10200$), grouped together in three $l=\pm 15^{\circ}$ latitude strips. If the metal-rich and metal-poor bulge populations have some intrinsic gradients with spatial location, like those reported here, a MDF assembled from fields spanning a large region might have components slightly smeared out in metallicity, thus displaying a more complex MDF. Upon completion of the next GES releases, including a larger number of observed fields, we will be able to test this possibility further. In addition, the ARGOS analysis (based on data of lower resolution and mean S/N than our GES data) is not fully spectroscopic since it made use of photometric constraints to estimate fundamental parameters and metallicity. A number of recent studies examining specific locations in the bulge region (Uttenthaler et al. 2012; Gonzalez et al. 2015; Schultheis et al. 2017) find a bimodal MDF. Instead, ARGOS MDFs are trimodal even in individual fields, where the sample size is comparable to that of the above-mentioned studies. In this sense, as discussed in Schultheis et al. (2017), it might be a possibility that the trimodal nature of the ARGOS MDFs is a feature that arises from their parametrization and the spatial distribution of their fields rather than as an effect of larger number statistics.

Two different origins can be proposed for the stars belonging to the metal-rich and metal-poor MDF components. Metal-rich RC stars participate in the B/P bulge, present bar-like kinematics, and are chemically comparable to metal-rich thin disk stars. We associate them, in agreement with the literature, with a population formed by the classical mechanism of secular evolution of the disk via bar formation and buckling into an X-shape structure (Combes & Sanders 1981; Raha et al. 1991; Athanassoula 2005; Martinez-Valpuesta et al. 2006). An internal vertical metallicity gradient, like the one reported here, is predicted by N-body simulations of secular bulge formation as an effect of the mapping in the vertical direction of horizontal (Martinez-Valpuesta & Gerhard 2013) or vertical (Bekki & Tsujimoto 2011) metallicity gradients initially present in the disk, or a combination of the two (Di Matteo et al. 2015). Thin disk stars in the bulge region can explain the presence of young stars (Bensby et al. 2013), and also the existence of the thin star forming inner disk in the bulge traced by classical Cepheids identified from VVV photometry (Dékány et al. 2015).

On the other hand, metal-poor RC stars do not participate in the B/P bulge, dominate in number density far from the plane, and display isotropic kinematics. These stars might be associated with a classical spheroid component formed at early times from the dissipative collapse of a primordial cloud accompanied by a strong burst of stellar formation. A radial internal metallicity gradient, like the one reported here, and a high-velocity dispersion are expected features from a dissipative collapse. The enhanced levels of $\alpha$-elements and the metallicity at which the knee of the sequence in the [Mg/Fe] vs. [Fe/H] plane takes place are both interpreted here with a chemical evolution model as signatures of a fast chemical evolution ($t\leq 1$ Gyr), dominated by massive stars and characterized by a high star formation efficiency. Shen et al. (2010) suitably reproduced the bulge rotation and dispersion curves from the BRAVA project with an N-body simulation which limits the mass contribution of a possible dissipative collapse-made bulge to be less than 8%. This scenario is incompatible with our results. The mean proportion of stars belonging to each of the MDF components is weighted toward metal-poor stars. Although our data set may not be fully adequate to make strong assessments on the mass contribution of both populations to the global bulge mass budget (the sample is not large enough and there is not enough spatial coverage, especially in the inner bulge), the data clearly hints at a bulge composed of a similar fraction of metal-poor and metal-rich stars. Moreover, from the observed mass-metallicity relation for galactic spheroids (Gallazzi et al. 2005), a system with mean metallicity around ${\rm[Fe/H]=-0.4\,dex}$, as our metal-poor component is, should have a mass of $\sim 10^{10}\,M_{\odot}$, which is comparable to the total bulge mass as estimated from observations (Valenti et al. 2016). Furthermore, RR Lyrae stars, as tracers of metal-poor old stellar populations, have been shown to display an axisymmetric spatial distribution, uncorrelated with the bar position angle (Dékány et al. 2013; Gran et al. 2016), and a high-velocity dispersion of around 130 km s{}^{-1}\leavevmode\nobreak\(Gratton 1987; Kunder et al. 2016). However, this is still under debate; a recent analysis of OGLE data suggests that the RR Lyrae distribution might be elongated with a pivot angle comparable to that of the main bar (Pietrukowicz et al. 2015). The fact that metal-poor bulge stars present cylindrical rotation has been taken as an argument for their secular origin (Ness et al. 2013b). However, recent N-body simulations have shown that an initially nonrotating classical bulge can spin-up into a bar-like structure by absorbing a significant fraction of the disk angular momenta emitted by the bar during its secular evolution (Saha et al. 2012). In this way, the classical bulge will become photometrically and kinematically indistinguishable from the B/P bulge. In the case of a massive initial classical bulge, its central parts might be less affected (Saha et al. 2016), providing a kinematic relic to be exploited by spectroscopic observations of the inner-bulge.

Despite the above discussion, the comparable $\alpha$-enhancement levels, along with the similar position of the knee of the bulge and thick disk sequences, argues for a possible common origin, or at least for a similar chemodynamical evolution of these populations. Moreover, it has been shown that it is possible to reproduce general kinematic and chemical bulge patterns from N-body models explaining its formation as the product of thin+thick disk evolution (Di Matteo et al. 2015). For the knee, we found a small difference of 0.06 dex and an upper limit of 0.24 dex. In the same vein, although similar $\alpha$-enhancement levels are found for the bulge and the thick disk, there is an indication of a larger dispersion in the bulge than in the thick disk. Moreover, enhancement differences in $r$- and $s$-process elemental abundances have been proposed as evidence of different formation timescales between the two structures (Johnson et al. 2012; Van der Swaelmen et al. 2016). Unfortunately, current efforts in this direction are based on the comparison of bulge giants and local dwarf samples, which might suffer from systematics given the different temperature and gravity regimes of the samples. All things considered, the origin and nature of the metal-poor bulge remains to be firmly defined on the basis of larger data sets and further detailed modeling.

The composite nature of galaxy bulges has been pointed out by several authors, from a theoretical point of view (Samland & Gerhard 2003; Athanassoula 2005; Obreja et al. 2013; Fiacconi et al. 2015) and from observational evidence in external galaxies (Gadotti 2009; Nowak et al. 2010; Williams et al. 2011; Erwin et al. 2015; Fisher & Drory 2016). Rather than an eccentricity of nature, the presence of composite bulges, with two or more structure types (disk, pseudobulge, classical bulge) coexisting in the same galaxy, appears as a common outcome of galaxy formation and evolution.

If we assume that the differences in the knee position and [Mg/Fe] dispersion between the bulge and thick disk sequences are real, we can draw a general picture of the bulge formation by interpreting the observational evidence in terms of two different formation episodes. On the one hand, the old bulge population that formed in situ is the product of a fast dissipative collapse in the early epochs of Milky Way evolution. As characterized by a strong SFR, the chemical enrichment of the gas may have reached super-solar metallicities before gas exhaustion, with the majority of stars produced around ${\rm[Fe/H]=-0.5\,dex}$. On the other hand, metal-rich stars in the X-shaped bulge are, as pointed out in the literature, the product of the secular evolution of the early inner disk. In this sense, the main epoch of chemical enrichment in the inner Galaxy occurred early, before the formation of the B/P bulge. A small fraction of the metal-rich stars, endemic to the central regions of the Galaxy, might be old, being currently outnumbered by stars with their origin in the early disk.

A semantic issue is then raised. In fact, we legitimately call “bulge” all stellar populations currently present in the central kiloparsecs of the Milky Way regardless of their origin and specific evolutionary histories.

The Gaia-ESO survey multi-method, model-driven, fully spectroscopic analysis of high-resolution high-S/N data provides a homogeneous self-consistent account of the main Galactic components. Using this exquisite data set, this is the first time that the bulge MDF has been characterized in a large spatial area, in individual fields, from a fully spectroscopic analysis. This is also the first attempt to compare the metallicity position of the bulge and thick disk knee based on a statistically significant sample of homogeneously analyzed stars. All in all, a composite picture of the Galactic bulge can be unambiguously established, with all the presented evidence pointing to the presence of two main components currently coexisting in the central regions of the Milky Way.

Acknowledgements.

This work was partly supported by the European Union FP7 programme through ERC grant number 320360 and by the Leverhulme Trust through grant RPG-2012-541. We acknowledge the support from INAF and the Ministero dell’Istruzione, dell’Università e della Ricerca (MIUR) in the form of the grant “Premiale VLT 2012”. The results presented here benefited from discussions held during the Gaia-ESO workshops and conferences supported by the ESF (European Science Foundation) through the GREAT Research Network Programme. A. Recio-Blanco, P. de Laverny, and V. Hill acknowledge the Programme National de Cosmologie et Galaxies (PNCG) of CNRS/INSU, France, for financial support. A. Recio-Blanco, P. de Laverny, and V. Hill acknowledge financial support form the ANR 14-CE33-014-01. M. Zoccali gratefully acknowledges support from the Ministry of Economy, Development, and Tourism’s Millenium Science Initiative through grant IC120009, awarded to the Millenium Institute of Astrophysics (MAS), by Fondecyt Regular 1150345 and by the BASAL-CATA Center for Astrophysics and Associated Technologies PFB-06.

Appendix A GMM output parameters

Individual GMM decompositions were performed on the 11 bulge fields studied in this paper. All except two MDFs were found to be better explained by a model consisting of two Gaussians, according to the Akaike Information Criterion. In Table. 3, we present the set of parameters characterizing these optimal models for comparison with future studies.