{'fraction_non_alphanumeric': 0.04385998222816304, 'fraction_numerical': 0.026118524481844624, 'mean_word_length': 4.533980582524272, 'pattern_counts': {'":': 0, '<': 54, '<?xml version=': 0, '>': 28, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Chemodynamical history of the Galactic Bulge

The Galactic Bulge can uniquely be studied from large samples of individual stars, and is therefore of prime importance for understanding the stellar population structure of bulges in general. Here the observational evidence on the kinematics, chemical composition, and ages of Bulge stellar populations based on photometric and spectroscopic data is reviewed. The bulk of Bulge stars are old and span a metallicity range -1.5<~[Fe/H]<~+0.5. Stellar populations and chemical properties suggest a star formation timescale below ~2 Gyr. The overall Bulge is barred and follows cylindrical rotation, and the more metal-rich stars trace a Box/Peanut (B/P) structure. Dynamical models demonstrate the different spatial and orbital distributions of metal-rich and metal-poor stars. We discuss current Bulge formation scenarios based on dynamical, chemical, chemodynamical and cosmological models. Despite impressive progress we do not yet have a successful fully self-consistent chemodynamical Bulge model in the cosmological framework, and we will also need more extensive chrono-chemical-kinematic 3D map of stars to better constrain such models.


INTRODUCTION
The Galactic Bulge was first identified as a distinct component of our Galaxy by Baade (1946) and Stebbins & Whitford (1947). In the same work, Baade also used for the first time the term stellar population of type II for the stars in the globular cluster NGC 6522 and its surrounding fields. The discovery of M giants by Nassau & Blanco (1958) first showed the presence of metal-rich old populations. Early historical work on the Galactic or Milky Way (MW) Bulge, hereafter the Bulge, is reviewed in Madore (2016).
Thorough studies of Bulge populations were carried out by Whitford & Rich (1983) and Rich (1988) who presented a metallicity distribution function (MDF) in the range −1<[Fe/H]<+1 1 . A Color-Magnitude Diagram (CMD) for Bulge field stars was obtained by Terndrup (1988), in parallel to the infrared survey by Frogel & Whitford (1987), the radial velocity study by Mould (1983), a velocity dispersion derivation from integrated spectra by Freeman et al. (1988), and the multi-fiber TiO band spectroscopic survey of 300 giants in Baade's Window (BW) by Sharples et al. (1990). McWilliam & Rich (1994) presented the first high-resolution abundance analysis of 11 Bulge giants that revealed their [α/Fe]-enhancement.
Photometric and spectroscopic surveys covering large regions of the Bulge together with advanced data interpretation tools have since led to a new challenge, of combining the multifaceted observations into a consistent chemodynamical evolution picture for our Bulge, and placing the Bulge into the more general context of bulges of spiral galaxies. Many issues are still open and debated. One of the critical issues concerns the age of the Bulge. Ages derived from CMDs suggest an uniform old age for the Bulge stellar populations (e.g. Ortolani et al. 1995;Clarkson et al. 2008;Gennaro et al. 2015), with little room for intermediate age stars, in contrast to findings by Bensby et al. (2017), based on spectroscopy, which suggest a fraction of at least 15% of metal-rich ([Fe/H] > ∼ −0.5) dwarf stars younger than 5-8 Gyr. A broad metallicity range of −1.5 < ∼ [Fe/H] < ∼ +0.5 has been confirmed with the MDF turning out to be dependent on the field location, which suggests the coexistence of different Bulge populations. Most Bulge stars with [Fe/H] < ∼ −0.4 show enhancements in α-elements, indicative of a fast chemical evolution, dominated by core collapse supernovae (CCSNe) chemical enrichment.
De Vaucouleurs (1964) first suggested a barred Bulge in the MW from the high gas velocities in the central Galaxy. The existence of the bar was put on firm grounds in the 1990's, from modeling the non-circular gas motions in the Bulge region (Binney et al. 1991;Englmaier & Gerhard 1999;Fux 1999),the near-infrared data in particular from the COBE satellite (Blitz & Spergel 1991;Weiland et al. 1994;Binney et al. 1997), and red clump giant star counts (RCG) (Stanek et al. 1994). Modern data have confirmed the Bulge to have a box/peanut (B/P) shape typical for barred galaxies. With star counts at high latitudes it was discovered that the Bulge has an X-shape component (McWilliam & Zoccali 2010;Nataf et al. 2010), which is a manifestation of a B/P bulge as later demonstrated by the reconstruction of the 3D Bulge density (Wegg & Gerhard 2013) from VVV data (see Table  1). Because B/P bulges are the inner 3D parts of a longer, flatter bar in the disk (e.g., Bureau & Freeman 1999;Athanassoula 2005;Erwin & Debattista 2013), the definition of 'the Bulge' has become somewhat blurred, affecting also estimates of the Bulge mass (see Section 2). General practice is to use a region on the sky with size ±10 • in both longitude and latitude (l, b), which however poses the problem of separating Bulge stars from foreground and background disk stars.
Bulges of spiral galaxies can be classified into two categories, according to their observational properties, often designated as classical bulges (ClBs) and pseudobulges, according to Kormendy & Kennicutt (2004); Athanassoula (2005) ;Fisher & Drory (2016): classical bulges have properties similar to elliptical galaxies, with steep surface brightness profiles with Sersic index n>2, pressure supported kinematics, generally rounder isophotes, steeper dispersion profiles, spheroidal rotation, anisotropy parameter (Vrot/σ) ≤ 1 2 (Kormendy & Illingworth 1982), bulge to total ratios B/T>1/3, and they follow elliptical galaxy scaling relations in the fundamental plane. Pseudobulges, including B/P bulges, have near exponential (n<2) surface brightness profiles, more flattened isophotes, approximately flat dispersion profiles, they are rotationally supported with near-cylindrical rotation (Vrot/σ) ≥ 1), B/T < 1/3, and are disjoint from the classical bulge scaling relations. In a number of galaxies both a classical bulge and a pseudobulge coexist (Erwin et al. 2015). Classical bulges are found in nearly half of early-type spirals and are absent in Sc and later-type systems. Pseudobulges spread over all types of spiral galaxies and their scalelengths are related to the disk scalelength.
The Milky Way galaxy is considered to be of SBbc type, but could be earlier (SBb) as indicated by its high bulge mass. It appears to fall into the so-called green valley region that separates the red and blue sequences in the galaxy colour-magnitude plot (Bland-Hawthorn & Gerhard 2016). Central velocity dispersions of σ 140 km s −1 measured for the Bulge (Blum et al. 1995;Ness et al. 2016;Zoccali et al. 2017) are intermediate between typical values for pseudobulges (70-150 km s −1 ) and for classical bulges (100-230 km s −1 ) as in Fabricius et al. (2012). Lick indices, e.g. Mgb∼3.0 (Puzia et al. 2002) compare well with estimated values of Mgb > 2.35 for classical bulges (Fisher & Drory 2016). The barred 3D density and cylindrical rotation of the Bulge suggest a box-peanut (B/P) pseudobulge (Shen et al. 2010;Wegg & Gerhard 2013).
The first detailed models focusing on the formation of the Bulge were chemical evolution models (Matteucci & Brocato 1990;Mollà et al. 1995;Ferreras et al. 2003), in which the common hypothesis was that the Bulge would form on short timescales, as inspired by early disssipative collapse models (e.g. Larson 1974). In modern hierarchical ΛCDM models, galaxies form upon the collapse of the baryonic matter within dark matter haloes. It was first found that these initial galaxies would then subsequently merge with other discs, thus spheroids would grow through early mergers (Steinmetz & Müller 1995;Springel & Hernquist 2005). However, when mergers are prevented, galaxies with predominant disks form from the inside out (e.g., Samland & Gerhard 2003). Nowadays, strong feedback is believed to decouple the gas from the merging dark matter halos, leading to the formation of realistic disk galaxies in cosmological simulations (see, e.g. Stinson et al. 2006;Aumer et al. 2013;Naab & Ostriker 2017). Models based on observations of nearby galaxies including the MW suggested another channel for the formation of spiral bulges, where through secular evolution by gas accretion and star formation in the disk first a bar forms and then the bulge is built by buckling of the bar (Combes et al. 1990;Raha et al. 1991;Norman et al. 1996). Rapid gas accretion in the early phases of the Universe, forming starbursts (Elmegreen 1999) and turbulent (Brook et al. 2004) clumpy disks have also been shown to produce bulges Bournaud et al. 2009).
These early different ideas, and sometimes conflicting scenarios for bulge formation, are discussed in excellent reviews such as Bouwens et al. (1999); Kormendy & Kennicutt (2004); Athanassoula (2005). IAU Symposia dedicated to galaxy bulges and to the Galactic Bulge were held by Dejonghe & Habing (1993), and Bureau et al. (2008). Annual Review articles dedicated to bulges and the Galactic Bulge were written by Frogel (1988), Wyse et al. (1997), and Kormendy & Kennicutt (2004). More recent reviews on bulges in general were presented by Laurikainen et al. (2016) and Naab & Ostriker (2017), whereas recent ones more focused on the MW Bulge were presented by Rich (2013), the PASA volume 33 (2016), Nataf (2017), and part of the contributions to IAU 334.
In this review, we describe the various tracers of the multiple populations in the Bulge, as well as their transitions into the present thin disk, thick disk, and extant early spheroid (old bulge/inner halo), and their chemodynamical properties. We begin in Section 2 with a review of the structure and dynamics of the Bulge. In Section 3 we describe our current knowledge of the Bulge stellar population properties (ages, metallicity distribution, chemical abundances, and initial mass function), discuss nucleosynthesis implications, and compare with halo and thick disk properties. In Section 4 we present various chemical and chemodynamical bulge models from the literature and discuss how the different approaches can lead to a more complete fully self-consistent chemodynamical bulge model in the cosmological framework. Conclusions and outlook are presented in Section 5. This review is being written to set the scene for several major surveys which will provide us with much-improved information and insights about the MW Bulge, such as the full Gaia astrometric catalog, the photometric surveys VVV/X, LSST, Plato, and several large spectroscopic surveys such as APOGEE-South, 4MOST, and MOONS (Table 1).

Structure and dynamics
In this section we describe the spatial distribution and kinematics of stars in the Bulge, and how these depend on the stellar metallicities. We review the evidence that the bulk of the Bulge stars constitute a barred bulge within which the metal-rich stellar populations have a pronounced box/peanut (B/P) shape similar as in other barred galaxies. Finally, we review constraints on the stellar mass and dark matter in the Bulge.

Bulge structure and morphology from star counts
The first secure morphological evidence for a triaxial bulge in our Galaxy came from COBE near-infrared photometry which unambiguously showed that the Bulge was both larger and brighter at l > 0. Several photometric models for the Bulge modeled this asymmetry (Dwek et al. 1995;Binney et al. 1997;Bissantz & Gerhard 2002) and showed that the Bulge had a barred shape, which extended to ∼ 2 kpc from the center and had its near side in the second Galactic quadrant, with a bar angle of φ bar =∼ 25 • ± 10 • between the bar's long axis and the Sun-Galactic center line. Early star counts using the Red Clump (RC) confirmed that indeed stars on the l > 0 near side of the Bulge are closer to us (Stanek et al. 1994).
Red clump giants (RCGs) are very useful tracers for the Bulge. They are He-core burning stars with standard candle properties except for young (ages < 2−4 Gyr) or metalpoor ([Fe/H]<−1) populations; see Salaris & Girardi (2002); Girardi (2016). For typical Bulge populations RCGs have a narrow range of absolute magnitudes and colors σ(Ks)  Rauer et al. (2014); Miglio et al. (2017) 0.17 and σ(J − K s ) 0.05, and individual distances can be determined to within ∼ 10%. Variations in Ks due to age are ∼ 0.03/Gyr at age 10Gyr, and due to metallicity lead to a spread σKs(Fe/H) ∼ 0.11 for the Bulge MDF. Relative numbers of RCGs are ∼ 10 −3 / M and are predicted to vary by less than 10% for old stars with metallicities in the range 0.02 to 1.5 solar. Among the 14,150 stars of the ARGOS survey within Galactic radius RG < 3.5 kpc, RC stars are prominent to [Fe/H]= −1.0, corresponding to ∼ 95% of the ARGOS sample (Ness et al. 2013a); i.e., RCGs are representative for most of the Bulge stellar mass.
RCGs have therefore played an important role in clarifying the structure of the Bulge (Babusiaux & Gilmore 2005;Rattenbury et al. 2007). RCG star counts in the OGLE, 2MASS, and VVV surveys (see Table 1) showed the presence of an X-shaped structure (Nataf et al. 2010;McWilliam & Zoccali 2010;Saito et al. 2012b) similar to B/P bulges in barred galaxies (Bureau et al. 2006) and N-body simulations (Li & Shen 2012). Ness et al. (2012) showed that only stars with [Fe/H]> −0.5 participate in the split RC (X-shape). Gonzalez et al. (2015b) showed that a scenario of one population significantly enhanced in helium does not explain the observed variation of the RC across different (l, b).

