© ESO, 2014

1. Introduction

Boxy and peanut-shaped bulges are present in about half of edge-on disk galaxies (Lütticke et al. 2000). The closest example of a boxy bulge can be found in our Galaxy (Okuda et al. 1977; Maihara et al. 1978; Weiland et al. 1994; Dwek et al. 1995). Even if some studies (Binney & Petrou 1985; Whitmore & Bell 1988) have proposed that these structures can be formed during accretion events, their high frequency and relation to the fraction of barred galaxies in disks (Eskridge et al. 2000; Menéndez-Delmestre et al. 2007; Marinova & Jogee 2007; Aguerri et al. 2009) suggest that a more common mechanism may be responsible for shaping the central regions of galaxies, giving them their boxy- or peanut-shaped morphology.

A number of numerical investigations (Combes & Sanders 1981; Pfenniger et Friedli 1991; Athanassoula 2005; Martinez-Valpuesta et al. 2006) have indeed shown that boxy bulges can be manifestations of secular processes that occur in disk galaxies such as thick stellar bars seen edge on. During their evolution, stellar bars can indeed go through one (or multiple) buckling phase(s), which are the consequences of vertical instabilities, and depending on the bar viewing angle, the resulting thick structure can appear boxy, if observed mainly along the bar major axis, or peanut-shaped, if observed mainly along the bar minor axis. It is also possible that a combination of these mechanisms, satellite accretion and bar instability, may be responsible for some of the observed bulge morphologies (Mihos et al. 1995).

Observations suggest that boxy bulges do not represent a homogeneous class of objects. Studies of galaxies indeed demonstrate that boxy bulges display a range of properties in their kinematics and stellar populations (Williams et al. 2011). Some of them show a constant rotation with height above the plane, while some others do not, thus indicating that cylindrical rotation does not necessarily characterize these structures; some of them show negative vertical metallicity gradients that can be accompanied by positive [α/Fe] gradients, while some are more homogenous, indicating that different stellar populations can dominate these structures at different vertical distances from the galaxy midplane. The complexity of the observed characteristics may also be due to the concomitant presence of a classical¹ bulge in the inner regions of some galaxy disks. It is not simple to identify the presence of such components in boxy or peanut-shaped structures by the characteristics of their stellar populations and/or kinematics. The existence of a vertical metallicity gradient, for example, does not necessarily imply the presence of a classical bulge (Bekki & Tsujimoto 2011; Martinez-Valpuesta & Gerhard 2013). Detection of cylindrical rotation does not necessarily imply that the boxy bulge is the result of pure bar instabilities (Saha et al. 2012, 2013).

Understanding how much of a classical bulge is present in structures otherwise mostly shaped by secular evolution processes is fundamental not only in interpreting observations of galaxies, but also for understanding the formation and evolution of the central regions of the Milky Way (MW). Indeed, over the last two decades, a number of studies have elucidated the complexity of the boxy, peanut-shaped structure at the center of our Galaxy (McWilliam & Rich 1994; Zoccali et al. 2006, 2008; Lecureur et al. 2007; Babusiaux et al. 2010; Shen et al. 2010; Ness et al. 2012, 2013a,b; Gonzalez et al. 2013; Bensby et al. 2013), but no consensus has yet been reached on how to interpret these important results.

A number of studies have shown evidence of the presence of a metal poor, α-enriched component, whose kinematic properties are significantly different from that of the metal-rich, nearly solar [α/Fe] bulge component (Babusiaux et al. 2010; Hill et al. 2011; Ness et al. 2013b). Babusiaux et al. (2010), for example, have pointed out the presence of two distinct populations along the bulge minor axis, with distinct kinematic properties: a metal-poor population ([Fe / H] ~ − 0.3 dex) whose radial velocity dispersion is constant with latitude, and a metal rich population ([Fe / H] > 0.1 dex) whose radial velocity dispersion decreases substantially with the distance from the Galactic mid-plane. The contribution of these two populations changes with latitude, the metal rich component disappearing when moving away from the plane, where the metal poor population is becoming dominant. By comparing the data with the N-body model of Fux (1999), Babusiaux et al. (2010) interpreted these two populations as being the signature in the inner Galactic disk of the simultaneous presence of a classical metal-poor bulge and a metal-rich population with bar-like kinematics. However, the presence of a classical bulge in the MW disk and even its role in explaining the characteristics of the observed metal poor, α-enhanced population is debated. Shen et al. (2010), for example, have questioned the existence of any classical bulge in the MW, ruling out the possibility that our Galaxy has a classical bulge with a mass greater than ~15% of the stellar disk mass.

In addition, recent studies have pointed out the similarities between the metal poor populations of the galactic bulge and the thick disk population at the solar vicinity (Meléndez et al. 2008; Ryde et al. 2010; Alves-Brito et al. 2010; Bensby et al. 2010; Gonzalez et al. 2011; Ness et al. 2013a). The ARGOS survey (Freeman et al. 2013), in particular, is currently mapping the Galactic bulge over a large extent of latitudes and longitudes, contributing significantly to our understanding of how the bulge populations differ in their spatial redistribution, chemical properties and kinematics (Ness et al. 2012, 2013a,b). The results of this survey mainly suggest the existence of at least three primary components in the MW bulge. Two components (defined respectively as component A and B in their paper), with [Fe / H] > − 0.5 dex, are part of the boxy/peanut shaped bulge, with component B ([Fe / H] ~ − 0.25 dex) being dynamically hotter than component A ([Fe / H] ~ 0.1 dex), rotating 20% faster than A, and being more prominent at high latitudes; component C ([Fe / H] < − 0.5 dex) has the highest radial velocity dispersions in the fields analyzed, nearly constant both in latitude and longitude, and has been explained by Ness and collaborators as part of the inner thick disk.

Is it possible to explain the complexity and richness of these data with a “simple” scenario for the formation of the MW bulge? How many of the observed characteristics, for example, are due to a complex accretion history for the Galaxy and how much is due to secular evolution of the disk, with the possible contribution of an older and more metal-poor thick-disk component and/or a classical spheroid?

In this paper, we try to answer to some of these questions by showing that at least part of the characteristics observed among the stellar populations of the Galactic bulge, such as the angular momentum support, radial velocity dispersions, and dependence of the fractional contribution of different stellar populations on latitude, can be explained as the simple result of the mapping of a stellar disk onto a boxy/peanut-shaped bulge. We will show that a large portion of the stellar disk, from the innermost regions to the outer Lindblad resonance (OLR) of the bar, is involved in the formation of a boxy/peanut structure, as a result of the radial migration initiated before the buckling instability of the formation of the bar.

In particular, we will show that the two populations contributing to the boxy structure (component A and B in Ness et al. papers) have kinematic and spatial characteristics compatible with a origin in different regions of the disk. Specifically, we show that, on average, component B formed from stars with initial radii larger than the stars that comprise component A. On the basis of the observed characteristics, the spatial distribution, chemistry and kinematics, we propose that component B is mostly made of the young MW thick disk (stellar ages between 8 and 10 Gyr, as observed at the solar neighborhood; Haywood et al. 2013), while A is mostly made of stars that originated in the inner thin disk. The presence of a small (B/D = 0.1), classical bulge is not excluded and we propose a possible signature of its presence should be searched for in observations of the bulge. We will discuss the possibility that the old thick disk (ages greater than 10 Gyr; Haywood et al. 2013) is the main contributor to component C, which is not part of the boxy/peanut structure.

