Was the Milky Way Bulge Formed from the Buckling Disk Instability, Hierarchical Collapse, Accretion of Clumps, or All of the Above?

Abstract

The assembly of the Milky Way bulge is an old topic in astronomy, one now in a period of renewed and rapid development. That is due to tremendous advances in observations of bulge stars, motivating observations of both local and high-redshift galaxies, and increasingly sophisticated simulations. The dominant scenario for bulge formation is that of the Milky Way as a nearly pure disk galaxy, with the inner disk having formed a bar and buckled. This can potentially explain virtually all bulge stars with [Fe/H] ≳ −1.0, which comprise 95% of the stellar population. The evidence is the incredible success in N-body models of this type in making non-trivial, non-generic predictions, such as the rotation curve and velocity dispersion measured from radial velocities, and the spatial morphologies of the peanut/X-shape and the long bar. The classical bulge scenario, whereby the bulge formed from early dissipative collapse and mergers, remains viable for stars with [Fe/H] ≲ −1.0 and potentially a minority of the other stars. A classical bulge is expected from Λ-CDM cosmological simulations, can accentuate the properties of an existing bar in a hybrid system, and is most consistent with the bulge abundance trends such as [Mg/Fe], which are elevated relative to both the thin and thick disks. Finally, the clumpy-galaxy scenario is considered, as it is the correct description of most Milky Way precursors given observations of high-redshift galaxies. Simulations predict that these star-forming clumps will sometimes migrate to the centres of galaxies where they may form a bulge, and galaxies often include a bulge clump as well. They will possibly form a bar with properties consistent with those of the Milky Way, such as the exponential profile and metallicity gradient. Given the relative successes of these scenarios, the Milky Way bulge is plausibly of composite origin, with a classical bulge and/or inner halo numerically dominant for stars with [Fe/H] ≲ −1.0, a buckling thick disk for stars with − 1.0 ≲ [Fe/H]] ≲ -0.50 perhaps descended from the clumpy-galaxy phase, and a buckling thin disk for stars with [Fe/H] ≳ −0.50. Overlaps from these scenarios are uncertain throughout.

1 INTRODUCTION

Formation scenarios for the Galactic bulge (and bulges in general) have been around for some time.

Eggen, Lynden-Bell, & Sandage (1962) suggested the scenario of early, rapid, monolithic formation by dissipative collapse, this is now referred to as the ‘classical bulge’ scenario. The hierarchical merging of smaller objects, typical of the earliest phases of galaxy formation as predicted by λ-CDM cosmological simulations (Tumlinson 2010; Kobayashi & Nakasato 2011), is now often included as part of the definition of ‘classical bulges’ (Athanassoula 2005). It is also known that bars can form from the dynamical evolution of a disk galaxy via the ‘buckling instability’ scenario, where the bulge would be predominantly composed of disk stars now on bar orbits (Miller 1978; Miller & Smith 1979; Hohl & Zang 1979), for example, the ×1 orbital family (Skokos, Patsis, & Athanassoula 2002). Finally, the migration of star-forming clumps towards the centres of disk galaxies due to dynamical friction emerges naturally from simulations, this is the ‘clump-origin’ scenario (Noguchi 1998, 1999).

Given the wide availability of plausible models, it is no surprise that it has long been acknowledged that the formation of the bulge may be a composite process. For example, in a previous review (of bulges in general), Combes (2000) wrote:

A fraction of bulges could have formed early (at first collapse); then secular dynamical evolution enrich them; in parallel, according to environment, accretion and minor mergers contribute to raise their mass.

The first and third process mentioned are the two components of the classical bulge scenario, the second is the buckling instability scenario. Relative contributions of the different scenarios of bulge formation were not assigned precise bounds, given that this is an intrinsically difficult problem.

Astronomers need an accurate, precise, and thorough census of the current Milky Way bulge before we can seriously consider disentangling its history.

That census is now slowly, but surely, becoming available. These include ground-based photometric surveys such as OGLE-IV in the optical (Udalski, Szymański, & Szymański 2015), VVV in the near-infrared (Saito et al. 2012), and Wide-field Infrared Survey Explorer (WISE) in the mid-infrared (Wright et al. 2010); Deep, multi-wavelength Hubble Space Telescope (HST) photometry which measures the main sequence and thus constrains ages (Brown et al. 2010). At the spectroscopic end, surveys such as BRAVA (Kunder et al. 2012), ARGOS (Freeman et al. 2013), Gaia-ESO (Rojas-Arriagada et al. 2014a), EMBLA (Howes et al. 2014), APOGEE (Ness et al. 2016), BRAVA-RR (Kunder et al. 2016), and GIBS (Zoccali et al. 2017), collectively provide kinematics and sometimes detailed chemistries for tens of thousands of bulge stars.

Concurrent with this, there have been vast improvements in the theory and simulations as well. For example, the N-body simulation of Athanassoula, Rodionov, & Prantzos (2017) contains 17.5 million particles, whereas that of Miller & Smith (1979) had 96,017 particles. The level of detail and sophistication of current models allow more robust analysis. Further, more questions are being asked of models, such as how a classical bulge will evolve when embedded within a massive disk (Saha, Martinez-Valpuesta, & Gerhard 2012), where younger stars might be distributed within a bar (Debattista et al. 2017), or how the separation between the two arms of an X-shape will appear as a function of direction (Nataf et al. 2015). These questions, largely motivated by the observations, should yield more predictive and discriminatory power in evaluating models, with the caveat that it is scientifically misguided to expect a perfect quantitative match between a simulation and an actual galaxy.

