Insight into the Galactic Bulge Chemodynamical Properties from Gaia DR3

The physical processes involved in the formation of the Galactic bulge encode important information about the formation and evolution history of the Milky Way (MW) (Barbuy et al., 2018, for review). In the early stage of galaxy formation, the bulges form through merging and dissipative collapse (Eggen et al., 1962) and later evolve through secular evolution. Severe dust extinction strongly affects photometric observations of our galaxy’s bulge region. However, infrared observations reveal that the MW hosts an elongated triaxial structure with the near-side in the first quadrant (Dwek et al., 1995). Early evidence for the existence of the bulge also exists from gas kinematics (Binney et al., 1991) with later observations revealing that the bulge rotates cylindrically (Howard et al., 2009; Kunder et al., 2012; Zoccali et al., 2014), a signature of a rotating rigid body. Our current view is that the MW hosts a boxy/peanut-shaped bulge or pseudo-bulge (Maihara et al., 1978; Weiland et al., 1994; Dwek et al., 1995; López-Corredoira et al., 2005), consistent with a buckled bar formed from the secular evolution of a disk. The presence of a small classical bulge in the MW, i.e., a merger-induced spherical structure similar in many aspects to a mini-elliptical galaxy, has not been excluded but its mass is likely less than 8% of the disk mass (Shen et al., 2010).

The Galactic bulge exhibits complex chemodynamical trends that have been revealed by spectroscopic surveys such as the BRAVA (Rich et al., 2007; Kunder et al., 2012), ARGOS (Freeman et al., 2013), GIBS (Zoccali et al., 2014), GES (Rojas-Arriagada et al., 2014, 2017), APOGEE (Majewski et al., 2017; Jönsson et al., 2020), the PIGS (Arentsen et al., 2020) and the BDBS survey (Lim et al., 2021). The metallicity distribution function in the bulge region spans a wide range of metallicities and appears to have multiple peaks (Rich, 1988; McWilliam & Rich, 1994; Hill et al., 2011; Fulbright et al., 2006; Zoccali et al., 2008; Ness et al., 2013a; Gonzalez et al., 2015; Ness & Freeman, 2016; Zoccali et al., 2017; Rojas-Arriagada et al., 2014, 2017, 2020), indicating the existence of different stellar populations in the bulge. Several works have confirmed a negative vertical metallicity gradient in the bulge (Minniti et al., 1995; Zoccali et al., 2008; Ness et al., 2013b; Hayden et al., 2015). Ness et al. (2013b) attributes this vertical metallicity gradient to the varying fraction of different stellar populations, which have almost fixed mean metallicity and metallicity dispersion, with distance from the galactic plane.

Galactic archaeology studies the kinematical and chemical properties of stars in order to reconstruct the history of formation of our Galaxy since the cosmic dawn. In the Galactic bulge region, chemo-kinematic studies have revealed that the different stellar populations, separated by their metallicity, have different kinematics (Ness et al., 2013b; Clarkson et al., 2018; Arentsen et al., 2020; Wylie et al., 2021). Stars with [Fe/H] $>-1$ dex are believed to be originated from disk material as they exhibit comparable cylindrical rotation and low velocity dispersion (Kunder et al., 2012; Ness et al., 2013b; Di Matteo, 2016; Fragkoudi et al., 2018) and there seems to be a chemical similarity between these bulge stars and thick disk stars (Rojas-Arriagada et al., 2017). Metal-poor stars with [Fe/H] $\lesssim-1$ dex, which exhibit high velocity dispersion (Arentsen et al., 2020; Dékány et al., 2013), have an unclear origin: they may be either formed in situ, e.g. ancient in-situ stars (White & Springel, 2000), halo interlopers (Lucey et al., 2020), Aurora (Belokurov & Kravtsov, 2022, i.e., ancient in-situ stars) or ex-situ e.g., remnants of early accretion events (Horta et al., 2020). More recently, Marchetti et al. (2023) used 2.3 million red clump stars in the Blanco DECam Bulge survey (BDBS) to study the chemo-kinematic properties of the Galactic bulge. They found that metal-rich stars ([Fe/H] $>-0.5$ dex) show bar signatures, with the metal-poor stars exhibiting isotropic motions.