The paper is organized as follows: after describing the N-body models used for the analysis (Sect. 2), we will present the main results from an analysis of these simulations in Sects. 3 and 4; then in Sect. 5, we state the main conclusions of our work.

2. Initial conditions and numerical method

The set of three simulations analyzed in this paper are the same already presented in Di Matteo et al. (2013): they consist of an isolated disk, with a varying bulge-to-disk ratio (B/D = 0., 0.1 and 0.25, respectively), and containing no gas. The dark halo and the optional bulge are modeled as a Plummer sphere (Binney & Tremaine 1987). The dark halo has a mass M_H = 1.02 × 10¹¹M_⊙ and a characteristic radius r_H = 10 kpc. The bulge, when present, does not rotate initially (but see Sect. 3.4 for its final rotational content), has a mass M_B = 9 × 10⁹M_⊙ and characteristic radius r_B = 1.3 kpc, for the case with B/D = 0.1, and M_B = 2.2 × 10¹⁰M_⊙ and r_B = 2 kpc, for the case with B/D = 0.25. The stellar disk follows a Miyamoto-Nagai density profile (Binney & Tremaine 1987), with mass M_∗ = 9 × 10¹⁰M_⊙ and vertical and radial scale lengths given by h_∗ = 0.5 kpc and a_∗ = 4 kpc, respectively. The initial disk size is 13 kpc (that is, initial stellar positions are generated between r = 0 and r = 13 kpc; see Appendix A for models employing an initially more extended stellar disk), and the Toomre parameter is set equal to Q = 1.8. The galaxy is represented by N_tot = 30 720 000 particles redistributed among dark matter (N_H = 10 240 000) and stars (N_stars = 20 480 000). To initialize particle velocities, we adopted the method described in Hernquist (1993). The amplitude of the initial fluctuations around the equilibrium virial ratio Q_vir in the resulting models is ΔQ_vir/Q_vir ~ 1% over the first Gyr of evolution.

To model galaxy evolution, we employed a Tree-SPH code, in which gravitational forces are calculated using a hierarchical tree method (Barnes & Hut 1986). The code has been presented in Semelin & Combes (2002) and we refer the reader to this paper for a full description. For the dissipationless simulations analyzed in this paper, the SPH part of the code has been switched off, and gravitational forces are calculated using a tolerance parameter θ = 0.7 and include terms up to the quadrupole order in the multiple expansion. A Plummer potential is used to soften gravity at scales smaller than ϵ = 50 pc. With this spatial resolution, it is possible to resolve the vertical structure of thin disks, and follow small scale inhomogeneities.

The equations of motion are integrated over 4 Gyr, using a leapfrog algorithm with a fixed time step of Δt = 0.25 × 10⁵ yr. Some of the main characteristics of the simulated models and their evolution are shown in Fig. 1.

3. Results

Several studies have pointed out the possibility that stars observed today at a given location of a galaxy disk may be born in very different regions from those where they are currently observed (Sellwood & Binney 2002; Debattista et al. 2006; Ros˘kar et al. 2008a,b; Minchev & Famaey 2010; Brunetti et al. 2011; Minchev et al. 2011, 2012). All these investigations have explored the impact this redistribution has on the observed properties of stellar disks, in particular on their external regions. None has yet investigated in detail the effect of this redistribution on the properties of the bar and bulge regions of disk galaxies. In the following sections, we show that this impact is significant and cannot be neglected when interpreting the observed properties of boxy/peanut shaped structures in the MW and other galaxies.

3.1. Birth radii and radial migration: face-on view

Our modeled bulgeless disk initially has an azimuthally symmetric stellar distribution, with no sign of asymmetries for the first 0.8 Gyr of evolution, as shown by the Fourier analysis of the surface density distribution of the face-on stellar disk (see Fig. 2, top panel). At about t = 0.8 Gyr, a bar and spiral arms start to develop, grow rapidly and stay strong for the following Gyr, until about t = 2 Gyr. At this time the stellar bar undergoes a vertical instability, as a consequence, its strength diminishes and its scale height increases considerably. This is the epoch when the bar changes from a thin structure within the plane to a thick structure which may appear either boxy or peanut-shaped depending on the angle between the observer’s line-of-sight and the bar major axis (see Fig. 2, bottom panels). In this simulation, the bar’s resonances are located at the corotation radius, r_CR = 8 − 9 kpc, the inner Lindblad resonance, r_ILR = 2 − 3 kpc and outer Lindblad resonance, r_OLR = 13 kpc for simulation ages between 1 and 2 Gyr. These resonances move outward slightly as the bar slows down and decreases in strength. Since the OLR in these models coincides with the initial extent of the disk (see Sect. 2), the mapping that we will discuss in this section is limited to the disk inside r_OLR. We will show in Appendix A that these results are also valid for simulated disks whose initial extent and resonances location are different from those presented in this section, showing indeed that the disk up to the OLR participates in the formation of the boxy-peanut shaped bulge structure. During the whole simulation, the bar semi-major axis length, r_bar, is approximately 7−8 kpc, a factor of about two greater than the length of the bar observed in the MW. To present our results in a more general way, independent on the bar size in the following analysis all spatial scales are given in units of r_bar. In this unit, between t = 1 and t = 2 Gyr, r_ILR = 0.3 − 0.4r_bar,r_CR = 1.1 − 1.2r_bar,r_OLR = 1.9r_bar.

How does the formation and presence of asymmetries in the disk, and in particular, the presence of the stellar bar – the strongest asymmetry in our simulations – affect the spatial redistribution of stars in the disk? To answer this question, similar to what was done in Di Matteo et al. (2013), we have selected stars according to their initial radii², defining six different regions in the disk, as follows: r_ini ≤ 0.4r_bar, 0.4r_bar ≤ r_ini ≤ 0.7r_bar, 0.7r_bar ≤ r_ini ≤ r_bar, r_bar ≤ r_ini ≤ 1.3r_bar, 1.3r_bar ≤ r_ini ≤ 1.6r_bar, 1.6r_bar ≤ r_ini ≤ 1.9r_bar, with r_ini being the distance, in the disk plane, of stars from the galaxy center. This selection is shown in Fig. 3.

As a consequence of the angular momentum (AM) redistribution initiated by the bar and spiral arms, stars tend to diffuse in the disk as soon as stellar asymmetries start to develop. However, while stars inside the inner Lindblad resonance stay mostly confined in the inner bar region (see also Martinez-Valpuesta & Gerhard 2013; Pfenniger et Friedli 1991), outer disk stars, in particular those at and beyond corotation, migrate both outward and inward, reaching both the edges and the center of the disk. Within a few rotational periods at the epoch of formation of the stellar asymmetries, outer disk stars are able to reach the inner disk, contributing to populating the bar: their distribution shows a clear m = 2 asymmetry elongated with the bar major axis and tends to accumulate in two stellar over-densities at the edges of this structure (see Fig. 3). At the onset of the bar vertical instability, those stars that are close to the vertical inner Lindblad resonance (VILR), which is at about 0.8r_bar from the center, are scattered to greater heights becoming part of the boxy/peanut-shaped structure. Since stars from a large range of initial birth radii are able to reach the bar region before its vertical buckling, as a consequence, the resulting bulge is populated by a mixture of populations, from stars born in situ (i.e. in the inner disk) to stars coming from all of the outer radii, from those just outside the bar to the outermost extent of the disk, at about 2r_bar.