In this review, three scenarios for bulge formation are discussed as well as the current evidence in their favour.

2 THE MILKY WAY AS A NEARLY PURE DISK GALAXY

This is the scenario where the Milky Way bulge is largely or nearly entirely a bar. This bar would have first formed from a disk, and then extended vertically by one or more buckling instability episodes. This scenario arguably has the most evidence in its favour.

The evidence presented is that

• • pure disk galaxies are observed to exist in the local universe, and the bar buckling process occurs naturally in N-body simulations;

• • the radial velocity measurements from large spectroscopic surveys are consistent with the theoretical predictions from buckling disk models;

• • the peanut/X-shaped morphology for the distribution of bulge stars is a prediction of these models, and is constrained to represent a large fraction of the bulge mass;

• • the morphology of the long bar is another specific, precise prediction of these models;

• • an old argument against this theory, that of the metallicity gradient, has been shown to be invalid as it can also be produced by the buckling instability model.

2.1. Pure disk galaxies exist, and simulations show that they can naturally evolve to be barred galaxies

Pure disk galaxies exist. Kormendy et al. (2010) obtained HST photometry of six nearby galaxies and measured their surface brightness profiles. They found an upperbound on the total (classical bulge + pseudobulge) stellar mass fraction of ~3% of the total disk mass. They also took an inventory of galaxies within 8 Mpc with $v_{\text{circ}} > 150\,\, \rm {km\, s}^{-1}$ , and found that 11 of 19 showed no evidence for a classical bulge, and four of the remaining eight may contain classical bulges contributing 5–12% of the stellar mass. It is the case, however surprisingly, that pure disk galaxies are common in the low-redshift universe, and thus it is plausible for the Milky Way to be one as well.

Simulations of disk galaxies consistently show that bars and subsequently buckling bars can occur spontaneously in disk galaxies (e.g. Combes & Sanders 1981; Raha et al. 1991; Athanassoula 2005; Martinez-Valpuesta, Shlosman, & Heller 2006; Quillen et al. 2014; Spinoso et al. 2017). The rate of evolution of the bar is sensitive to various variables, such as the ratio of the disk to halo mass (Combes & Sanders 1981). The bar may have a north–south asymmetry at first, but this rapidly dissipates (Raha et al. 1991). Overtime, the bars tend to grow in vertical extent, for example, via multiple, recurrent buckling episodes (Martinez-Valpuesta et al. 2006, see Figure 3).

Separately, a bar does form in the gas-rich cosmological “zoom-in” simulation of Spinoso et al. (2017). The bar's properties are a sensitive function of redshift and related to both gas transport and AGN feedback, suggesting the need for diligence in associating bar properties to Galactic archaeology.

The combination of these two facts, that pure disk galaxies are common in the local universe, and that simulations predict that pure disk galaxies can spontaneously evolve to have bars, render it a plausible model for the Milky Way’s bulge as well. The necessary initial conditions—a disk galaxy—are ubiquitous, and the necessary evolution is natural.

2.2. Radial velocity measurements from large spectroscopic surveys

The Milky Way is a barred galaxy (Dwek et al. 1995), however it has not been clear until recently how quantitatively dominant the bar is relative to the hypothetical classical bulge contribution. The advent of large spectroscopic surveys has allowed us to disentangle the bulge into plausible subcomponents.

Shen et al. (2010) fit a suite of N-body models to bulge radial velocity data from the BRAVA survey, specifically the mean radial velocity (and thus rotation curve) and radial velocity dispersion as a function of longitude and latitude. They recover a best-fit viewing angle, between the bar’s major axis and the line of sight between the Sun and the Galactic centre, of α_Bar = 20°, consistent with more recent measurements. More significantly, they show that models with significant classical bulges, with parameters defined to lie on the fundamental plane of ellipticals and bulges (Kormendy et al. 2009), are not well fit by the data. They constrain the classical bulge mass of the Milky Way to be no more than 8% of the total mass of the disk.

Ness et al. (2013) analysed the data from the larger ARGOS survey, specifically 17 400 red giants with best-fit parameters including R _GC ⩽ 3.5 Kpc, which were spread over 30° of longitude and 20° of latitude. They find several features in the data consistent with N-body expectations for disk origin to the bulge. Bulge stars located between 5 and 10° from the plane ( $0.7\, \rm {Kpc} \lesssim |{\it z}| \lesssim 1.5\, \rm {Kpc}$ ) are cylindrically rotating, with the exception of the 5% of stars with [Fe/H] ⩽−1.0. They suggest that those stars are the inner extensions of the halo and the metal-weak thick disk. It is also the case that the thin disk includes very, very few stars with [Fe/H] ⩽−0.50 (Delgado Mena et al. 2017), so it was always an implausible origin source for the metal-poor component of the bulge.

