Insights into the Galactic Bulge Chemodynamical Properties from Gaia Data Release 3

We explore the chemodynamical properties of the Galaxy in the azimuthal velocity V ϕ and metallicity [Fe/H] space using red giant stars from Gaia Data Release 3. The row-normalized V ϕ –[Fe/H] maps form a coherent sequence from the bulge to the outer disk, clearly revealing the thin/thick disk and the Splash. The metal-rich stars display bar-like kinematics, while the metal-poor stars show dispersion-dominated kinematics. The intermediate-metallicity population (−1 < [Fe/H]< − 0.4) can be separated into two populations, one that is bar-like, i.e., dynamically cold ( σVR∼80 km s−1) and fast-rotating (V ϕ ≳ 100 km s−1), and the Splash, which is dynamically hot ( σVR∼110 km s−1) and slow-rotating (V ϕ ≲ 100 km s−1). We compare the observations in the bulge region with an Auriga simulation where the last major merger event occurred ∼10 Gyr ago: only stars born around the time of the merger reveal a Splash-like feature in the V ϕ –[Fe/H] space, suggesting that the Splash is likely merger-induced, predominantly made up of heated disk stars and the starburst associated with the last major merger. Since the Splash formed from the proto-disk, its lower metallicity limit coincides with that of the thick disk. The bar formed later from the dynamically hot disk with [Fe/H] > − 1 dex, with the Splash not participating in the bar formation and growth. Moreover, with a set of isolated evolving N-body disk simulations, we confirm that a nonrotating classical bulge can be spun up by the bar and develop cylindrical rotation, consistent with the observations for the metal-poor stars.


Introduction
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; see Barbuy et al. 2018 for a 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 nearside 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).
Galactic archeology studies the kinematical and chemical properties of stars in order to reconstruct the history of the formation of our Galaxy since the cosmic dawn.In the Galactic bulge region, chemokinematic 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]  − 1 dex, which exhibit high velocity dispersion (Dékány et al. 2013;Arentsen et al. 2020), have an unclear origin: they may be formed in situ, e.g., ancient in situ stars (White & Springel 2000), be halo interlopers (Lucey et al. 2020), be formed in the Aurora (Belokurov & Kravtsov 2022, i.e., ancient in situ stars), or be ex situ stars, 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 BDBS to study the chemokinematic 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 also form from a hot thick disk (Ghosh et al. 2023).Disks at high redshifts (z ∼ 4-5) have been observed by the Atacama Large Millimeter/submillimeter Array (Parlanti et al. 2023;Roman-Oliveira et al. 2023) and JWST (Ferreira et al. 2023) and even bars have been clearly observed at z > 2 (Costantin et al. 2023;Guo et al. 2023;Le Conte et al. 2023).At early times, the emergence of stellar bars can be rapid when the dark matter fractions in the centers of disk galaxies are low (Bland-Hawthorn et al. 2023).In this work, we will further show, using high-redshift snapshots of an MW-like simulation, that it is possible to form a bar immediately after the last major merger ∼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 (DR2 and DR3) revealing plenty of previously unknown galactic features (Antoja et al. 2018;Gaia Collaboration et al. 2018, 2021).Gaia's astrometric measurements combined with radial velocities from the Radial Velocity Spectrometer or ground-based spectroscopic surveys have enabled the study of individual stars in the orbital-chemical space, toward 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 f -[Fe/H] space.In the aluminumselected "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 "spinup" 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 km s −1 ), likely originating from the protodisk, heated by the Gaia Sausage Enceladus (GSE; Belokurov et al. 2018;Helmi et al. 2018) merger event.The presence of a kinematically heated disk 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 f -[Fe/H] space.In Sections 4 and 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.

Data
In this work, we use a catalog of 13.3 million red giants with stellar metallicity from Andrae et al. (2023).The metallicities are derived from the low-resolution XP spectra in Gaia DR3 (Collaboration et al. 2022a) using the data-driven method XGBoost, which 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 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 6D Gaia observables are transformed to Galactocentric cylindrical coordinates using Astropy8 (Robitaille et al. 2013;Collaboration et al. 2018Collaboration et al. , 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) km s −1 (Schönrich 2012) in a right-handed frame.The spatial distribution of this sample is shown in Figure 1.From the figure, it is obvious that the Gaia data only cover the nearside of the bulge (distance to the Sun less than 8.3 kpc) and almost avoid 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 Figure 2 we compare the metallicity and 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 underestimated 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 with 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 Galactic center (GC; see the red circle in Figure 1), decreasing the sample A size to 330,414 stars (sample B).