It is interesting to note that the orbital characteristics of stars that end up in the bar region depend on where they originated in the disk. Figure 4, for example, shows the distribution of orbital pericenter (r_min) and apocenter (r_max) radii for stars originating in five different disk regions: the median of both of these radial distributions increases with increasing originating radii of the stars. This effect is such that outer disk stars which end up in the bar/bulge region tend to orbit over a larger portion of the disk than stars born in the inner disk. This memory of their initial location in the disk translates directly into only limited variations (see bottom panel in Fig. 4) in their Jacoby energy (see Binney & Tremaine 1987, for a definition), as already pointed out by Martinez-Valpuesta & Gerhard (2013). From Fig. 4 and Table 1 one can also see that while stars born inside the VILR are trapped in the inner disk (for example, 90% of stars with 0.4r_bar ≤ r_ini ≤ 0.7r_bar have apocenters inside the bar), outer disk stars can have pericenters penetrating in the bar region, and apocenters outside it; in other words not all the outer disk stars penetrating the bar region are confined within the bar region over their whole orbit. As an example, among the stars with r_bar ≤ r_ini ≤ 1.3r_bar, 90% of them have pericenters inside the bar (r_min<r_bar), while only 30% have their whole orbit confined inside the bar (r_max<r_bar). There is thus a high fraction – about 60% – that have pericenters penetrating inside the bar region, but apocenters outside which lie outside of it. It is also interesting to note that, among stars with r_bar ≤ r_ini ≤ 1.3r_bar, there is a significant fraction (~20%) whose apocenters are inside the bar region and are also smaller than their initial birth radii (i.e. both the conditions r_max<r_bar and r_max<r_ini are satisfied, see Table 1): those stars may be considered as the true migrators, since their orbital radii have significantly changed (this is guaranteed by the condition r_max<r_ini). However, since all the stars whose orbit penetrates, at least partially, into the bar contribute to its properties, in the following, we adopt a more general definition: “outside-in migrators” are stars that enter the bar region and spend part of their orbit at distances from the galaxy center smaller than their initial birth radii.

3.1.1. Angular momentum

The spatial redistribution of stars in the disk initiated at the epoch of bar formation is a consequence of AM redistribution, as previously shown in several papers (e.g. Minchev & Famaey 2010; Brunetti et al. 2011; Minchev et al. 2011, among others). Figure 5 (left column) shows maps of face-on projections in the variation of angular momenta, ΔL, of stars born at different radii. This variation is evaluated at times between t = 0 and t = 2.75 Gyr. At this epoch, t = 2.75 Gyr, the bar has already acquired its boxy/peanut shaped morphology. As expected, these maps clearly show that for outside-in migrators trapped in the bar region, those with larger birth radii experience larger AM changes. However, it is interesting to note that, even if outside-in migrators have experienced the largest changes, their final AM still retains the memory of their initial birth radii: Fig. 5, middle column, shows projected face-on maps of L_norm = L/L_circ, that is of the AM, at time t, normalized by the AM of the corresponding circular orbit at that radius, at the same time t. As usual, stars have been selected with respect to their initial provenance in the disk. From these plots, it is evident that, in the bar region, the AM content depends on the birth radius, and that the AM content increases for stars with larger birth radii. For example, among the stars ultimately found in the bar, those born inside 0.4r_bar have, on average, an AM content 6−7 lower than those of stars born beyond 1.6r_bar that have subsequently migrated in the central kpcs. In other words, in the bar region, the larger the initial birth radius of a star, the greater its final AM. As we will discuss in the following sections, this finding has important consequences for the spatial redistribution and for the kinematics of stars in the boxy/peanut-shaped structure. Note also that stars born inside the VILR show also an interesting trend in the vertical velocity dispersion, σ_z (Fig. 5): inside the bar region, stars which have smaller birth radii also have the lower σ_z.

3.2. Dissecting the structure of boxy/peanut bulges

If the stellar bar is the result of the mixing of various stellar populations with different AM contents and vertical velocity dispersions, it is natural to investigate how these populations are mapped into the vertical structure of the boxy/peanut-shaped bulge.

In Fig. 6, we show face-on and edge-on maps of the final evolved stat of our modeled bulgeless disk (i.e. t = 3.95 Gyr). Stars have been classified according to their birth radii. To facilitate the comparison to the MW bar, the final configuration has been rotated so that the bar viewing angle is about 20 degrees (Bissantz & Gerhard 2002; Shen et al. 2010, see also Wegg & Gerhard 2013, for a more recent estimate of the bar inclination). Figure 6 shows the resulting face-on maps of stars which originate from six different annuli (defined in Sect. 3.1). Stars born in the outer regions of the simulated disk that are close to the OLR, are found in the central kpcs of the galaxy and are part of the bar structure. From the corresponding edge-on maps (Fig. 6) one sees that stars with different birth radii do not redistribute in the same way in the boxy/peanut shaped structure. Stars born in the innermost disk regions tend to make up populations which exhibit rounder isophotes in the projected maps than the population made up of outside-in migrators. The shape of the stellar distribution indeed becomes increasingly boxy/peanut-like for stars which migrated inward from the outer disk or were born around the VILR. Even stars which migrated from the edge of the disk participate in the formation of the thick bar (most keenly visible in the edge-on view of their stellar distribution which shows the characteristic X-shape profile). This is a consequence of the AM redistribution initiated by the formation of the bar: these stars migrate from the outer to the inner disk at the epoch of bar growth (see Figs. 2 and 3), and are, at the time of bar buckling, captured at the VILR and thus participate in the orbital families supporting the boxy/peanut shaped structure.

The dependence of the isophotal shape on individual populations with different birth radii implies that, at any given location in the bulge, stars of different provenance have different relative contributions to the local stellar density. This is clearly shown in the edge-on maps of Fig. 6, where the fractional contribution to the local stellar density of stars of different provenance are elucidated. Stars born in the inner disk (r_ini ≤ 0.4r_bar) dominate the local stellar density only in the very inner regions of the bulge, inside ~0.4r_bar, and their contribution to the local density decreases very quickly when moving horizontally or vertically out from the center. The relative contributions of stars of increasing birth radii reaches a maximum in progressively more distant regions of the bulge. Also stars born outside the VILR are found in the bulge and their contribution is most significant at its outer edges.

Another way of showing this inhomogeneous mapping of a stellar disk into a boxy bulge is shown in Fig. 7, where we have divided the edge-on (x − z) projection of our modeled bulge into a 5 × 5 grid, with each panel in the grid having a size of 0.3r_bar × 0.3r_bar, divided in such a way as to cover the whole extent of the modeled bulge (that is a [− 0.7r_bar,0.7r_bar] × [ − 0.7r_bar,0.7r_bar] region; see the top-right panel for a x − z view of the grid). To avoid contamination from foreground and background stars, in analogy with what commonly is done in bulge surveys, we have selected stars only inside the bar region, which in our model corresponds to the region defined by | y | ≤ 0.8r_bar (see the top-left panel in Fig. 6). The resulting panels in Fig. 7 show the fractional contribution to the local stellar density of stars born in different regions of the disk.

From this analysis, we can deduce that:

• the fractional contribution to the boxy bulge of stars of different provenance changes both with latitude (i.e. vertically) and with longitude (i.e. horizontally);