Di Matteo et al. (2015) compared the ARGOS data to three different N-body models and found that the bulge could not be explained as having resulted from the evolution of a pure thin disk galaxy, and thus a joint origin with the thick disk is likely required as well. This may appear as requiring too many free parameters from models, but this is a parameter that exists in nature, the Milky Way is not a pure thin disk galaxy, it also has a thick disk. Di Matteo et al. (2015) show that if the bulge was purely due to a thin disk buckling event, the curve of mean radial velocity versus longitude would decrease in amplitude, and the velocity dispersions would decrease, with increasing metallicity. The decrease in the amplitude of the rotation curve is seen for stars satisfying [Fe/H] ≳ −0.50, but not for lower metallicities. The velocity dispersions do decrease, but not as quickly as predicted from a pure thin disk model. Di Matteo et al. (2015) suggest that the solution to this issue is that bulge stars with [Fe/H] ≲ −0.50 formed from the thick disk, rather than from the more metal-poor component of the thin disk.

Zasowski et al. (2016) investigated the skewness (third moment) and kurtosis (fourth moments) of the radial velocity distribution functions in APOGEE data towards the bulge. The correlation between skewness and mean velocity (first moment), a known diagnostic of bars, is only observed for the more metal-rich fraction of the stars, defined in that paper as stars with [Fe/H] ≳ −0.40. The data have a flat kurtotsis, Kurt (V) ≈0, consistent with that expected from models. The data have slightly higher velocity dispersion and slightly lower skewness than expected from N-body models of a simple single disk galaxy undergoing buckling.

The recent review of Di Matteo (2016) goes over these lines of evidence in greater detail than possible here, and concludes that bulge stars with [Fe/H] ≳ −1.0 formed from buckling thick and thin disks.

2.3. The peanut/X-shaped spatial distribution of bulge stars

The sightline-dependent bifurcation of the apparent magnitude distribution of red clump stars (comprising ~99% of bulge horizontal branch stars, see Nataf et al. 2013) was first reported by Rattenbury et al. (2007). They speculated that it might be due to a distinct population lying at a different distance. Nataf et al. (2010) and McWilliam & Zoccali (2010) showed that the double peak was spread across a large swath of the bulge, at large separations (|b| ≳ 5°) from the plane. McWilliam & Zoccali (2010) argued convincingly, that the feature had to be due to an X-shaped bulge (a very strong peanut shape), a feature of strong bars in N-body models.

Subsequent analyses confirmed the claim. Ness et al. (2012) found that the bright and faint red clumps had the same difference in mean radial velocity as expected from N-body models. Li & Shen (2012) showed that the N-body model of Shen et al. (2010), already demonstrated to be a good match to radial velocity observations, also includes an X-shape. The relative brightness, and the dependence of the relative number counts on direction, was qualitatively similar to that seen in the data. Vásquez et al. (2013) showed that the proper motion distributions were consistent. Nataf et al. (2015) showed that precision measurements of the brightness difference, relative number counts, and total number density of stars in the data from OGLE-III (Udalski et al. 2008) were qualitatively matched by N-body models from Athanassoula (2003) and Shen et al. (2010). An unambiguous image of the X-shaped bulge in integrated mid-infrared photometry of the Milky Way is shown in Figure 1, and the development of one from an N-body model is shown in Figure 2.

Figure 1. Figure 1 from Ness & Lang (2016), the X-shaped bulge is unambiguous in the integrated mid-infrared photometry of the Milky Way.

Figure 2. Figure 3 from Martinez-Valpuesta et al. (2006). An X-shaped bulge is a natural outcome of N-body models of disk galaxies undergoing the buckling instability.

Precision modelling of the photometric data by Wegg & Gerhard (2013) has been followed by detailed dynamical modelling by Portail, Wegg, & Gerhard (2015a) and Portail et al. (2015b, 2017). A suite of N-body models coarsely consistent with Milky Way constraints was adjusted to be more consistent using the made-to-measure method (Syer & Tremaine 1996), whereby weights of the particles in the model are shifted up or down to match the observational data. The end models always contain a significant X-shape, comprising 40–50% of the stellar mass of the bulge, but only dominant for stars with [Fe/H] ≳ −0.50. That is not a repetition of the results mentioned in the previous subsection, as Portail et al. (2015a, 2015b, 2017) also constrained their N-body models to match the spectroscopic data from the APOGEE survey, which goes much closer to the plane (Ness et al. 2016), and the global photometric parameters of the bulge which trace the number density and distance distribution function of stars (Wegg & Gerhard 2013).

Reiterating, this is an impressive series of observational tests for models to have passed. First, the alternative hypotheses (RGB bump, metallicity differences, etc.) do not work. Second, N-body model predictions qualitatively match the radial velocity offsets, the proper motion measurements, the mean brightness difference, and the dependence of the relative and total number of counts as a function of direction. This picture survives, and is in fact dominant for stars with [Fe/H] ≳ −0.50, when models are required to match observations from several different surveys.

2.4. The morphology of the long bar is another specific, precise prediction of these models

An issue of Galactic modelling in the last two decades is that of the long bar of the Milky Way (Benjamin et al. 2005; Cabrera-Lavers et al. 2007, 2008). Evidence has accumulated from multiple investigations of a long, in-plane bar, with a half-length of ~4 Kpc, and an orientation angle of α_Bar = 45°. This represented a challenge to Galactic structure studies, as the orientation angle is far larger than the value found for the triaxial bulge towards the inner few Kpc, and the morphology is not consistent with the triaxial ellipsoid models used to fit for the bulge (Dwek et al. 1995; Stanek et al. 1997). It appears as a double-barred system where the two bars are of comparable length and are not aligned, which is an unstable configuration (Athanassoula 2012).