V f -[Fe/H] Maps from Disk to Bulge Region
The morphological, kinematical, and chemical properties of the bulge are mainly shaped via its formation and evolution processes; the merger-induced classical bulge is spheroidal-like and dispersion-dominated, while the disk-instability-induced pseudo-bulge is disk-like and rotation-dominated (Kormendy & Kennicutt 2004;Kormendy & Fisher 2008).Shen et al. (2010) demonstrated that the Galactic bulge is purely diskoriginated via a simple but realistic N-body model that reproduces the stellar kinematics of the BRAVA survey (Rich et al. 2007) and the photometric asymmetric boxy shape of the Galactic bulge (Weiland et al. 1994;Dwek et al. 1995;López-Corredoira et al. 2005) remarkably well.Therefore, the chemodynamical properties of many bulge stars should closely resemble those of disk stars.
To explore the possible chemodynamical connection between the bulge and disk stars, we follow the approach of Belokurov et al. (2020) and build row-normalized number density maps of our sample A (see Section 2) in the V f -[Fe/H] space for consecutive radial annuli, as shown in Figure 3.In order to better quantify the general trend of the V f -[Fe/H] distribution across the Galactic disk, for each radial bin we split the subsample into small V f bins and fit a Gaussian-Hermite function to their metallicity distribution profiles in order to obtain a sequence of best-fit mean [Fe/H] values and dispersions, which we overlay in the white curve in Figure 3.In this procedure, we do not include stars with [Fe/H] < − 1 dex that are believed to be mostly accreted stars, halo stars, or in situ Aurora stars.We observe a similar trend in all panels, from the outer disk to the inner few kiloparsecs: a chevron-like pattern at positive azimuthal velocities with a vertical extension toward small and negative azimuthal velocities.The upper and lower branches of the chevron pattern and the vertical extension have been suggested to trace the thin disk, the thick disk, and the Splash, respectively (Belokurov et al. 2020;Belokurov & Kravtsov 2022;Lee et al. 2023).The upper branch of the chevron, with high rotational velocities, is overall more metal-rich, while the lower branch, with intermediate rotational velocity, is more metal-poor, in agreement with our knowledge of the thin and thick disks.As for the vertical extension toward lower v f , it is consistent with the Splash, i.e., a heated disk population or the in situ halo.Note that the density map in Figure 3 is row-normalized: in the Splash region, the number of prograde stars actually outweighs that of retrograde stars, resulting in an apparent small net rotation.
Although the chevron pattern looks similar at different radii, there is a systematic variation across the disk in the slopes and conjunction points of the thin-and thick-disk branches.To better illustrate it, we draw the sequences (the white curves in Figure 3) together in Figure 4, where each color now indicates the radial bin.First, we notice that the vertical extensions (corresponding to the Splash) at different radii overlap and all have similar peak metallicity, indicating that they have a common origin in a global event, i.e., the GSE merger.Second, the chevron patterns form a coherent sequence with the peak [Fe/H] decreasing with radius, except for the innermost radial bin (R < 5 kpc; see the purple line in Figure 4).By extrapolating the trends revealed by the profiles of the outer radial bins, we would expect the chevron pattern of the  innermost radial bin (sample B) to be located at the highest metallicities.However, it is located in the middle of the other profiles, with a conjunction point at [Fe/H] ≈−0.25 dex.This could be partially caused by the observational bias: there are almost no stars in our sample close to the Galactic plane due to survey-specific limitations, in particular the strong dust extinction.Since stars closer to the midplane are more metalrich, the absence of disk stars would shift the overall metallicity trend toward the more metal-poor region.Moreover, those relatively metal-poor stars with metallicities larger than −1 dex of different origin (accreted stars, halo interlopers, metal-poor pre-existing in situ stars, etc.) are also present, further shifting the overall metallicity trend to lower [Fe/H] values.
For the thin-disk population, the rotational velocity V f is anticorrelated with metallicity.In a scenario in which the MW forms inside out (e.g., Frankel et al. 2019), this trend is mainly due to radial migration through churning and blurring (Schönrich & McMillan 2017;Vickers et al. 2021) and the negative metallicity gradient of the thin-disk stars: at a specific radius, stars with smaller V f likely migrated from the inner region that is more likely metal-rich, giving rise to the negative correlation between V f and [Fe/H].On the other hand, for the thick-disk population, V f is positively correlated with [Fe/H].Thick-disk stars are older and as such they have had more time to experience perturbations, resulting in larger velocity dispersion and slower rotation; because of their less circular orbits and higher distances from the Galactic plane z, they have also experienced less radial migration (Vera-Ciro et al. 2014).
From Figure 4, it is also apparent that the chevron patterns shift toward lower metallicities as the radius increases.The intervals between the upper branches (the thin disk) of the chevron-like pattern at different radii are larger than the intervals of the lower branches (the thick disk), which are almost aligned with each other; the intervals reflect the steeper negative metallicity radial gradient of the thin disk, and the weaker metallicity gradient in the thick disk, seemingly confirming that the thick disk is more disturbed and well mixed than the thin disk, with a weaker metallicity gradient (Vickers et al. 2021).The negative metallicity gradient of the thin disk can be reproduced in an inside-out disk formation scenario in cosmological simulations (Valentini et al. 2019).In addition, the slopes of the upper branches (thin disk) increase with radius R.This may also be explained with an inside-out formation scenario, where the inner disk has more time to enrich itself, thus becoming more metal-rich and with a broader metallicity distribution than the outer regions.Considering the flat rotation curve and the asymmetric drift effect, the slope of the V f -[Fe/H] profile of the outer disk is expected to be steeper than the inner disk.The slopes of the lower branches (thick disk) are shallower and almost constant with radius, implying that the thick disk could have formed in one single episode within a short time frame, an event that likely occurred in the turbulent early epoch when the gas supply from mergers was abundant and well mixed.

