Metallicity Properties of the Galactic Bulge Stars Near and Far: Expectations from the Auriga Simulation

Using the high resolution Milky Way-like model from Auriga simulation we study the chemical properties of the Galactic bulge, focusing on the metallicity difference between stars on the near side (in front of the Galactic center) and the far side (behind the Galactic center). In general, along certain sight lines the near side is more metal-rich than the far side, consistent with the negative vertical metallicity gradient of the disk, since the far side is located higher above the disk plane than the near side. However, at the region $l<0^\circ$ and $|b|\le6^\circ$, the near side is even more metal-poor than the far side, and their difference changes with the Galactic longitude. This is mainly due to the fact that stars around the minor axis of the bar are more metal-poor than those around the major axis. Since the bar is tilted, in the negative longitude region, the near side is mainly contributed by stars close to the minor axis region than the far side to result in such metallicity difference. We extract stars in the X-shape structure by identifying the overdensities in the near and far sides. Their metallicity properties are consistent with the results of the whole Galactic bulge. The boxy/peanut-shaped bulge can naturally explain the metallicity difference of the double red clump stars in observation. There is no need to involve a classical bulge component with different stellar populations.

Observations of the red clump (RC) stars revealed an Xshaped structure in the Galactic bulge (McWilliam & Zoccali 2010;Nataf et al. 2010) which was later associated with the boxy/peanut-shaped bar after the buckling instability (Li & Shen 2012;Ness et al. 2012;Li & Shen 2015). The X-shaped structure reveals itself as a bimodal distribution of the appar- ent magnitude of the red clump stars, which are commonly used as standard candles to trace the bulge structure. Gonzalez et al. (2015) summarized the main observational properties of the RC double peaks: 1. The double peaks are clearly visible in |b| 5 • region and they merge to a single peak in |b| 5 • (P1).
2. The magnitude difference between the two peaks increases with the increasing absolute value of the Galactic latitude (P2).
3. At a fixed Galactic latitude, the brighter RC is more dominant in the positive Galactic longitude region than the fainter one, which is more important in negative longitudes. In |l| > 4 • region only one of them is visible (P3).
4. The variation of the RC magnitudes as well as the magnitude difference between the two RCs are independent with the observational bands (P4).
All these properties are consistent with an X-shaped structure in the bulge (X-shape scenario). Li & Shen (2015) also pointed out that there is no simple letter "X"-like structure in the Galactic bulge, but a boxy/peanut-shaped component; the X-shape is likely a visual effect caused by the concave curved isodensity contours in the inner region of the peanut bulge.
Recently, Lee et al. (2015) and Joo et al. (2017) proposed a different scenario to use the multiple population classical bulge to explain the split of RC stars (MCB scenario). They argued that the magnitude difference between the two RCs reflects the intrinsic luminosity difference of two stellar populations in a classical bulge, with an additional unbuckled bar embedded in; the intrinsic magnitude difference of the two RC populations originates from the helium abundance difference between the first generation (G1) stars and later generations (G2+). Such phenomena have been observed in globular clusters (see Lee et al. 2015 and references there in). In Lee et al. (2015), their Fig. 4 showed a superposition of the classical bulge and a thin bar in the MCB scenario, with the two different RC populations in the classical bulge only. The bar component blurs the magnitude difference of the double RCs. Under this scenario, only one peak is observable near the midplane of the MW (P1). With larger |b|, such blurring effect becomes weaker to make the double peaks gradually visible, with their magnitude difference increasing with the Galactic latitude (P2). Since the bar is tilted, the blurring effect is different in different longitudinal regions, leading to the spatial variation of the relative importance between the double peaks (P3). Condition P4 is also accounted for in this scenario.
Metallicity information could be a key factor in discriminating the two scenarios. The MCB scenario predicts an intrinsic metallicity difference between the two RCs' stars (Lee et al. 2015;Joo et al. 2017), which was supported by following observations (Lee et al. 2019;Lim et al. 2021). However, Gonzalez et al. (2015) argued that the X-shape scenario could also produce a similar but smaller metallicity difference between the two RCs. In the X-shape scenario, the MW has a vertically thickened boxy/peanut-shaped bar, indicating that the Galaxy evolution is dominated by the secular process. On the other hand, the MCB scenario has a dominant classical bulge component corresponding to a more violent early merger history. Therefore, towards the Galactic bulge region, in a given line-of-sight, the double peaks in the apparent magnitude distribution will have different metallicity properties in the two scenarios. In the X-shape scenario at intermediate Galactic latitudes, stars are bar dominated so the metallicity difference between the bright (near side) and faint (far side) of the RC stars in the Galactic bulge is expected to vary in different (l, b). For MCB scenario, a dominant classical bulge component will result in similar metallicity difference between the bright and faint RC stars across the Galactic bulge region.
Aiming to better understand the Galactic bulge structure and chemical properties, with the high resolution MW-like Nbody model Auriga halo 23 (AU23) in Grand et al. (2017), we investigate chemical properties of the Galactic bulge stars, particularly focusing on the difference between stars in the near and far sides of the Galactic center. In Section 2, we briefly describe the model used here. In Section 3, we explore the metallicity properties of the Galactic bulge of AU23. In Section 4, we compare the simulation results with observations and discuss the influence of the bar angle on the metallicity difference between the near and far sides. The results are summarized in Section 5.
2. NUMERICAL MODEL OF THE GALACTIC BULGE Auriga project introduces 30 cosmological magnetohydrodynamical zoom-in simulations of galaxies in isolated MW-like mass dark halos (Grand et al. 2017), namely AU1−AU30. Fragkoudi et al. (2020) systematically analyzed the chemodynamics of the simulations in the Auriga project. Among all the galaxies in the Auriga project, AU23 exhibits a clear boxy/peanut-shaped bar similar to the Galactic bulge, as shown in Fig. 2 of Grand et al. (2017) and Fig. 4 of Fragkoudi et al. (2020). Mass of the baryons (stars and gas) in AU23 is about 9 × 10 10 M in a dark matter halo of ∼ 1.6 × 10 12 M . The disk scale length R d is about 4.99 kpc and the disk to total stellar mass ratio (D/T ) is 0.63, with an accreted stellar mass fraction within 0.1R 200 at ∼ 0.1. The metallicity of the single stellar population (stellar particles in simulation) is determined by the initial metallicity of the parent gas cell, and the metal return for SN II, SN Ia and asymptotic giant branch (AGB) stars are concerned in each time-step 4 . In this paper, we set the solar position at (X, Y , Z) = (-8, 0, 0.02) kpc with the Galactic center (GC) at (0, 0, 0) kpc and the angle between the bar major axis and Sun-GC line at 30 • . Fig. 1 shows the face-on and edge-on density maps of the stellar disk in AU23, which exhibits a clear boxy/peanutshaped bulge. Grand et al. (2018) showed that the stars in AU23 reproduced a chemical thin/thick disk dichotomy as seen in the MW. The thick α-rich disk was produced via a merger-induced starburst event 10-11 Gyr ago, after that the thin low-α disk and bar began to form. Therefore, the stars younger than 11 Gyr in AU23 are more bar-like with strong boxy/peanut-shaped bulge. The older stars older than 11 Gyr resembles a spheroidal structure similar to a possible classical bulge or inner halo structure, contributing to ∼ 15% of the total stellar mass. Similar tendency has been found in both observations (e.g. Babusiaux et al. 2010;Ness et al. 2012;Catchpole et al. 2016) and other simulations (e.g. Debattista et al. 2017;Portail et al. 2017). In the following analysis, stars younger than 11 Gyr are used to trace the double peak positions in the Galactic bulge region, while the full sample (including stars older than 11 Gyr) is used to estimate the metallicity difference between the two peaks.