A suggestion to this resolution was independently proposed by Martinez-Valpuesta & Gerhard (2011) and Athanassoula (2012). Their idea was the analytical triaxial ellipsoid models then widely used to model the bulge were only (somewhat) suitable for the peanut/X-shaped component, whereas bars in simulations often have long and thin extensions. This hypothesis made a specific prediction, that the difference in measured angle Δα_Bar, then around 20°, was due to observational errors. A smaller difference might remain if the ends of the bar develop interactions with spiral arms Martinez-Valpuesta & Gerhard (2011).

This prediction was confirmed by Wegg, Gerhard, & Portail (2015), who combined photometry from Spitzer-GLIMPSE (Benjamin et al. 2005), 2MASS (Skrutskie et al. 2006), VVV (Inoue & Saitoh 2012), and UKIDSS (Lucas et al. 2008) to make a global map of the bulge. The long bar had an angle of α_Bar = 28–33°, consistent with the values of α_Bar = 27 ± 2° and α_Bar ≈ 29° then recently measured by Wegg & Gerhard (2013) and Cao et al. (2013), respectively, for the inner peanut/X-shaped component.

Thus, the issue of the long bar, which was previously a challenge to Galactic structure models, ended up being a triumph. The long bar is not only not surprising, but expected from simulations of buckling disk galaxies. Models predicted that the difference in orientation angle should shrink with better analysis, and it did.

2.5. The metallicity gradient can also be reproduced by buckling disk

The Galactic bulge has a vertical metallicity gradient, with stars further from the plane having a lower mean metallicity. Tiede, Frogel, & Terndrup (1995) used near-IR photometry and estimated a gradient in ∇[Fe/H] =−0.06 ± 0.03 dex deg⁻¹ or −0.43 ± 0.21 dex kpc⁻¹ between b = −3° and b = −12°. The metallicity gradient was independently and concurrently confirmed by Minniti et al. (1995), and since confirmed numerous times with spectroscopic data (e.g. Zoccali et al. 2008; Johnson et al. 2011).

This was argued by Minniti et al. (1995) to be evidence for a dissipative collapse, as models of spheroid formation predicted a radial metallicity gradient (Carlberg 1984). However, Minniti et al. (1995) acknowledged that ‘the alternative interpretation that the gradient itself is caused by the mixing of different components in the inner Galaxy cannot be ruled out’.

That second point, that the bulge has different components with both different metallicity distribution functions and scale heights, is very much the case, as should be clear from the literature evidence summarised in this review. Further, it also turns out that a metallicity gradient for the bulge can emerge from a pure disk galaxy. Martinez-Valpuesta & Gerhard (2013) evolved an N-body model of a pure disk galaxy which included an initial radial metallicity gradient. They chose [M/H](R) = +0.60 − 0.40(R) Kpc⁻¹, where ‘R’ is the galactocentric radius of particles at the start of the simulation. Given that particles from initial radii are scattered to different orbits, a metallicity gradient emerges which is qualitatively consistent with that observed for the whole bulge—there is impressive agreement with the global photometric metallicity maps from Gonzalez et al. (2013), which are derived from the morphology of the red giant branch. That said, it is worth noting that the metallicity gradient required by Martinez-Valpuesta & Gerhard (2013) is extremely large.

This result is not altogether surprising as it has been known for a while that dynamical mixing processes do not erase gradients. For example, White (1980) demonstrated that simulations of mergers were expected to reduce, but not erase, gradients in metallicity. Sánchez et al. (2014) find that the metallicity gradients in interacting galaxies, measured in terms of oxygen abundance of H ii regions versus effective radius, is reduced by ~1/2 relative to non-interacting galaxies.

The metallicity gradient of the bulge, first largely argued to be evidence for dissipative collapse and thus evidence against a disk–galaxy origin, turns out not to be discriminating. A metallicity gradient can be matched by models of dissipative collapse, pure disk galaxies undergoing the buckling instability, a combination of the two, or as we will see, the clump-origin bulge scenario (Inoue & Saitoh 2012). The existence of a metallicity gradient, by itself, yields no significant constraints bulge formation.

2.6. Caveat against the bar scenario: Gas fractions and the challenge of initial conditions

The properties and evolution of bars in N-body simulations of pure exponential disks are a well-researched subject. However, real galaxies (such as the Milky Way) contain gas, and generally contained more gas in the past.

Athanassoula, Machado, & Rodionov (2013) investigated the relative predicted properties of bars in N-body simulations with and without gas. They found that the gas-rich galaxies remain axisymmetric for longer. When they do develop bars, they do so at a slower rate, and end up much weaker. Their Figure 7 shows that bar strengths in gas-rich galaxies are predicted to end up lower even after the gas has been depleted. In contrast, as discussed in this section, the Milky Way has a very strong bar.

Indeed, the bar fraction is considerably lower at high redshift. Melvin et al. (2014) find that it drops to 11% at z = 1, corresponding to a lookback time of 7.8 Gyr. Simmons et al. (2014) also estimate a bar fraction of 11%, in the redshift range 0.50 < z < 2.0. By that time, most Milky Way bulge stars were already formed (Clarkson et al. 2011; Bensby et al. 2017) and thus they already have a kinematic distribution.