The classic picture of bulge formation requires a dynamically cold disk in which a bar instability can be triggered. Recent numerical simulations were able to demonstrate that a bulge can form also from a hot thick disk (Ghosh et al., 2023). Disks at high redshifts ($z\sim 4-5$) have been observed by ALMA (Roman-Oliveira et al., 2023; Parlanti et al., 2023) and JWST (Ferreira et al., 2023) and even bars have been clearly observed at $z>2$ (Guo et al., 2023; Le Conte et al., 2023; Costantin et al., 2023). At early times, the emergence of stellar bars can be rapid when dark matter fraction in the centres of disk galaxies is low (Bland-Hawthorn et al., 2023). In this work we will further show, using high redshift snapshots of a Milky-Way like simulation, that it is possible to form a bar immediately after the last major merger $\sim 10$ Gyr ago when the disk was still dynamically hot.

The European Space Agency’s Gaia satellite has revolutionized our understanding of the MW’s past, with the second and third data releases (DR3) revealing plenty of previously unknown galactic features (Gaia Collaboration et al., 2018, 2021; Antoja et al., 2018). Gaia’s astrometric measurements combined with radial velocities from the Radial Velocity Spectrometer (RVS) or ground-based spectroscopic surveys enabled the study of individual stars in the orbital–chemical space, towards the Galactic bulge (Queiroz et al., 2021). Using Gaia and APOGEE chemical abundance measurements, Belokurov & Kravtsov (2022) conducted a deep analysis of the Galactic disk in the $V_{\phi}-{\rm[Fe/H]}$ space. In the aluminium-selected ‘in-situ’ population, they identified an ancient metal-poor population that is kinematically hot with an isotropic velocity ellipsoid and modest net rotation (dubbed ‘Aurora’), and a subsequent ‘spin-up’ phase of the disk stars which culminates with the formation of the Splash. The Splash (Belokurov et al., 2020) is a relatively metal-rich population ([Fe/H] $>$ -0.7 dex) with low azimuthal velocity dispersion (50-60 $\rm km\ s^{-1}$), likely originating from the proto-disk, heated by the Gaia Sausage Enceladus (GSE, Helmi et al., 2018; Belokurov et al., 2018) merger event. The presence of a kinematically heated disc component in the halo as an imprint left by an accretion event had already been suggested by Di Matteo et al. (2019).

In this work, we aim to look for evidence of the Splash and the bar/bulge in the Galactic bulge region, in an attempt to decipher the origin of the chemodynamical properties of the Galactic bulge. For this goal, we use a recently released catalog based on Gaia-XP spectra that reaches the bulge (Andrae et al., 2023). We describe the catalog in Section 2 and analyze the observations in Section 3, with a focus on the $V_{\phi}-$[Fe/H] space. In Section 4,5, we employ the Auriga (Grand et al., 2017) cosmological simulation and three $N-$body simulations to better understand the chemodynamical properties of the bulge region as revealed by our sample. In Section 6 we discuss our results and in Section 7 we summarize our findings.

In this work, we use a catalog of 13.3 million red giants with stellar metallicity from Andrae et al. (2023). The metallicity is derived from the low-resolution XP spectra in Gaia DR3 (Collaboration et al., 2022a) using the data-driven method XGBoost, that is trained on 510,413 APOGEE stars (Majewski et al., 2017) augmented by 291 ultra-metal-poor stars from LAMOST (Li et al., 2022). According to Andrae et al. (2023), the typical uncertainty on metallicity is 0.1 dex. We select stars with radial velocity errors smaller than 2.0 $\rm km\ s^{-1}$ and parallax/parallax_error > 5, which results in a sample of 5,497,737 stars (sample A). We use the photogeometric distances calculated by Bailer-Jones et al. (2021) to complete the 6D phase-space information. The six-dimensional Gaia observables are transformed to Galactocentric cylindrical coordinates using Astropy¹¹1https://www.astropy.org/index.html (Robitaille et al., 2013; Collaboration et al., 2018, 2022b) with the Sun position at ($X$, $Y$, $Z$) $=$ (-8.34, 0, 0.027) kpc (Reid et al., 2014) and peculiar motion ($V_{X}$, $V_{Y}$, $V_{Z}$) $=$ (11.1, 12.24, 7.25) $\rm km\ s^{-1}$ (Schönrich, 2012) in a right-handed frame. The spatial distribution of this sample is shown in Fig. 1. From the figure it is obvious that the Gaia data only covers the near side of the bulge (distance to the Sun less than 8.3 kpc) and almost avoids the low latitude region. We do see an overdensity in the bulge region within the red circle, confirming the coverage of the bulge stars in our sample.