RESULTS
In this section, we first show the metallicity distribution of stars in each latitude slice. Then we explore the metallicity difference between the stars in the near and far sides of the bulge in details. Finally we focus on the metallicity difference between the stars in the double peak position in the Galactic bulge. . We exclude the stars <4 kpc or >12 kpc that are not located in the bulge. At larger Galactic latitude regions (|b| > 6 • ), the near side is more metal-rich than the far side at all longitudes. However, closer to the mid plane(|b| ≤ 6 • ), at positive Galactic longitude region the near side is more metal-rich than the far side, and vice versa in the negative longitude side.

Metallicity Distributions in Latitudinal Slices of the
Galactic Bulge In Fig. 2, we divide the stars in AU23 into different Galactic latitude slices, from b = −9 • to b = 9 • with 1 • height for each slice. We further superpose the opposite latitudinal slices (±b) to increase the number statistics. In each panel of Fig. 2 ([Fe/H] color-coded), the contour represents the number density distribution of stars younger than 11 Gyr to highlight the position of the double peaks, useful to identify the stars belonging to the X-shape (similar to observations). In order to better visualize the density contours, larger bins are used in the slices at |b| > 6 • . Including all the stars will produce a smooth distribution. In the slices at |b| 5 • , there are two clear over-density regions, corresponding to the double peaks in the apparent magnitude distributions of the RC stars as in observations (McWilliam & Zoccali 2010;Nataf et al. 2010) and previous simulation works (e.g. Li & Shen 2012;Ness et al. 2012).
From the metallicity distribution in Fig. 2, we can see that in slices |b| ≥ 4 • , the two over-density regions are more metalrich than the region between them. In slices |b| ≤ 3 • , the stellar metallicity is higher close to GC and along the major axis of the bar. The slices closer to the mid-plane are more metal-rich. In slices |b| > 4 • , the near side (bright RCs) seems slightly more metal-rich than the far-side (faint RCs).