The combination of these issues should give pause to the notion that a bar + buckling instability from a pure exponential disk is independently sufficient to explain the inner Milky Way’s dynamics. Separately from the issue of gas weakening bars, subsequent sections of this review will discuss how a minor classical bulge can actually strengthen the bar, and why the Milky Way is expected to be a former ‘clumpy galaxy’ which has implications for bulge formation.

3 THE CLASSICAL BULGE SCENARIO

The classical bulge scenario is one where the bulge is formed in a combination of early, dissipative collapse, and accretion of objects via minor or major mergers.

The evidence we present is that

• • it is expected from theory;

• • the metal-poor bulge stars are kinematically consistent with a classical bulge behaviour;

• • the theoretical interaction between classical bulges and buckling disks is consistent with observations;

• • the bulge trends in the α-elements are not consistent with those of the disk.

3.1. A classical bulge is expected from theory

The most obvious advantages of this scenario are that it is predicted by straightforward theory of gravitational, dissipative collapse (Eggen et al. 1962), and hierarchical clustering in a cold dark matter universe (White & Rees 1978; Kauffmann, White, & Guiderdoni 1993).

White & Rees (1978) proposed that most of the matter in the universe condensed early into small ‘dark’ objects. They suggested a model with Ω_m = 0.20 and dark matter making up 80% of the matter, impressively similar to the modern values (Planck Collaboration et al. 2014). Within this picture, the ‘pure disk galaxy’ that has particles and no gas is not a viable initial condition, as the universe and star formation begin with larger number of small haloes that coalesce and accrete additional small haloes.

As far as current, more up-to-date models are concerned, a classical bulge can be considered a requirement. Kobayashi & Nakasato (2011) simulate 150 galaxies from various cosmologically motivated initial conditions. They formed disk structures in 48 of their galaxies, including five galaxies that had masses comparable to the Milky Way. Though some of their simulated galaxies had small classical bulges, none completely lacked a classical bulge. It is a concern that many local galaxies are found not to have classical bulges locally (Kormendy et al. 2010), but one can argue against that by saying that information is lost when looking at those galaxies in integrated light. That cannot be argued for the Milky Way.

3.2. The kinematics of metal-poor bulge stars

There has been much written in this review of the spectacular consistency between N-body models of buckling disks and the dynamics of bulge stars with [Fe/H] ≳ −1.0 or ≳ −0.50, depending on the study. These consistencies are not found for the metal-poor bulge stars. For example, the Galactic bar is either null (Dékány et al. 2013) or weak (Pietrukowicz et al. 2012, 2015) in the RR Lyrae stars, variable stars of standardisable distance that can be used to directly probe the metal-poor bulge. Dékány et al. (2013) investigated 7 663 fundamental-mode RR Lyrae in I and K _s bands, and found that they do not trace a strong bar, but rather a more spheroidal and centrally concentrated distribution. No correlation is observed between dereddened distance and longitude, a correlation detected at high significance in the more metal-rich red clump stars which trace a bar (Stanek et al. 1997).

Pietrukowicz et al. (2012) and Pietrukowicz et al. (2015) obtained different results in their analysis of 16 836 and 27 258 RR Lyrae in the OGLE-III and OGLE-IV surveys, respectively. They do find a bar, but it is kinematically hotter than the bar measured in red clump stars, and does not show a peanut/X-shape at large separations from the plane.

Analysis of RR Lyrae kinematics further confirms these findings. Kunder et al. (2016) analyse spectroscopic data for 947 RR Lyrae as part of the ongoing BRAVA-RR survey. These RR Lyrae, measured towards |l| ≲ 4° and −6° ≲ b ≲ −3° show higher velocity dispersions and weaker rotation than the metal-rich M-giants studied as part of the BRAVA survey. The velocity dispersion is ~15% higher, and the rotation is null or negligible. When the RR Lyrae are split into two metallicity bins with [Fe/H]=−0.75 marking the bifurcation point, no difference in rotation is measured, though the metal-poor RR Lyrae do have a higher velocity dispersion. The RR Lyrae radial velocity measurements are shown in Figure 3.

Figure 3. Left panel of Figure 2 from Kunder et al. (2016). RR Lyrae stars show null or negligible Galactic rotation, as well as very high velocity dispersion, in contrast to the majority of bulge stars.

A further issue with the RR Lyrae is the very fact that they are RR Lyrae. The bulge RR Lyrae are measured to have a mean metallicity of [Fe/H]≈−1.0, whereas metallicities of [Fe/H]≈−1.60 are more typical of RR Lyrae in the globular clusters and in the halo (Dotter et al. 2010). More metal-rich stars have a higher turnoff mass at fixed age, and thus for these horizontal branch stars to be on the instability strip at a higher metallicity, they need to have a lower turnoff mass by other means, and thus likely a greater age. Lee (1992) estimated that the bulge RR Lyrae had to be Δt ~ 1.3 ± 0.30 Gyr older than halo stars, for stars formed with metallicities − 1.5 ≲ [Fe/H] ≲ -1.0. A bulge which is older than the halo is more consistent with a classical bulge then say, a buckling disk origin.

This argument also applies to bulge globular clusters. For example, NGC 6522, has a blue horizontal branch with a mean metallicity of [Fe/H]=−1.0 (Barbuy et al. 2009), necessitating an extremely old age. This is a consistent pattern of bulge globular clusters (Barbuy et al. 2006, 2007). This information, combined with that of the RR Lyrae, strongly suggests that bulge stars of metallicity [Fe/H]=−1.0 are the oldest or among the oldest stellar populations in the Galaxy, and thus were around prior to there being a massive disk that could form a bar.