Chemodynamical Properties of the Galactic Bulge
In the last subsection, similar general trends in the V f -[Fe/H] space have been revealed for both the bulge and disk stars.In this section, our focus turns to the bulge region (sample B, introduced in Section 2).
In the disk region (R > 5 kpc), the V f -[Fe/H] pattern can be attributed to the thin and thick disks and the Splash.However, in the bulge region (R < 5 kpc), it is more complicated to attribute this pattern to distinct stellar populations, e.g., the bar, Splash, Aurora, or inner disks.The bar formation is a key event in the history of the MW, which requires a pre-existing stellar disk.Evidence of the bar in the V f -[Fe/H] space could help reveal the possible formation epoch of the Galactic disk.
We show the row-normalized number density map in the V f -[Fe/H] space for the bulge sample in the top panel of Figure 5 (note that it is just an enlarged version of the top left panel in Figure 3) and the radial velocity dispersion s V R in the bottom panel.The radial velocity dispersion increases from the upper right corner to the lower left corner monotonically.
In previous works, the metallicity is usually the only criterion used to separate different stellar populations in the Galactic bulge region (Ness et al. 2013a(Ness et al. , 2013b;;Rojas-Arriagada et al. 2020;Wylie et al. 2021).However, in the Splash region (−1 < [Fe/H] < − 0.5), the dynamically cold (V f > 100) and hot (V f < 100) populations co-exist in the same metallicity range: the traditional metallicity cut is not accurate enough for dissecting the different structures-kinematics plays a crucial role.Nonetheless, for comparison with previous works, we first examine the kinematic properties of populations with different metallicities.The results are shown in Figure 6, with the metallicity specified in the upper left corner of the third column.The left two columns show the average radial velocity á ñ V R and radial velocity dispersion s V R maps in the X − Y plane, respectively.The more metal-rich populations, with [Fe/H] > − 1 dex, all exhibit "butterfly"-like patterns in the average radial velocity maps, a signature of bar-like kinematics (Bovy et al. 2019;Fragkoudi et al. 2020;Queiroz et al. 2021;Drimmel et al. 2023).They are caused by the positive and negative radial velocities of stars on bar-supporting orbits in adjacent quadrants (Zoccali et al. 2014;Bovy et al. 2019;Drimmel et al. 2023).The corresponding s V R maps show larger amplitudes closer to the bar major axis (zero average radial velocity), where the radial velocity V R changes its sign.On the other hand, the metal-poor populations with [Fe/H] < − 1 dex do not show bar-like kinematics in the á ñ V R maps, and their s V R dispersion maps are quite isotropic, suggesting that they are not part of the bar structure but represent a dispersion-supported structure.Notice that the zero-radial-velocity line in the á ñ V R map should align with the major axis of the bar rather than the X-axis (i.e., the  Sun-GC line).This is mainly caused by the distance errors, which will be discussed in Section 6.
The rotation profiles of the bulge sample at various latitude slices are shown in the right column of Figure 6, where the errors were obtained via bootstrapping.We overlay data from the PIGS survey (black points; Arentsen et al. 2020) and the BRAVA survey (blue points; Shen et al. 2010) at the corresponding metallicity.There is good agreement between Gaia and the previous surveys, highlighting the quality of the Gaia DR3 data even in the bulge region.Except for the most metal-poor population, all stars in sample B display cylindrical rotation, with the more metal-rich stars, [Fe/H] > − 1 dex, rotating faster than the metal-poor stars, with [Fe/H] < − 1 dex.It is interesting that the metal-poor stars (−1.5 < [Fe/H] < − 1) display cylindrical rotation, albeit slower than the metal-rich ones and without a "butterfly"-like pattern: Shen et al. (2010) has suggested that a pre-existing classical bulge could be spun up by a rotating bar.In this case, the spun-up classical bulge can also develop cylindrical rotation, as an evidence of secular evolution of the bar (Saha et al., 2012).This effect will be investigated in more detail in Section 5.
Recent works have shown that the spinup phase of the MW disk ends at [Fe/H] ∼ − 1 dex and the bulk of the disk stars form afterward, with the azimuthal velocity continuously increasing toward higher metallicity (Belokurov & Kravtsov 2022;Semenov et al. 2023).Figure 6 shows that stars in the Galactic bulge region with [Fe/H]  − 1 dex exhibit bar-like kinematics, implying that the primordial inner disk at the time when the bar formed should have similar metallicity properties (i.e., [Fe/H]  − 1 dex).The result in the Galactic bulge is consistent with previous works focusing on the Galactic disk (Belokurov & Kravtsov 2022).
To better understand the chemodynamical properties of the Galactic bulge and to disentangle the different stellar populations, we divide the V f -[Fe/H] space into a grid (dashed black lines in Figure 5) of small cells.In this way, we do not assign the different loci of the V f -[Fe/H] diagram to different stellar populations, but systematically investigate each cell in order to avoid a priori bias; in Figure 7, we show the mean radial velocity maps in the X − Y plane for the stars in each grid cell.Neglecting the panels with low number statistics (the number of stars in each cell is specified in the bottom left corner), we can see that only some of the remaining panels show bar-like kinematics (enclosed within the blue frame).This result is consistent with the radial velocity dispersion map (bottom panel of Figure 5), where the stars with bar-like kinematics have lower radial velocity dispersion (s  90 V R km s −1 ), i.e., are dynamically cold.The panels in the middle column all have the same metallicity range of [−1, −0.4] dex: the top two panels at higher V f and lower s V R display bar-like kinematics, while the bottom three panels corresponding to the Splash region (enclosed with the red line) display dispersiondominated kinematics.At even lower metallicities (the left two columns of Figure 7 with [Fe/H]  −1 dex), the kinematics is consistent with a dispersion-dominated population (i.e., accreted stars, the classical bulge, halo stars, or the Aurora).The kinematics is therefore crucial for disentangling the bar from other populations (e.g., the Splash) in the bulge region.

Auriga Cosmological Simulation
In this section, we make use of a cosmological zoom-in simulation of an MW-like galaxy from Auriga to understand the observed pattern in the V f -[Fe/H] space and gain insight into the formation and evolution history of the Galaxy.
Auriga is a suite of 30 magnetohydrodynamical cosmological zoom-in simulations of MW-mass dark matter halos, with varying mass from ∼1 × 10 12 to 2 × 10 12 M e and evolving from redshift 127 to the present day, incorporating active galactic nuclei feedback, star formation, and magnetic fields (see Grand et al. 2017 for details).The Auriga simulations (Grand et al. 2024) are publicly available. 9We choose Auriga 23 (Au23), which resembles many of the MW properties.For example, it displays clear [α/Fe] bimodality for the thin/thick populations (Grand et al. 2018) and its last major-merger event happened at a lookback time t lb ≈ 10 Gyr, 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 f is crucial to the identification of the nonbar population with negative V f 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 Gyr), around the time of the merger (9.6-10.4Gyr), and after the merger (8.5-8.9Gyr).The last major merger of Au23 happened at around 10 Gyr (see Fragkoudi et al. 2020).Thus, stars of ages 9.6-10.4Gyr are mainly contributed by the starburst population and dynamically hot disk.
The face-on and edge-on distributions of the stars in the inner 5 kpc of Au23 are shown in Figure 8 (top and bottom rows, respectively) at t lb = 8.38 Gyr (left column) and 0 Gyr (right column).The 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 the results in the literature that stars with hotter kinematics will have thicker and rounder bars than stars with colder kinematics (Athanassoula et al. 2017;Debattista et al. 2017;Fragkoudi et al. 2017).According to Figure 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, resembling an oblate structure.The snapshot at present (t lb = 0) is shown in the right column.The stars born before the merger maintain a dynamically hot spherical structure, while the stars born around and after the merger resemble the thick and thin disks (both components showing the bar structure), respectively.The stars are selected in the bulge region with R < 5 kpc and |Z| < 3 kpc.
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 f − [Fe/H] space are shown in Figure 9 for the same lookback times as Figure 8.The stars born around the merger (i.e., the middle column) exhibit the chevron-like pattern, with a vertical extension toward lower V f , very similar to the one in observations (Figure 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 Figure 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 Figure 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 the thin-disk component.The stars born after the merger (top right panel of Figure 9) distribute within the upper right region at high metallicities and high V f , which is expected since they are born in the disk (see also the red contours of the left column of Figure 8).In Figure 9, comparing the V f -[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 a significant difference in the V f -[Fe/H] pattern during the same time period: the pattern shifts to lower azimuthal velocities and has a much larger velocity dispersion.Both external and internal mechanisms could have contributed to the disk heating.After the last major merger occurred about 10 Gyr ago, the galaxy experienced subsequent minor mergers (see Figure 11 in Fragkoudi et al. 2020), which can also contribute to the disk heating.In addition, as shown in Figure 8, the stars born after the merger host a strong bar.The formation of the bar requires a loss of angular momentum, which leads to a larger velocity dispersion.
We now investigate the V f − [Fe/H] space color-coded by the radial velocity dispersion s V R in Figure 10 and compare it to the data (bottom panel of Figure 5): we notice that stars born around the merger (the middle column in Figure 10) have a distribution similar to the observations.The velocity dispersion decreases monotonically from the lower left corner to the top right corner, consistent with the observed pattern in Figure 5.For the stars born before the merger (left column), s V R is relatively uniform at any given [Fe/H], which is consistent with the behavior of a dispersion-dominated structure.There is not much difference between the V f -[Fe/H] s -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 Figure 9: 8 Gyr of evolution (e.g., the accretion of cold gas, the following disk formation, and the 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 the bottom right panel).
We now focus on the stars born around the time of the merger in Au23.We separate them into three subpopulations-MPD (metal-poor disk), MRD (metal-rich disk), and the Splash (see the three boxes in Figure 9)-to mimic the observed transition from the "disk" (the top three panels of the fourth column in Figure 7) to the "disk coexisting with the Splash" (the middle column in Figure 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 Figure 10).MPD and MRD have similar V f distributions, but MRD is more metal-rich.The face-on and edge-on distributions of these three subpopulations are presented in Figure 11 at t lb = 8.38 Gyr (blue isodensity contour maps) and 0 Gyr (green).The corresponding velocity maps are shown in Figure 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 (∼2 kpc half-length) following the bar orientation in MPD and MRD.However, as shown in the á ñ V R color-coded maps in Figure 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 Figure 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 only looking at the mean V R maps alone, we cannot exclude the possibility that a few of the Splash stars indeed move on barsupporting orbits.Beyond 2 kpc, the Splash morphology becomes more spheroidal, resembling a "fluffy classical bulge" (see Figure 11 in Belokurov et al. 2020).   5. Interaction between the Bar and a Pre-existing Bulge:

N-body Simulation Test
To investigate the influence of a rotating bar on a preexisting nonrotating spheroidal structure (i.e., a classical bulge) and the origin of the observed cylindrical rotation for the metalpoor stars in the Galactic bulge region, we make use of three isolated N-body simulations of MW-like galaxies that selfdevelop 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 ∼ 1.2 × 10 12 M e ), a stellar disk (M d = 4.3 × 10 10 M e ) with a scale length of 2.5 kpc and a scale height of 0.3 kpc (resembling a thin disk), and a nonrotating stellar spheroidal structure-the pre-existing classical bulge.We build three models-namely M1, M2, and M3 -each with a 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 × 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 modeling modules 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 ∼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 Figure 13, for the disk (left two columns) and nonrotating 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 Gyr (top row) have an axisymmetric distribution without any velocity pattern.Once a bar forms (middle row), the butterfly pattern appears in the á ñ V R map due to the systematic motion of the bar orbits in different quadrants.As the bar grows, the á ñ V R butterfly pattern becomes larger (bottom row).On the other hand, the spheroidal nonrotating component does not show the á ñ V R 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 Galactic bulge from the simulations by setting the observer at R = 8 kpc and a bar viewing angle of 20°.The longitudinal rotational profiles of the particles are shown in Figure 14 for the disk (blue curves) and spheroidal nonrotating component (red curves) at different times (increasing from the left to right columns).M1, M2, and M3 are shown in the top, middle, and bottom rows, respectively.The rotation profiles of the nonrotating 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 Figure 13) from spheroidal to ellipsoidal, with the major axis aligned with the bar major axis.In Figure 13, we have seen that classical bulges may not display the á ñ V R 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 a 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| ∼ 10°, the rotation curves in the 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; see also Figure 6).This result could put constraints 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 nonrotating 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) metalpoor population ([Fe/H] − 1 dex).In our bulge sample B, stars with −1.5  [Fe/H]  − 1 dex do not exhibit the á ñ V R butterfly pattern (Figure 7) but show cylindrical rotation (Figure 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  [Fe/H]  − 1 dex) could represent the 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.