Metallicity Difference Between the Near and Far Sides
of the Galactic Bulge According to the stellar distance from the GC, we can divide the stars into two groups: 4 kpc< d <8 kpc as the near side, and 8 kpc< d <12 kpc as the far side. By analyzing the metallicity properties, we can get a qualitative picture of the metallicity difference between the near and far sides of the bulge. (the difference between previous two panels). In the high latitude region (|b| > 6 • ) the near side is more metalrich than the farther one at all longitudes. This is due to the negative vertical metallicity gradient of the MW; along certain line-of-sight towards the Galactic bulge region, the fainter/farther peak is located at larger vertical height above the disk plane than brighter/closer one. In addition, in the right panel we can see a horizontal [Fe/H] gradient in the intermediate latitude region. Close to the mid-plane (|b| ≤ 6 • ) at negative longitude side (l < 0 • ), the metallicity difference even becomes negative. This is unexpected, since the vertical metallicity gradient of the disk should indicate the near side to be more metal-rich than the far side.
To better understand the variation of ∆[Fe/H] across the bulge region, we show the faced-on and edge-on maps of AU23 color-coded with the mean [Fe/H] in the left panel of Fig. 4. The metallicity is larger along the major axis of the bar than that along the minor axis in the face-on map. The more metal-rich (younger) stars are more concentrated around the bar major axis with higher metallicity, while the metal-poor (older) stars show rounder distribution with smaller ellipticity, as shown in Fig. A1 of Fragkoudi et al. 2020. It can help to understand the phenomena in Fig. 3. The right panels in Fig. 4 show the latitude-distance distributions of stars in two vertical slices at l = ±3 • . The black solid line marks the position of the GC distance at 8 kpc. For the sight line at l = 3 • , close to the mid-plane (|b| ≤ 6 • ), the near side is dominated by the stars in the major axis of the bar (metal-rich), while the far side is populated by the stars in the minor axis, which is more metal-poor. On the contrary, for the sight line at l = −3 • and |b| ≤ 6 • , the near side is mainly populated by stars in the minor axis region (metal-poor) and the far side is mainly dominated by the major axis region of the bar (metal-rich). This difference further enhances the horizontal metallicity gradient close to the mid-plane. Interestingly, in the left panels of Fig. 4 the most metal-rich region of the bar seems to coincide with the central boxy core identified in Li & Shen (2015) for numerical buckled bars.
In the region away from the mid-plane, the metallicities of both the near and far sides decrease, with the far side showing a steeper latitudinal gradient. Therefore, at both positive and negative longitude regions, the near side gradually becomes more metal-rich than the far side. This is demonstrated in the lower panel of an edge-on projection of the model with two lines of sight at b = ±5 • . The vertical metallicity gradient is more important to result in a positive metallicity difference between the near and far sides. At such latitude range we can still find a longitudinal gradient of the metallicity difference. The negative longitude side generally shows smaller (or even negative) metallicity difference. This difference arises from the fact that at negative longitude region, the near side is mainly contributed by the more metal-poor minor axis of the bar than the more metal-rich far side around the major axis. The latitudinal metallicity gradient decreases faster in the far side, to gradually result in a positive metallicity difference at larger latitude regions.

