The chemistry of stars in the bar of the Milky Way

We use a sample of 938 red clump giant stars located in the direction of the galactic long bar to study the chemistry of Milky Way bar stars. Kinematically separating stars on bar orbits from stars with inner disc orbits, we find that stars on bar-like orbits are more metal rich with a mean iron abundance of<[Fe/H]>=+0.30 compared to<[Fe/H]>=+0.03 for the inner disc. Spatially selecting bar stars is complicated by a strong vertical metallicity gradient of -1.1dex/kpc, but we find the metallicity distribution varies in a manner consistent with our orbital selection. Our results have two possible interpretations. The first is that the most metal rich stars in the inner Galaxy pre-existed the bar, but were kinematically cold at the time of bar formation and therefore more easily captured onto bar orbits when the bar formed. The second is that the most metal rich stars formed after the bar, either directly onto the bar following orbits or were captured by the bar after their formation.


Introduction
The vertically extended bulge region of the Milky Way is shaped as a box/peanut bulge (B/P bulge; Lopez-Corredoira et al. 2005;Saito et al. 2011;Wegg & Gerhard 2013). Such shapes naturally arise in simulations of barred galaxies (e.g. Athanassoula 2005;Martinez-Valpuesta et al. 2006), and they are commonly seen in external galaxies (Bureau et al. 2006;Laurikainen et al. 2014). In these simulations and in external galaxies the bar becomes thinner outside of the central B/P bulge (Erwin & Debattista 2013). An equivalent longer, thinner bar structure also exists in the Milky Way outside the bulge; this structure is termed the long bar (Hammersley et al. 1994(Hammersley et al. , 2000. However, because the long bar has comparable thickness to the surrounding disc, with consequently higher extinction, fundamental parameters such as its length and orientation and even relationship to the barred bulge have been difficult to determine (Lopez-Corredoira et al. 2006;Cabrera-Lavers et al. 2008;Wegg et al. 2015).
In this work, we study the chemistry of stars in the long bar of the Milky Way from a spectroscopic sample of ∼2500 stars. From these we select ∼1000 red clump giants (RCGs) whose standard candle nature provides ∼ 10% accurate distances. When combined with Gaia proper motions, this provides a homogeneous sample of long bar stars with full 6D kinematic information. We use this, together with spectroscopic iron abundances, to study the long bar and compare it to the surrounding inner disc. The recent work of Bovy et al. (2019) found, using APO Galactic Evolution Experiment (APOGEE) data on the inner Galaxy, that the bar was more metal poor than the surrounding disc. We study a smaller region of sky, which has a more homogeneous sample, and find contrasting results.
We first describe our sample of stars and their spectroscopic analysis in section 2. In section 3 we describe the construction of our sample of RCGs. In section 4 we make spatial and kine-matic selections of stars that are likely to be bar stars and show the chemistry of these bar stars compared to the inner disc of the Milky Way. In section 5 we discuss the implications of our results.

Target selection
We selected fields in the bar region of the Milky Way, targeting low extinction regions close to the galactic plane, where the bar is most prominent in star counts (see Figure 1). Field centres were placed at longitudes of l = 16.2, 18, 21 deg, each with latitudes of b = −2.0, −2.33, −2.66, −3.0 deg, for a total of 12 fields (see Figure 2).
In order to construct a clean sample of bar stars with accurate distances, we targeted standard candle RCGs for which ∼ 10% accurate distances can be calculated (Bovy et al. 2014). In each field we selected stars whose extinction free K s band magnitude would place them in the bar if they were RCGs. Specifically we targeted stars whose (H − K s ) extinction free magnitude would place then at ∼ 4−8 kpc (i.e. we used eqn 1 of Wegg et al. 2015).
To avoid biasing the sample we did not make a strict colour cut, removing only the small fraction of stars with (H − K s ) < 0.06. These stars are too blue to be RCGs even without extinction. We note that at our chosen latitudes almost all our target RCGs have Gaia astrometry available (Figures 1 and 2   employed, providing a spectral coverage spanning from 8484 to 9001Å with a resolving power of R ∼ 18000. This set-up covers the calcium triplet region and was chosen as it is also a set-up used by the Gaia-ESO Survey (GES) (Gilmore et al. 2012) and therefore facilitates comparison with other regions of the Galaxy.
The sample comprises 2456 spectra of unique stars that have a mean signal-to-noise ratio (S/N) of 40. The raw data were reduced through flat-fielding, extraction, wavelength calibration, and correction to heliocentric reference system by ESO as part of the GIRAFFE stream release.
3 VISTA (Visible and Infrared Survey Telescope for Astronomy) Variables in the Via Lactea 3 United Kingdom Infrared Deep Sky Survey -Galactic Plane Survey 3 Two Micron All-Sky Survey We used the IRAF 4 task skytweak to remove the sky emission lines from the individual spectra by adopting a median sky spectrum computed from individual spectra of dedicated fibres at each pointing. After this, we applied the continuum task to normalise the stellar continuum to the unit through a cubic spline fit. The radial velocity of each individual spectrum was derived by cross correlation against a coarse grid of synthetic templates covering the parameter space of FGK type stars. The radial velocity pipeline selects in a first iteration the best template to the observed spectrum in a χ 2 sense, which is then used to re-compute the final radial velocity.
The set of fundamental parameters (T eff , log(g), [M/H]) and global alpha-elements enhancement [α/Fe] were obtained by comparing the observed spectra against a grid of synthetic models. We adopted the synthetic spectral library of the GES (de Laverny et al. 2012) and the code FERRE. This code performs a global fit search for the set of fundamental parameters providing the best agreement (by minimisation of χ 2 ) between the observed spectrum and the synthetic grid, allowing for interpolation between the models. During the χ 2 minimisation search, both the models and observed spectra are re-normalised by applying a fourth-degree polynomial. This is intended to be a linear and homogeneous process (i.e. without applying any kind of iterative pixel rejection process), thereby decreasing normalisation systematic differences between observed and model spectra. We tested a couple of different strategies for re-normalisation (using a sample of GES spectra; see below) and we found that our results do not change significantly with the adopted strategy, and are of better quality than those obtained without any re-normalisation. A pixel weight mask was constructed to avoid fitting the core of Ca II triplet (CaT) lines whose modelling suffer from larger uncertainties. A number of other small spectral regions were discarded because of their systematically larger residuals, many of which are attributable to residuals of sky line subtraction.
In order to verify the quality of the fitted stellar parameters and abundances obtained with our procedure, we processed a sample of spectra from the GES in the same way as our dataset. We considered stars whose set of recommended fundamental parameters (T eff , log(g), [M/H], [α/Fe]), and metallicity [Fe/H] were available from the third public data release (DR3). From this, we selected a random sample of ∼ 600 bulge giant (log(g) < 3.0) stars observed in the same instrumental set-up, in order to be representative of those we are studying in this work. We then ran our complete pipeline on the retrieved spectra and obtained estimates of T eff , log(g), [M/H], and [α/Fe]. In the upper panels of Figure 3 we compare our set of parameters with those of GES. There is good agreement between both sets of estimates; comparisons follow linear trends reassuringly near to the 1:1 line with small dispersion and systematic scale differences. These plots show the ability of our pipeline to produce results of comparable quality as those of the GES. We also used these comparisons to estimate the small zero points in the stellar parameters compared to GES and these are shown in the figure.   We applied these zero points to calibrate the parameters of our programme stars into the astrophysical scale of the GES survey.
Once calibrated fundamental parameters were obtained for our programme stars, we adopted these parameters to estimate elemental abundances through line-by-line spectrum synthesis calculations using the BACCHUS code (Masseron et al. 2016). In the case of the [Fe/H] iron abundances adopted in this letter, they were computed from the mean of individual estimates over five clean Fe I lines 5 . In the large lower panel of Fig. 3, we compare the results of this procedure applied on the ∼ 600 test stars from the GES. The comparison between our [Fe/H] estimates and those from the GES is linear and essentially unbiased. In consideration of that, we do not apply any zero point calibration to the [Fe/H] estimates of our programme stars. In Fig. 4 we show the distribution of our programme stars in the Hertzsprung-Russell diagram, colour coded by iron abundance.
In a forthcoming paper (Rojas-Arriagada et al. in prep.) we will measure further elemental abundances beyond [Fe/H], and the complete sample and further details on our processing pipeline will be presented.

Red clump giant sample selection
In this paper we focus our analysis on RCGs. We use RCGs because first, they are standard candles with distances accurate to ≈ 10% (Bovy et al. 2014) and second, they form a very homogeneous sample, which makes differential analysis such as whether the bar is more metal rich than the disc, more straightforward. Therefore, from the sample of 2456 stars we select a clean sample of RCGs by making very similar cuts in the (T eff , log(g)) plane to those described by Bovy et al. (2014)   Iron abundance distribution of stars within 0.5 kpc of the bar major axis (blue) compared to those in front and behind the bar (red and green, respectively). Middle panel: Iron abundance distribution of stars in three vertical slices. The slices are chosen so each contains one-third of the sample. Lower panel: Iron abundance distribution of stars spatially in the bar, separated kinematically into bar stars (red) and inner disc stars (blue; see subsection 4.2 for details). Although these stars are in the same region of space, the stars on bar orbits are on average more metal rich than their inner disc counterparts. In black we show the iron abundance distribution of Baade's window measured using the same instrument set-up (Schultheis et al. 2017). The highest metallicities in the bar are similar to those in the bulge, but there is a much smaller fraction of metal poor stars. second removes first ascent giants. We also remove stars a small number of sample stars flagged as having problems with the NIR photometry or spectroscopic fitting.
We cross match this sample of 944 stars to Gaia DR2 in which every sample star has a companion within 0.35" (Figure 1). We remove the 6 clear foreground stars whose parallax places them closer than 3 kpc at 2σ, i.e. those with − 2σ( ) > 1/3 kpc.
We do not use the Gaia parallaxes beyond this foreground cut because they provide limited information at the distance of the bar for our sample. Instead we compute distances to our stars using the distance modulus where we adopt M K s ,RC = −1.61 (Hawkins et al. 2017), (J − K s ) RC = 0.68 and A Ks E(J−K s ) = 0.528 (Gonzalez et al. 2011). An alternative approach would be to use spectro-photometric distances, however for RCGs the population corrections predicted by isochrones are highly uncertain (Girardi 2016;Chan & Bovy 2019). Fortunately, since RCGs are standard candles, these are at the ∼ 0.05 − 0.1 level and so we adopt the more straightforward approach of Equation 2. We also checked that we find similar results using distances computed using spectro-photometric distances derived using the method of Rojas-Arriagada et al. (2017). We also checked that selecting RCGs from APOGEE DR14 and applying Equation 2 gives distances that are consistent with the spectro-photometric distances computed by Queiroz et al. (2018) with a mean absolute deviation of less than 9%. Our final sample comprises 938 RCGs stars with full 6D kinematic information together with spectroscopic iron abundance measurements.

Spatial selection of bar stars
To investigate how the chemistry of the bar compares to that of the inner disc we use the 3D position of our stars to divide the sample into three parts: (i) Those in the bar region, defined to be within ±0.5 kpc of the major axis of the bar. In making this selection we assume R 0 = 8.2 kpc (The GRAVITY Collaboration et al. 2019) and that the bar major axis lies at 27 deg (Wegg & Gerhard 2013) to the Sun-galactic centre line of sight.
(ii) Those in front of this bar region, but still more than 3 kpc distant from the Sun. (iii) Those behind this bar region, but still less than 10 kpc distant from the Sun. For reference, at a longitude of l ∼ 18 deg typical for our sample, these correspond to (i) bar stars with distance 4.6 kpc < d < 5.8 kpc, (ii) stars in front of the bar with 3 kpc < d < 4.6 kpc, and (iii) stars behind the bar with 5.8 kpc < d < 10 kpc.
In the upper panel of Figure 5 we show the iron-abundance distribution of each of these sub-samples. We find that the metallicity distribution function (MDF) of the stars in the bar and in front of the bar are very similar, while those behind the bar are more metal poor. Quantitatively the mean [Fe/H] of the bar stars is 0.16 ± 0.01 dex compared to 0.18 ± 0.02 dex for those in front of the bar and 0.06 ± 0.02 dex for those behind.
However, we also find a very strong vertical metallicity gradient. We show the distribution of [Fe/H] in three vertical slices in Figure 5. The population of stars changes as a function of distance to the galactic plane; a higher fraction of more metal poor stars is found in the slice furthest from the galactic plane. Fitting a linear relation between [Fe/H] and distance from the galactic plane across the entire sample, the vertical gradient is −1.1 ± 0.2 dex kpc −1 .
Because we use pencil beam fields just away from the galactic plane, our more distant stars are also on average further from the galactic plane. We must therefore account for this vertical gradient to compare the [Fe/H] distribution in the bar to the surrounding inner disc. Using the measured vertical gradient of −1.1 dex kpc −1 to correct our entire sample to 200 pc from the galactic plane, we find the mean [Fe/H]  Lower panels: The mean metallicity (left) and fraction of super solar metallicity stars (right). In the middle and lower panels we also plot the model circular velocity at these galactocentric distances (red) and orbits with zero angular momentum in a frame rotating at Ω = 35 km s −1 kpc −1 (orange). These zero angular momentum orbits correspond to an pencil thin bar and have galactocentric velocity -80 km s −1 to 80 km s −1 in steps of 20 km s −1 . In the model, bar orbits have kinematics near these low angular momentum orbits, and the stars in this region of kinematic space are preferentially metal rich. region is 0.18 ± 0.01 dex compared to 0.15 ± 0.02 dex for those in front of the bar and 0.15 ± 0.02 dex for those behind.
We conclude that, at the same distance from the galactic plane, the bar region in our sample is slightly more metal rich than the surrounding disc. We note that this conclusion is the opposite from that found by Bovy et al. (2019), a disagreement we discuss further in section 5.

Kinematic selection of bar stars
In this section we present the 2D kinematics of bar stars and make a kinematic selection of stars whose orbits indicate they follow the bar. To guide our selection and methods we use a barred made-to-measure Milky Way model of the Galaxy, which was fitted to a range of data on the inner Galaxy by Portail et al. (2017, hereafter P17). We note however, that this model was fit only to star count data in the long bar (the region studied in this work) because of the lack of kinematic data at that time. We choose to show a model with a pattern speed of Ω = 35 km s −1 kpc because its kinematics agrees qualitatively slightly better than the fiducial model (Ω = 40 km s −1 kpc), however our conclusions would remain unchanged with any model from P17. To compare this model to our data we use the Galaxia code (Sharma et al. 2011) to generate mock catalogues of stars generated from the Padova isochrones with a 10 Gyr old population. For our fiducial MDF we use a Gaussian with [Fe/H] = 0, σ([Fe/H]) = 0.3. To these mock catalogues we apply our selection criteria, both spatially and in the (T eff , log(g)) plane, to generate a mock sample of RCG stars. We also compute distances to this mock catalogue using Equation 2 in the same manner as the real sample. In the mock sample the mean absolute deviation of this RCG distance to the true distance is 6.1%. To replicate the treatment of the real sample we use these RCG distances when using the mock catalogue in what follows.
In the upper left panel of Figure 6 we show the kinematics of the stars selected to be in the bar region in the longitudinal proper motion versus radial velocity plane (i.e. µ l , V r plane). This 2D plane encodes similar information as the well-known (U, V) plane locally. We compare this to the mock sample finding qualitatively very good agreement. The sample is centred at µ l ≈ −6 mas yr −1 similar to the µ l = −6.24 mas yr −1 of Sgr A* (Reid & Brunthaler 2004). This is largely due to the motion of the Sun with respect to the galactic centre and it is reassuring for the purity of our sample that there is a complete lack of stars at µ l > 0 mas yr −1 where nearby foreground contaminants would lie.
Using our mock sample of stars we can investigate the kinematic signatures of the bar in this 2D kinematic plane. We integrate the model stars and define long bar stars as those which spend more than 40% of their time in the long bar region; that is within ±0.5 kpc of the bar major axis and more than 2 kpc from the galactic centre. We also used the orbital frequency based definition from Portail et al. (2015) (specifically f r / f x ∈ 2±0.1) finding very similar results. However we use the orbital time definition because it both better selects particles which were humanclassified to be on bar orbits, and better photometrically separates bar stars from inner disc stars in face-on projections.
We split the model based on this orbital selection of bar stars in the middle panels of Figure 6. In these panels we also plot the position of the circular velocity for our l = 18 deg longitudinal field (the red line), together with the orbits of stars travelling directly along the major axis of the bar, i.e. with zero angular momentum in the co-rotating frame (orange lines). These pencil thin bar orbits cover a galactocentric radial velocity range -80 km s −1 to 80 km s −1 in steps of 20 km s −1 . As expected, the model bar orbits lie in the region of these pencil thin orbits with a scatter due to the width of the bar and its underlying orbits (Valluri et al. 2016).
In the lower panel of Figure 6 we plot the iron abundance of our sample in the same plane. In the region of kinematic space where the bar stars lie, both the mean metallicity and the fraction of stars with super-solar abundance ([Fe/H] > 0) appear significantly higher.
To quantify these differences we use the µ l , V r plane shown in Figure 6 to assign each of our sample stars a probability, p bar , that it is on a bar following orbit. We do this by for each sample star computing the fraction of bar following orbits among the 100 nearest neighbour mock sample stars in this plane. This provides a sample of stars with ([Fe/H],p bar ). We use this to plot in the lower panel of Figure 5 the MDF of stars with bar orbits, P([Fe/H]|p bar = 1), and inner disc orbits, P([Fe/H]|p bar = 0). We compute this using a conditional kernel density estimate (KDE) P([Fe/H]|p bar ) = P([Fe/H], p bar )/P(p bar ), where both probability distributions on the right-hand side are estimated us-