For consistency checks, in Fig. 2 we compare the metallicity, azimuthal and radial velocities of stars common between the Gaia and APOGEE DR17 surveys. We also compare the APOGEE-based StarHorse distances (Queiroz et al., 2020) with Gaia-based Bayesian distances (Bailer-Jones et al., 2021) in the upper right panel of the figure. The radial velocities, azimuthal velocities and metallicities are in good agreement while the Gaia distances are slightly under-estimated compared to the StarHorse distances. As mentioned earlier, Andrae et al. (2023) primarily trained the Gaia XP spectra on APOGEE, therefore we expect good agreement between the two surveys for the common stars. The cylindrical azimuthal velocities are computed using Gaia proper motions for both surveys, but Gaia and StarHorse-derived distances respectively. The bulge sample we analyze in the following sections is selected to be within 5 kpc of Galactocentric radius from the GC (see red circle in Fig. 1), decreasing the sample A size to 330,414 stars (sample B).

In this section we make use of a cosmological zoom-in simulation of a MW-like galaxy from Auriga to understand the observed pattern in the $V_{\phi}$-[Fe/H] space and gain insight into the formation and evolution history of the Galaxy.

Auriga is a suite of 30 magneto-hydrodynamical cosmological zoom-in simulations of MW mass dark matter halos with varying mass from $\sim 1\times 10^{12}$ to $2\times 10^{12}$ $M_{\odot}$ and evolving from redshift 127 to the present day, incorporating active galactic nuclei feedback, star formation, magnetic fields (see Grand et al. 2017 for details). The Auriga simulations (Grand et al., 2024) are publicly available²²2https://wwwmpa.mpa-garching.mpg.de/auriga/data. We choose Auriga 23 (Au23) which resembles many of the MW properties. For example, it displays clear $[\alpha$/Fe] bimodality for the thin/thick populations (Grand et al., 2018) and its last major merger event happened at a lookback time $t_{\mathrm{lb}}\approx$10 Gyrs, similar to the GSE merger event (Montalbán et al., 2021; Xiang & Rix, 2022). Additional information about Au23 regarding its star formation history, bar morphology, rotation curve, and merger history, can be found in Fragkoudi et al. (2020)

As we have shown in the previous section, chemical information is not enough to distinguish the bulge stellar populations; the addition of $V_{\phi}$ is crucial to the identification of the non-bar population with negative $V_{\phi}$ and high radial velocity dispersion in the metallicity range $-1<$[Fe/H]$<-0.5$, which we associate with the Splash. Belokurov et al. (2020) suggested that it is constituted by disk stars heated by the GSE merger in the early epoch. To study the effects of the merger on the formation of the Splash and the bar in the Au23 simulation, we select stars of different ages in the bulge region, dividing them into three epochs: before the last major merger (12 $-$ 13 Gyrs), around the time of the merger (9.6 $-$ 10.4 Gyr) and after the merger (8.5 $-$ 8.9 Gyrs). The last major merger of Au23 happened at around 10 Gyrs (see Fragkoudi et al., 2020). Thus, stars of ages at 9.6 $-$ 10.4 Gyr are mainly contributed by the starburst population and dynamically hot disk.