Figure 1
Top left: WISE 3.4 µm/4.6 µm contrast-enhanced image of the Bulge region, showing the X-structure. Bottom left: surface density of the Galactic bulge viewed side-on, obtained by projecting the 3D density of red clump giant stars measured by Wegg & Gerhard (2013) from VVV data. Between the dashed lines (|z| = 160 pc, |b| = 1 • ), the density has been extrapolated down from higher latitudes, due to crowding uncertainties. Contours are plotted every 0.33 mag and show the exponential vertical density profile. Right: stellar surface density for the entire bulge/bar region, according to best-matching model of Portail et al. (2017a). This shows how the B/P bulge transits continuously into the outer planar bar which extends to ∼ 5 kpc from the Galactic Center. Also shown are lines-of-sight from the Sun towards longitudes l=−10 • , −5 • , 0, +5 • , +10 • . Adapted respectively from Ness & Lang (2016), Portail et al. (2015a), Portail et al. (2017a), with permission.
most crowded region |b| < 1 • , and with typical (systematic) error of ∼ 10%. They found a strongly barred and peanut-shaped density with bar angle φ bar = 27 • ± 2 • , and face-on projected axis ratio (1 : 2). Viewed from the Sun the RCG Bulge appears boxy, slightly peanut-shaped, consistent with the earlier COBE and 2MASS data and with recent WISE photometry (Ness & Lang 2016). Figure 1 illustrates these results. Unsharp masking results in a strong off-centered X-shape structure (Portail et al. 2015a), similar to some galaxies in the sample of Bureau et al. (2006). Along the bar axes, the central (< 1 kpc) density distributions were all found to be nearly exponential. In particular the minor axis profile is exponential in the range 500 pc < z < 1.2 kpc, with short scalelength z0 = 180 pc, and shows no indication of a central R 1/4 component. Subsequently Wegg et al. (2015) showed, using the combined 2MASS, UKIDDS and VVV surveys, that the MW's B/P Bulge continuously transits outwards into a planar bar with half-length of 5.0 ± 0.2 kpc and bar angle φ bar =∼ 28 • − 33 • (assuming R0 = 8.3 kpc). This work superseded earlier work on the long bar which had suggested different orientations of the Bulge bar and long bar (Hammersley et al. 2000;Benjamin et al. 2005;Cabrera-Lavers et al. 2008); see also Bland-Hawthorn & Gerhard (2016). The MW is thus similar to external barred galaxies that often contain a B/P bulge embedded in a longer, thinner bar (Kuijken & Merrifield 1995;Bureau & Freeman 1999;Erwin & Debattista 2013). N-body models for bar-unstable disk evolution regularly find that B/P bulges are the inner 3D parts of a longer, planar bar and form through a buckling instability and/or through orbits in vertical resonance (Combes et al. 1990;Raha et al. 1991;Athanassoula 2005). Fig. 1 shows the surface density in the MW bulge/bar and inner disk region of the best dynamical model of Portail et al. (2017a) fitted, amongst others, to the RCG density of Wegg & Gerhard (2013) and the RC magnitude distributions of Wegg et al. (2015).
Recently, large samples of Bulge RR Lyrae (RRL) stars have been identified in the OGLE and VVV surveys. RRL are pulsating, low-metallicity core-helium burning horizontal branch giants (Marconi et al. 2015) with ages > 11 Gyr (Walker 1989)   . They do not participate in the B/P-Bulge (Dékány et al. 2013;Pietrukowicz et al. 2015;Gran et al. 2016), consistent with the ARGOS result that only stars with [Fe/H] > ∼ − 0.5 participate in the split RC (Ness et al. 2012). While Dékány et al. (2013) find weak if any triaxiality in the central kpc, Pietrukowicz et al. (2015) report strong triaxiality out to ∼ 1.5 kpc, approximately aligned with the main Bulge/bar. This controversy still needs to be resolved. In the outer Bulge (Gran et al. 2016, |b| = 8 • −10 • ) and near the Galactic plane for 350 • > l > 295 • ) the RRL distribution do not trace the bar seen in the RCG. Recent spectroscopy shows that the RRL have no significant rotation for |l| < ∼ 3.5 • , rising to vrot 50 km s −1 at l = 9 • , b = 3 • (Kunder et al. 2016, and in preparation), compared to 90 km s −1 for RGB and RCG stars (Kunder et al. 2012).
These properties suggest that the Bulge RRL trace a different component from the RCG stars. From the number ratio of RRL to RCG this very old population is estimated to contribute of order 1% of the Bulge stellar mass (Kunder et al. 2016). A natural question is whether it is related to the inner stellar halo, which must invariably be present in the Bulge. There is currently no evidence for a break in the mean metallicity (Pietrukowicz 2016) or density (Pérez-Villegas et al. 2017) of the RRL from the Bulge into the halo.  suggest from the RRL MDF that their sample is dominated by halo RRL. More detailed density maps, kinematics, and metallicity maps must show whether there exists a separate Bulge RRL component. Pérez-Villegas et al. (2017) illustrate that in regions where the main B/P Bulge dominates the gravitational potential, also the RRL orbits are generally non-axisymmetric. They show that an initially spheroidal halo-bulge perturbed by the evolving bar and B/P Bulge would at the end show some triaxiality in the vicinity of the main stellar Bulge, due to a fraction of orbits trapped by the bar. This process (see also Saha et al. 2012) induces rotation velocities consistent with Kunder et al. (2016) if the initial rotation is small. At larger distances from the B/P Bulge, the perturbing effects on the initial halo distribution decrease rapidly, with only mild effects expected in the outer Bulge.
In conclusion, RCG star counts show that the Bulge is a B/P bulge similar to bulges in many barred galaxies. RCG trace ∼ 95% of the Bulge stars over a wide range of metallicities, and the combined distribution of these populations is strongly barred with a nearly exponential minor-axis density profile. However, only stars with [Fe/H]> −0.5 participate in the B/P Bulge feature. A small fraction of the Bulge traced by very old, metal-poor RRL stars (∼ 1%) appears to have a distinct morphological structure.

Kinematics
Early kinematic studies of the Bulge concentrated on BW, measuring a line-of-sight (LOS) velocity dispersion σr = 110 km s −1 for M and K giants in this low-extinction field (Mould 1983;Sharples et al. 1990;Rich 1990). With the recent large spectroscopic surveys we now have radial velocity information for many 1000s of stars spread over most of the Bulge, leading to a greatly improved understanding of Bulge kinematics and, ultimately, dynamics.
2.2.1. The velocity dispersion profile. Studying K-giants in two fields, Minniti et al. (1992); Minniti (1996) found that the velocity dispersion for Bulge giants steeply decreases outwards. This was confirmed in much more detail first using Planetary Nebulas (PNs) by Beaulieu et al. (2000), and then by the recent large surveys, the BRAVA survey for M-giants (Rich et al. 2007b), and the ARGOS, GES, and GIBs surveys for (mostly) RC giants. Using NIR spectroscopy of the 2.3 µm CO bandhead, Blum et al. (1995) found that the dispersion increases inwards to σr 140 km s −1 at 1.2 • from the center. This is consistent with the minor axis profile from BRAVA (Kunder et al. 2012) and the σ = 140 km s −1 peak at b = 1-2 • found by GIBS using the infrared Ca triplet . Combining data from the ARGOS and APOGEE surveys, Ness et al. (2016) also found that the Bulge dispersion reaches the highest values (σ ∼ 130 km s −1 ) within 2-3 • of the center. The last two papers give approximate velocity dispersion maps across the Bulge, which show the overall decrease of the dispersion from the center outwards.
2.2.2. Cylindrical rotation. The first Bulge rotation curve was measured by Minniti et al. (1992), combining K-giant velocities in two outer-Bulge fields at l, b = (8 • , 7 • ), (12 • , 3 • ) with results from BW. They found an increase of vrot to 80 km s −1 in their field at l = 12 • which was confirmed by Beaulieu et al. (2000). Much more detailed studies of the Bulge rotation field became possible with modern surveys. The BRAVA survey measured velocities for ∼8,600 M giants in sequences of Bulge fields at different longitudes for b = −4 • , −6 • , −8 • . These data showed that the Bulge rotation curve turns away from solid body rotation at l ∼ 5 • and that at any given l the rotation velocity is nearly independent of b for b < −4 • . This so-called cylindrical rotation was suggested earlier by the PN survey of Beaulieu et al. (2000) and is confirmed by the ARGOS survey for (mostly RC) stars with [Fe/H]> −1.0. The combined ARGOS and APOGEE surveys show that the Bulge rotates cylindrically across all latitudes up to b < ∼ 10 • . Both the rotational properties and the velocity dispersion show a smooth transition from the Bulge into the disk and bar region further out.
These kinematic properties are similar to the cylindrical rotation seen in edge-on barred galaxies (Shaw 1993;Molaeinezhad et al. 2016) and in N-body models of B/P bulges that originated from dynamical instabilities of their initial disks. Shen et al. (2010) compared such models to the BRAVA data, and found a very good match to the cylindrical rotation in the Bulge with a pure B/P bulge model while models containing an initial classical bulge (ClB; see Introduction and Kormendy & Kennicutt 2004) gave less good fits to the data. Saha et al. (2012) and Saha & Gerhard (2013) studied the spin-up of ClB by evolving bars and found that in some models the ClB may acquire cylindrical rotation as well, but that the induced rotation profile rises slowly and approximately linearly with radius. Thus while the cylindrical rotation per se cannot be used to rule out a ClB component in the MW Bulge, more detailed studies of the rotation properties of different stellar populations in the Bulge may set more stringent constraints on a possible ClB.
2.2.3. Streaming velocities along the line-of-sight. Several studies have obtained distance estimates from the magnitudes of their target stars, thereby measuring variations of the radial streaming velocities with distance. For a bar at the orientation of the Galactic bar, one expects LOS streaming velocities to be more negative on the near side of the Bulge where stars approach the observer, and more positive on the far side. While at low latitudes (b < ∼ 6 • ) such streaming velocity differences have been seen (Rangwala et al. 2009;Babusiaux et al. 2010Babusiaux et al. , 2014Vásquez et al. 2013), studies at higher latitudes have not found a similar signal (De Propris et al. 2011;Uttenthaler et al. 2012;Rojas-Arriagada et al. 2014). The ARGOS data were recently reanalyzed by Portail et al. (2017a) including the full survey selection function. Their Fig. 17 illustrates the dependence of the streaming motions on distance modulus for b = −5 • and b = −10 • , both in the ARGOS data and in their overall best B/P bulge model. A related variation with distance has been found in some proper motion studies, due to stars streaming towards positive (negative) l on the near (far) side of the bar (e.g., Kuijken & Rich 2002;Poleski et al. 2013, and below). These streaming velocity patterns provide clear evidence for the barred nature of the Bulge.
2.2.4. Proper motions, vertex deviation and anisotropy. In a barred bulge tilted with respect to our LOS, the principal directions of the velocity ellipsoid are also tilted. To see this effect requires measuring line-of-sight velocities and proper motions (PM). In a pioneering paper, Spaenhauer et al. (1992) measured PMs for ∼ 400 stars in BW from photographic plates taken in 1950 and 1983. An analysis of 61 of these stars with 3D kinematics by Zhao et al. (1994) indeed showed velocity ellipsoid tilts suggesting bar dynamics. Another important application of PM studies is to separate foreground disk stars from Bulge stars (Kuijken & Rich 2002;Vieira et al. 2007). The disk stars have an offset centered around µ l = 4 mas yr −1 relative to the mean (used as zero point in the absence of an absolute reference frame). For population studies a pure Bulge sample can then be isolated (e.g., requiring µ l < 0), but this will be kinematically biased.
Selecting instead Bulge stars in the CMD and measuring their PMs and radial velocities enables kinematically unbiased studies of the 3D velocity distribution. This can generally be characterized by the velocity dispersions σr, σ l , σ b and the cross terms σ rl , σ rb , and σ lb . Soto et al. (2012) showed that the vertex deviation lv, where tan 2lv = 2σ 2 rl /[σ 2 r − σ 2 l ], decreases in a field at b = −6 • relative to BW and the SgrI field at b = −2.5 • . The actual value of the vertex deviation is not easily interpreted -it depends on the distance distribution of the stellar sample and must be related to 3D information through a good dynamical model for the Bulge. Babusiaux et al. (2010) and Soto et al. (2012) give measurements for the anisotropy in their Bulge fields, typically σ b /σr = 0.8-0.9, and Koz lowski et al. (2006) and Rattenbury et al. (2007) measure anisotropies σ l /σ b = 1.1-1.2 and a significant negative cross term σ 2 lb /(σ l σ b ) −0.10. Altogether these measurements imply a barred potential. PMs are also very useful for more targeted investigations into the properties of the stars in the Bulge X-shape. Poleski et al. (2013) and Vásquez et al. (2013) studied kinematics of the nearby and distant arms of the X defined through RC stars and detected clear differences in the radial velocities and PMs of the two arms. A quantitative interpretation of these data with Bulge dynamical models is still pending.

Kinematics with Metallicity
First detections of kinematic differences with metallicities go back to at least Sharples et al (1990) and Rich (1990). Rich found a lower σr = 92 ± 14 km s −1 for a metal-rich subsample of K giants with [Fe/H]≥ 0.3, than σr = 126 ± 22 for the metal-poor subsample with [Fe/H]< −0.3. Minniti (1996) observed that halo giants with [Fe/H]< −1.5 in a field at l, b = (8 • , 7 • ) have higher dispersion and lower rotation than the Bulge giants ([Fe/H]> −1.0), and also that the Bulge giants themselves have kinematics correlated with metallicity. Combining with BW data he found that for the Bulge giants σr decreases steeply with distance from the Galactic Center, whereas for the halo giants it hardly decreases. Zhao et al. (1994); Soto et al. (2007); Babusiaux et al. (2010) and Hill et al. (2011) found differences between the vertex deviations of metal-rich and metal-poor Bulge stars. Fig.2 gives the vertex deviation as a function of metallicity in BW, together with the chemodynamical model of Portail et al. (2017b) discussed in Section 4.2. The more metalrich ([Fe/H] > ∼ − 0.5) stellar population is found to have a clear tilt (σ rl = 0), providing evidence that its spatial distribution (and the potential) is barred.
Recent spectroscopic surveys agree on distinct kinematic properties between metal-rich and metal-poor Bulge stars. The ARGOS survey observed the Bulge in multiple fields with |l| ≤ 20 • , |b| = 5−10 • , selecting equal numbers of stars in 3 magnitude bins based on 2MASS photometry. This strategy includes Bulge stars along the entire LOS, giving high weight to the outer Bulge (Freeman et al. 2013;Portail et al. 2017a). For LOS distances within ±3.5 kpc around the MW's center, Ness et al. (2013b) (Ness et al. 2012). For stars in C which do not participate in the B/P Bulge, they suggest a thick-disk interpretation, consistent with cylindrical rotation, while stars in D may be part of the stellar halo. Dynamical models for these different stellar populations are in broad agreement with this interpretation, as discussed in Section 4.2.
The APOGEE survey has observed a large number of 2MASS giants in the NIR, thereby also covering the low-latitude inner Galaxy. Combining the APOGEE and ARGOS data, Ness et al. (2016) constructed latitude-combined rotation and dispersion profiles for the same metallicity bins, showing a smooth transition from the Bulge to the thin bar and disk. In these average profiles, A and B are very similar, with C showing higher dispersion and lower rotation.
GIBS surveyed fields in the inner Bulge (|l| < ∼ 8.5 • , |b| = 1.4−8.5 • ), targeting RC and underlying RGB stars selected from VVV, which are mostly located in a slice of ±600 pc width along the LOS around the Bulge center ). The GIBS MDF divides into two components, separated approximately at [Fe/H]= 0 (see §3), for which Zoccali et al. (2017) show separate surface density maps and kinematic profiles. While for the metal-poor stars σr varies mildly with latitude, consistent with ARGOS, the metal-rich component has a steep gradient with |b|, so that its σr becomes roughly equal at |b| = 2 • and higher than that of the metal-poor component at |b| = 1.4 • .