3.3. The theoretical interaction between classical bulges and buckling disks is consistent with observations

There is relatively sparse research on the theoretically predicted interaction of a central concentration of stars on a buckling disk, but what research is available now turns out to be consistent with the data. This is true both of the predictions of the classical bulge behaviour and the bar behaviour.

Athanassoula & Misiriotis (2002) studied three different N-body models of disk galaxies with variable initial central concentrations. The disk mass and disk-to-halo mass ratios are fixed. Two models have a non-centrally concentrated halo, one of those two also has a bulge. A third model has a centrally concentrated halo, but no bulge. Within this suite of models, the galaxy with a non-centrally concentrated halo and no bulge ends up forming the weakest bar, it does not develop cylindrical rotation, and its boxy shape does not evolve to an X-shape. Interestingly, the model with a centrally concentrated halo but no bulge ended up with the strongest bar. This is suggestive that greater central concentration was needed given the exceptionally strong bar of the Milky Way, though in and of itself, it is not sufficient as more models would be needed for a more convincing picture.

Conversely, Saha et al. (2012), Saha & Gerhard (2013), and Saha, Gerhard, & Martinez-Valpuesta (2016) have investigated how a joint classical bulge/buckling disk origin impacts the development of the bulge by means of N-body models. Saha et al. (2012) find that a small (7% of the total disk mass) classical bulge can pick up angular momentum from the larger rotating bar, and thus even develop into a triaxial object with cylindrical rotation. Saha & Gerhard (2013) show that the composite bulge always ends up rotating cylindrically, but may have deviations from cylindrical rotation at specific moments in its evolution. Curiously, they find that the final size of the composite bulges are reduced if the initial classical bulge has its own angular momentum. Saha et al. (2016) extend their prior results to show that even massive classical bulges can pick up as much specific angular momentum as low-mass classical bulges, but that the resulting rotation is non-cylindrical. All composite systems eventually form a boxy/peanut bulge.

Pérez-Villegas, Portail, & Gerhard (2017) model the evolution of a disk galaxy with a halo via an N-body model and find that the properties of the RR Lyrae are consistent with that of an inner halo, specifically their number density, their slow rotation, the lack of a peanut/X-shape in their spatial distribution, and their higher velocity dispersion. Within their model, only 12% of RR Lyrae end trapped on bar-like orbits, which can be considered as a prediction that the fraction will be very small once proper motions are available and once astronomers can compute orbits for these stars.

It is been established in the previous sections of this review that observations and analysis thereof have ruled out a predominantly classical bulge for the Milky Way bulge. However, the N-body models are also consistent with either a centrally concentrated inner halo, or an additional classical bulge in addition to that. Such a feature could help explain the Milky Way’s very strong bar as well as the behaviour of low-metallicity bulge stars. Further analysis and more data of metal-poor bulge stars are needed.

3.4. The bulge trends in the α-elements are not consistent with those of the thin and thick disks

The abundances of the α-elements relative to iron, [α/Fe], are a historic argument for a distinct formation to the bulge and the disk. That is because, the [α/Fe] abundance ratios trace the efficiency of star formation and possibly the initial mass function of a stellar population (Matteucci & Greggio 1986). A higher level of [α/Fe] can be due to a lower contribution to the chemical enrichment of interstellar gas from type Ia SNe (which take longer to form), and thus star formation would need to be faster.

The landmark study of McWilliam & Rich (1994) found that the trends for magnesium and titanium were enhanced by ≈0.30 dex relative to the solar neighbourhood trends over the full range of [Fe/H], whereas the trends for calcium and silicon were consistent with those of the disk. These different ratios suggested a different origin, and McWilliam & Rich (1994) said it may reflect a common enrichment process between bulges and ellipticals. Fulbright, McWilliam, & Rich (2007) compiled a more sophisticated analysis of a higher resolution, higher signal-to-noise sample. They found that the bulge has a magnesium trend elevated by ~0.30 dex relative to the disk, whereas the trends of oxygen, silicon, calcium, and titanium are slightly elevated relative to the disk. The bulge also has higher aluminum abundances at fixed iron abundance. Fulbright et al. (2007) conclude that the relative abundance offsets between the bulge and the disk are inconsistent with models where the bulge forms from the buckling of the disk. Separately, Fulbright et al. (2007) found that the metal-poor bulge stars also show higher mean abundances of silicon, calcium, and titanium than the halo. They concluded that the metal-poor bulge stars could not have formed from gas with the present-day halo composition.

The discussion shifted with the work of Alves-Brito et al. (2010). That investigation used high resolution of optical spectra of 25 bulge giants and 55 comparison giants and analysed abundances in a homogeneous manner to minimise systematic offsets. Their results were that metal-poor bulge stars ([Fe/H] ≲ −0.50) have the same abundances as the thick disk, and more metal-rich stars have the same abundances as the thin disk.

The relative abundance offsets between the bulge and the disk have not converged to zero as more data have come in. Bensby et al. (2017) compiled what is among the best datasets of bulge abundances, as they have high-resolution, high signal-to-noise abundances for 90 bulge stars located on the main-sequence turnoff and subgiant branch, analysed using the same methods as their comparison disk sample, which is composed of stars in the solar neighbourhood. Some elements trace the same abundance trends as the disk, but magnesium, titanium, and aluminum do not. They typically trace the upper end of the larger distributions spanned by the thick and thin disks. The abundance trend of [Mg/Fe] vs. [Fe/H] is shown in Figure 4.