The face-on and edge-on distributions of stars in the inner 5 kpc of Au23 are shown in Fig. 8 (top and bottom row respectively) at $t_{lb}=$ 8.38 Gyr (left column) and 0 Gyr (right column). Stars born before the merger (filled blue contours) resemble an oblate spheroid appearing more circular in the face-on view and elliptical in the edge-on view. The stars born around the merger (black contour lines) and after the merger (red contour lines) contribute to the bar structure. The bar containing stars formed around the merger is slightly shorter and thicker compared to the bar containing stars formed after the merger, consistent with results in literature that stars with hotter kinematics will have thicker and rounder bars than stars with colder kinematics (Fragkoudi et al., 2017; Debattista et al., 2017; Athanassoula et al., 2017). According to Fig. 8, in the Au23 simulation, the bar forms quickly after the last major merger when the disk is still dynamically hot. The stars born around and after the merger constitute the thicker and thinner bar, respectively, while the stars born before the merger did not participate in the bar formation process to resemble an oblate structure.

To compare the data (sample B) with the simulation, we narrow down the selection of simulation particles to $2<R<5$ kpc and $|Z|<3$ kpc, to avoid the innermost region with complex kinematics and halo stars. The row-normalized density maps in the $V_{\phi}-$[Fe/H] space are shown in Fig. 9 for the same lookback times as Fig. 8. The stars born around the merger (i.e., the middle column) exhibit the chevron-like pattern with a vertical extension towards lower $V_{\phi}$, very similar to the one in observations (Fig. 5). In the MW, the top and bottom branches of this chevron-like pattern are attributed to the thin and thick disks, respectively (Belokurov & Kravtsov, 2022; Lee et al., 2023). However, as shown in Fig. 9, in the Auriga simulation, such a chevron-like pattern has already emerged 8.38 Gyr ago, after the last major merger. In the left column of Fig. 8, the spatial distribution of these stars (black contours) resembles a thick disk with a clear bar structure. This seems to indicate that the chevron-like pattern can be produced even without a significant fraction of thin disk component. The stars born after the merger (top right panel of Fig. 9) distribute within the upper right region at high metallicities and high $V_{\phi}$, which is expected since they are born in the disk (see also red contours of Fig. 8, left column). In Fig. 9 comparing the $V_{\phi}$-[Fe/H] distributions of stars born before (left column) the merger at $t_{lb}=8.38$ Gyr (top panels) and $t_{lb}=0$ Gyr (bottom panels), there is no clear variation after 8 Gyr of evolution. However, the stars born after the merger (right column) show significant difference in the $V_{\phi}$-[Fe/H] pattern during the same time period: the pattern shifts to lower azimuthal velocities and has a much larger velocity dispersion. Both the external and internal mechanisms could have contributed to the disk heating. After the last major merger occurred about 10 Gyrs ago, the galaxy experienced subsequent minor mergers (see Fig. 11 in Fragkoudi et al., 2020) which can also contribute to disk heating. In addition, as shown in Fig. 8, the stars born after the merger, host a strong bar. The formation of the bar requires loss of angular momentum which leads to larger velocity dispersion.

We now investigate the $V_{\phi}-$[Fe/H] space color$-$coded by the radial velocity dispersion $\sigma_{V_{R}}$ in Fig. 10 and we compare it to the data (bottom panel of Fig. 5): we notice that stars born around the merger (middle column in Fig. 10) have a distribution similar to the observations. The velocity dispersion decreases monotonically from lower left corner to the top right corner, consistent with the observed pattern in Fig. 5. For the stars born before the merger (left column), $\sigma_{V_{R}}$ is relatively uniform at any given [Fe/H], which is consistent with the behaviour of a dispersion$-$dominated structure. There is not a much difference between the $V_{\phi}-$[Fe/H]$-\sigma_{V_{R}}$ patterns at $t_{lb}=8$ (top row) and $t_{lb}=0$ (bottom row) for the stars born before (left panels) and around the merger (middle panels), confirming the findings in Fig. 9: 8 Gyrs of evolution (e.g. the accretion of cold gas, the following disk formation and dynamical heating from secular evolution) do not greatly influence the morphology and kinematics of the stars formed early-on. On the other hand, the younger stars born just after the merger (right column) which are kinematically cold, experience significant heating over the same time period (see bottom right panel).