• for increasing longitude at any given latitude, the contribution of outside-in migrators increases, in agreement with the dependence of the AM content on the star birth radius (Fig. 5). For example, for the region, − 0.1r_bar<x< 0.1r_bar and − 0.1r_bar<z< 0.1r_bar, only 5% of stars have birth radii greater than r_bar (i.e., external to the bar region itself). This fraction increases to more than 25% at the edge of the bar, close to the galaxy midplane (− 0.7r_bar<x< − 0.4r_bar and − 0.1r_bar<z< 0.1r_bar; 0.4r_bar<x< 0.7r_bar and − 0.1r_bar<z< 0.1r_bar);

• the redistribution depends also on latitude, as expected from the trends found for σ_z (Fig. 5). Along the bar minor axis, for example, at low latitudes about 50% of the stars have r_ini ≤ 0.4r_bar, while, at high latitudes, half of the stars have r_ini ≥ 0.7r_bar, with more than 20% of the stars having an outer disk origin (r_ini>r_bar);

• at the edges of the X-shaped bulge structure, stars whose origin is external to the bar (i.e. with r_ini>r_bar) represent about 30% of the total local density.

We emphasize that it is stars from all of the disk, from the center to the OLR, that contribute to the formation and structure of a boxy bulge, and that stars formed outside the bar region and subsequently by bar instabilities can represent a significant fraction of the stellar density even at high latitudes. In the next section, we show how the birth radius of a star reflects into its subsequent kinematics as a contributor to the bulge morphology.

3.3. Imprints of the stars birth radii on the bulge kinematics

In Fig. 8, we sort, as previously done, stars according to their birth radius. The resulting line-of-sight velocity and velocity dispersion maps for these different regions are shown in Fig. 8. Note that in this plot, both quantities are line-of-sight velocities, that is we are considering only the velocity component parallel to the y-axis, that is perpendicular to the x − z plane. As for Figs. 6 and 7, the bar is inclined at an angle of 20 degrees with respect to the line-of-sight. We will fix the inclination of the bar for all the following analysis to facilitate direct comparison with the properties of the MW. From the maps in Fig. 8, we can deduce the following trends:

• in the bulge region, the larger the birth radius of a star, the higher its line-of-sight velocity. As shown in Fig. 8, moving from smaller to larger birth radii, the average absolute line-of-sight velocity in the galaxy mid-plane increases from ~100 km s^-1 (r_ini<r_bar) to more then 220 km s^-1 (1.6r_bar<r_ini< 1.9r_bar). This result is clearly in agreement with the results shown in Fig. 5, that is the larger the birth radius of a star, the higher its AM content. As a consequence of this higher AM, while stars born in the inner-disk regions stay confined there, outer disk stars which have reached the central kpcs, and which have a higher AM, span a larger volume of the bar and disk than in situ, inner-disk stars (as already shown in Fig. 4 and Table 1);

• in the bulge region, the larger the birth radius of a star, the greater its line-of-sight velocity dispersion. For example, in the central regions of the bulge, the velocity dispersions increases from about 150 km s^-1 (for stars with r_ini< 0.4r_bar) to about 200 km s^-1 (for stars with r_ini>r_bar). This trend is a consequence of the previous point and of the orientation angle of the bar. Indeed because the AM content of outside-in migrators is higher than those of in situ stars, these stars mostly contribute to populating the edges of the bar, thus their orbit span a larger radial extent (i.e. a larger distance from the galaxy center) than those of in-situ, inner-disk stars. The bar being oriented at 20 degrees, the velocity component parallel to the line-of-sight is determined mainly by motions along the bar major axis. As a consequence, outer disk stars populating the bar span the larger extent along the bar, and thus for an orientation of 20 degrees, show the highest line-of-sight velocity dispersions. Note that this trend of increasing velocity dispersion with increasing star birth radius is mostly evident inside the VILR. Outside this radius, stars in the boxy bulge seem to show a line-of-sight velocity dispersion which is independent on their birth radius and similar to that of stars born near VILR.

To elucidate the dependence of the stellar kinematics on the birth radius of the stars and on latitude and to facilitate comparison with observations, we have plotted the line-of-sight velocity and the velocity dispersion of stars in the bulge region, as a function of the projected distance x from the center, for different vertical distances z from the galaxy mid-plane. This is similar to the analysis done for observations of stars in the Galactic bulge by Ness et al. (2013b), except that, to avoid having to define the Sun position in our model – which is always somewhat arbitrary – we evaluate the line-of-sight integrated quantities instead of radial velocities, i.e. instead of velocities evaluated with respect to the Sun-star direction, as done in their investigation. Since the relation between the two is v_los = v_rcos (l), where l is the longitude of a star, it means that at most, at the edge of the bulge (l ~ 15°, see for example Freeman et al. 2013), we underestimate the radial velocities of stars in our model by ~4%. Figure 9 shows the result of this analysis. Line-of-sight velocities (upper panels) and velocity dispersions (lower panels) are shown as a function of x, for three different regions of provenance of the stars (r_ini ≤ 0.7r_bar, 0.7r_bar<r_ini ≤ 1.3r_bar, and r_ini> 1.3r_bar) and for six different vertical distances from the galaxy mid-plane. It is clear that:

• 1.

  stars born in the inner disk show an almost cylindrical rotation, while stars originating in the outer disk do not. The resulting total (i.e. taking into account all the stars, independent of their birth radii) line-of-sight velocity curve (Fig. 10) still shows a cylindrical rotation, with values similar to those observed in the Galactic bulge. For example, comparing the vertical distances 0.14r_bar ≤ z ≤ 0.35r_bar with the ARGOS fields between b = − 5° and b = − 10°: at x/r_bar = ± 0.7, which corresponds to approximately l = ± 20°, the differences in the velocity curves with latitude are ~20 km s^-1, similar to what observed.

• 2.

  at any given latitude, stars that have their origins further out in the disk have higher line-of-sight velocities. For example (top panels, Fig. 9), stars born between 0.7r_bar<r_ini ≤ 1.3r_bar have a v_los which is about 30% higher than that of stars born inside 0.7r_bar. This value increases up to about 40% if one compares the v_los of stars born inside 0.7r_bar with that of stars born outside 1.3r_bar.

• 3.

  Whilst the total line-of-sight velocity curve is dominated by stars originating between 0.7r_bar<r_ini ≤ 1.3r_bar (cf. Figs. 9 and 10), innermost (r_ini< 0.7r_bar) and outermost (r_ini> 1.3r_bar) stars represent, respectively, the low and high velocity tails of the line-of-sight velocity distribution at different latitudes and longitudes. This is a consequence of point 2.

• 4.

  the line-of-sight velocity dispersions show a trend with latitude. This trend is in the sense that the velocity dispersion decreases with increasing latitude (bottom panels, Fig. 9)

• 5.

  the line-of-sight velocity trends are the same for outer and inner disk stars, but their absolute values are larger for the former compared to the latter (Fig. 9).

The trends listed in points 3 and 4 are exactly those found by Ness et al. (2013b) for populations A and B in the galactic bulge. Moreover, population B rotates 20% faster than A, as would be the case if the average birth radius of the stars that make up this component was larger on average than that of population A (point 2). We will discuss on this important point in Sect. 4.

3.4. Classical bulges