Figure 2
Vertex deviation θ l,r in BW as a function of metallicity. Data are from Babusiaux et al. (2010) and Soto et al. (2007). Overplotted is the fiducial chemod-dynamical model of Portail et al. (2017b, see Section 4.2) including all stellar particles between 6 kpc and 10 kpc along the line-of-sight from the Sun. The shaded area indicates the one-sigma statistical error on the model vertex deviation. The large errors in the model predictions at low metallicity arise from the fact that σ 2 r and σ 2 l are increasingly similar to each other with decreasing [Fe/H], and therefore the angle of the near-circular velocity ellipsoid is poorly defined. Adapted from Portail et al. (2017b) with permission.
targeting also RC and RGB stars selected from VVV in a distance range ±2. In summary, the Bulge kinematics suggest a complex multi-component chemodynamical structure of the Bulge. Differences between recent surveys are likely to be related, at least in part, to different survey selection functions coupled with LOS metallicity gradients. The metal-rich Bulge stars are clearly part of the B/P structure, they show cylindrical rotation, a clear vertex deviation, and their velocity dispersion profile decreases steeply with latitude. The more metal-poor stars, which do not participate in the B/P structure, also show nearcylindrical rotation but no vertex deviation, and their dispersion profile is approximately constant with latitude. The most metal-poor have lower, non-cylindrical rotation, but these stars are rare in the Bulge and not well-sampled in the current data.

Bulge Dynamics and Mass
After the COBE data had established the non-axisymmetric, boxy shape of the Bulge, several barred dynamical models were built (Zhao et al. 1994;Fux 1997). Little stellar kinematical data was available at the time to constrain these models, but a barred bulge formed through disk instability appeared to be a promising explanation for the observed kinematics (Beaulieu et al. 2000). More detailed modeling became possible with the modern survey data. Shen et al. (2010); Li & Shen (2012); Ness et al. (2013b) showed that the observed X-shape and cylindrical rotation could be reproduced with an N-body B/P Bulge model. Gardner et al. (2014); Qin et al. (2015) explained the different kinematics of stars on the near and far parts of the X in the Bulge. Portail et al. (2015b); Williams et al. (2016); Abbott et al. (2017) showed that the B/P Bulge is maintained by a wide range of orbits, both resonant and non-resonant.
The mass of the Galactic Bulge has been estimated many times; see Licquia & Newman (2015); Bland-Hawthorn & Gerhard (2016). The most accurate determinations of the dynamical mass are from recent made-to-measure (M2M) dynamical models which combine well-determined tracer densities with multiple field kinematics (Portail et al. 2015a(Portail et al. , 2017a. The stellar mass in the Bulge is best estimated from RCG star counts containing information on the LOS density distribution, combined with an empirical estimate of the RCG abundance per unit mass Portail et al. 2017b). The main systematic uncertainty is that, because of the continuous transition between the 3D inner B/P Bulge bar and the flatter long bar, quantities like the total Bulge mass (vs. the mass within a fixed volume) require a somewhat arbitrary definition of what stars constitute the Bulge. Portail et al. (2015a) constructed M2M dynamical models of the Bulge matching the 3D density of RCG stars measured by Wegg & Gerhard (2013) from VVV data together with the BRAVA kinematics. They obtained a total dynamical mass of the Bulge in the 3D VVV box, (±2.2 × ±1.4 × ±1.2 kpc), of 1.84 ± 0.07 × 10 10 M , and comparing to stellar population models, were able for the first time to estimate the fraction of dark matter mass in the Bulge. Portail et al. (2017a) subsequently included additional kinematic data from the ARGOS and OGLE surveys, as well as the planar long bar based on its RCG density from Wegg et al. (2015). This leads to better constraints on the Bulge as well, because long-bar stars moving through the inner Bulge influence the projected velocity dispersions. With their new models, Portail et al. (2017a) obtained an improved measurement of the pattern speed of the MW's B/P Bulge and bar, Ωp = 39.0 ± 3.5 km s −1 kpc −1 , putting corotation at 6.1 ± 0.5 kpc. They redetermined the total dynamical mass in the VVV box to 1.85±0.05×10 10 M , in excellent agreement with Portail et al. (2015a), and found a dark matter fraction there of 17 ± 2%. Together with the Galactic rotation curve in R = 6−8 kpc, this low dark matter fraction in the Bulge can only be accommodated if the MW's dark matter halo has a core or mild cusp. Portail et al. (2017a) determined photometric stellar masses for the Bulge, bar, and disk, based on their calibration of the number of RCG stars per unit stellar mass. Their method is analogous to that of Valenti et al. (2016), but identifies RCGs statistically as the excess above the smooth background of stars in the extinction-corrected magnitude distribution. This results in 1000 ± 100 M /RCG, in good agreement with the predicted value based on the initial stellar mass function (Kroupa 2001;Calamida et al. 2015;Wegg et al. 2017, see Section 3). Then defining the Bulge and long bar as 'photometric' density excess over the constant-density inner disk, Portail et al. (2017a) obtain 1.34±0.04×10 10 M , 0.54±0.04× 10 10 M , and 1.29 ± 0.12 × 10 10 M for the mass of the barred Bulge, long bar, and inner disk around the long bar. Added to this is the mass of the nuclear disk, 2.0 ± 0.5 × 10 9 M . In total, the region inside radius 5.3 kpc thus contains ∼ 65% of the MW's stellar mass (∼ 5 × 10 10 M , Bland-Hawthorn & Gerhard 2016). Valenti et al. (2016) estimated the stellar mass from counting the number of RCG+RGB stars near the peak of the RCG in the Bulge fields, and calibrating on a NICMOS field at (l, b) = (0, −6 • ) where the stellar mass function was earlier estimated from HST data (Zoccali et al. 2000). Scaling to the entire projected (|l| ≤ 10 • , |b| ≤ 9.5 • ) Bulge region, they obtained 2.0 ± 0.3 × 10 10 M which is higher than the result from Portail et al. (2017a).
In determining this number an 18% correction for disk-star contamination is applied which is determined at (l, b) = (0, −6 • ) and is taken to be constant in all fields. This is likely to lead to an overestimate of the Bulge stellar mass because the fraction of foreground disk RGB stars increases towards low latitudes . In addition, the projected mass of Valenti et al. (2016) is likely to include a fraction of the long-bar mass of Portail et al. (2017a), reflecting the above-mentioned problem of Bulge definition.

Summary
The evidence from the structural and kinematic properties of the metal-rich stars suggests that a major part of the Bulge was built from a thin disk through evolutionary processes leading to bar and B/P Bulge formation, as also inferred in many external galaxies (socalled secular evolution, Kormendy 2013). There is still debate about the part of the Bulge in the thick-disk metallicity range. As we will see in Section 4.2, a fraction of these stars is in a barred thick disk in the Bulge, but there is a central high-density part whose nature is currently open.
The exponential vertical density of the Bulge is typical of pseudobulges in the nomenclature of Kormendy & Kennicutt (2004) ;Fisher & Drory (2016), however, a detailed analysis is required to quantify a possible ClB component. Howard et al. (2008) estimate a global vmax/σ for the Bulge and find that in the vmax/σ − diagram of Kormendy & Kennicutt (2004) it is located below the oblate-isotropic line, near ClBs and the B/P bulge of NGC 4565, but less rotationally supported than other pseudobulges.

STELLAR POPULATIONS
Stellar populations in the Bulge are characterized by their star formation history (SFH), initial mass function (IMF), and dynamical history, and can be studied in terms of their ages, radial velocities and proper motions, metallicities and chemical abundances. Observational studies use tracers in many stages of stellar evolution, from dwarf to red giants, red clump, blue horizontal branch, RR Lyrae, Cepheids, Asymptotic Giant Branch (AGB) stars, and Planetary Nebulae, besides AGB stages as masers, OH/IR stars among others, and these stellar populations in the Bulge are described. The observed chemical evolution of selected elements in low mass stars is presented (for a general view on abundances and chemical evolution see e.g. Pagel 1997; Matteucci 2001).

Stellar population components and Ages
3.1.1. Color-Magnitude Diagrams. Age derivation from CMDs in the Bulge requires accurate photometric data. A first modern attempt to dating the Bulge was carried out by Ortolani et al. (1995), based on a combination of HST and ground-based data, indicating that the main sequence turn-off (MSTO) is faint, and therefore older than ∼10 Gyr, as later confirmed by several authors (Kuijken & Rich 2002;Zoccali et al. 2003;Brown et al. 2010;Clarkson et al. 2008;Valenti et al. 2013;Gennaro et al. 2015).
Even so, one of the still debated issues in the literature is which fraction of young stars ( < ∼ 5 Gyr) would still be compatible with the observed Bulge CMDs. In the past, one of the main challenges in obtaining a precise answer to this question was the CMD contamination from disk stars. More recently, thanks to proper-motion cleaned CMDs the estimates of the contribution from young stars became more robust. Clarkson et al. (2008) presented a proper-motion cleaned CMD of field stars at (l,b)=(1.25 • , −2.65 • ) within SWEEPS (field indicated in Fig. 3), based on the Advanced Camera for Surveys -ACS/HST data. They found the CMD to be dominated by an old population with an age of 11±3 Gyr and a mean metallicity of [Fe/H]=0.0±0.4. They estimate a maximum fraction of 3.5% of young stars. Clarkson et al. (2011) reported the first detection of Blue Straggler Stars (BSS) in a SWEEPS fields. These are old hydrogen-burning stars, that due to accretion, have larger luminosities and temperatures which are typical of younger stars. Therefore, to better quantify the fraction of genuine young (≤ 5 Gyr) stars, BSS have to be accounted for. Clarkson et al. (2011) were able to estimate the fraction of young-non BSS population, made possible with accurate proper-motion cleaning based on extensive HST data, finding it to be below 3.4%.
More recently, Gennaro et al. (2015) presented the first results of the multi-epoch, multi-band GBTP Program (fields indicated in Fig. 3). They build proper-motion cleaned CMDs in four low-extinction Bulge fields, with a set of filters suitable to derive reliable estimates of metallicities and temperatures of individual stars (Brown et al. 2009;Brown et al. 2010). Their preliminary results show that there is no evidence for young stars (5-10 Gyr), suggesting that most Bulge stars were formed in a short burst at early times. Consistently with the findings of Clarkson et al. (2008), the authors claim that their observed CMDs of Bulge stars rule out an extended star formation period for the Bulge population.
On the other hand, as we will see in §3.1.2, Bulge microlensed dwarfs suggest a different picture, with a much larger fraction of young stars (Bensby et al. 2017). The apparent discrepancy between the results obtained from the CMDs and microlensed dwarfs led Haywood et al. (2016) to suggest that age-metallicity degeneracies could lead to a young population not detectable in the observed CMDs of the Bulge. One of the focuses of the latter work is to show that an old population with a broad metallicity range would not produce the tight TO as observed in the SWEEPS HST CMD, and that by assuming a large proportion of younger stars this claimed discrepancy disappears. Although the effect is clear when plotting simple isochrones, the simulation of the observed CMD is more complex as it depends strongly on the number density of stars corresponding to each isochrone. The authors' analysis focused only on colors (horizontally) to extract constraints in age. However, in their simulated CMD (their Fig. 8) a plume of MS is seen, clearly higher than the number of RGBs at the same magnitude, whereas this is not observed in the data (their Fig. 6, specially if one accounts for BSS). If one takes these extra constraints into account (not only focusing on the TO) it becomes more difficult to believe that a young population hides in the CMD. Further more detailed analysis of the Bulge CMD are needed before firm conclusions can be drawn.

Microlensed Dwarfs.
Microlensing occurs when a star, or a compact object act as a lens, magnifying a more distant star behind it, Bensby et al. (2017, and references therein) analysed 90 dwarf stars in the Bulge. Ages were derived based on spectroscopic determinations of effective temperature (Te) and gravity log g: for [Fe/H]< −0.5 stars are 10 Gyr or older, 35% of metal-rich stars with [Fe/H]>0.0 are younger than 8 Gyr, and 26% stars of the overall sample younger than 5 Gyr. Bensby et al. (2013Bensby et al. ( , 2017 analysed the reliability of their findings, and pointed out biases inherent to their sample, namely, a) the observing strategy to compose the microlensed dwarf sample tends to be biased toward metal-rich and younger stars, and b) age derivations from a probability distribution have a large uncertainty. Bensby et al. (2017) propose a conservative fraction of genuinely young   Table 1). The background is the region as seen by 2MASS.
stars of 15%, when taking into consideration their observational bias. As concerns the He abundance effect on ages (e.g. Nataf & Gould 2012), Valle et al. (2015) pointed out that a He enrichment as high as ∆Y/∆Z=5 implies an age difference of ∼0.6 Gyr. Even adopting the upper values of age by Bensby et al. (2013), there would remain 4 stars younger than 4 Gyr, and 15 younger than 9 Gyr. The location of microlensed dwarfs by Bensby et al. (2017) Figure 3 are well-distributed in (l, b), whereas in depth terms, Bensby et al. (2013) argue that only <3% of their stars might correspond to foreground/far-side populations.

indicated in
3.1.3. Globular clusters and RR Lyrae. Globular clusters located in the Bulge probably have formed in situ, therefore their ages should be indicators of an early Bulge formation.
The term Bulge globular cluster comes from the landmark studies by Minniti (1995) and Côté (1999), based on metallicity, spatial distribution, rotation with respect to field stars, and radial velocities have shown that metal-rich globular clusters within 3 kpc from the Galactic center are associated with the Bulge/bar, and not with the thick disk. However, §2 reports that presently a fraction of the Bulge is dynamically a thick disk.
Among Bulge globular clusters, a prototypical moderately metal-poor case is NGC 6522, estimated to be the oldest globular cluster in the Galaxy by Barbuy et al. (2009). Kerber et al. (2018) found that constraints from RR Lyrae stars require a normal helium abundance Y∼0.26, therefore taking into account [Fe/H]∼ −1.0, [α/Fe]∼+0.3 , and Y=0.26, NGC 6522 appears to have an age close to the Universe. The prototypical metal-rich globular clusters NGC 6528 and NGC 6553 with [Fe/H]∼−0.2, were shown to have similar age and metallicity to BW field stars, and as well being coeval with the inner halo/thick disk reference cluster 47 Tucanae by Ortolani et al. (1995). Therefore, Bulge globular clusters are old and α-enhanced (e.g. Renzini 2006), similarly to most Bulge field stars. Larger samples of Bulge globular clusters with accurate age and metallicity are needed, in order to further confirm the bulge clusters as representative of the old Bulge population, being probably the oldest extant objects in the Galaxy.
Another important indicator of the age and metallicity of stellar populations in the Bulge are the RR Lyrae (RRL). Walker & Terndrup (1991)  for Bulge RRL has been found consistently over the years, and this old stellar population most probably has a spheroidal distribution , and references therein).