Effect of Distance Error on the 〈V R 〉 Map
Bulge stars have large distance measurement uncertainties due to the their larger distance, strong dust extinction, and crowding.To test the influence of the distance uncertainty on the á ñ V R maps projected onto the X − Y plane (i.e., the butterfly pattern), we first transform the Cartesian coordinates (X, Y, Z, V X , V Y , V Z ) of the simulation into observables (ra, dec, d, pmra, pmdec, rv), then perturb the heliocentric distance d with various levels of uncertainty, and finally transform the observational quantities into Galactic cylindrical coordinates.Figure 15 shows the original á ñ V R map (top left panel) and distanceperturbed maps, with the corresponding relative distance uncertainty (d err /d) listed in the upper left corner.The butterfly pattern becomes more significant as the relative uncertainties increase and the zero-velocity line gradually aligns with the Sun-GC line (the X-axis) as they reach d err /d ≈25%.Beyond 25%, the zero-velocity line remains aligned with the X-axis and the butterfly pattern keeps becoming more significant.The á ñ V R butterfly pattern is therefore a robust bar signature: in the case of bar presence, large distance uncertainties will not cause the pattern to disappear; rather, the pattern will merely twist such that the zero-velocity line becomes more aligned with the Sun-GC line of sight (see also Vislosky et al. 2024).