All our previous analysis has been restricted to the case of a boxy/peanut shaped bulge formed in a pure stellar disk (B/D = 0.). It is natural to ask whether these results are sensitive to the inclusion of a classical spheroid and what signature would be left by the presence of such a component underlying a boxy/peanut-shaped structure. With the goal of answering these questions, we have analyzed two simulations, with an initial bulge-to-disk ratio respectively equal to B/D = 0.1 and 0.25. In this section we present the main results of this analysis.

To mimic the likely configuration of the bulge/bar of the MW, we have rotated the stellar bar so that it has an inclination of about 20 degrees with respect to our line-of-sight in the analysis. After removing foreground and background stars as we have done previously (Fig. 7), we discuss the mapping of the stellar disk into the boxy bulge for these two models (Fig. 11)³. Concerning the birth radii of stars⁴ in the boxy bulge, we find the following results:

• as was the case for B/D = 0, the fraction of inner stars (r_ini< 0.4r_bar) decreases with increasing longitude and/or latitude. For the case with B/D = 0.1, their fractional contribution changes from more than 50%, in the galaxy central regions, to about 20% when moving vertically along the bar minor axis. For the case with B/D = 0.25, at high latitudes the contribution of inner disk stars is higher, about 35% of the total, because at these distances from the plane the contribution of the classical spheroid increases (Fig. 11).

• The relative contribution of inner stars to the local density is very sensitive to changes in longitude: these stars constitute the majority of the stellar mass in the inner bulge region (− 0.1r_bar<x< 0.1r_bar and − 0.1r_bar<z< 0.1r_bar), while their contribution is ~10% in the outer bulge region, for a B/D = 0.1. For a B/D = 0.25, their fractional contribution in the outer regions of the bulge is higher (at most ~20%), and this is a consequence of the presence of stars from the classical bulge at these longitudes, as would be expected as the classical bulge contribution increases overall.

• as was the case for B/D = 0, outer disk stars (r_ini>r_bar, the bar scale length in both cases) are part of the bulge structure. Within the outer regions of the bar, indeed, they constitute between 20 and 30% of the local stellar content.

Because of the existence of a mass-size relation for bulges (see, for example, Gadotti 2009), the fractional contribution of a pre-existing spheroid to the local stellar density depends on the mass of the spheroid:

• a B/D = 0.1 classical bulge is mostly confined to the inner bulge regions (mostly inside 0.4r_bar, for a bulge core radius equal to 0.18r_bar), and, in our models, its contribution to the total stellar mass never exceeds 20%, regardless of the longitude or latitude;

• a B/D = 0.25 classical bulge, which is more extended (its initial core radius is equal to 2 kpc) and more massive, contributes stars all along the boxy bulge structure. Moreover, its contribution increases with latitude: as the thick bar stellar density diminishes, the classical bulge stellar contribution increases (up to about 50% of the total, see the behavior along the bulge minor axis in Fig. 11).

In our models, classical spheroids hidden in a boxy bulge thus show trends which depend on their mass. To summarize:

• while a B/D = 0.1 bulge nowhere in the boxy structure contributes for more than 20% to the local density, a B/D = 0.25 bulge can contribute as much mass as that of the bar at high latitudes (if no other structure, like a thick disk, is present, of course).

• While the fractional contribution of a B/D = 0.1 bulge is constant or slightly decreases with latitude, a B/D = 0.25 bulge shows a trend such that its contribution increases with latitude, with trends similar to those observed for outside-in migrators. However, while outside-in migrators are part of the boxy/peanut shaped structure (see Fig. 6), classical bulges are not, because of their higher velocity dispersions. In other words, any stellar population whose orbits support the boxy structure (as populations A and B in Ness et al. 2013b) cannot be dominated by a classical component. Thus even if the metallicities of the stars are compatible with those expected from a classical bulge, on the basis of a mass-metallicity relation (Gallazzi et al. 2005; Thomas et al. 2010), the morphology and kinematics of these stars will not be consistent with those expected for a classical bulge. We discuss the evidence for a classical bulge in the MW in Sect. 4.2.

A classical component leaves signatures in the kinematics of stars in the bulge. As found for the case with B/D = 0, we confirm that stars that originate further out the disk have higher rotational support and velocity dispersions. That is, outside-in migrators in the bulge are still recognizable because of their higher AM and because they are a dynamically warmer replica of populations born in the inner disk (Fig. 12). The trends found in Sect. 3.3 are thus robust and independent of the mass of a possible (hidden) spheroid. At the same time, stars originally in the classical bulge:

• are characterized by a much lower rotational support than stars that originate in the disk. Indeed even if stars that comprise the classical bulge acquire some AM during the disk and bar formation and evolution (see also Saha et al. 2012, 2013, and Fig. 12), these stars never reach rotational velocities similar to those of disk stars, whatever their location in the boxy structure.

• have line-of-sight velocity dispersions which diminish with latitude, but in such a way that the steepness of the relation is lower than that observed for stars originating in the disk (moving vertically from low to high latitudes, the central velocity dispersions decrease by ~20%)

• have line-of-sight velocity dispersions which also diminish with increasing longitude, but the gradient of the relation is flatter than that measured for disk stars, regardless of their birth radii.

4. Discussion: the MW bulge

The models presented in the previous sections are intended to simulate the formation of the bar and its subsequent buckling. We found that the spatial redistribution the stars in the disk undergo during this period generates a structure that has properties similar to that observed in the bulge of the MW. However, these simulations have their limitations. These simulations were not intended to model all the components present in a galaxy at the epoch of bar formation, but modeled only those that are sufficiently kinematically cold to be susceptible to having their orbital parameters changed significantly by bar formation, evolution, and the instabilities and asymmetries that this evolution generates.

The fossil chemical and dynamical record of stars in the solar vicinity suggest that the thick disk of the MW formed during a period where the gas out of which stars formed was well-mixed chemically and highly turbulent (Haywood et al. 2013; Snaith et al. 2014, and references therein). If the MW formed its bar at the end of its turbulent phase, about 8−9 Gyr ago, more than 50% of its current stellar mass was probably already in place (Snaith et al. 2014). This value is also in agreement with the extrapolation of the MW thick and thin disk masses (e.g. Fuhrmann et al. 2012) and with the estimates of the stellar mass growth of MW-type galaxies as a function of redshift (van Dokkum et al. 2013; Leitner 2012). The oldest part of this stellar disk (old thick disk in the nomenclature of Haywood et al. 2013) had probably too high vertical velocity dispersions to be involved in the bar instability process current values of the vertical velocity dispersions of old thick disk stars at the solar neighborhood are around 40 km s^-1, see Haywood et al. (2013). However, the youngest disk (young thick disk in the nomenclature of Haywood et al. 2013) had much lower velocity dispersions – current values at the solar neighborhood are ~25−30 km s^-1, see Haywood et al. (2013), similar to the values characterizing the stellar disks modeled in this paper. In other words, the youngest and thinnest component of the MW thick disk, which formed between 8 and 10 Gyr ago, must have been much more susceptible to the influences of bar instabilities than the oldest component of the disk, the old thick disk, which formed more than 10 Gyr ago (Haywood et al. 2013).

With these considerations in mind, we now investigate how much of the general scheme presented in the previous section is applicable to our Galaxy and its stellar populations. We will mainly compare our results to the set of ARGOS papers, which provide a detailed mapping of the MW bulge, and a “dissection” of its stellar populations that allows us to make a direct comparison to the results of our simulations.

4.1. The disk origin of the MW bulge