We now focus on the stars born around the time of the merger in Au23. We separate them into three sub-populations,MPD (metal-poor disk), MRD (metal-rich disk) and the Splash (see the three boxes in Fig. 9) to mimic the observed transition from ‘disk’ (top three panels of fourth column in Fig. 7) to ‘disk coexisting with Splash’ (the middle column in Fig 7) seen in the bulge Gaia sample. The Splash and MPD have the same metallicity range but the Splash is kinematically hotter with weaker rotation than MPD (see Fig 10). MPD and MRD have similar $V_{\phi}$ distributions but MRD is more metal-rich. The face-on and edge-on distributions of these three sub-populations are presented in Fig. 11 at $t_{lb}=8.38$ Gyr (blue iso-density contour maps) and 0 Gyr (green). The corresponding velocity maps are shown in Fig. 12 for the face-on projection. At $t_{lb}=8.38$ Gyr, after the last major merger when the bar just formed, the Splash is nearly spherical with the inner 1 kpc showing a bar-like morphology. MPD and MRD both contain a bar structure within 2 kpc from the GC, with the bar in MPD being slightly thicker than the bar in MRD. However at redshift 0, both MPD and MRD have evolved to have almost the same morphology, including a similar-size bar/bulge. The Splash component also develops a weak bar ($\sim 2$ kpc half length) following the bar orientation in MPD and MRD. However, as shown in the $\left\langle V_{R}\right\rangle$ color-coded maps in Fig. 12, the Splash does not display the typical bar-like butterfly pattern seen in MPD and MRD with a very prominent bar in the inner 2 kpc – this is consistent with the data shown in Fig 7. Therefore, the Splash stars likely do not contribute to the bar supporting orbits. The bar-like appearance of the Splash in the inner 2 kpc at $t_{lb}=0$ is likely due to the gravitational attraction from the bar. Nonetheless, by looking at the mean $V_{R}$ maps alone, we cannot exclude the possibility that a few of the Splash stars indeed move on bar supporting orbits. Beyond 2 kpc, the Splash morphology becomes more spheroidal, resembling a ‘fluffy classical bulge’ (see Fig. 11 in Belokurov et al. 2020).

To investigate the influence of a rotating bar on a pre-existing non-rotating spheroidal structure (i.e., a classical bulge) and the origin of the observed cylindrical rotation for the metal-poor stars in the Galactic bulge region, we make use of three isolated $N$-body simulations of MW-like galaxies that self-develop a bar out of disk particles. We set up the initial condition (IC) following Tepper-Garcia et al. (2021), with a dark matter halo ($M_{h}\sim 1.2\times 10^{12}$ $M_{\odot}$), a stellar disk ($M_{d}=4.3\times 10^{10}$ $M_{\odot}$) with a scale length of 2.5 kpc and a scale height of 0.3 kpc (resembling a thin disk) and a non-rotating stellar spheroidal structure - the pre-existing classical bulge. We build three models, namely, M1, M2 and M3, each with different classical bulge mass $M_{b}=0.1M_{d},0.2M_{d}$ and $0.3M_{d}$. The dark matter halo, disk, and bulge in all models are approximated by $10^{6}$, $10^{6}$, and $6\times 10^{5}$ particles, respectively. More simulation details and the final resemblance with the MW can be found in Tepper-Garcia et al. (2021). The ICs are generated using iterative self-consistent modelling module in AGAMA (Vasiliev, 2019), which are then evolved with GADGET4 (Springel et al., 2021) for 5 Gyrs. A bar is fully formed and buckles into a peanut-shaped bulge within $\sim$1.5 Gyr for M1, the model with the smallest bulge mass. The formation time of the bar is delayed for M2 and M3 with larger bulge mass. This is consistent with the findings of Li et al. (2023).