Asymptotic Giant Branch Stars.
The Bulge contains a large number of AGB stars, supplying gas to the central parts of the Galaxy (Habing 1996). They are useful tracers of the structure and history of the Bulge, as they cover a large age-baseline. Even if the agedating of AGBs, Miras and Planetary Nebulae (PNe) is not well-constrained, as exemplified in Guarnieri et al. (1997), it is important to note that the fraction of young microlensed dwarfs ( §3.1.2) have to show a counterpart in more advanced stages of evolution. Uttenthaler et al. (2015) report that among 45 AGBs, those in the Bulge are oxygenrich, therefore old, and those in the inner Bulge locations are metal-rich with kinematics compatible with belonging to the bar. Outer-bulge ones, located in Plaut's field at b∼−10 • are metal-poor ([Fe/H]∼ −0.6) with radial velocity dispersion compatible with a metalpoor spheroidal component of the Bulge. Previously, Uttenthaler et al. (2007) identified short-lived technetium in 3 outer-Bulge AGBs.
Mira variables are AGBs in their last stage of evolution, to end later as PNe. Catchpole et al. (2016, and references therein) observed 648 bulge oxygen-rich Mira variables combined with another 8057 inner-Bulge ones from the OGLE, MACHO and 2MASS surveys. The 5 Gyr and younger long period Miras, including the OH/IR stars show a triaxial bar-like structure. The shorter period older ones appear spheroidal or near-spheroidal in the outer Bulge, making evident that the Bulge shape depends on the tracer used. Matsunaga et al. (2017) confirm 4 C-rich Miras located in the Bulge (indicated in Fig. 3). They would be of intermediate age, but could still be old if formed from binary evolution. Old type II Cepheids are a magnitude more numerous than younger classical Cepheids in the Bulge. In a reanalysis of extinction and distances of young classical Cepheids by Matsunaga (2015), and Dékány et al. (2015), Matsunaga et al. (2016) verified that no classical Cepheids are found in the Bulge volume within 2.5 kpc of the Galactic center, except for four of them in the Nuclear Stellar Disk (NSD, §3.1.5). Bulge masers, that are single stars of intermediate mass (1-5M ), and ages of 0.8 to 5 Gyr, in a last-episode of AGB stage, in a phase slightly less evolved than PNe, are found in the bar (Habing 2016). Gesicki et al. (2014) derived masses and ages of central stars of PNe and find ages of 3-10 Gyr. They derive ages scaling them to those by Bensby et al. (2013), therefore it is not an independent determination, and further investigation is needed.
In conclusion, there might be presence of intermediate-age objects in the Bulge, although in small numbers. It is important to point out that the fraction of intermediate age/young dwarf stars in the Bulge, suggested by results from Bensby et al. (2013Bensby et al. ( , 2017, has to show a counterpart as C-rich Miras, and PNe. The current data discussed in this section shows that there is still no clear evidence for that, and that even in the case of the 4 C-rich Miras recently confirmed, there is the possibility that these evolved objects are old if belonging to binary systems. 3.1.5. Detected young populations in the Nuclear Disk. In the innermost parts of the Bulge resides the Nuclear Stellar Disk (NSD) (Launhardt et al. 2002), that has a vertical extent of b<0.4 • or ∼50 pc and Galactocentric radius of The gas to form these young stars could have been lost from AGB (and RGB) Bulge giants, and/or be due to gas flow from the dynamical evolution of the bar, or else to accretion of a dwarf galaxy.

Metallicity Distribution Function
3.2.1. Metallicity Distribution Function traced by giants and dwarfs. Metallicity Distribution Function (MDF) surveys of Bulge stellar populations have used as tracers Red Giant Branch (RGB), Red Clump (RC), and M giant stars. RCG tracers ( §2) are useful tracers also in terms of spectroscopy because their gravity is known, with a typical value of log g=2.5±0.5 for metal-rich and log g=2.5±0.2-0.3 for metal-poor stars of all ages. In the case of a single-aged population, log g is known within 0.1 dex (Girardi 2016). A limitation of RCGs is that they characterize metal-rich populations, with [Fe/H] > ∼ −1.5. RGB stars have the advantage of including all metallicities, but require limits in magnitude range. M giants are bright, and as RCGs, they might miss the most metal-poor component of the bulge population (Rich et al. 2007b).
Early MDF diagrams conclusively indicated that the Bulge contains a metal-rich population ( §1). More recently, Zoccali et al. (2008) presented MDFs for over 500 RGB stars in the Bulge, in three fields along the minor axis, showing that metal-rich stars were present in inner regions, gradually disappearing with increasing latitudes. Babusiaux et al. (2010) further inspected these data in terms of kinematics, confirming that there are two stellar populations: a metal-rich one showing a clear vertex deviation compatible with the bar, and a metal-poor population compatible with an spheroid or a thick disk (see §2). The MDF studies by Hill et al. (2011), Johnson et al. (2013), Bensby et al. (2013Bensby et al. ( , 2017, Rich et al. (2012), Schultheis et al. (2017), and the larger surveys ARGOS, GES, and GIBS are summarized in Table 2. In this table are listed the main MDF studies, indicating the location (l,b), tracer, number of stars, (J-K) color selection, the metallicity peaks reported, and corresponding references. A few more details on the MDF of these authors, including selection functions, can be found in §2.3. Figure 4 shows the MDF from the ARGOS, GIBS and GES surveys, along the minor axis. It shows the probability density distribution of [Fe/H], using a kernel function with a smoothing parameter of 0.17 dex, which is close to the typical accuracy in metallicity. In the top left panel corresponding to the BW at (l,b) = (1 • ,−3.9 • ), the data from Hill et al. (2011) are also included. The top right and bottom left panels show similar comparison for fields located at (0.2 • ,-6 • ), and (359.6 • ,-8.5 • ), respectively. The central field coordinates are from GIBS, and they are compared to GES and to the several small ARGOS fields located nearby (within 1.0 • of the GIBS field centers). In the case of the GIBS survey, the two first panels have data coming from high-resolution spectra, whereas the right panel comes from moderate-resolution spectra (see Zoccali et al. 2017, their Fig. 4).
One of the caveats of comparing MDFs from different surveys is that their shapes may be affected by selection effects as illustrated in the last panel of Fig. 4. Because of the survey selection function designed to give high weight to the outer bulge (see Section 2), the ARGOS MDF in this (l,b)=(0 • ,-5 • ) field has a larger number of stars around [Fe/H]= −0.3 than the MDF reconstructed by Portail et al. 2017b when including the selection function. Their relative number further decreases when only stars within 1 kpc around R0 are included, which is closer to the selection function of GIBS. Finally, this panel also illustrates the effect of survey statistics. In a coarse histogram as would be appropriate for ∼ 200 stars, the peak around [Fe/H]=-0.3 seen in the original data would be very hard to see.
For this reason, we here do not concentrate on the detailed peak structure differences among the samples, but on more global aspects such as their metallicity range. In this framework, the main discrepancy seen between the MDFs traced by ARGOS with respect to GIBS, GES or Hill et al. (2011) in BW, is that it finds fewer metal-rich ([Fe/H] > 0.3 dex) stars. Note that the better agreement between GIBS and Hill et al. (2011) is most probably due to their similar sample selection, and that ARGOS used a different calibrator (the open cluster Mellott 66) for the metal-rich end. Also GES finds a large number of metal-rich stars in BW. Another caveat is that the ARGOS fields close to BW are are not exactly on BW, located around 0.5-1 • below. However it is difficult to explain such a large discrepancy only in terms of field location or any rigid metallicity shift between the ARGOS and the other metallicity scales. One explanation could be that the ARGOS field close to BW (but not on BW) goes less deep into the Bulge, due to the large extinction, and might miss the metal-rich Bulge population at lower latitudes. Samples of RCGs in low extinction regions such as the BW (Hill et al. 2011;Schultheis et al. 2017) or samples including more luminous RGB stars, as for instance the one of Babusiaux et al. (2014), including high extinction fields, are able to detect the Bulge metal-rich component. Schultheis et al. (2017) suggest that the additional peak found in the ARGOS survey might be an artifact from the derivation of temperatures from the (J-K)• colors, and the adoption of samples over a large longitudinal area, where Otherwise, despite the differences in temperature and metallicity, the [α/Fe] ratios from ARGOS and those from UVES are compatible for 4 out of 5 stars. Cross-checks of larger samples with the ARGOS metallicity scale are still missing towards the more metal-rich end. Moreover, as discussed before, the ARGOS selection effects can also play a role on the location of the peaks. In the second panel the MDF of GIBS and ARGOS are similar, whereas GES still sees a large number of metal-rich stars. When combining GIBS fields at different longitudes, but with constant latitute |b| = −6 • , the over-enhancement in metal-rich stars reappears in GIBS (see Zoccali et al. 2017, their Fig. 7), being again compatible with GES. This indicates that at that latitude, a large number of metal-rich stars are present, which is again not seen by ARGOS. In the third pannel, at latitute |b| = −8.5 • , both GIBS and GES give consistent results towards the metal-rich end.
In summary, this figure illustratres how complex still is the analysis of the Bulge MDF, and that it is strongly dependent on target selection, tracers, observed fields, and anaysis methods. This prevents to clearly disentangle the mixed populations in the Bulge. Therefore, it is still unclear if there is or not a genuine chemical gradient in the Bulge, or if it is only the result of the gradient of stellar populations. Because the MDF represents a crucial  constraint to models of bulge formation ( §4), more complete and homogenous surveys of the inner Galaxy region are urgently needed.

Other tracers of the metallicity distribution.
It is important to discuss how other tracers of the Bulge MDF compare with those for giant and dwarf stars. Planetary nebulae (PNe) are the representatives of the end of AGB stellar evolution for stars with initial mass between ∼0.8 and ∼8 M . The currently known number of Bulge PNe is reported to be of ∼800, although estimates as high as 3500 are given in the literature, with abundances derived for about 300 of these according to Cavichia et al. (2017). They analysed 17 high-extinction PNe located in the inner 2 • of the Bulge, detecting considerable difference between inner-and outer-Bulge PNe, with inner-Bulge ones indicating higher progenitor masses, with enhanced N/O, being probably younger objects, compatible with their location closer to the bar. A large sample of 166 Bulge, and 60 inner-disk PNe, were analysed by Chiappini et al. (2009), carrying out a detailed comparison between oxygen abundances in these samples and Bulge field RGBs. PNe have a great advantage with respect to other tracers, namely, the fact that in this case it is possible to resort different selection criteria for assigning them to the Bulge, and hence more easily identify foreground (inner-disk) objects. Chiappini et al. (2009) adopted several independent selection criteria to define their Bulge PNe sample, based on their diameter and radio fluxes, as first proposed by Stasińska & Tylenda (1994), and objects not fitting these criteria, but in the Bulge direction, were identified as inner-disk PNe. They found that abundances of Bulge and inner-disk PN samples coincide, whereas they show a shift of 0.3 dex lower relative to RGB stars, as recently confirmed by Pottasch & Bernard-Salas (2015). To explain this difference, apart from several hypotheses on the analyses, a possible one in terms of stellar evolution is that very metal-rich progenitors would not reach the PNe stage, or that PNe had slightly more massive progenitors, and in this case different populations are being compared.

Abundances in the Galactic Bulge
Chemical abundances provide signatures of the formation and evolution of different stellar populations, and their studies are interpreted as near-field cosmology (Freeman & Bland-Hawthorn 2002). This task, which can be complex even in a very local volume near the Sun (Fuhrmann 2011), turns out to be extremely challenging in the Bulge. As discussed in §3.2, the different stellar populations that co-exist in the Bulge comprise not only the innerhalo, and traces of younger stars associated with the bar and NSD, but also other two very important unknowns, namely, the inner thick and thin disks. Therefore important caveats need to be kept in mind when trying to infer any conclusion from the chemical abundances of the Bulge and how it compares to that of other Galactic components.  3.3.1.1. Oxygen. Oxygen abundances as a function of metallicity are one of the most robust indicators on the process of Bulge star formation rate (hereafter SFR) and chemical evolution, clearly more so than the other α-elements, especially because in this case no contribution from SNIa is expected (e.g. Friaça & Barbuy 2017). Woosley et al. (2002) described the nucleosynthesis processes in massive stars. Zoccali et al. (2006), Lecureur et al. (2007), and Fulbright et al. (2007), found Bulge oxygen-to-iron abundances slightly higher than those found for thick-disk stars by Bensby et al. (2004). This would suggest that the thick disk and Bulge components did not share the same chemical-enrichment history (or IMF). On the other hand, Meléndez et al. (2008) and Alves-Brito et al. (2010) examined a sample of thick-disk and Bulge giants, covering the same stellar parameter space, and found no difference between the two samples. So far it is not clear if there is a systematic shift between the oxygen-to-iron abundances in the Bulge and thick disk, as discussed below.  Zoccali et al. (2006) and Lecureur et al. (2007) are replaced by revised calculations carried out by Friaça & Barbuy (2017), where both the C2(0,1) bandhead at 5635.2Å and the CI 5380.3 A atomic line were taken into account. The Ryde et al. (2010) oxygen abundances, where a few of the same stars from the Zoccali et al. (2006) sample were observed, also take C into account, given that it is derived from OH and CN lines in the H-band. Jönsson et al. (2017) reanalysed 23 stars from Zoccali et al. (2006), and another 10 red giants. There is agreement between the results for the larger samples by Friaça & Barbuy (2017) and Jönsson et al. (2017) in the optical, and Ryde et al. (2010) and Cunha & Smith (2006) in the near-infrared. All NIR data used the H-band, where molecular lines allow a proper derivation of C, N, and O altogether. These data tend to show lower oxygen abundances relative to other work, and should be more accurate. In Fig. 5 the results in grey correspond to oxygen abundance derivations that did not take into account C and N abundances in red giants, or that did not employ high-resolution data. Other issues with data on oxygen are detailed in Friaça & Barbuy (2017). The dwarfs by Bensby et al. (2013) (we plot those with ages >11Gyr) have higher temperatures and are not affected by CNO dissociative equilibrium, but the results might instead be affected by NLTE effects.
Overplotted in Fig Table 3 reports metallicities at which the knee is located in different samples.  Friaça & Barbuy (2017) indicated that νSF = 0.5 Gyr −1 best fits the knee corresponding to the start of SNIa enrichment, and νSF = 3 Gyr −1 fits-well the level of the plateau. Given that chemical evolution models critically depend on the nucleosynthesis prescriptions, their study is more useful in terms of trends as a function of metallicity, in particular indicating the location of the knee and as constraints on nucleosynthesis (Friaça & Barbuy 2017).
3.3.2. Odd-Z elements Na, Al. Both 23 Na and 27 Al are produced, together with Mg, mainly in carbon and neon burning during hydrostatic phases of massive stars (WW95, Sukhbold  Also, about 10% is produced in the hydrogen envelope, in the neon-sodium cycle, and some 23 Na is made from neutron capture on 22 Ne in He burning, in the AGB phase of intermediate-mass stars, through Hot-Bottom Burning processes. In this phase Na is produced, but it is also destroyed , and references therein). At the highmetallicity end, Na appears to increase with metallicity. Cunha & Smith (2006) were the first to point out a high Na abundance in the most metal-rich Bulge stars. Johnson et al. (2014) confirmed this result, considered as a behavior typical of a secondary-vs.-primary element. This can be due to metallicity-dependent yields from massive stars, or contribution from AGBs (Ventura & D'Antona 2009 (2016). Massive AGB stars can also produce some smaller amount s of Al through the MgAl cycle (Ventura et al. 2013). Kobayashi et al. (2006) presented the unique available models for the Bulge, well-reproducing the behavior of Al-to-Fe, and that of Na-to-Fe but with overenhanced Na with respect to observations, and with no upturn of Na at the very metal-rich end.
It is interesting that Rich (2013) and Johnson et al. (2014) noted that in disk samples, [Na/Fe]∼+0.1 is constant for all metallicities, thus differing from the Bulge samples.  3.3.3. Iron-peak elements Sc, Ti, V, Cr, Mn, Co, Ni, Cu, Zn. The lower iron-peak element group includes Sc, Ti, V, Cr, Mn, Fe with 21 ≤ Z ≤ 26. Depending on temperatures and densities they are produced in explosive oxygen burning, explosive Si burning (WW95, Nomoto et al. 2013). The upper iron-group elements Co, Ni, Cu, Zn, Ga and Ge, with 27 ≤ Z ≤ 32, are produced in mainly two processes, namely, neutron capture on iron-group nuclei during He burning and later burning stages, also called weak s-component and the α-rich freezeout in the deepest layers (WW95, Limongi et al. 2003, hereafter LC03, Woosley et al. 2002 (2016) showed that [Cu/O] has much less spread than [Cu/Fe] data, indicating a production of Cu and O in the same massive stars. Both Co and Cu are produced in the alpha-rich freezeout as primary elements (Sukhbold et al. 2016). They are also produced in the weak-s process in massive stars, and in this case as secondary elements. It is interesting to note that their abundances can constrain the relative efficiencies of these two contributions (Woosley, private commun.). The data appear to indicate that the weak s-process is dominant for Cu, and alpha-rich freeze-out for Co.
Zinc Available data are gathered in Barbuy et al. (2015), and more recent data are given in Bensby et al. (2017), and for metal-poor Bulge stars (same references as in Fig. 5).  (Nomoto et al. 2013). At the high-metallicity end, a constant [Zn/Fe] would imply that the ISM is being enriched in both Zn and Fe on similar timescales (either both are coming from SNIa, which is unlikely, or the Zn yields from CCSNe vary with metallicity such that it mimics the timescales of Fe enrichment by SNIa), whereas a decrease of [Zn/Fe] would correspond to the enrichment in Fe (and none in Zn) by SNIa, likewise for the α-elements.
In summary, Ti and Zn are primary elements, produced in massive stars, and have an α-like behavior, similar to O, Mg, Si, Ca. Similarly to halo and thick-disk stars (e.g. Sneden et al. 1991;Nissen & Schuster 2011;Reddy et al. 2006), Mn is deficient in metal-poor stars, confirming that it is underproduced in massive stars; Mn-over-Fe steadily increases at the higher metallicities due to a metallicity-dependent enrichment by SN Ia. Co vary in lockstep with Fe, compatible with a dominant production in alpha-rich freeze-out. Cu behaves as a secondary element, indicating its production in the weak-s process in massive stars. V, Ni, Cr, and possibly Sc and Co tend to vary in lockstep with iron. The models well-reproduce their behavior in the Bulge, except for yields of Sc and V. It is important to stress that [Mn/O], and maybe also [Zn/Fe] at the high-metallicity end, appear to be discriminators between Bulge and thick disk.
3.3.4. Heavy elements Y, Sr, Zr, Ba, La, Eu. Heavy-element abundances have different abundance patterns and scatter for different stellar populations, given their complex nucleosynthesis (Cescutti et al. 2018, in preparation). The slow neutron-capture, s-process, is subdivided in weak, main and strong components, according to their production sites. The weak s-process takes place in He core and C-shell-burning phases of massive stars, producing heavy elements up to A∼ 90. The main s-process takes place in the He burning layers of low-mass AGB stars, during thermal pulses, producing elements with 90≤A≤208. AGBs of low metallicity are responsible for a strong s-component (Bisterzo et al. 2014;Käppeler et al. 2011). The production site of rapid neutron-capture r-process elements has not yet been firmly established. The main candidates are: a) Neutrino-driven winds arising in CCSNe (?) were reviewed extensivey by Arcones & Thielemann (2013), and they might account for the weak r-process, or the 'lighter element primary process (LEPP)' (mainly Sr, Y, Zr) process (Bisterzo et al. 2017, and references therein); b) Neutron star mergers (NSM), where nucleosynthesis occurring during the merging and milliseconds afterwards (Freiburghaus et al. 1999), might naturally provide the neutron-rich environment needed for the r-process. Wanajo et al. (2014) and Goriely et al. (2015) relativistic models produce a strong r-process, including r-elements in the full-mass range (90<A<240). The recent detection of r-process elements in the neutron star merger GW170817 revealed by LIGO and VIRGO experiments greatly reinforces this as a confirmed r-process site (Pian et al. 2017). Other possible sites are, among others: c) neutron-rich material ejection via magnetic turbulence in magnetohydrodynamically driven jets (MHDJ) from CCSNe, resulting from massive stars characterized by a high rotation rate and a large magnetic field necessary for the formation of bipolar jets (Nishimura et al. 2015, and references therein), d) Electron-capture supernovae could produce elements of the first neutron-capture peak (Wanajo et al. 2011).
A new interpretation of neutron-capture heavy-element enrichment takes into account rotation in massive stars. An efficient s-process can be triggered in rapidly-rotating massive low-metallicity stars ([Fe/H]<-1.0) stars (Frischknecht et al. 2012). This is due to rotational mixing allowing large amounts of 14 N from the H-burning shell to migrate to the He-burning core, being transformed to 22 Ne, and then favoring neutron production via the 22 Ne(α,n) reaction (Meynet et al. 2006;Frischknecht et al. 2016). These authors are able to show that large amounts of s-process elements can be produced. This interpretation is in contrast with the well-accepted view that very early in the Galaxy, the typical s-elements (such as Ba, La) abundances are due to their r-process fraction only, as first suggested by Truran    Most probably more than one process is needed to explain the complex observed abundances of neutron-capture elements in the Bulge, as it seems to be the case in the halo (Spite et al. 2018). Reviews on abundances of neutron-capture elements in halo stars can be found in Beers & Christlieb (2005), Sneden et al. (2008). In this review we will describe abundances of heavy elements in Bulge stars only. Fig. 8 shows [Zr,La,Ba,Eu/Fe] from the literature, together with inhomogeneous models by Cescutti et al. (2018).
Europium is an r-process element at a 94% level (Bisterzo et al. 2011), and behaves as an α-element. Similarly to α-elements, it is enhanced in metal-poor stars with  terzo et al. 2014;van der Swaelmen et al. 2016). The available samples correspond to a large spread in these ratios in between these two values, indicating variable contributions of s-and r-processes. As discussed above, s-process elements can be due to transfer from a companion AGB star, or else to spinstars. A first-to-second neutron-capture peak ratio [Zr/La] has a large spread at all metallicities, with the spread increasing for [Fe/H] < ∼ −1.0 that might be explained with the inhomogeneous models.
3.3.5. Summary. As a general conclusion, abundances of the α-elements are powerful indicators of the SFH in the early Galaxy, and of the time of start of enrichment by SNIa. The enhanced [α/Fe] measured in Bulge stars demands that an early-fast chemical enrichment has occurred. The odd-Z elements Na and Al are little studied in terms of chemical evolution, and these elements could represent the bridge between the field and globular cluster stars through Na-O and Mg-Al anticorrelations. Iron-peak elements can distinguish between nucleosynthesis processes, and some ratios (e.g. [Mn/O]) can be discriminants of stellar populations. Nucleosynthesis of heavy elements is the most complex, given the various paths possible. Inhomogeneous models might explain their spread in abundances at lower metallicities.

Comparison of Bulge and thick-disk chemical properties
The identification of Bulge stellar populations with thick-disk ones, or else an older Bulge with its own identity is an issue of intense debate in the literature. As discussed in §4 several mechanisms are at play, such as high-redshift clumpy formation of bulges and thick  (2017) disk, mergers and secular evolution. The debate now concentrates on which one of these are dominant in the evolution of spiral galaxies and bulges. It is difficult to disentangle these processes by studying external galaxies (see §4.4), and the hope is that for the MW it would be possible to offer new insights by obtaining accurate abundance ratios in samples for bona-fide Bulge and thick-disk stars. Expectations are that differences on a) the level of the abundances plateau (in the case of α-elements), and b) the metallicity at which the so-called knee occurs, which reflects the enrichment by SNIa, could be used to assess differences in the SFR (and IMF) of the two components, i.e. their SFH. The comparisons are usually done between Bulge and local thick-disk stars. The controversy on oxygen abundances of Bulge and thick-disk populations being similar has so far no conclusion (see §3.3.1.1) -(see also discussion in Friaça & Barbuy 2017). For Mg, Ca, Si, Ti, Alves-Brito et al. (2010) concluded that they are similar. As for Mg, in a reanalysis of Mg abundances from Fulbright et al. (2007 has shown that the Bulge is more Mg-enhanced than the thick-disk data from Bensby et al. (2005), and the same is found by Schultheis et al. (2017) from APOGEE data, confirming previous findings by Anders et al. (2014), and Bensby et al. (2017), suggesting that SFR was faster in the inner Galaxy. Rojas-Arriagada et al. (2017) found their Bulge stars only slightly more enhanced than thick disk ones, in terms of plateau level, but detected a difference of ∆[Fe/H]∼0.06 in the turnover between the two populations.  concluded that the Bulge [Mg/Fe] appears shifted by ∼+0.16 dex compared to the solar-neighborhood thick-disk stars by Bensby et al. (2005). Such behavior is also confirmed by Bensby et al. (2017). Table 3 gathers some of these results from the literature. For the odd-Z element Na, Rich (2013) and Johnson et al. (2014)

Initial Mass Function (IMF)
The present-day stellar mass function (PDMF) is needed to infer the total stellar mass of a stellar population, including their fainter dwarfs and stellar remnants. The IMF is even more fundamental because of its implications on the chemical enrichment of subsequent stellar generations. Several methods have been employed to measure the PDMF and infer the IMF. In the solar-neighborhood, counting of stars in clusters and the field has given fairly robust measurements except for the lowest mass brown dwarfs (see Bastian et al. 2010;Krumholz 2014, for recent reviews).
In the Bulge, two methods have been used for estimating the PDMF. The first is by counting stars in deep HST observations of small Bulge fields and converting the luminosity function to a stellar mass function using a stellar mass-luminosity relation (Holtzman et al 1998). Zoccali et al. (2000) used NIR HST data in a field at b = −6 • to estimate the PDMF (and for these masses, IMF) in the range 0.15−1 M . They found dN/dM ∝ M −α with a best main sequence slope αms = −1.33±0.07, and a hint of a steepening for masses >0.5 M . This agrees with the IMF obtained by Holtzman et al (1998) in BW over the common mass range, and is similar to the disk IMF obtained by Kroupa (2001). Calamida et al. (2015) derived the IMF in the SWEEPS field at (l = 1 • .35, b = −2 • .65), using proper motions to separate Bulge stars from foreground disk stars, and also correcting for unresolved binary stars. They found that for the same mass range their IMF could be fitted with two power laws with a break at 0.56 M , a slope of −2.41±0.50 for higher masses, and a shallower slope −1.25±0.19 for the lower masses, for 50% binary fraction. Alternatively, Calamida et al. (2015) obtained an excellent log-normal fit to their IMF with a central mass Mc = 0.25 ± 0.07 M and dispersion σ = 0.50 ± 0.01 M . These results are in good agreement with the nearby disk IMF as in Kroupa (2001) andChabrier (2003).
A second method makes use of the observed timescale distribution of microlensing events (MLTD). The timescale of each event depends on the square root of the lens mass, as well as on the relative distances and velocities of the lens and the source. Then, with a dynamical model providing these distances and velocities statistically, the lens mass distribution and then IMF can be inferred from the MLTD (Han & Gould 1996;Calchi Novati et al. 2008). Using the MLTD of ∼3000 microlensing events from the OGLE-III survey (Wyrzykowski et al. 2015) and an accurate dynamical model from Portail et al. (2017a), Wegg et al. (2017) measured the distribution of lens masses in the Bulge and inner disk and from this inferred the IMF down to brown dwarf masses. For a broken power-law parametrisation they obtained αms=1.31±0.10|stat ± 0.10|sys and α bd = −0.7 ± 0.9|stat ± 0.8|sys in the main sequence and brown dwarf regimes, where the systematic errors cover the range 0− 100% in the binary fraction for unresolved lenses. These values, and equivalent results obtained for a log-normal IMF, are again similar within uncertainties to local disk values from Kroupa (2001) andChabrier (2003). In summary, the IMF in the Bulge and inner disk is indistinguishable from that measured locally, despite these being regions where the stars are predominantly 10 Gyr old and formed on a fast α-element enhanced timescale.

Models for the chemodynamical properties of the Bulge
The modeling of the Bulge has seen an enormous progress in recent years, mainly due to the large and complex observational data (both in the local and high-redshift Universe) combined with the evolution in the tools used for N-body and hydrodynamical simulations. The data and models have suggested six (most probably mixed) scenarios for Bulge formation (see recent reviews Somerville & Davé 2015;Babusiaux 2016;Bournaud 2016;Brooks & Christensen 2016;Shen & Li 2016;Naab & Ostriker 2017;Nataf 2017), namely: • Primordial collapse where Bulge and the thick disk form early on simultaneously (also referred to as in situ scenario) via strong gas accretion; • Bulge formation by hierarchical merging of sub-clumps, prior to the disk formation; • Merging of early thick-disk sub-clumps, migrating to the center and forming the Bulge; • Major merger hypothesis, where disk galaxies result from major mergers of gas-rich galaxies; • Formation of the Bulge from the disk through a bar instability ( "secular evolution"); • Formation of a Bulge component triggered by the accretion of dwarf galaxies.
Ideally one would like to predict the dynamics and chemical abundances of Bulge stars within all of these scenarios and compare them with observations in the Bulge. In practice we still lack such complete models and large unbiased samples of genuine Bulge stars. Despite these difficulties, the theoretical descriptions of Bulge formation and evolution have seen advances thanks to the parallel efforts made by using classes of models that do not treat the full problem but concentrate on specific aspects (just as it happens with the data).
Historically the above listed scenarios for bulge formation have been grouped in essentially two main families (see §1 for a brief description on early theoretical work on bulge formation), namely, those forming classical bulges and those involving disk evolution forming pseudobulges (disky and B/P bulges). As it will become clear along §4, early bulge formation scenarios that once were considered to form only classical bulges, can also lead to the formation of pseudobulges on short timescales. In the first group we find scenarios where the Bulge was formed early, on short timescales, and in the second group those where most of the stars now in the Bulge were formed from gas (and stars) transferred from the disk. The bar-like nature of the MW bulge seems to rule out most of the ClB models for the bulk of the stars, but these models might still apply to a minor component of the present Bulge. What are the most common expectations in each of the two groups? Expectations are that in the first group the stellar populations would be dominated by old stars (ages ≥10 Gyr), thus showing enhanced α-elements (with respect to Solar) for most of the metallicity range, and a more isotropic kinematics. Typical expectations in disk evolution scenarios are that the Bulge will have a B/P structure, with stars showing anisotropic kinematics. The expectations of scenarios in this latter group with respect to chemistry are less clear as they can differ from each other depending on the assumptions on when/how was the pseudobulge formed. If very early on, then most of the stars would be thick-disk like (i.e. old), whereas in a more continuous secular scenario, an important part of the Bulge will be made by younger (inner thin-disk like) populations. In the former case, the chemistry as well as the dominance of old stars could be similar to that expected by early bulge formation scenarios. Therefore, the expected fraction of young stars (<5Gyr) and the detailed chemical-patterns expected from the second scenario will be dependent on the timescales of the secular processes.
The B/P shape of the Bulge (see §2) is not incompatible with the properties expected from scenarios where most of the Bulge stars formed on short timescales. In addition, the observational evidence of coexisting probably merger-built and secular-built bulges reinforces this idea (Prugniel et al. 2001;Méndez-Abreu et al. 2014;Fisher & Drory 2016). These analyses, based on morphological, photometric, and kinematic properties of bulges have found that many of them have a composite nature, showing pseudobulges, central disks or B/P bulges (Erwin et al. 2015), and conclude that the high frequency of composite bulges in barred galaxies points toward a complex formation and evolutionary scenario.
Finally, a minor contribution of young populations seem to exist in the innermost regions of the Bulge (see §3.1.5), where the current star formation could be powered by radial gas flows coming from the possible following sources: stellar mass-loss (AGBs), gas infall, barinduced inflow, and bar-spiral arms interactions. This gas is most probably enriched, and expectations are that the metallicity of these young stars would be at least Solar.
Observational constraints such as the Bulge MDF (and its variation along different fields - §3.2), the fraction of young and intermediate age stars (as a whole, as well as its variations within the Bulge region - §3.1), and a detailed map of the chemo-kinematic relations of Bulge stars ( §2.2 and §2.3) are key to constrain the best scenario for Bulge formation. The challenge in the observational side has been to extract these constraints, and at the same time account for selection biases of the studied samples (using different tracers) which potentially affect ages, metallicities and kinematic distributions (see §2, §3).

Chemical evolution models
Three basic observational early results for the Bulge have led to the development of the first Bulge chemical evolution models in the early 90's. These were: a) the fact that the bulk of the Bulge seemed to be composed of old stars, as indicated by the pioneering CMDs for field stars by Terndrup (1988), b) the broad metallicity distribution of Bulge stars reaching oversolar metallicities, as shown by Rich (1988) and c) the [α/Fe] decline observed in Bulge stars appeared to occur at higher metallicities with respect to thick-disk and thin-disk stars of the Solar-vicinity (see §3.3.1, and Table 3).
Triggered by these early results, Matteucci & Brocato (1990) developed one of the first Bulge chemical evolution models, aiming to explain Rich (1988)'s result which suggested a metallicity peak in the Bulge at roughly twice the Solar value (a value that was later revised downwards with better data by McWilliam & Rich 1994, see §3.2), and a metallicity distribution that resembled much more to a closed box model than the one traced by solarneighborhood G-dwarfs, where a deficiency of metal-poor stars with respect to predictions of a closed box model is seen (e.g. Gilmore et al. 1989;Pagel 1997).
Differently from previous simpler approaches, where only the global metal content had been predicted by models assuming the instantaneous recycling approximation, Matteucci & Brocato (1990) traced the chemical evolution of Fe, Si, Mg and O individually, including the contribution of SNIa, important for Fe yields. After exploring the effects of different assumptions on the star formation efficiencies, IMFs, and collapse timescales, the authors concluded that the Bulge should have formed on short timescales (< 1Gyr), with an efficient SFR and an IMF with a larger proportion of massive stars (Salpeter 1955, or a flatter IMF) than the one in the solar-neighborhood (e.g. Scalo 1986;Kroupa et al. 1993;Kroupa 2001;Chabrier 2003). A solid prediction of these models was that almost all Bulge stars would be α-enhanced with respect to the Sun, which is now an observationally established fact (see §3.4). It became clear that the [α/Fe] vs. [Fe/H] diagram observed in the Solar Vicinity was not universal, but should vary for different Galactic components and different galaxies (nowadays known as the time-delay model), and proven very useful to studies of Galactic Archaeology.
With the more detailed chemical information for Bulge stars that have started to appear thanks to high-resolution spectrographs installed on 8-m class telescopes, chemical evolution predictions for a large number of individual chemical species appeared in the literature (e.g. Matteucci et al. 1999;Ballero et al. 2007;Cescutti & Matteucci 2011). However, the assumption that in the Bulge the IMF is flatter than locally seems now ruled out by more recent data (see §3.5). By simultaneously predicting different chemical species, important constraints on the nucleosynthetic sites of the different elements can be used a) to test stellar yields from the different stellar models, at different metallicities, and b) help constraining the SFH of different stellar populations.
Confirming earlier findings by Zhao et al. (1994) and Soto et al. (2007), in 2010-2014 observations of larger samples towards the Bulge have demonstrated this Galactic component to host at least two main stellar populations (e.g. Babusiaux et al. 2010;Hill et al. 2011;Babusiaux et al. 2014, and more recent ones), with different chemistry, ages and kinematics (as discussed in §2 and §3). As a response to the new data, several groups started modelling the chemical evolution of the Bulge as two distinct populations, an older metalpoor and α-enhanced population, and a more metal-rich with disk-like abundance ratios (e.g. Bekki & Tsujimoto 2011;Grieco et al. 2012;Tsujimoto & Bekki 2012). Interestingly, the first chemodynamical models of a MW-like galaxy had already suggested a double (or even multi-) component Bulge stellar populations ( §4.4).
Models forming the Bulge in a two-step process like the ones quoted above are able to reproduce the bulge MDF by combining a metal-poor and a metal-rich population, with different star-formation histories. In particular, the work of Grieco et al. (2015) models a large large variety of chemical species. The authors assume that the metal-poor population formed first and on a short timescale, similarly to previous models, and that later on, a more metal-rich population formed from enriched gas whose origin can be either from the inner-disk gas flows or from material left by the first Bulge population. The authors also suggest the oldest Bulge population to present chemical gradients similar to the ones found in Pipino et al. (2008Pipino et al. ( , 2010 chemodynamical models ( §4.4), adding one more reason for variation of the MDF at different Bulge locations, i.e. those would not only be due to an overlap of the (at least) two different stellar populations co-existing in the Bulge, but also due to an intrinsic abundance gradient of the oldest Bulge populations created by dissipative collapse (see discussion on vertical-gradient in §3.2.3).
Tsujimoto & Bekki (2012), on the other hand, suggested that a significant fraction of the present Bulge is composed of stars initially in the inner part of the thick disc, introducing a radially varying IMF, which tends to be more top-heavy in the inner parts. However, a radially varying IMF seems not to be supported by the recent APOGEE data (e.g. Hayden et al. 2015) which show a thick-disk [α/Fe]-[Fe/H] sequence that is independent both from the Galactocentric distance and the height from the Galaxy's mid-plane. On the other hand, the APOGEE data do not rule out that part of the Bulge could be composed of thick-disk stars. Finally, Cavichia et al. (2014) presented a chemical evolution model for the Bulge and disk, which includes radial flows. In this case an IMF similar to the local one is assumed, together with a collapse timescale of 2 Gyr (slightly longer than most of the models for the Bulge presented above). The authors have explored different gas flow pattern velocities, also investigating the role of a Galactic bar of different ages. According to the authors, a radial inflow of gas with velocities lower than 1 km/s already increases the SFR in the Bulge region by ∼19%, which helps shifting the MDF to larger values, without the need of a flatter IMF or a second infall episode.

Chemodynamical equilibrium models: the dynamical structure of stellar populations in the Bulge
Another important piece of information for understanding the formation of the Bulge is provided by dynamical models addressing how the different stellar populations are distributed in the Bulge. The dynamical models reviewed in § 2 constrain the mass distribution and gravitational potential of the inner Galaxy, which combines contributions from the Bulge and long bar, and the NSD, main disk, and dark matter halo. Different stellar populations 3 will generally have different orbit distributions in this joint potential, resulting in different spatial distributions and kinematics for each population. From the observed positions and velocities of a sufficient number of stars of that population, its specific orbit distribution can be reconstructed using chemo-dynamical equilibrium (CDE) models. These models assume that the stars have reached dynamical equilibrium in the overall gravitational potential. They do not consider how the chemical properties of the gas and stars have evolved over time. A notable property of these models, resulting from Jeans' Theorem (e.g. Binney & Tremaine 2008), is that they constrain also the stellar density and kinematics in regions where no data were available for the orbit reconstruction. Portail et al. (2017b) recently presented the first self-consistent CDE models for the Bulge and entire inner MW. They reconstructed the 3D density, kinematics, and orbital structure of stellar populations in different metallicity bins, using the overall barred gravitational potential from Portail et al. (2017a), distance-velocity-metallicity data in multiple fields from the ARGOS and APOGEE (DR12) surveys, and correcting for the survey selection functions.  Portail et al. (2017b), such that bin C' includes the metal-poor side of the local thick disk MDF but most of the Bulge RRL MDF is contained in bin D'. In order of decreasing metallicity, the photometric selection of the ARGOS sample leads to (23%, 43%, 29%, 5%) of ARGOS stars in the bins (A, B, C, D), while the final CDE model has (52%, 34%, 12%, 2%) in bar-supporting and (38%, 47%, 14%, 1%) in non-bar supporting orbits in these bins, respectively. Since estimated metallicity errors are ∼ 0.1 dex (Ness et al. 2013a), much smaller than the total range, current uncertainties in the MDF (Section 3.2) can only have small effects on the metallicity ordering. Fig. 9 illustrates that stars in all metallicity bins are significantly barred, but in terms of mass, most of the support to the bar is provided by metal-rich stars. Bin A stars contribute most to the Galactic bar and B/P Bulge; they have dynamical properties consistent with a disk origin. Stars in bin B are hotter and rotate slightly faster than stars in A, they are more extended vertically, and contribute somewhat less to the bar and B/P shape. They are consistent with a disk origin formed from stars located initially at larger radii (Di Matteo et al. 2014). Metal-poor stars in bin C' ([Fe/H]≤−0.5) rotate slower and have higher dispersion than the more metal-rich stars. They indicate a weaker support to the bar and do not support the B/P shape. Outside the central kpc, these metal-poor stars are found to have the density distribution of a thick disk; in these regions their vertical profile is exponential with scale-height ∼ 500 pc (Portail et al. 2017b). These stars also show cylindrical rotation (Ness et al. 2013b;Portail et al. 2017b), suggesting a thick-disk nature.
In the inner Bulge (x < ∼ 1 kpc, z < ∼ 0.6 kpc), however, Fig. 9 shows evidence for an extra component of these metal-poor stars with short scale-height rising towards the Galactic center; this component is seen also for stars in the even more metal-poor bin D'. A central concentration of metal-poor stars was also reported by Dékány et al. (2013);Pietrukowicz et al. (2015) in RRL stars, by Schultheis et al. (2015) in metal-poor APOGEE stars, and by Zoccali et al. (2017) from the GIBS survey. These could be thick-disk stars on orbits compressed by the deep gravitational potential of the NSD, stars from the inner halo-Bulge, or stars from a ClB formed by early mergers or strong gas accretion at early times but whose existence has not yet been established. Finally, the very metal-poor stars in bin D' constitute a thick, high-dispersion, more slowly and not cylindrically rotating, but still barred population, which is only weakly constrained by the ARGOS data (see Portail et al. 2017b). This metallicity range includes the RRL population but probably also the metal-poor tail of the bin C' populations; clarifying the nature of this population requires further data and analysis.
The combined orbit distributions of all metallicities in the model of Portail et al. (2017b) naturally reproduce the observed vertex deviations in the Bulge, as seen in Fig. 2. The absence of a significant vertex deviation for stars in this plot with [Fe/H]< −0.5 has been interpreted as a signature of a ClB by Soto et al. (2007). However, the stars with [Fe/H]< −0.5 in BW shown in Fig. 2 are predominantly from the thick-disk component of bin C, which therefore must be the main cause for this decrease in the vertex deviation. The fraction of stars in a ClB defined structurally as in Section 2, is small according to these models. This implies that the dominant population of the old Bulge stars is likely to include the barred thick disk and many of the metal-rich barred Bulge stars.
In summary, the fraction of stars in a ClB defined structurally as in §1 is small according to these models, less than of order 10 percent for the whole bulge (see also §4.4). This implies that the dominant population of old Bulge stars is likely to include the barred thick disk and many of the metal-rich barred Bulge stars (see §3.1).

Inhomogenous models and the first stellar generations
Although the previously discussed dynamical models cannot constrain the timescales of enrichment of the different coexisting populations in the Bulge volume, previous results, combined with those discussed in Sect. 3.1 imply that many of the old stars in the Bulge are in thick disk-like or even bar-like configuration (and that the fraction of not very old stars so far found in the Bulge seems to be larger in the strongly barred populations). Independently from its present kinematic and structural configurations, one can still constrain its origin by studying the detailed chemistry of the old stars. The advantage is that chemistry is preserved from birth (for most elements), and that abundance ratios would hopefully be much more sensitive to the different SFH of the different components (which could instead overlap both in metallicity and ages). In this Section we now discuss some constraints added by the analysis of the detailed chemical patterns and scatter of the oldest bulge stars. The nature of the first stellar generations remains a hot debated topic in the literature, as models of their formation are complex, and have difficulties in predicting from first principles their properties, among which their mass distribution which can range from topheavy IMFs, to more normal ones also including low-and intermediate-mass stars (Bromm 2013;Karlsson et al. 2013, and references therein). Studies of detailed chemistry in the oldest stars found in the MW offers extra constraints to those models (e.g. Frebel & Norris 2015, and references therein). Of relevance to the topic of the present review is that some of these fossil records should be found not only in the Galactic halo but also in the central regions of the MW (see §3.2.4). Among the old Bulge stellar population, the oldest stars would have formed soon after the Big Bang, carrying in their atmospheres the fingerprints left by the first stellar generations on the gas from which these fossil objects formed. There has been an effort to obtain detailed chemical abundances for very metal-poor stars first in the halo (Beers & Christlieb 2005), and more recently in the Bulge (see §3.2.4), where so far, very few stars with [Fe/H] < ∼ −2.0 have been found ( §3.2.4). This is expected if, in the old Bulge the chemical enrichment proceeded in a faster pace than in the MW halo. This would imply that stars as old as the ones found in the halo at metallicities [Fe/H]<−2.5 could be found in the Bulge at larger metallicities (−2.5<[Fe/H]<−1) as proposed by Chiappini et al. (2011, see also Wise et al. 2012. Inhomogeneous chemical evolution models are a useful tool to test the above ideas and help to constrain both the nature and the nucleosynthesis of the first stellar generations contributing to the chemical enrichment, while also giving insights on the enrichment timescale of different stellar populations, such as the halo and the old Bulge. These models have been very useful for interpreting the intrinsic scatter observed in several neutron-capture elements in very-metal poor stars ([Fe/H]<−2.5), while preserving the small scatter in other abundance ratios such as the α-elements (e.g. Argast & Samlamd 2004;Font et al. 2011;Cescutti et al. 2008Cescutti et al. , 2013Chiappini 2013), and are the appropriate tool for studying the chemical enrichment of the earliest phases of the Galaxy assembly, when the interstellar medium was less well-mixed (Thornton et al. 1998). Such models take into account the stochasticity of the IMF and SFR, and provide an important test bench for studies of the production sites of chemical elements, such as r-elements, helping one identifying key abundance ratios, as well as the most promising stellar yields, before implementation in self-consistent cosmological simulations (Scannapieco 2017).
Although the first inhomogenous models have been developed for the halo, they can also be useful to interpret the metal-poor data in the Bulge (Chiappini 2013;Barbuy et al. 2014). Via the study of the detailed chemical patterns in the most metal poor Bulge stars (−2.5<[Fe/H]< -1), one can also address the question of how metal-poor are the oldest Bulge stars, and if the old Bulge is just an extrapolation of the halo towards the inner regions, or if its enrichment has proceeded in a different pace. The results of Cescutti et al. (2018) are promising as they can potentially be used to chemically disentangle bonafide old Bulge from halo stars populating the central regions of the Galaxy. Their results suggest that some of the most metal-poor stars found in the Bulge region by e.g. Howes et al. (2016) could actually be halo stars trapped in the bar potential (see §8).

Evolutionary chemodynamical models
Chemodynamical models not only describe the chemical evolution of the gas and stars, but also follow the physical processes related to the assembly of the Galactic components, which in pure chemical evolution models are introduced by hand (gas accretion, mergers, radial gas flows and Galactic winds). Because of their complexity, various simplified approaches have also been used. Here we describe these different approaches in the following subsections discussing chemodynamical models that consider evolutionary aspects on a) dynamics, but assume a simple "tagging" for the chemical properties ( §4.4.1), b) dynamical and global metallicity properties, but do not consider aspects related to the evolution of individual chemical elements ( §4.4.2, §4.4.3) and c) hydrodynamical and chemical ( §4.4.4). Models in fully cosmological framework are described in §4.5. 4.4.1. Dynamical models with chemical tagging. B/P bulges originate through the buckling instability of bars in disks (Raha et al. 1991;Debattista 2006;Martinez-Valpuesta et al. 2006) or through orbit trapping in the rotating barred potential (Combes et al. 1990;Quillen 2002). These systems are characterized by exponential density profiles and cylindrical rotation. An initial ClB can be spun up by the rotating bar during the growth of the bar and B/P bulge which can lead to cylindrical rotation for the ClB stars in some models (Saha & Gerhard 2013). However, the models predict a gently outwards-increasing rotation profile which is different from the typical rotation properties of B/P bulges. Moreover, N-body simulations also show that, if sufficiently massive, pre-existing ClBs lead to final composite bulges with Sersic index n>1 (Blana et al. 2017) and non-cylindrical rotation (Saha et al. 2016). Therefore, limits to the mass of the ClB can be obtained both from structural and kinematic data. Comparing N-body simulation models with BRAVA data, Shen et al. (2010) constrained the mass of a possible pre-bar Bulge in the MW to < 25% of the total Bulge mass (see also Di Matteo et al. 2014).
Because bars have the tendency to mix stars from different radii and thus to erase population gradients (Friedli et al. 1994), the observed metallicity gradient in the Bulge had long been taken as evidence for a ClB formed through dissipative collapse early in the history of the MW. Dynamical models indicated that this picture needed to be revised. There are two key points: (a) The reshuffling of stellar matter during the bar and buckling instabilities believed to be responsible for creating B/P bulges is too weak to erase population gradients present in the prior galaxy disk, and (b) depending on the spatial distribution and kine-matics of the stars in the disk, the instabilities may spread out stars in the final system, creating or enhancing vertical gradients. Martinez-Valpuesta & Gerhard (2013) showed that violent relaxation during the bar and buckling instabilities in their N-body disk-galaxy-evolution model changed the Jacobi energies EJ of the particles only slightly, compared to the range of energies in the part of the disk participating in the evolution. This implies that population gradients with EJ in the prior system can largely be preserved even if the stars are scattered widely in angles. They illustrated this by showing that a radial metallicity gradient in the disk can be mapped into outward and vertical metallicity gradients in the final B/P Bulge. Bekki & Tsujimoto (2011) had earlier concluded that vertical metallicity gradients in a thick disk would survive through the bar instability, thus explaining the final Bulge vertical metallicity gradient.
Di Matteo et al. (2014) analysed three simulations with different ClB mass fractions. Because the stars migrate in the disk due to angular momentum transfer around the time of bar formation, stars from up to the radius of the outer Lindblad resonance of the bar are mapped into the final B/P Bulge. Particle EJ were again approximately conserved, and stars born at small and large radii in the disk end up preferentially at small and large radii in the B/P Bulge, respectively. Di Matteo et al. (2014) also found that Bulge stars coming from larger birth radii have greater rotational support and velocity dispersion. ClB stars extend to larger heights, have lower rotation, and less steeply decreasing dispersion profiles. Comparing their results to Bulge chemo-kinematics they limit the mass of a ClB in the MW to < 10% of the disk mass.
While these thin-disk evolution models can reproduce the global kinematic and chemical properties of the Bulge, Di  argued that to reproduce the chemokinematic and structural properties at different metallicities an additional pre-cursor thick disk is needed. The primary arguments are that for pure thin-disk initial models, the B/P Bulge should show a split red clump for all stars, and that the metal-poor population (−1.0 < [Fe/H] < −0.5) should be a kinematically hotter and faster-rotating replica of the more metal-rich populations. Both are not observed in the ARGOS data. To test this scenario, Di  investigated N-body models in which the B/P bulge results from the evolution of a galaxy with thin, intermediate, and thick disks with inital velocity dispersions adjusted to those measured near the Sun, and Fragkoudi et al. (2017) analysed models with thin and thick disks. In these models, the B/P Bulge is found to be more pronounced in the kinematically cold populations, and the peanut shape is visible in the thick-disk component only at heights > 1 kpc from the Galactic plane. Also, the thicker and kinematically warmer disks lead to weaker and rounder stellar bars that still rotate cylindrically, as observed in the Bulge. Debattista et al. (2017) describe a process termed kinematic fractionation whereby cospatial stellar populations with different in-plane random motions separate when the bar forms and evolves, in that the radially cooler populations form a stronger bar and a peanutshaped Bulge, while the radially hotter populations form a weaker bar and a box-shaped Bulge. In their experiments an initially radially hot thin disk can become thicker than an initially cool thick disk, provided the ratio σR/σz is about three times greater in the initial thin disk. In this regime it is mainly σR that drives the evolution, not σz (i.e., the vertical thickness). In their star-forming simulation (see next subsection), the younger populations are indeed similarly anisotropic at early times, although the physical origin for this remains to be clarified. In conclusion, according to these models, chemical gradients in the B/P bulge reflect at least in part the chemical and dynamical structure of the thin and thick disks, from which it might have grown.