Briefly recapitulating the Introduction, the current picture of the MW bulge suggests that it is made of several distinct stellar populations, which have different chemical and kinematic characteristics, and whose fractional contribution to the bulge stellar density depends on latitude (Ness et al. 2012, 2013a,b). Interestingly, as shown by Ness et al. (2013a), the three main components found in the bulge – A, B and C – form a sequence in chemical characteristics, from population A with a high mean abundance and a low α enrichment (mean [Fe/H] and [α/Fe] ~ 0.1 dex), to population B with a intermediate mean metal abundance and moderate α enhancement (mean [Fe/H] ~ −0.25 dex and [α/Fe] ~ 0.2 dex), to population C, the most metal-poor (mean [Fe/H] ~ −0.7 dex) and α-enriched (~0.3 dex) of the 3 populations. This sequence is reminiscent of that found at the solar vicinity (see, e.g. Gonzalez et al. 2011). This is also confirmed by the findings of Bensby et al. (2013): the chemical characteristics of the bulge almost perfectly tracks that of the local disk⁵. The similarity in the elemental abundance distributions may imply, for stars in the bulge, the existence of an age sequence similar to that found within the solar vicinity (Haywood et al. 2013). In such a scheme, stars comprising component B would be associated with the stars that were originally part of the young thick disk, and thus would be kinematically colder than its more metal-poor counterpart in the old thick disk which may be the stars that comprise component C.

Therefore, the chemical properties of the Galactic bulge suggest a possible disk origin for all its main components. In the following, we will discuss how such an evolutionary path is supported by our modeling, by showing in particular that: (a) the kinematic characteristics of stars that comprise components A and B are indicative of their different provenance in the disk, with stars that makeup component B formed at larger distances from the galactic center on average than stars that comprise component A; (b) a significant classical spheroid can be excluded in the MW, thus leaving a disk origin as the only possibility for most of the stellar mass present in the MW bulge.

4.1.1. Components A and B

Of the three main components, A and B are the only ones that participate in supporting the boxy peanut-shaped morphology of the bulge (Ness et al. 2012). This indicates that they must have been sufficiently kinematically cold originally to be influenced by bar instabilities. Components A and B are detected throughout the bulge region, but in different relative proportions: B becomes dominant at intermediate (b = − 7.5 deg) and high latitudes (b = − 10 deg), while its relative contribution is comparable to that of A at lower heights (b = − 5 deg) above the plane (Ness et al. 2013a). Component B rotates 20% faster than A, and has radial velocity dispersions that, at all latitudes and longitudes, are a dynamically warmer replica of those of component A. Our models suggest that all these characteristics can be explained if the stars that makeup component B formed, on average, further out in galactic disk than the stars of component A⁶. We have seen, indeed, in Sect. 3.3 that that stars with larger birth radii that end up in the boxy bulge have higher final rotational support and radial velocity dispersions (this last trend is valid for the small bar viewing angle of the MW bar). Moreover, the fractional contribution to the peanut-shaped density distribution of component B tends to increase with latitude. This is another argument in favor of the scenario where the stars that comprise component B formed, on average, further away from the Galaxy center than the stars that makeup component A.

When did these two components become part of the boxy structure? A first possibility is that, components A and B were both present in the MW disk at the time of bar formation. Because of the chemical characteristics of these two components, and because of the fact that A must have been, on average, more centrally concentrated than B, this would imply the existence of a negative [Fe/H] gradient and of a positive [α/Fe] gradient in the disk at the onset of the bar formation and evolution phase of the Galaxy. Note that, with a difference Δ_{[Fe / H]} = ⟨ [Fe / H] ⟩ _A − ⟨ [Fe / H] ⟩ _B ~ − 0.4 dex over most of the extent of the disk, the resulting radial gradient would be significantly flatter than the −0.4 dex/kpc needed in the work of Martinez-Valpuesta & Gerhard (2013) to reproduce the vertical bulge abundance gradient observed in the MW. In our opinion, a radial gradient as steep as that proposed by Martinez-Valpuesta & Gerhard (2013) is unrealistic: it would imply that the typical metallicity of the MW disk at 4–5 kpc from the Galaxy center was between −1 and −1.4 dex. Such low metallicities are never reached in the MW thin disk (see, for example Fuhrmann 1999; Bovy et al. 2012; Haywood et al. 2013), but are typical of the metal-poor tail of the thick disk and of the MW halo (see, among others, Beers & Sommer-Larsen 1995; Reddy & Lambert 2008; Nissen & Schuster 2010).

The second possibility is that, when the bar formed, the young thick disk was already in place, while the thin disk had not yet formed. This second scenario is compatible with a formation epoch for the MW bar around z ~ 1, the epoch of bar growth in external galaxies (Sheth et al. 2008; Melvin et al. 2013). This age also corresponds to the transition epoch between the thick and thin disk formation at the solar vicinity and a significant dip in the star formation rate (Haywood et al. 2013; Snaith et al. 2014). In this scenario, the young thick disk was captured first in the bar instability, forming component B of the bulge. Subsequent star formation formed the thin disk, which was then also involved in the formation of the peanut structure: stars close to the vertical resonance were scattered to great heights from the plane, contributing to the thick part of the peanut (in agreement with stars of component A being detected at high latitudes as well), while the stars from the inner disk contributed to the formation of population A close to the mid-plane. To differentiate between these scenarios will require N-body models which include realistic treatments of gas dynamics, star formation and chemical enrichment. We note, however, that stars younger than 8−9 Gyr seem to be present in the Galactic bulge (Bensby et al. 2013). If confirmed, this result will favor this scenario of continuous enrichment of the Galactic bulge by younger disk populations.

4.1.2. Component C

Of the three main components, C is the only one in the bulge that does not participate in contributing to the boxy/peanut structure. Its kinematics is indeed compatible with a dynamically hotter population, characterized by radial velocity dispersions nearly constant with longitude and latitude (Ness et al. 2013b). As we will discuss in the following section, it is difficult to reconcile the chemical and kinematic properties of component C with a classical bulge. The other possibility is that stars that comprise component C are stars which are part of the old thick disk – a population that was not significantly influenced by the buckling of the bar because of its high velocity dispersions. Consistent with this picture is the fact that its relative contribution to the bulge increases with latitude as would be expected if it originated from a component with a larger scale height than components A and B (see the discussion in Ness et al. 2013b). Even if the models discussed in this paper do not include this component, such a scenario seems compatible with the N-body simulations which do include a thick disk and which show that indeed such a kinematically hot component is much less sensitive to being influenced by bar instabilities (Bekki & Tsujimoto 2011). While further studies taking into account both thick and thin disk populations are needed to quantify the evolution of these thicker components in the presence of a stellar bar, and thus to validate the suggested scenario, we emphasize that the chemical properties of population C (mean [Fe/H] and [α/Fe]) are also characteristics of old thick disk stars within the solar vicinity (Haywood et al. 2013). Similarly, the high rotational support of component C, which is similar to that of component A, also supports a (in situ) disk origin for component C.

4.2. Is there any classical bulge in the MW?