Splash Formation Scenario and the Early Formation of the Disk
The Splash in the V f -[Fe/H] space appears as a smooth transition from the thick disk to a dynamically hot component with small or retrograde rotation velocity (Figures 3 and 4).A natural explanation for the formation of the Splash is that it consists of thick-disk stars heated by a major merger (Belokurov et al. 2020).The Auriga cosmological simulation (Grand et al. 2020) also contains a Splash-like population.This scenario assumes that a thick disk was already in place before the major merger.Another mechanism that could be responsible for the formation of the Splash is the clumpy formation scenario; clumps form in the disk due to the gas fragmentation in the early gas-rich stage (Clarke et al. 2019).The clumps have high a star formation rate density to self-enrich in α elements rapidly, which are then destroyed by supernovae to dissolve into a high-α thick disk.Meanwhile, across the disk, the continuous long-term star formation gradually gives rise to a low-α thin disk.In the clumpy formation scenario, it is possible to form a Splash-like population, which corresponds to the low-angular-momentum tail of the thick disk (Amarante et al. 2020).However, the clump formation scenario could not produce a significant fraction of retrograde stars as a merger does.
We have shown in Section 3.1 and Figures 3 and 4 that the Splash exists across the disk as a uniform population with similar metallicities at different Galactic radii, which can be explained by both scenarios above.In the first scenario, there was already a thick disk in place at the time of the last major merger, which had formed in the early chaotic epoch of the galaxy.The chemistry of such a disk would be spatially well mixed, implying that stars kicked out of the thick disk into halo-like orbits by the major merger would have similar metallicities at different radii.In the clumpy formation mechanism, the proto-thick-disk originates from the clumps formed in the first Gyr during the gas-rich phase.These clumps were formed out of the gas with similar chemistry, experiencing a similar chemical enrichment process.In this scenario, the metallicity of the thick disk is also well mixed.In addition, Lee et al. (2023) reported that there is a relatively more metalpoor and α-depleted component in the Splash region.The authors attributed approximately half of the stars in this component to the accreted GSE and the other half to the starburst population formed during the merger.However, the low-α Splash component can also be explained by the clumpy formation scenario, if the early clumps scattered thin-disk stars with low-α abundances to hotter orbits.The current observations of the Splash stars cannot discriminate between the two scenarios.It is also possible that both scenarios could have played a role in the formation of the Splash.