The temporal evolution of the morphology and kinematics of M1 is shown in Fig. 13, for the disk (left two columns) and non-rotating spheroidal component (right two columns). The number density and mean $V_{R}$ maps in the $X-Y$ plane reveal that the disk particles at 0 Gyrs (top row) have an axisymmetric distribution without any velocity pattern. Once a bar forms (middle row), the butterfly pattern appears in the $\left\langle V_{R}\right\rangle$ map due to the systematic motion of the bar orbits in different quadrants. As the bar grows, the $\left\langle V_{R}\right\rangle$ butterfly pattern becomes larger (bottom row). On the other hand, the spheroidal non-rotating component does not show the $\left\langle V_{R}\right\rangle$ butterfly pattern even at the end of the simulation, indicating that the pre-existing spheroidal bulge component is not involved in the bar formation. Instead, it morphs into a slightly elliptical structure aligning with the bar major axis, due to the gravitational attraction from the bar.

To further investigate the evolution of the bulge in a more quantitative way, we generate mock observations of the Galactic bulge from the simulations by setting the observer at $R=8$ kpc and a bar viewing angle of $20^{\circ}$. The longitudinal rotational profiles of the particles are shown in Fig. 14 for the disk (blue curves) and spheroidal non-rotating component (red curves) at different times (increasing from left to right columns). M1, M2 and M3 are shown in the top, middle, and bottom rows, respectively. The rotation profiles of the non-rotating component are initially flat, with no net rotation, but become steeper with time as it gains angular momentum from the rotating bar (Saha et al., 2012) for all three models; its morphology also changes (see Fig 13) from spheroidal to ellipsoidal, with the major axis aligned with the bar major axis. In Fig. 13, we have seen that classical bulges may not display the $\left\langle V_{R}\right\rangle$ butterfly pattern, suggesting that they are not building blocks for the bar. However, the classical bulges in all three models evolve to exhibit cylindrical rotation, which is a known property of bars and boxy/peanut-shaped bulges in $N$-body simulations (Athanassoula & Misiriotis, 2002; Saha & Gerhard, 2013). This work confirms the result of Saha et al. (2012) that cylindrical rotation is not an unique property of bars, e.g. a pre-existing classical bulge can be spun up by the bar to result in cylindrical rotation. However, the disk particles (the bar) do not strictly follow the cylindrical rotation. Although the rotational velocities at different latitude bins almost converge at $|l|\sim 10^{\circ}$, the rotation curves in inner region change with latitude and the discrepancy is more significant for models with larger classical bulge mass. In the MW, the rotation curves in the bulge region barely change with latitude (Howard et al., 2009; Kunder et al., 2012; Zoccali et al., 2014, also see Fig. 6). This result could put constrains on the mass of the MW classical bulge, which should be less than 10% of the disk mass, consistent with Shen et al. (2010).

Although our simulations are pure $N$-body simulations with no chemical information, the initially non-rotating spheroidal component with high velocity dispersion could be considered a classical bulge, formed in the early chaotic stages of the MW. In observations, it would likely correspond to the (old) metal-poor population ([Fe/H]$\lesssim-1$ dex). In our bulge sample B, stars with $-1.5\lesssim$[Fe/H] $\lesssim-1$ dex do not exhibit the $\left\langle V_{R}\right\rangle$ butterfly-pattern (Fig. 7) but show cylindrical rotation (Fig. 6), while the more metal-poor stars ($-2<$[Fe/H]$<-1.5$ dex) do not display cylindrical rotation or the butterfly-pattern. Our interpretation is that the more metal-rich stars ($-1.5\lesssim[Fe/H]\lesssim-1$ dex) could represent classical bulge spun up by the bar while the more metal-poor ones ($-2<$[Fe/H]$<-1.5$ dex) could be halo stars, which spend more time outside of the bar’s influence. Therefore, their kinematics is less easily affected by the bar presence.

Recent Solar neighbourhood studies have found that the thin disk, thick disk and the Splash appear as distinct features of the $V_{\phi}$-[Fe/H] diagram (Belokurov et al., 2020; Lee et al., 2023). In this work, we used data from Andrae et al. (2023), a catalog with metallicities extracted from the Gaia-XP spectra, to study the $V_{\phi}$-[Fe/H] diagram over a large range of radii, from the Galactic bulge to the outer disk, up to radius of 14 kpc. We find that the characteristic patterns of the disks and Spash in the $V_{\phi}$-[Fe/H], exist across the disk. The $V_{\phi}$-[Fe/H] profiles of the thin disk show systematic variation with radius, which is in agreement with the inside-out formation scenario. The $V_{\phi}$-[Fe/H] profiles at various radii for the thick disk and the Splash are quite consistent, implying their early formation timescale when the chemistry was spatially well mixed.