Figure 4. From Figure 21 of Bensby et al. (2017). The [Mg/Fe] vs. [Fe/H] abundance trend for bulge stars (coloured points) are shown superimposed on the disk trends (grey points). The magnesium abundances are elevated with respect to the disk abundances at all [Fe/H] values.

Relative abundance offsets between the bulge and the disk remain, which is a challenge to models of the Galaxy where the bulge is simply due to buckling thin and thick disks. As the abundance ratios are higher in magnesium in particular, suggesting a more rapid star formation, the suggestion is that the bulge formed faster than the disk, characteristic of an early, dissipative collapse.

4 THE CLUMP-ORIGIN BULGE SCENARIO

Before proceeding, a brief description of star-forming clumps will be given, though the description is itself a matter of active research. Whereas ‘smooth exponential disks’ are decent approximations to local disk galaxies and exact descriptions of many N-body models, star-forming galaxies at high redshift are generally clumpy and gas rich (Cowie, Hu, & Songaila 1995; van den Bergh et al. 1996; Genzel et al. 2006), with off-centre clumps accounting for 7% of the stellar mass and 20% of the star formation in massive, star-forming galaxies (Wuyts et al. 2012).

The clumps are regions of excess star formation seen in the disks and proto-disks of high-redshift galaxies. Wisnioski et al. (2012) studied the properties of eight clumps in three redshift z ~ 1.3 observed as part of the WiggleZ Dark Energy Survey. They had an average size of 1.5 Kpc and average Jeans mass of 4.2 × 10⁹M_⊙, and accounted for roughly half the stellar mass of the disks.

Within this scenario, motivated by predictions from simulations (Noguchi 1998, 1999), some of these clumps will migrate to the centre of their galaxies due to dynamical friction, and thus form a bulge.

The evidence presented here is that

• • high-redshift galaxies largely appear as clumpy galaxies, and thus this is the plausible set of ‘initial’ conditions for the Milky Way;

• • the simulations succeed at predicting many of the observations.

This is less evidence than for the other two scenarios, but that is plausibly simply due to this being less researched topic. It is hoped that there will be further research testing whether or not the Milky Way bulge may be a clump-origin bulge. Zoccali et al. (2014) also discussed the issue. They pointed out that the mean age of bulge stars corresponds to an epoch of gas-rich disks, with gas fractions sometimes exceeding 50% (Tacconi et al. 2010; Daddi et al. 2010), which are more consistent with simulations of clumpy galaxies than the usually gas-free N-body simulations of bar formation in disks.

Readers interested in a more thorough review of bulge growth in high-redshift galaxies are referred to the excellent review by Bournaud (2016), in particular, Section 3, ‘Mechanisms of bulge growth through high-redshift disk instabilities’.

4.1. This is what high-redshift galaxies actually look like

Guo et al. (2015) analysed 3 239 high-redshift, star-forming galaxies studied as part of the Cosmic Assembly Near-Infrared Deep Extragalactic Legacy Survey (CANDELS). They use a conservative definition for clumps—a clump has to contribute at least 8% of the UV light of a galaxy, which excludes smaller clumps. They also require clumps to be off-centre, which excludes central clumps and is thus limiting in our context as a central clump could obviously contribute to bulge formation. One of their mosaics of clumpy galaxies is shown in Figure 5.

Figure 5. Figure 4 from Guo et al. (2015) showing HST images of high-redshift galaxies. The norm is that of massive, star-forming clumps within disks, rather than simple and smooth exponential disks.

They find that galaxies with log(M _⋆./M_⊙) > 9.8 have a 55% probability of being clumpy at redshift z ~ 3, down to 15% at z ~ 0.5, which is largely due to the fact star formation declines with decreasing redshift. At all redshifts, the clump contribution to rest-frame UV light peaks at log(M _⋆./M_⊙) > 10.5—the current stellar mass of the Milky Way. Integrating over both clumpy and non-clumpy galaxies, they find that 4–10% of the star formation takes place within these massive clumps.

In other words, the clumpy galaxy is a very plausible assumption for the initial conditions of the Milky Way. Much has been said in this review that the pure disk galaxy can work, as Kormendy et al. (2010) has pointed out. However, the smooth and massive exponential disks are widely seen in observations of the local universe. That clumpy galaxies are the norm for high-redshift observations suggests that most massive galaxies have passed through a clumpy phase.

4.2. Promising insights from simulations

The prevalence of clumps in high-redshift galaxies renders them a legitimate point of discussion for the origin of the Milky Way.

A way to test this is with comparison to the Milky Way, which Inoue & Saitoh (2012) did. They used an N-body/SPH model to study the evolution of an isolated disk galaxy where the clumps migrate to the centre via dynamical friction and form a clump-origin bulge. The final bulge resembles what they call a pseudo bulge, and is referred to as a bar throughout this paper and most of the Milky Way literature. The surface density profile is nearly exponential, the final shape of the bar is boxy, and the rotation is significant. The resulting stars are old and metal rich, with a flat star-formation history in an interval of Δt ~ 2 Gyr followed by a rapid decline in star formation. They obtain a metallicity gradient that stretches across the full vertical extent of the bulge. All of these properties are qualitatively consistent with what is observed for the Milky Way bulge. However, these properties can be matched by other models, and further this is a comparison of a single clump-origin bulge model to those of the bulge. It is worthy of consideration, but it is far too premature to declare victory.