The possibility that the MW hosts a classical bulge, hidden within the overall boxy/peanut-shaped structure, has been debated for at least the last two decades (Zhao et al. 1994; Soto et al. 2007; Zoccali et al. 2008; Babusiaux et al. 2010). Some general arguments should be considered in this debate. First, if a classical spheroid is underlying the MW bulge, how massive can it be? For a MW-type galaxy, the average bulge-to-total light ratio, B/T, is ~0.1 (Laurikainen et al. 2007; Weinzirl et al. 2009), with typical values between 0 and 0.2 (Weinzirl et al. 2009). This implies bulge-to-disk ratios between 0. and 0.25, with an average value of about 0.11. However, the bulge-to-disk ratio depends on the environment as well. Kormendy et al. (2010), for example, remarked that out of the 17 giant disk galaxies found within 8 Mpc from the MW, 15 are compatible with being bulgeless or having very small bulges. Extending the sample to 11 Mpc, Fisher & Drory (2011) confirm these results, with 80% of disk galaxies being either bulgeless or hosting a pseudobulge. This suggests that, for the MW, a B/D greater than 0.1 is perhaps unlikely.

In any event, if the Galaxy hosts a classical bulge, the expected metallicity of this component can be estimated from the results of bulges in other galaxies. Assuming a MW stellar mass of about 5 × 10¹⁰M_⊙ (McMillan 2011), the mass-metallicity relation suggests typical [Fe/H] values around −0.3−0.4 dex, for a spheroid ~10¹⁰M_⊙ (Gallazzi et al. 2005; Thomas et al. 2010). Specifically, from Gallazzi et al. (2005), one sees that for a ~10¹⁰M_⊙ spheroid, a metallicity below − 0.7–− 0.8 dex is unlikely, with about only 15% of galaxies with those masses having such a low metal content. Another independent approach is to use the i-band absolute magnitude versus metallicity relation of classical bulges given by Zhao (2012). For an absolute B-band of the MW bulge M_B = − 17.65 (Kormendy & Gebhardt 2001), and thus an i-band absolute magnitude necessarily lower than this value, one obtains [Fe/H] values generally not below − 0.35–− 0.4 dex. Both arguments suggest that the classical bulge, if present, likely would have metallicities compatible with the component B of Ness et al. According to our models, at B/D = 0.25, a spheroid this massive would imply that a substantial fraction of stars (between 30% and 50%) at high latitudes should have kinematic characteristics typical of a spheroid⁷ (Fig. 11) but the kinematics of component B excludes this possibility.

Could the classical spheroidal population be masquerading as component C or hidden among its stars? For a number of good reasons, this is an appealing idea. For example, component C shows the trends expected for a rather massive bulge: its contribution increases with latitude and it does not contribute to the boxy/peanut-shaped morphology of the bulge. However, the mass-metallicity relation of bulges does not support this hypothesis: a bulge metallicity around −0.7 dex is typical of 10⁸M_⊙ spheroids (Lee et al. 2008; McConnachie 2012; Kirby et al. 2013). Such a mass would be an absolutely negligible contribution to the bulge for a galaxy with mass and size of the MW, or conversely, in order to maintain the necessary B/D ratio of the MW, the MW would have to have an unrealistically low mass of only about a few ×10⁹M_⊙.

From these considerations, we favor a scenario where a classical bulge, if present, is hidden in component B, and it is small (B/D ≲ 0.1). The results presented in Sect. 3.4 suggest that a possible signature of this classical component may be left in the radial velocity distribution of the stars. We have shown, indeed, that even if the classical component acquires some rotation during the secular evolution, its rotational support is always lower than that of the disk counterparts. Moreover, the difference in the rotational support of these two components increases with longitude. As a result, the presence of a low velocity tail in the radial velocity distribution whose strength increases with longitude, may be indicative of the presence of a small classical bulge component.

But is a classical spheroid really needed to explain the characteristics of the MW bulge? The works cited at the beginning of this section (Zoccali et al. 2008; Babusiaux et al. 2010; Hill et al. 2011) all recognized that the complex kinematics and chemistry found in the MW bulge is at odds with the presence of a single stellar population with bar-like kinematics. They suggested that the observed trends may be explained by a double-component bulge, with an inner boxy-peanut structure formed by vertical instabilities, and an outer classical bulge. However, is this component really a classical spheroid? Babusiaux et al. (2010) noted that, at high latitudes, the total radial velocity dispersion of stars along the bulge minor axis is dominated by metal-poor stars that are dynamically hot. In other words, it is the metal poor component that dominates the observed trend of the radial velocity dispersion at high vertical distances from the Galaxy plane (e.g. Fig. 2 in Babusiaux et al. 2010). Similar trends have been found in models where the classical bulge becomes the dominant component at high latitudes and determines the global trend of the velocity dispersion (Fux 1999). However, in the Fux model, the classical bulge is very massive – in fact, too massive for a MW-type galaxy – the bulge-to-disk ratio being of the order of 0.5. With a lower bulge mass, around 10% of that of the disk, consistent with that of the MW, it would not have been possible to explain the kinematic trends observed in the Galactic bulge: in the absence of a thick disk component, indeed, at all latitudes, disk stars with bar-like kinematics would have been the dominant component to the local stellar mass content, rather than stars associated with a classical spheroid (Fig. 11).

This implies that the metal-poor, α −enhanced, dynamically hot population dominating the MW bulge at high latitudes must have a different origin. The arguments presented at the beginning of Sect. 4 point indeed to a different scenario. There is a massive stellar component (that could be about half of the present-day MW stellar mass) already in place by z ~ 1. Even assuming that an old (i.e. ages >8 Gyr) classical spheroid about one tenth of the current MW disk mass exists in the Galaxy, this would imply that we are still left with about 80% of the stellar mass formed at z ≥1⁸ that is not in a classical spheroid. As a consequence, it is this massive, disk-like, metal-poor, α-enriched component – the thick disk – the one that is impossible to exclude from consideration when interpreting the observed trends in bulge populations with vertical distances within the MW bulge. We know the thick disk exists and it must be contributing to the characteristics of the bulge.

5. Conclusions

By means of idealized, dissipationless N-body simulations that follow dynamical influence of the formation and subsequent buckling of a stellar bar on disk stars, we have studied the formation and characteristics of boxy/peanut-shaped bulges and compared them with the properties of the stellar populations of the MW’s bulge. The main general results of our modeling, valid for the general family of boxy/peanut shaped bulges, are the following:

• 1.

  Because of the radial migration initiated at the time of the bar formation, boxy bulges are populated by stars born both in the inner disk and in the outermost regions of the disk, up to the OLR. That is, it is essentially the entirety of stellar disk that is mapped into the stellar populations that comprise the boxy bulge.

• 2.

  Stars formed outside the bar radius can constitute a non-negligible fraction of the stars that currently comprise the boxy bulge, in our modeling up to 30% of the stellar mass at high latitudes.

• 3.

  The fraction of outside-in migrators in a boxy bulge increases both with galactic latitude and longitude.

• 4.

  The contribution of stars to the local bulge density depends on their birth radius: stars born in the inner disk tend to stay confined in the innermost regions of the boxy bulge, while stars born close to or beyond the VILR tend to populate the more extended regions of the boxy/peanut-shaped structure.

• 5.

  The stellar birth radii are imprinted on the stellar kinematics of the bulge stars: bulge stars with larger birth radii have higher levels of rotational support and line-of-sight velocity dispersions compared with stars with smaller birth radii (but note that the trends in the line-of-sight velocity dispersion depend on the angle at which the bar is viewed).