4.4.2.
Chemodynamical models with simplified chemistry. One of the first chemodynamical models to appear in the literature was the model of Samland et al. (1997). This was a two-dimensional chemodynamical model forming a large spiral galaxy, which the authors suggested to be representative of the MW. The model starts from an isolated, rotating protogalactic gas cloud in virial equilibrium, which collapses forming a halo, disk and bulge components. We here concentrate on what the authors have defined as Bulge, particles for which [R 2 + Z 2 ] < 2 kpc. As clearly shown in their Fig. 3, the predicted SFH for this inner-bulge shows multi components. Two of them seem to be predominant, resembling the modern data we discussed in §3 suggesting at least two components, but details on their chemistry were not given.
This model was then extended to a 3-D model and placed into the cosmological framework by Samland & Gerhard (2003) by building in a more realistic cosmological accretion history and angular momentum distribution for the gas. The accretion history was found to be closely linked to the SFH, which together with dynamics then determined the enrichment history of the different components of the Galaxy. According to this simulation the halo forms first, followed by the bulge and and a thick disk before finally the thin disk forms as the SFR goes down. Around redshift one, a bar forms and produces a triaxial bulge. The final bulge was found to be composed of at least two stellar subpopulations, namely: an early collapse population and a population that formed later in the bar. 4.4.3. Classical and pseudobulges formed from early major mergers. Athanassoula et al. (2017) presented a coupled chemical-kinematical-morphological approach, consisting of Nbody simulations, where gas and star formation are included, coupled with a chemical evolution code. The galaxy formation model is the result of a major merger of two spherical protogalaxies initially composed of dark matter and gas, taking place at 8-10 Gyr ago, with an initial metallicity of [Fe/H]=−1.0. Stars formed before the merging form a spheroidal centrally concentrated population, whereas those formed during the merger form mainly a thick disk and/or an extended stellar halo. This simulation is able to form a classical bulge with 9-12% of the mass, and a bar containing a boxy/peanut inner part.
This formation scenario is in agreement with the star formation history prediction by the so-called two-infall model of Chiappini et al. (1997). The latter authors argue for a break in star formation history as the cause for the [α/Fe] structure in the MW local disk (Fuhrmann 2011), soon after the first 1-2 Gyrs (by redshift z∼2). Interestingly, during the recent IAU Symposium 334 devoted to the Milky Way, Wetzel et al. (in preparation) presented FIRE simulations where a bimodal [α/Fe] distribution appeared to result from a major merger event early on in the history of the MW. Moreover, as shown by Minchev et al. (2013), radial migration alone is unable to produce the observed discontinuity (now seen beyond the local volume Anders et al. 2014;Hayden et al. 2015;Mikolaitis et al. 2014). In this scenario, the MW as observed today would have resulted from an early collision of two galaxies that gave the origin to the bulge. Thereafter the merging the evolution is secular, starting the formation of a thin disk.
Gas-rich two-spiral major mergers are the starting point also for Sauvaget et al. (2018). In their simulations, the mergers take place at high-redshifts, starting 12 Gyr ago, with final coalescence occurring at 1.5 to 3 Gyr from the initial time. Their motivation is the fact that two thirds of nearby spiral galaxies have pseudobulges or are bulgeless, whereas it has been shown that more than 50% of present-day galaxies have experienced a major merger in the past (Hammer et al. 2009). Their simulations show that most of the early gas-rich major mergers could lead to the formation of disks and bars, identified as pseudobulges.
These results indicate that the observed properties of the MW may be consistent with an early merger. Then the evolution after the merging would be secular, starting with the formation of a thin disk.
4.4.4. Chemodynamical models with predictions for individual chemical elements. Chemodynamical models with more detailed chemical information are discussed next. Nakasato & Nomoto (2003) have built a 3D hydrodynamical N-body model for the formation of the Galaxy, using the GRAPE (Sugimoto et al. 1990) smoothed particle hydrodynamics (SPH). This is combined with chemical evolution of the Bulge. In their model, most Bulge stars form during the subgalactic merger, and another group forms later in the inner-disk region. This again leads to two chemically different groups, in terms of [α/Fe] ratios, the first characterized by SN II and the second by SNIa enrichment. In their model, there is a merging of sub-clumps occurring at t=0.4-0.5 Gyr, causing a starburst, when 60% of Bulge stars are formed already at 0.5 Gyr.
Kobayashi & Nakasato (2011) used similar GRAPE simulations but with different prescriptions for feedback and chemical enrichment: the models now include not only SNII and SNIa, but also hypernovae with the goal to explain the large [Zn/Fe] chemical abundances observed in the MW halo most metal-poor stars. The initial conditions are the same as in Nakasato & Nomoto (2003). In the Kobayashi & Nakasato (2011) simulations, the Bulge again forms from the assembly of small galaxies early on (redshifts z > ∼ 3) and as a consequence, the bulge-particles are generally older than around 10 Gyr, and around 60% of them do show an important [α/Fe] enhancement ([O/Fe]>+0.3). As typical in the chemical evolution models discussed in previous Sections, due to the short timescale of the Bulge formation, this simulation also predicts an [α/Fe] plateau extending to above Solar metallicities ([Fe/H]∼+0.3 dex), creating an excess of very metal-rich stars reaching [Fe/H]∼1, compared with observations ( §2 and §3). Predictions for other abundance ratios, including iron-peak and odd-Z elements are presented in their models ( §3.3).
Chemodynamical models adapted from hydrodynamical models of elliptical galaxies, not in the cosmological framework, are presently available. Pipino et al. (2008Pipino et al. ( , 2010) developed a 1-D hydrodynamical model of the collapse of an spherically symmetric gas cloud, following the evolution of the gas density, momentum and internal energy, metallicity and [α/Fe] ratio (in particular oxygen). For a MW-like bulge model (the collapse of a gas cloud on a very short timescale due to fast gas consumption -below 0.5 Gyr -resulting on a stellar mass of 2×10 10 M , and an effective radius of ∼1 kpc) the authors find that a stellar abundance gradient can be created (especially in its central regions), which would have an impact on the final MDFs at different locations in the Bulge. The model also predicts integrated Lick-indices Hβ, M g2, Mgb and < F e > in good agreement with the integrated light derived Bulge values of Hβ (∼1.45), M g2 (∼0.25), Mg b (∼3.0) and < F e > (∼2.0) given in Puzia et al. (2002). This model is employed by Grieco et al. (2012) in order to provide a comparison of the expected Bulge chemical gradients from pure chemical evolution models with those of the chemodynamical model (their Table 2). Friaça & Terlevich (1998) presented a 1D chemodynamical evolution model for elliptical galaxies. This code was later modified to model the Bulge, available in 1D and 2D (Friaça & Barbuy 2017), where a single massive dark halo hosts baryonic gas which is turned into stars, by assuming different star formation rates, a baryonic mass of 2×10 9 M , corresponding to a small classical spheroid, and a dark halo mass MH = 1.3×10 10 M . The cooling function of the gas Λ(T ) is calculated consistently with the chemical abundances in the ISM, taking into account separately the abundances of the three main coolant elements, oxygen, carbon and iron. Feedback from heating, ionisation, mechanical pressure and chemical enrichment is taken into account.

Recent bulge formation models in the cosmological framework
In galaxy formation models within the ΛCDM scenario bulges would form via mergers. In the last years, high-redshift observations have contributed to better understand bulge formation.
From the observations discussed in §3.1, the conclusion that most of the Bulge stars are ∼10 Gyr old (still currently debated, but essentially confirmed by the most robust age indicators) would imply that at redshift z∼2 galaxies would be building their bulges. Indeed, these galaxies have shown to contain around 50% of gas (Daddi et al. 2010;Genzel et al. 2015;Tacconi et al. 2017), are mostly clumpy (e.g. Förster Schreiber et al. 2009) and show bulges in formation even in absence of any merging (Tacchella et al. 2015). The latter authors have shown evidence of mature bulges at redshifts of z∼2.2 by measurements of stellar mass and SFR surface-density distributions in massive star-forming galaxies. This study has shown that in most massive galaxies (with stellar masses ≥ 10 11 M ), although a high star-formation activity is sustained at large radii, mature bulges can reside in their central parts due to a quenching of the star formation already on timescales shorter than ∼1 Gyr. This important observational result provides insights on the mechanisms leading to bulge formation, and show that high central stellar densities are the result of gas-rich dissipative processes acting at very early epochs (and not by dissipationless mergers, nor by slow secular evolution). These observations go in the direction of models such as the one presented by Dekel & Burkert (2014), where violent disk instability in high-redshift gas rich galaxies also lead to bulges fully formed by gas dissipation.
In situ fast formation of bulges has also been suggested by several simulations of the so-called clumpy disks, where the main formation process relies on long lived clumps which migrate into the central parts of galaxies, early on (e.g. Noguchi 1999;Immeli et al. 2004a). Note that the formation of clumps in simulations of disks with typical masses and gas fractions as observed at z>1 does not lead to clump formation, unless turbulence is also powered by high rates of cosmological infall and stellar feedback (Bournaud 2016, and references therein). Elmegreen et al. (2008) and Bournaud et al. (2009) have presented simulations of clump migration resulting in the formation of both classical and pseudobulges. The authors show that long-lived clumps migrate and coalesce into the central bulge, but the final fate of the bulge formed by this channel will be very dependent on assumptions on feedback. Weak feedback simulations cannot prevent the clumps from becoming very massive, resulting in the formation of large bulge-to-disk B/D ratios and in the formation of ClBs (B/D∼0.2-0.3). On the other hand, when strong radiative feedback is assumed, the clump mass is regulated and they remain gas-rich, which can lead to rotating pseudobulges Bournaud 2016). This means that the results are still strongly dependent on sub-grid physics such as supernova feedback, thresholds in star formation, AGN feedback, among others (see Scannapieco et al. 2009Genel et al. 2012), but bulge growth through high-redshift disk instabilities leading to migrating dense gas clumps could have been a common process in the early Universe.
An encouraging result comes from the hydrodynamical cosmological simulations of two MW mass galaxies in which bars could form naturally within the ΛCDM scenario (Scannapieco & Athanassoula 2012; Martig et al. 2012). Ideally, one would like to study the bar formation process in simulations where the galaxy itself is growing and accreting mass both by infall of pristine gas and mergers. In the fully cosmological simulations of Scannapieco & Athanassoula (2012) for two MW-like galaxies one of them resulted in a system with a well-defined disk, bulge and halo components and a strong and long bar, whereas in the second case a smaller and weaker bar formed in an almost bulgeless galaxy. Also very important, these simulations show that bars can start forming early on, without w aiting for the galaxy to be fully formed before starting their own formation, as in idealised dynamical simulations discussed in §4.2 whose focus is the study of dynamical mechanisms, the role of resonances, and the formation of buckled bars.
Furthermore, bulges and elliptical galaxies, although similar in several aspects, do differ because in bulges the contribution of secular evolution during disk instability is important (Kormendy & Kennicutt 2004, and references therein). Smaller bulges of late-type galaxies should also form on longer timescales from bar-driven flows of disk gas towards the galaxy center. Secular evolution through disk instability can produce a pseudobulge via bar formation and buckling (e.g. Samland & Gerhard 2003;Athanassoula 2005;Obreja et al. 2013). These processes have also been modeled before (see §4.2) but have mainly focussed on observations in the nearby Universe, and on the more recent Bulge data (discussed in §3). The latter weak non-axisymmetric instabilities of nearby galaxies are substantially different from the violent clump instability process taking place in high-redshift galaxies, as discussed above. Interestingly, wide and deep IR surveys of galaxies in the 1-3 redshift range do not show barred disks but massive young clumps (e.g. Sheth et al. 2008;Simmons et al. 2014). The latter data confirm the lack of regular barred disks, giving support to the idea that these star forming clumps are bound to their host galaxies, as suggested already in O'Neil et al. (2000), Immeli et al. (2004b), Elmegreen et al. (2005Elmegreen et al. ( , 2007, and are a common property of the progenitors of nowadays massive spirals such as the MW and one of the main channels for an early (mostly in situ) bulge formation.
Higher resolution fully cosmological simulations forming bulges were presented by Guedes et al. (2011), Obreja et al. (2013), and Hopkins et al (2014) (see Brooks & Christensen 2016), and these usually provide a fraction of old-to-total components in the bulge of 30-50%. More recent chemodynamical models, in the cosmological framework, including chemical predictions have been presented by Grand et al. (2017Grand et al. ( , 2018, although in this case the focus is more on the formation mechanisms of disk galaxies. These very complex simulations are based on the magneto-hydrodynamical simulation code AREPO (Springel 2010), that solves fluid equations on a moving mesh. Each particle represents an single stellar population (SSP), where chemical evolution and stellar evolution are included, and gas enriched feedback from SNI, SNII and AGBs are considered (although here too, sub-grid physics related to feedback still prevent a more robust prediction for the mixing/distribution of the chemical elements along the several stellar generations). These models produce a red spheroidal bulge, but subdominant in many of the haloes. Moreover, in these simulations massive bulges are produced mainly as a result of major mergers.
In conclusion, the theoretical picture emerging is one in which different processes are at work during Bulge formation, including disk instabilities due to strong gas accretion, major and minor mergers, and secular processes. This complexity is confirmed by the multi-stellar population observed in the Bulge (as discussed in §3). Although the last few years have brought an enormous development in cosmological simulations of the formation and evolution of MW-like galaxies, a successful self-consistent fully cosmological chemodynamical evolutionary model of the MW-like bulges is still missing.

Conclusions
The Galactic Bulge hosts one of the oldest stellar populations in the Galaxy, carrying the imprints of the earliest Galaxy formation phases. Because in the Galactic Bulge we can measure kinematic and chemical properties for individual stars, unique constraints on bulge formation models can be obtained. In this review we summarize recent results on the structure, kinematics, and dynamics of the Bulge, and on the physical properties of its stellar populations, including their ages, metallicities and chemical abundances, considering elements produced through different nucleosynthetic channels. These results have established a composite nature of the Bulge where the spatial distributions, kinematics, and dynamics of Bulge stellar populations depend on their population parameters.
In order to provide a link between these observations and the formation of the Bulge, we also review possible bulge formation scenarios and models. At present, a successful fully self-consistent model for the Bulge in the cosmological framework is not yet available (see § 4). Different modeling approaches have, so far, focused on reproducing different aspects of the data, investigating several formation scenarios in parallel. Therefore, we first provide a list of the main observational facts that have emerged regarding kinematics and dynamics ( § 2) and stellar populations ( § 3), and then proceed to their implications for the several interpretations of how the Bulge was formed. We also indicate which of these resulsts are still debated: • Most of the stellar mass of the Bulge is barred (overall ∼ 90% based on RCG star samples). Bimodal magnitude distributions at high-latitudes and reconstruction of the stellar density from RCG show a B/P Bulge (X-shape), similar to other barred galaxies. Vertical density profiles are exponential above |z| >500 pc where they are most reliably measured. At 2 − 3 kpc along the long axis, the B/P Bulge transits into a planar long bar, similar to other barred galaxies. The RRL population, tracing ∼ 1% of the stellar mass, is mostly unbarred but a fraction of orbits may have been trapped by the bar. • Cylindrical rotation across all latitudes (|b| < 10 • ) is seen in the BRAVA, ARGOS and APOGEE surveys, characteristic of stellar bars. The cylindrical rotation is observed for most Bulge stars over a broad metallicity range. However, metal-poor and metalrich stars differ in their vertical velocity dispersion profiles, and only the metal-rich stellar component shows a kinematic vertex deviation. • Spectroscopic measurements show a broad MDF ranging from below 1/10 solar to supersolar metallicities. Metal-rich stars are relatively less numerous at high |b|. This result has to be further confirmed, and studies of innerthick-disk samples would be needed for a more conclusive comparison. The α-element enhancements with respect to Fe extend to higher metallicities, therefore the turnover or knee in the Bulge occurs at higher metallicities than in the thick disk. This result is however dependent on sample selections details. These differences are likely not due to IMF variations, because current data have shown that the IMF in the Bulge is similar to that in the solar-vicinity.
This above list of observational facts constitutes an impressive achievement of the Galactic community in the past decade. These data come together with the information obtained from high-redshift (z∼2) external bulges and thick disks (briefly summarised in §4). It is this body of data that must eventually be accounted for by galaxy formation models in the cosmological framework, which then should be able to predict chrono-chemical-kinematic maps of the MW Bulge/inner-disk region. Below we now list some of the more straightforward implications of the observational facts, based on the various modelling approaches reviewed in Section 4.
• Detailed dynamical models of RCG starcounts and kinematics in multiple fields in the inner Galaxy confirm the barred structure of the Bulge and place the corotation radius of the Bar at 6.1 ± 0.5 kpc. These models find that the Bulge and Bar region out to ∼ 5.3 kpc radius contains about 65% of the MW's total stellar mass, and the Bulge alone about 25%. The relative fraction of dark matter in the Bulge is found to be small, ∼15-20%, implying a low dark matter density in the Bulge and a core or mild cusp in the density of the dark matter halo. • Best-fit chemo-dynamical equilibrium models of distance-kinematics-metallicity data find that the most metal-rich populations in the Bulge and Bar are in a strong bar concentrated towards the Galactic plane, while the more metal-poor stars are in a weaker, thicker bar. The superposition of the different populations causes the apparent vertical metallicity gradient in the Bulge. Stars in the thick disk metallicity range are mostly found in a thick disk bar of vertical scale-length ∼ 500 pc, but a fraction of these stars is in a centrally concentrated component which could be a small classical bulge but may also have several other possible origins.
• Bars and B/P Bulges in N-body dynamical models of unstable, pre-assembled disks qualitatively explain many of the structural and dynamical properties of the Galactic B/P Bulge and long bar. Differences between metal-rich and metal-poor populations can largely be explained if the metal-poor stars are in a thicker and/or hotter disk when the bar and buckling instabilities set in. Chemodynamical collapse models leading to bar-unstable disks can quantitatively reproduce much of the chemo-kinematic data in the Bulge. • The oversolar abundance ratios with respect to iron in most Bulge stars more metalpoor than [Fe/H] < ∼ −0.5 can be explained by chemical evolution models where most of the Bulge was assembled on short timescales of around ∼2 Gyr, where the enrichment was dominated by SNII and SNIa did not have time to contribute much to the interstellar medium from which these stars were born. The oversolar [Zn/Fe] ratios in metal-poor Bulge stars suggest that hypernovae contributed to the early chemical enrichment of the Bulge. The presence of rapidly rotating massive stars might contribute to explain the intrinsic scatter observed in the neutron-capture abundance ratios of the oldest Bulge stars. • The chemical differences observed in the comparison of stars in the local thick disk with Bulge stars suggest that that thick disk and Bulge could be two genuinely different components, and the Bulge is not just a superposition of halo, thick disk and bar components. The caveat here is that still very few data are available for thick disk stars towards the inner regions of the MW. • Some chemical differences between Bulge and thin-disk stars at higher metallicities are useful constraints on stellar yields of intermediate mass-stars and SNIa, as well as on their metallicity dependency (also in massive stars). • Although the star formation timescale seems constrained by the data, the efficiency of the star formation is still poorly constrained. Pure chemical evolution models seem to require larger star formation efficiencies than those assumed in chemodynamical (hydrodynamical) models; this parameter affects not only the predictions for the knee (Table 3), but also for the MDFs ( §3.2). • The dominant old ages of Bulge stars is compatible with models where the Bulge is assembled on short timescales, which includes chemical evolution models, as well as simulations where the Bulge forms from strong gas accretion, early mergers, or gasrich merging clumps from the thick disk migrating towards the inner regions. The fact that most of the Bulge mass is barred would be accounted for by models where the Bulge stars formed at high redshifts, and are able to form a bar sufficiently early on depending on the assumptions on feedback in the simulations.
Ongoing and future survey data will allow us to test and refine our current understanding. On the dynamical side, proper motions from VVV and Gaia are expected to greatly improve mass models of the Bulge. Stellar ages obtained by combining spectroscopic and asteroseismological data will be revolutionary in allowing us to dissect the Bulge popula-tions as a function of their age. Combined with the metallicity structure and abundances this will give very strong constraints on formation models, which themselves are expected to improve greatly in the coming years.