The bulge stars share a similar pattern in the $V_{\phi}$-[Fe/H] space with disk stars, implying that it originated from the disk. By investigating the bar signatures of stellar populations in $V_{\phi}$-[Fe/H] space and only considering the populations with reliable number statistics, we found that:

• •

  Stellar populations with [Fe/H] $\lesssim-1$ dex are dispersion-dominanted.

• •

  Stellar populations with $-1\lesssim$ [Fe/H] $\lesssim-0.4$, show bimodal distribution. The stars with $V_{\phi}\gtrsim 100$ $\rm km\ s^{-1}$ follow the bar kinematics while the stars with $V_{\phi}\lesssim 100$ $\rm km\ s^{-1}$ (Splash) have dispersion dominated kinematics similar to those more metal-poor populations.

• •

  Stellar populations with [Fe/H] $\gtrsim-0.4$ dex all have bar-like kinematics.

By analyzing a MW-like Auriga simulation Au23, we found that only the stars born around the time of the last major merger ($\sim 10$ Gyr ago), which is mainly contributed by the starburst population and the dynamically hot disk, have a $V_{\phi}$-[Fe/H] pattern similar to the observation in the Galactic bulge. A possible picture for the Galactic evolution emerges based on the great resemblance between the simulation and observations, that is, a preliminary thick disk was already in place not long before the last major merger. Afterwards, the occurrence of the last major merger heated some of the thick disk stars and triggered starburst formation that contribute to the Splash, whose lower metallicity limit is $\sim-1$ dex that is inherited from the thick disk. Thus the lower metallicity limit of the thick disk should be $\sim-1$ dex at the time of the last major merger. Subsequently, the thick disk stars with [Fe/H] $\gtrsim-1$ dex form the bar structure, their kinematics becomes bar-like and detach from those more metal-poor stars. In another aspect, the merger kicks some of thick disk stars to halo-like orbits and induces starburst formation, mainly giving rise to the Splash population, which does not participate in the bar formation. Thus its kinematics is more dispersion-dominated.

Moreover, the observed metal-poor stars ($-1.5<$ [Fe/H] $<-1$ dex) show cylindrical rotation without butterfly pattern in their mean $V_{R}$ map, possibly as a consequence of the bar losing angular momentum to a pre-existing classical-bulge like structure (Saha et al., 2012) during secular evolution. We also perform three $N$-body simulations to study the interaction between an initially non-rotating spheroidal component (a classical bulge) and a later formed bar and confirm the result of Saha et al. (2012) that an initially non-rotating spheroidal component can be spun up to develop cylindrical rotation under the bar influence without following the bar orbits. In addition, we found that the bulge mass affects the characteristic rotation profiles of the bar. The three $N$-body models are initiated with different bulge masses. Although the rotation profiles ($l-V_{gsr}$) of the bars in the three models almost converge at $l\sim 10^{\circ}$ for the different Galactic latitudes we consider, in the inner region they show relatively large discrepancies:

• •

  For the model with the least massive spheroid component ($M_{1}$), the rotation profiles of the bar at different latitudes are similar to each other, consistent with the cylindrical rotation pattern.

• •

  For the model with the most massive spheroid component ($M_{3}$), the rotation profiles of the bar at higher latitudes have much shallower slopes at $|l|\lesssim 4^{\circ}$ compared to the lower latitudes.

In the Milky Way, the rotation profiles at different latitudes are almost the same (see Fig. 6), which is most close to the model with the least massive bulge, implying that the classical bulge in MW, if there is one, should not be too large (less than 10% disk mass). In the future, we plan to consider other parameters of the classical bulge such as spin, velocity dispersion, and velocity anisotropy, etc., to quantitatively match to observations and further constrain the mass and other properties of the MW classical bulge.