Why is [Fe/H] ∼ − 1 Special?
The multiple stellar populations and structures in the Galactic bulge make its chemodynamics complex.Most Galactic bulge studies have used metallicity as a rough tracer of different stellar populations, to explore its structure and kinematics.Many surveys have revealed that the stellar populations with [Fe/H]  − 1 dex show bar-like kinematics (i.e., have a disk origin; Kunder et al. 2012;Ness et al. 2013b;Di Matteo 2016;Fragkoudi et al. 2018), while stars with [Fe/H]  − 1 dex are consistent with a dispersion-dominant spheroidal structure (Dékány et al. 2013;Arentsen et al. 2020).An important question naturally arises: what is the reason for this transition at [Fe/H] = − 1 dex?
As shown in Figure 3, the metallicity distribution of the Splash at different radii is within −1 < [Fe/H] < − 0.4 dex; the lower-metallicity limit of the Splash is approximately [Fe/H] ∼ − 1 dex, where the kinematics changes from bar-like to dispersion-dominated.Under the assumption that the formation of the Splash is merger-induced and mainly composed of heated protodisk stars and the subsequent starburst population that is more metal-rich than the protodisk (with partly overlapping metallicities), then the most metal-poor stars in the Splash cannot be more metal-poor than the protodisk itself.Therefore, the [Fe/H] ∼ −1 dex limit should be the lower limit of the protodisk at the time of the major-merger event.Afterward, the bar forms out of the dynamically hot disk stars.This is consistent with recent simulations showing that a bar can form from the thick disk (Ghosh et al. 2023).Thus, the disk stars with [Fe/H]  −1 dex become the building blocks of the bar and kinematically detach from the more metal-poor stars.However, most of the hotter Splash population with [Fe/H]  −1 dex does not participate in the formation of the bar; instead, the kinematics of the Splash is similar to that of stars with [Fe/H] < −1 dex.However, in Section 4, we have shown that the morphology of the Splash could become bar-like in the inner 2 kpc, which could be due to the long-term gravitational attraction of the bar.Beyond 2 kpc, the Splash becomes more spheroidal-like (e.g., Belokurov et al. 2020).
In addition, other works also suggest that [Fe/H] ∼ − 1 dex is a turning point in the galaxy chemical evolution.Based on the Gaia data with the metallicity derived from XGBoost, Rix et al. (2022) Figure 14.Evolution of the rotation profiles of the disk and spheroidal components in the mock observation of the bulge region for the three models.The blue and red curves represent the disk (bar) and bulge rotation profiles, respectively.The mass fraction of the pre-existing classical bulge increases from the top to bottom panels.Apparently, the classical bulge would be spun up by the bar.And in a more massive bulge, the bar deviates from a perfectly cylindrical rotation.
reported that the slope of the metallicity distribution function is Using a large sample of subgiants with precise age measurements, Xiang & Rix (2022) separated the sample into early-phase (low angular momentum and high α; i.e., thickdisk and halo) and late-phase (high angular momentum and low α; i.e., thin-disk) stars; the thick-disk stars in the early phase have [Fe/H]  − 1 dex.This metallicity boundary agrees with our inference from the Galactic bulge observations that the low-metallicity limit of the thick disk is [Fe/H] ≈ − 1 dex.