• 6.

  If a classical spheroid is hidden among the stellar components of a boxy bulge, then our modeling and empirical relationships between spheroidal mass and size and metallicity indicate the following:

  • Because of the existence of a mass-size relation for spheroids, its contribution to the local density of the boxy bulge depends on its mass: the larger the classical bulge-to-disk ratio, the greater is its fractional contribution to the stellar densities at high vertical distances from the galaxy mid-plane.

  • A boxy bulge which contains a small classical spheroid (B/D = 0.1) is dominated everywhere by stars originating in the disk. In our models, the classical spheroid indeed contributes, at most, about 20% of the local mass density regardless of galactic longitude or latitude.

  • For a more massive classical spheroid (B/D = 0.25), its contribution to the boxy bulge mass density stays constant or increases with latitude. In the absence of a thick disk component, our models predict that the ratio between the classical spheroid and disk stars mass can become equal to 1 at high latitudes, that is the classical spheroid becomes a non-negligible fraction of the stellar content of the boxy structure far from the mid-plane.

  • Even if classical spheroids acquires some rotational angular momentum during the secular evolution of the bar (see also Saha et al. 2012, 2013), their rotational support is generally smaller than that of disk stars at the same location: classical spheroids may be thus be revealed as a low velocity tail in the line-of-sight velocity distribution of stars in the boxy bulge, whose strength increases with longitude.

Comparing these results with the properties of the stellar populations in the MW bulge obtained by the ARGOS survey, we conclude that:

• 1.

  The two most metal-rich components of the MW bulge have a disk origin, with component B formed, on average, at larger radial distances than component A.

• 2.

  Because of their chemical and kinematic characteristics, we suggest that component A and B are dominated by the Galactic inner thin disk and by the young (ages ≤10 Gyr; see Haywood et al. 2013) thick disk, respectively.

• 3.

  Because of the existence of a mass-metallicity relation for spheroids, and because of the properties of extra-galactic classical bulges, it is difficult to associate component C with a classical spheroid. Its high level of rotational support, as deduced by radial velocity measurements, suggests a disk origin for this component as well. On the basis of its chemical characteristics, we suggest that it is indeed associated to the old galactic thick disk (ages in the range 8−10 Gyr; see Haywood et al. 2013).

• 4.

  The presence of a massive classical spheroid, with B/D ~ 0.2 can be excluded for the MW. If present, such a massive component would indeed be significant at high latitudes (30−50% of the local stellar density) and should have metallicities comparable to those of population B in the bulge. The bar-like kinematic and morphology of population B excludes this possibility. As a result, if a classical bulge is hidden in the populations of the boxy-peanut structure, it cannot be massive (B/D ≤ 0.1). This result is in agreement with those found by Shen et al. (2010); Kunder et al. (2012), on the basis of kinematic arguments, thus supporting the scenario that most of the mass of the MW’s bulge has a disk origin.

Many observational arguments discussed in this paper support the presence of a massive thick disk in the inner Galactic regions. Such a component cannot be neglected when tallying the mass budget of the MW bulge and inner disk. Modeling its interplay with and contribution to the bar and the thin disk will be the focus of future studies.

Online material

Appendix A: Exploring different galaxy models

The galaxy simulations presented in Sects. 2 and 3 develop a bar whose length (r_bar = 7 kpc) is about half of the initial disk size (13 kpc). In these simulations, the corotation radius and the OLR are located at 7 − 8 kpc, and 13 kpc, respectively. Since the initial disk does not extend further, it is natural to investigate what is the response and redistribution of stars to bar formation and evolution that are initially at even larger radii than in our original simulations. In this Appendix, we have analyzed two additional simulations which are part of the sample of simulations described in Hallé & Combes (2013). These simulations include gas physics, implementing a cooling prescription which allows the gas to cool low temperatures and the formation of a molecular hydrogen component, together with star formation and several different stellar feedback efficiencies (see Hallé & Combes 2013, for details). For the purpose of the comparison to the models presented in this paper, we have chosen two simulations from this suite without molecular hydrogen and alternatively, with and without including stellar feedback. We have chosen these two simulations because they develop stellar bars with different characteristics, thus allowing us to investigate the impact of those characteristics, and specifically, the impact of the location of the associated resonances, on the overall disk evolution. We summarize the main features of the simulated galaxies for these two simulations from Hallé & Combes (2013): they consist of a gSb-like galaxy, composed of a stellar disk whose mass is M_d = 4.5 × 10¹⁰M_⊙, and whose scale length is a_d = 5 kpc; a classical bulge whose mass is 0.25 M_d and characteristic radius a_b = 1 kpc; and a gaseous disk whose mass is 0.2 M_d, and whose scale length is 11.8 kpc. The gaseous and stellar disk components follow a Miyamoto-Nagai density profile while the stellar bulge and dark halo components have Plummer density profiles (see Hallé & Combes 2013, for a complete description of the models). In particular, in these models the initial stellar disk extends to 36 kpc from the galaxy center, thus allowing us to trace

the evolution of stars up to distances of ~8 disk scale lengths. The total number of particles employed is 1 200 000, equally distributed in number as gas, stars and dark matter particles.

Both models develop prominent stellar bars, whose length is about 7 kpc, at 1 Gyr after the start of the simulation. At about t = 5 Gyr, the corotation and OLR are located at r = 10 kpc and r = 17 kpc respectively for the simulation with stellar feedback (Fig. A.1), while they are located slightly further out in the simulation without stellar feedback (at r = 12 kpc and r = 20 kpc; Fig. A.2). This difference is a result of the fact that the simulation without stellar feedback has a lower bar pattern speed than the simulation with stellar feedback. As was the case in our original simulations (Sect. 3), the formation of the bar is accompanied by a significant redistribution of the stars in the disk, with stars initially as far as 15−20 kpc, reaching the inner regions in less than 2 Gyr. These outside-in migrators become part of the stellar bar, as is evident in their asymmetric distribution in the inner 5−10 kpc which tracks the orientation of the bar. However, it is also interesting to note that the outermost regions of the simulated disks, those at distances greater than 15 kpc (Fig. A.1) and 20 kpc (Fig. A.2) in the two simulated galaxies, stay mostly confined to the outer disk, with only less than 1/1000 of those stars reaching the bar region, after 9 Gyr in the simulations. The transition radius where the different responses of the disk to the spatial redistribution seems to occur is approximately the radius of the OLR. In particular, comparing Figs. A.1 and A.2, one see that the greater the radius of the OLR, the larger the disk region involved in contributing stars to the outside-in migration. As a corollary to this, it also means that the regions where stars do not migrate inward also occurs beyond a larger radius. While the study of this behavior will be the subject of subsequent studies, the point we wish to emphasize here is that in the MW all kinematically cold populations inside the solar circle, where the MW bar OLR is probably located (see, for example, Gerhard 2010), may have contributed to make up the Galactic boxy/peanut-shaped bulge.

---

Acknowledgments

All the simulations have been performed on CURIE at CCRT, CEA, within the framework of the “Projet Grand Challenge” and with support provided by the GENCI grant x2012040507. The authors acknowledge the support of the French Agence Nationale de la Recherche (ANR) under contract ANR-10-BLAN-0508 (GalHis project). PDM warmly thanks M. Ness, O. Gerhard, and L. Origlia for comments, suggestions and criticisms on the content of this paper, and C. Babusiaux for remarks on a first version of this manuscript. We are grateful to the anonymous referee for a careful report which helped us to improve the presentation and clarity of our results.

References

All Tables

All Figures