One concern is that of whether or not the clumps actually do migrate to the centre, as they do the simulations previously discussed in this review. This is largely a question for theory, as the migration duration is too long for the baseline of observations. In the simulations of Hopkins et al. (2012), the inclusion of their prescription for stellar feedback disrupted the clumps, and prevented them from migrating to the bulge where they can coalesce. This prediction is not reproduced by the simulations of Bournaud et al. (2014), who find that the ejection of stars in clumps due to stellar feedback is compensated by their accretion of gas from the gas-rich disks in which they are contained.

Mandelker et al. (2014) studied 770 snapshots of 20 simulated galaxies using adaptive-mesh refinement cosmological models of Galaxy evolution. The global number of clumps were consistent with those in observations, an important check of the violent-disk instability hypothesis of clump formation. They did not study the properties of the final resulting bulge, but did say that if the clumps can survive accretion onto the centre of the galaxy, they are expected to accrete gas from the surrounding interstellar medium (similarly to Bournaud et al. 2014), and will thus show gradients in their mean properties with respect to separation from the centre of their galaxy, such as those measured by Förster Schreiber et al. (2011) that clumps closer to the centres of their disks are redder, older, and more massive. Further, they found in their simulations was that a full 91% of galaxies develop a bulge clump. These are massive, typically equivalent to 40% of the disk mass, with 20% of the star formation, and gas fractions of less than 1%.

Mandelker et al. (2017) study 34 galaxies with more sophisticated prescriptions. Among their findings, they find that the inclusion of radiation pressure disrupts the smaller clumps, reducing their lifetimes to a few free-fall times, but that the more massive and dense clumps still nevertheless survive and migrate to the centre. The inclusion of radiation pressure reduces the number of long-lived clumps by 81%. Radiation pressure has little to no effect on the bulge clumps, with ~83% of simulated galaxies hosting a bulge clump.

5 A CLUE AS TO THE ORIGINS OF THE BULGE FROM APOGEE

A significant clue as to the origin of the metal-poor stars in the bulge has been identified by the APOGEE collaboration. Previously, the detailed chemical abundance trends for the bulge have only been interpreted in their mean, due to the large observational error. The mean [X/Fe] vs. [Fe/H] can and has be compared to the thin disk, thick disk, and halo, but the scatter has not yet been of particular insight.

Schiavon et al. (2017) found a population of nitrogen-rich stars in the bulge, predominantly at [Fe/H] ⩽−1.0. These stars have enhanced nitrogen, aluminum, and depleted carbon, characteristic of the ‘second-generation’ stars in globular clusters (Carretta et al. 2009). Given that even surviving globular clusters must have been far more massive at birth (Conroy 2012) to produce their second generation, this suggests that between 50 and 100% of bulge stars with [Fe/H] ⩽−1.0 formed in disassociated globular clusters.

This is a clue to the origin of the bulge, but it is not clear which line of evidence it can be used to support. That is why it is left as a separate section.

An additional clue, also from the APOGEE dataset, may have been identified by Fernández et al. (2017), who found plausible second generation globular clusters in the Milky Way with [Fe/H] > − 1.0, enhanced Na, Al, and deficient C, Mg relative to Galactic abundance trends. This is noteworthy in part as the surviving Galactic globular clusters do not show the Al-Mg anti-correlation at such high metallicities.

6 DISCUSSION AND CONCLUSION

The level of inputs new to the last decade, both observational and theoretical, that can inform and constrain bulge formation scenarios is truly spectacular. Global photometric maps are now available from the optical through to the mid-infrared, with substantial coverage in the variability domaine. Spectroscopic data sets are now available towards a large fraction of bulge stars, and towards the full metallicity range. Knowledge of what high-redshift galaxies look like, including plausible Milky Way precursors, is greater than it is ever been. The breadth and depth of models are constantly increasing.

It would be tempting to say that the situation remains one of uncertainty between different scenarios, but that would be so limiting as to be inaccurate. In a competition between the buckling disk and the classical bulge, the buckling disk is winning. The peanut/X-shape, the long bar, the correlation between the mean and skewness of the velocity distribution functions, and so on are non-trivial predictions that are required of any theory of bulge formation and at this time require a buckling disk. The classical bulge may dominate for the 5% of stars with [Fe/H] ≲ −1.0 and a minority of more metal-rich stars, but that is an upper limit on its contribution.

The clump-origin bulge scenario may prove to be a viable alternative. It cannot be ignored given the ubiquity of star-forming, gas-rich clumps in high-redshift galaxies. It is more likely than not that the Milky Way was, at one time, a clumpy galaxy. More simulations and comparisons of said simulations to observations are needed to ascertain whether or not this describes the bulge assembly history. One plausible hybrid scenario, is if the clumpy phase of the Milky Way led to the thick disk (Inoue & Saitoh 2014), with the thick disk predominantly responsible for bulge stars with − 1.0 ≲ [Fe/H] ≲ -0.50 (Di Matteo et al. 2015).

One thing is certain, the Milky Way bulge is a sensitive probe of Galactic assembly history, and research of its properties will continue yielding insights thereof.