Conclusion
Recent Solar neighborhood studies have found that the thin disk, thick disk, and Splash appear as distinct features of the V f -[Fe/H] diagram (Belokurov et al. 2020;Lee et al. 2023).In this work, we have used data from Andrae et al. (2023), a catalog with metallicities extracted from the Gaia-XP spectra, to study the V f -[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 Splash in the V f -[Fe/H] space exist across the disk.The V f -[Fe/H] profiles of the thin disk show systematic variation with radius, which is in agreement with the inside-out formation scenario.The V f -[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 f -[Fe/H] space with disk stars, implying that it originated from the disk.By investigating the bar signatures of the stellar populations in the V f -[Fe/H] space and only considering the populations with reliable number statistics, we found that: By analyzing an MW-like Auriga simulation, Au23, we found that only the stars born around the time of the last major merger (∼10 Gyr ago), which mainly resulted from the starburst population and the dynamically hot disk, have a V f -[Fe/H] pattern similar to the observations 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.Afterward, the occurrence of the last major merger heated some of the thick-disk stars and triggered the starburst formation that contributed to the Splash, whose lowermetallicity limit is ∼ − 1 dex, which is inherited from the thick disk.Thus, the lower-metallicity limit of the thick disk should be ∼ − 1 dex at the time of the last major merger.Subsequently, the thick-disk stars with [Fe/H]  − 1 dex form the bar structure, their kinematics becomes bar-like and detaches from those of more metal-poor stars.In another aspect, the merger kicks some of V R butterfly pattern becomes more significant, with the zero-velocity line more closely aligned with the X-axis rather than the bar major axis (20°tilt from the X-axis).
the thick-disk stars to halo-like orbits and induces the 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 a 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 nonrotating spheroidal component (a classical bulge) and a later-formed bar and confirm the result of Saha et al. (2012) that an initially nonrotating spheroidal component can be spun up to develop cylindrical rotation under the bar's 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 ∼ 10°for the different Galactic latitudes we consider, in the inner region they show relatively large discrepancies: 1.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.2. 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|  4°compared to the lower latitudes.
In the MW, the rotation profiles at different latitudes are almost the same (see Figure 6), being most close to the model with the least massive bulge, implying that the classical bulge in the 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.

Figure 1 .
Figure 1.Spatial distribution of the full Gaia-XP-spectra-derived sample (5,497,737 stars; sample A), projected onto the X − Y (top) and X − Z (bottom) planes.The red cross and dot denote the Sun position and the GC, respectively.The red circle of 5 kpc radius encloses the bulge stars (330,414 stars; sample B).

Figure 2 .
Figure 2. Comparison of the metallicity, distance, azimuthal velocity, and radial velocity for the common stars in the APOGEE DR17 and Gaia samples.The two samples are in good agreement.The distances adopted for Gaia and APOGEE are Bayesian distances (Bailer-Jones et al. 2021) and StarHorse distances (Queiroz et al. 2020), respectively.

Figure 3 .
Figure 3. Row-normalized number density maps of the V f -[Fe/H] space of stars at various radial bins, indicated in the top left corner of each panel.The white curve in each panel represents the best-fit mean [Fe/H] values from a Gaussian-Hermite function at a sequence of small V f bins (see the text for details), with the horizontal white lines indicating the corresponding 1σ metallicity dispersion within each V f bin (stars with [Fe/H]< − 1 dex are not included).The number in the bottom right corner of each panel denotes the number of stars in that radial bin.The horizontal dashed black line indicates V f = 0.

Figure 4 .
Figure 4. V f -[Fe/H] patterns at different radii from Figure 1 (white curves) are shown here together to better illustrate their systematic variation across the Galactic disk.The color indicates the radial bin, with the bulge region (R < 5 kpc) shown in purple.

Figure 5 .
Figure 5. Top: the row-normalized number density of bulge stars (R < 5 kpc; sample B) in the V f -[Fe/H] space, with the dashed lines indicating the division of the stars into different grids.Bottom: similar to the top panel, but color-coded with the radial velocity dispersion.A clear pattern emerges with higher dispersion at the bottom left corner and lower dispersion at the top right corner.

Figure 6 .
Figure 6.Kinematic properties of the Galactic bulge stars (sample B) within different metallicity bins sliced as the dashed vertical lines marked in Figure 5.The mean radial velocity á ñ V R (left column) and radial velocity dispersion s VR (middle column) are projected onto the X − Y plane.The position of the GC is marked with a purple cross.The rotation profiles (l − V gsr ) at different latitudes of 0°<b< 5°(red), 5°< b <8°(green), and 8°< b <10°(yellow) are shown in the right column.The PIGS survey (black points; Arentsen et al. 2020) and the BRAVA survey (blue points; Kunder et al. 2012) show good agreement with the Gaia sample.

Figure 7 .
Figure 7.The á ñ V R maps in the X − Y plane for the grid defined in the V f -[Fe/H] space in Figure 5.The positions of the panels in this figure correspond to the grid in Figure 5.The number in the bottom left corner of each panel is the number of stars.The panels with good number statistics showing bar-like kinematics are enclosed within the blue frame, consistent with the dynamically cold region in the s VR map in the bottom panel of Figure 5.The panels denoting the Splash region are delimited by the red box.No bar-like kinematics can be seen in the Splash and the other more metal-poor stars.

Figure 8 .
Figure8.Evolution of the spatial distribution of stars born before, during, and after the last major merger about 10 Gyr ago in the Auriga simulation Au23.The top and bottom rows show the number density maps in face-on (X − Y plane) and edge-on (X − Z plane) views, respectively.The number density contours of the three populations are shown in different colors, which are indicated in the lower left panel.The left column shows a snapshot with lookback time t lb = 8.38 Gyr, which is ∼1.5 Gyr after the last major merger.The snapshot at present (t lb = 0) is shown in the right column.The stars born before the merger maintain a dynamically hot spherical structure, while the stars born around and after the merger resemble the thick and thin disks (both components showing the bar structure), respectively.The stars are selected in the bulge region with R < 5 kpc and |Z| < 3 kpc.

Figure 9 .
Figure 9.The row-normalized number density maps in the V f -[Fe/H] space for stars in the Auriga simulation that were born before the merger (left column), during the merger (middle column), and after the merger (right column).The top and bottom rows show two different lookback times, at t lb = 8.38 Gyr and t lb = 0 Gyr, respectively.Note the three populations show dramatically different patterns.The three boxes in the middle column represent three subpopulations that are analyzed in detail later.

Figure 10 .
Figure 10.The V f -[Fe/H] space color-coded with the radial velocity dispersion for stars born before (left column), during (middle column), and after the merger (right column) in the Auriga simulation Au23.The layout of this figure is the same as Figure 9.The top and bottom rows represent t lb = 8.38 Gyr and t lb = 0 Gyr, respectively.

Figure 11 .
Figure11.The face-on (top) and edge-on (bottom) images for the Splash, MPD, and MRD components at t lb = 8.38 Gyr (left three columns) and t lb = 0 Gyr (right three columns).MPD and MRD all show the same bar structure.The Splash is overall spherical, with a small bar early on and a relatively weak bar in the late stage.

Figure 12 .
Figure 12.The á ñV R color-coded face-on distribution of the Splash, MPD, and MRD components at t lb = 8.38 Gyr (left three columns) and t lb = 0 Gyr (right three columns).Consistent with the results in Figure7, both MPD and MRD show the bar-like butterfly pattern, while the Splash is random-motion-dominated.

Figure 13 .
Figure 13.Evolution of the morphological and kinematical properties of the disk and the spheroidal component (classical bulge) of the model M1.From top to bottom, the rows show snapshots at 0.0, 2.4, and 4.5 Gyr, respectively.The first two columns are the number density map and median v R map in the X − Y plane for disk particles, while the right two columns are the number density and velocity maps for the spheroidal component.
bulge region for the bulk of stars with [Fe/H]  − 1.In the one-zone chemical evolution model(Weinberg et al. 2017), consequence of a self-enriching in situ system.The slope then becomes steeper at [Fe/H]  − 1 dex, which implies an event with rapid gas accretion (or a merger) corresponding to the onset of the thick disk(Belokurov & Kravtsov 2022;Semenov et al. 2023).Chandra et al. (2023) also used red giants with Gaia-XP spectra to study the MW evolution in a similar parameter space (L z /L c -[Fe/H]) as this work (V f -[Fe/H]).They found that the old high-α disk starts at [Fe/H] ∼ − 1 dex, which agrees well with our results for the protodisk, for which we estimate a lowermetallicity limit of [Fe/H] ∼ − 1 dex.

1.
Stellar populations with [Fe/H]  − 1 dex are dispersiondominated. 2. Stellar populations with −1  [Fe/H]  − 0.4 show a bimodal distribution: the stars with V f 100 km s −1 follow the bar kinematics, while the stars with V f  100 km s −1 (Splash) have dispersion-dominated kinematics similar to those of more metal-poor populations.3. Stellar populations with [Fe/H]  − 0.4 dex all have barlike kinematics.

Figure 15 .
Figure 15.Effect of the relative distance error (d err /d) on the mean á ñ V R map in the X − Y plane of the Au23 simulation at t lb = 0 Gyr.The top left panel shows the original map and the other panels are the perturbed different levels of distance uncertainty, as indicated by the text in the upper left corner of each panel.The contours in each panel are the density contours of the original bar structure.With the increasing distance uncertainties, the á ñV R butterfly pattern becomes more significant, with the zero-velocity line more closely aligned with the X-axis rather than the bar major axis (20°tilt from the X-axis).