Metallicity Difference between the Bright and Faint
peaks of the X-shape The results in Section 3.2 hint that the bright (near) and faint (far) peaks of the X-shape may show similar patterns To quantify such phenomenon more precisely, we transform the stellar spatial distribution into distance histograms for different (l, b) fields, similar to Fig. 2 in Li & Shen (2012). To better visualize the double peaks in the distance distribution, we mainly focus on stars younger than 11 Gyr (as shown in Fig. 2) to determine the peak position. Note that the double peaks can be identified in the region |l| < 3 • and 3 • |b| 8 • . Then we apply a Gaussian mixture model provided by scikit-learn (Pedregosa et al. 2012), which implements an expectation-maximization algorithm, to identify the location of the double peaks. More details are shown in the Appendix. In the following analysis, we pick out stars (including all age populations) that lie within 0.25 kpc distance from the peaks to study their metallicity.
The left column of Fig. 5 shows the mean metallicity and metallicity differences as a function of the longitude in different latitude slices from |b| = 4 • to 8 • of the bright and faint double peaks. For lower Galactic latitude slices (|b| = 4 • , 5 • , 6 • ), with l increasing, [Fe/H] of the brighter (fainter) peak gradually increases (decreases). ∆[Fe/H] gradually increases from negative (at l < 0 • ) to positive values (at l > 0 • ) due to the metallicity difference between the major and minor axis of the bar demonstrated in Fig. 4. In slices |b| ≥ 6 • almost all the metallicity difference values are positive due to the vertical metallicity gradient.
The right column of Fig. 5 shows the similar result but in different Galactic longitude slices (l = −3 At each slice, the brighter peak always has a shallower latitudinal metallicity gradient than the fainter one, since at the same latitude the fainter/farther stars are located at larger vertical heights than the brighter/closer stars. In the l = −3 • slice, ∆[Fe/H] gradually increases from negative to positive values with increasing |b|. For l ≥ 0 • slices, at different latitudes, ∆[Fe/H] values are all positive. This has been demonstrated in Fig. 4: at l ≥ 0 • , the near side is mainly contributed by the major axis part of the bar, which is more metal-rich than the minor axis part to produce a positive metallicity difference. Therefore, according to the simulation, we can see that the near side is not always more metal-rich than the far side. The near side is even more metal-poor than the far side at −3 • < l < −1 • and |b| < 6 • . As we have explained in Section 3.2, it is the metallicity difference between the minor axis (metal-poor) and major axis (metal-rich) of the bar to result in negative ∆[Fe/H] values in such region. In the other fields, the near side is generally more metal-rich, which is mainly due to the vertical metallicity gradient in the bulge/bar structure. Although this model may not truly reflect the chemical properties of the Galactic bulge, it is still valuable to shed light on the spatial variation of the metallicity differences between the near and far sides across the bulge region. Our results confirm that the buckled bar can naturally explain the metallicity difference between the two RC peaks. A multi-population classical bulge component, i.e., the MCB scenario, is not necessarily needed to explain the observations. Most importantly, we find a gradient of the metallicity difference along the longitudinal direction at |b| ≤ 6 • in the simulation shown in Figs. 3 and 5. Such a horizontal gradient is mainly caused by two factors. The first one is that the metallicity is higher close to the major axis of the bar, and lower around the minor axis. The tilted bar angle is the second factor. However, such a phenomenon is not expected in the MCB scenario; a prominent spheroidal classical bulge tends to produce similar ∆[Fe/H] values between the double RCs at the opposite ±l regions. Particularly at |b| = 5 • region, we find a monotonically increasing pattern of ∆[Fe/H] with l. This trend should not be expected in the MCB scenario. In addition, along the latitudinal direction, we note that the vertical metallicity gradient in brighter peaks is shallower than fainter ones. In the MCB scenario, the bright and faint peaks should have no significant difference in the vertical metallicity gradient.
The metallicity difference map between the near and far sides in Fig. 3 shows similar pattern to the profiles in Fig. 5. We apply linear fitting to the result in Fig. 3 to get a more quantitative estimation of the longitudinal (horizontal) and latitudinal (vertical) gradients of ∆[Fe/H] with error estimated by Monte Carlo simulation. The results are listed in Table 1. We can see that both directions generally show positive gradients except for the latitudinal gradients at positive l. The longitudinal gradient is highest at smaller latitude values (|b| = 3 • ) and decreases towards higher latitudes. The latitudinal gradients are much higher at the negative longitudinal side than the positive longitudinal side with slightly negative gradients.
Our results show that a self-consistent buckled bar in AU23 can account for the metallicity difference between the bright and faint RCs in the Galactic bulge with a similar value to the observation at the bulge field (l, b) = (−1 • , −8.5 • ). Compared to the MCB scenario, a boxy/peanut bulge can naturally explain the structure, kinematics, and metallicity properties of the bulge stars. On the other hand, the MCB scenario has to include an unbuckled but very thick bar (extending to |b| ∼ 8 • ) to produce the Galactic longitude dependence of RCs at high |b| region (Joo et al. 2017). Such a bar configuration is actually quite unlikely to be found in the MW or other disk galaxies. Future observations of the longitudinal/latitudinal metallicity difference gradients will be the smoking gun evidence to discriminate the MCB scenario. seems sensitive to the bar tilting angle. Fig. 7 shows the gradient d∆[Fe/H]/d l at different |b| slices for different α. The gradient curves at |b| < 6 • show a sinusoidal pattern, with the steepest gradient at α ∼ 45 • . For higher |b| slices (|b| ≥ 7 • ), the gradient is much shallower and varies mildly at different bar tilting angles.
With the gradient curves in Fig. 7, future observations of the bulge stars could help to constrain the bar tilting angle. It is suggested to measure the gradients in multiple latitudinal slices lower than 6 • since the other curves at higher latitudes in Fig. 7 are relatively flat with low amplitudes (< 0.005 dex/deg). If the gradient in only one latitudinal slice is used, even a small error will result in a large uncertainty of the bar tilting angle. As listed in Table 1, the gradient at |b| = 5 • is 0.012 dex/deg with the error of 0.002, corresponding to α = 15 • ∼ 45 • . Therefore, to better constrain the bar angle, the error of the gradients should be small, and measurements at multiple latitudinal slices are also needed.

CONCLUSION
We analyze the Auriga simulation AU23 to investigate the chemical properties of the boxy/peanut-shaped bulge, especially focusing on the difference between the near and far sides of the Galactic bulge. By identifying the peak positions of the stellar distance distributions towards the bulge region, we can also utilize the metallicity difference properties to discriminate the X-shape scenario and the MCB scenario for the origin of the RC double peaks in the Galactic bulge.
Stars in the bulge are separated into the near side (d < 8 kpc) and far side (d > 8 kpc) based on their relative position with respect to the Galactic center. Generally speaking, the near side is more metal-rich than the far side, especially for regions at larger latitudes. However, we find an unexpected result that the far side could be more metal-rich than near side at l < 0 • and |b| 6 • . Such metallicity difference pattern is resulted from the tilted bar and the fact that the bar majoraxis region is more metal-rich than the minor-axis region. Such phenomenon can also be found for stars in the double peaks of the X-shaped structure. For the negative longitude region (l < 0 • ), at intermediate latitude (|b| ≤ 6 • ), the line-ofsight will first go across the metal-poor region at the near side (d < 8 kpc) and then the metal-rich major axis region at the far side (d > 8 kpc), resulting in the negative metallicity difference. At positive longitude region, the line of sight crosses the more metal-rich major axis for the near side and metalpoor minor axis for the far side to generate positive ∆  away from the mid-plane than the near side, the metallicity of the far side decreases faster with latitude compared to the near side. Therefore, for the negative longitude region the near side gradually becomes more metal-rich than the far side at high latitude region (|b| 7 • ), and it makes the longitudinal gradient shallower than the lower latitudinal region. This simulation may not be the real portrait of the Milky Way. However, the chemical patterns revealed in the simulation are still valuable to help us understand the evolution and formation history of the Galactic bulge. We found that a secular evolution produced boxy/peanutshaped bulge can account for the metallicity difference between the double RCs in observations. The MCB scenario is not necessary for such phenomenon. Further more, we note that the horizontal gradient is sensitive to the bar titling angle. Such relationship can help to constrain the bar tilting angle of the MW with future observations in the Galactic bulge (e.g., SDSS-V, BDBS, Gaia, CSST, etc). 6. ACKNOWLEDGEMENTS We thank the referee for helpful comments and suggestions. We also thank Robert Grand for providing the Auriga simulation, and Juntai Shen for helpful discussions.

APPENDIX IDENTIFYING THE DOUBLE PEAKS
In order to identify the double peaks in difference (l, b) windows, we apply the sklearn.mixture.GaussianMixture provided by scikit-learn (Pedregosa et al. 2012). The results are shown in Fig. 8. The cyan histograms are the stellar distance distributions for stars at 5 kpc < d < 11 kpc in AU23, which serve as the input for sklearn.mixture.GaussianMixture.