Life in the fast lane: a direct view of the dynamics, formation, and evolution of the Milky Way's bar

TL;DR

This study uses APOGEE and Gaia data to map the Milky Way's bar, revealing its formation about 8 billion years ago and providing detailed insights into its dynamics, composition, and evolution.

Contribution

It offers the first detailed kinematic and compositional maps of the Milky Way's bar, determining its pattern speed, shape, and age with unprecedented precision.

Findings

The bar's pattern speed is 41+/-3 km/s/kpc.

The bar formed approximately 8 billion years ago.

Stars in the bar share similar ages and metallicities with the oldest disk stars.

Abstract

Studies of the ages, abundances, and motions of individual stars in the Milky Way provide one of the best ways to study the evolution of disk galaxies over cosmic time. The formation of the Milky Way's barred inner region in particular is a crucial piece of the puzzle of disk galaxy evolution. Using data from APOGEE and Gaia, we present maps of the kinematics, elemental abundances, and age of the Milky Way bulge and disk that show the barred structure of the inner Milky Way in unprecedented detail. The kinematic maps allow a direct, purely kinematic determination of the bar's pattern speed of 41+/-3 km/s/kpc and of its shape and radial profile. We find the bar's age, metallicity, and abundance ratios to be the same as those of the oldest stars in the disk that are formed in its turbulent beginnings, while stars in the bulge outside of the bar are younger and more metal-rich. This implies that the bar likely formed ~8 Gyr ago, when the decrease in turbulence in the gas disk allowed a thin disk to form that quickly became bar-unstable. The bar's formation therefore stands as a crucial epoch in the evolution of the Milky Way, a picture that is in line with the evolutionary path that emerges from observations of the gas kinematics in external disk galaxies over the last ~10 Gyr.

Equations

$$\begin{array}{rl}0.99+\frac{0.060}{\pi }& \left[\frac{\pi }{2}+\arctan \left(\frac{\tilde{D}/\mathrm{kpc}-4.35}{0.30}\right)\right]\\ +\frac{0.055}{\pi }& \left[\frac{\pi }{2}+\arctan \left(\frac{\tilde{D}/\mathrm{kpc}-5.50}{0.15}\right)\right],\end{array}$$

$$-\frac{\partial \Sigma }{\partial t}={\Omega }_b\frac{\partial \tilde{\Sigma }}{\partial \varphi }=\frac{1}{R}\frac{\partial (R{v}_R\tilde{\Sigma })}{\partial R}+\frac{1}{R}\frac{\partial (v_T\tilde{\Sigma })}{\partial \varphi },$$

$$\frac{\partial \ln \tilde{\Sigma }}{\partial R}=-\frac{1}{R{v}_R}\left(v_R+R\frac{\partial v_R}{\partial R}+\frac{\partial v_T}{\partial \varphi }\right).$$

$$\frac{\partial \ln \tilde{\Sigma }}{\partial R}$$

$$\frac{\partial \ln \tilde{\Sigma }}{\partial \varphi }$$

$$\begin{array}{rl}{\Omega }_b=\frac{v_T}{R}& +\frac{2q^2{R}_a}{(1-q^2)R^2\sin 2\varphi \partial \ln \tilde{\Sigma }/\partial R_a}\left(\frac{1}{R}\frac{\partial v_T}{\partial \varphi }+\frac{v_R}{R}+\frac{\partial v_R}{\partial R}\right)\\ & +\frac{2v_R\left(q^2{\cos }^2\varphi +{\sin }^2\varphi \right)}{(1-q^2)R\sin 2\varphi }.\end{array}$$

Full text

††pubyear: 2019††pagerange: Life in the fast lane: a direct view of the dynamics, formation, and evolution of the Milky Way’s bar–A

1 Introduction

The Milky Way’s bar was discovered in the distribution of near-infrared emission (Blitz & Spergel, 1991) and in the kinematics of gas (Binney et al., 1991) in the central regions. However, the large and variable amount of interstellar extinction (Nataf et al., 2013) has made direct studies of the bar difficult in the 30 years since its discovery. Measurements of stars in fields typically away from the Galactic mid-plane to avoid the effects of interstellar extinction demonstrate that the bar dominates the kinematics in the inner few kpc (Shen et al., 2010; Ness et al., 2013b), leaving little room for the classical bulge component that is expected to result from the hierarchical galaxy formation paradigm (Kormendy & Kennicutt, 2004). Bars rotate as solid bodies and the bar’s angular pattern speed is its most important property aside from its mass for determining its influence on, e.g., the observed perturbed kinematics of stars in the Galactic disk (Antoja et al., 2018; Hunt & Bovy, 2018; Fragkoudi et al., 2019; Fujii et al., 2019) or on the structure of stellar streams in the halo and this structure’s interpretation in terms of dark-matter substructure (Pearson et al., 2017; Banik & Bovy, 2019).

The distribution of stellar ages and chemical abundances reveals the formation and evolution of the Milky Way’s stellar components. Originally thought to be primarily old (Zoccali et al., 2003), spectroscopic observations of stars in the bulge region have revealed a complex mix of populations spanning a wide range of ages, metallicities, and abundance ratios (Bensby et al., 2013; Ness et al., 2013a). However, directly separating the bar populations from the inner disk, inner halo, and spheroidal bulge components has been difficult using past observations and the age of the bar and the chemical evolution of its stars remains largely unknown. The chemical evolution of galactic components is traced by their distribution in the $([\mathrm{Fe/H}],[\mathrm{O/Fe}])$ plane, because oxygen is mainly produced through type II supernovae that occur soon after the commencement of star formation, while iron has a large contribution from type Ia supernovae, which lag star formation by typically $\approx 1\,\mathrm{Gyr}$. Determinations of detailed abundances of stars in the bulge show that the bulge follows a chemical evolution track similar to that of the local old, thicker disk component (e.g., Bensby et al., 2010), indicative of a high star-formation efficiency, but it remains uncertain how consistent the bulge and old disk chemical evolution really is.

In this paper, we combine two powerful new data sets—APOGEE (Majewski et al., 2017) and Gaia (Gaia Collaboration et al., 2018)—to for the first time directly reveal the bar in spatial maps of the kinematics, chemical abundances, and ages of stars in the Milky Way. We demonstrate that these maps allow a direct, purely kinematic determination of the bar’s pattern speed and its variation with radius, the bar’s radial profile, and the bar’s shape using a novel application of the continuity equation. Furthermore, the fact that we clearly see the bar in the abundances and ages of stars lets us unambiguously separate the bar from the inner bulge and disk, which reveals the bar’s chemical evolution history.

2 Data

2.1 APOGEE data

We use data from the 16th data release (DR16) of the Apache Point Observatory Galactic Evolution Experiment (APOGEE; Majewski et al. 2017). APOGEE is a high-resolution ($R\approx 22,500$), high signal-to-noise ratio (typically $\gtrsim 100$), near-infrared ($1.5$ to $1.7\mu\mathrm{m}$), large spectroscopic survey (Nidever et al., 2015; Blanton et al., 2017; Zasowski et al., 2017; Holtzman et al., 2018; Wilson et al., 2019). DR16 contains 473,038 stars observed using the APO 2.5m telescope (Gunn et al., 2006) and the DuPont telescope at Las Campanas Observatory (LCO). The dual-hemisphere, full-sky coverage is crucial for mapping the center of the Milky Way and for producing a consistent data set covering the inner bulge to the outer disk.

We analyze the data using a custom set of analysis pipelines that is based on, but separate from, the normal APOGEE pipeline (Shetrone et al., 2015; Zamora et al., 2015; García Pérez et al., 2016). We start from the one-dimensional, rest-frame corrected, combined-exposure spectra in the apStar files. We perform a (pseudo-)continuum normalization to remove the effect of the overall continuum on stellar abundance, age, and luminosity measurements, by fitting a second-order polynomial to a set of wavelength pixels previously determined to trace the continuum, for each of the three APOGEE detectors separately. To determine abundances, we then run these continuum-normalized spectra through the astroNN trained neural-network (Leung & Bovy, 2019a) to determine the stellar parameters effective temperature $T_{\mathrm{eff}}$ and surface gravity $\log g$ and to determine the elemental abundances of iron Fe, oxygen O, as well as other elements that are not used here. The DR16 spectral reduction imprints an overall difference between the DR16 apStar spectra that we use here and the DR14 spectra that the neural network was trained on. That is, even when the raw data underlying a given apStar spectrum is the same, differences between the DR14 and DR16 spectral reduction cause a small difference in the reduced apStar spectrum. Moreover, this difference is different for DR16 spectra observed at APO and at LCO (DR14 contained only APO observations). To apply the previously trained neural network, we determine the median difference between continuum normalized spectra in DR16 (separately for APO and LCO) using the set of spectra observed both from APO and from LCO to place all spectra onto the same footing. Comparing stellar parameters and Fe and O abundances for stars in the APO/LCO overlap sample shows results consistent between DR14 and DR16 APO (which are the same underlying spectra) to a level of $15\,\mathrm{K}$ in $T_{\mathrm{eff}}$, $0.035\,\mathrm{dex}$ in $\log g$, $0.015\,\mathrm{dex}$ in $[\mathrm{Fe/H}]$, and $0.018\,\mathrm{dex}$ in $[\mathrm{O/Fe}]$; between spectra of the same star taken at APO and at LCO, differences are at the level of $40\,\mathrm{K}$ in $T_{\mathrm{eff}}$, $0.1\,\mathrm{dex}$ in $\log g$, $0.04\,\mathrm{dex}$ in $[\mathrm{Fe/H}]$, and $0.05\,\mathrm{dex}$ in $[\mathrm{O/Fe}]$. These differences are smaller or similar to the typical uncertainties in the derived parameters of $60\,\mathrm{K}$ in $T_{\mathrm{eff}}$, $0.2\,\mathrm{dex}$ in $\log g$, $0.06\,\mathrm{dex}$ in $[\mathrm{Fe/H}]$, and $0.05\,\mathrm{dex}$ in $[\mathrm{O/Fe}]$.

We determine distances to the stars in APOGEE DR16 using a previously trained astroNN neural-network (Leung & Bovy, 2019b) that determines the luminosity of giants from continuum-normalized APOGEE spectra and derives distances from these by combining this with the observed apparent $K_{s}$ magnitude from 2MASS (Skrutskie et al., 2006) and extinction values determined using mid-infrared photometry with the RJCE method (Majewski et al., 2011). Below, we determine that our distances are systematically biased low by 5% to 10% near the mid-plane at distances of 4 kpc and beyond and we correct this bias using its determination below. We convert the 3D heliocentric positions to Galactocentric coordinates by setting the distance to the Galactic center to $8.125\,\mathrm{kpc}$ (Gravity Collaboration et al., 2018) and the Sun’s height above the plane to $20.8\,\mathrm{pc}$ (Bennett & Bovy, 2019).

To determine ages, we re-train the neural network that we used to determine ages using DR14 data (Mackereth et al., 2019) on the overlap between stars in the APOKASC (Pinsonneault et al., 2018) sample with asteroseismic ages and APOGEE DR16, using the DR16-corrected spectra as explained above. As in our previous work, limitations in the training data set make it such that we cannot obtain ages for stars with $[\mathrm{Fe/H}]<-0.5$, but only a small fraction of stars that we consider in the bulge and disk region have such low metallicities and this limitation does not impact our results. As in Mackereth et al. (2019), the raw ages determined by the neural network at ages $\gtrsim 6\,\mathrm{Gyr}$ are underpredicted, because the reduced sensitivity of carbon and nitrogen spectral features to age at high age causes the predictions from the regularized neural-network used to tend to the mean age of the sample. We measure this bias on the training sample, smoothly interpolate the bias using a Locally Weighted Scatterplot Smoothing (lowess), and apply the lowess bias correction to all predicted ages by the neural network. For our data set, this correction has the effect of shifting all ages towards larger ages, but it has little effect on the comparison of ages of stars in the bar to those of stars outside the bar, because the ages of stars in these different samples are corrected by a similar amount. The typical age precision is $\approx 30\,\%$.

For all kinematic, abundance, and age maps and distributions that we determine, we select stars that have (a) measured values of the relevant quantity, (b) astroNN uncertainty in $\log g<0.2$, which removes dwarf stars (Leung & Bovy, 2019a), and (c) relative astroNN distance uncertainty less than 20 %. For the $[\mathrm{Fe/H}]$ maps we further limit to stars with uncertainties $<0.05\,\mathrm{dex}$ and for the $[\mathrm{O/Fe}]$ maps we use the $[\mathrm{Fe/H}]$ cut as well as requiring the uncertainty in $[\mathrm{O/H}]$ to be $<0.05\,\mathrm{dex}$.

2.2 Gaia data

To determine the 3D kinematics of stars, we match the APOGEE data to the Gaia catalog (Gaia Collaboration et al., 2018) using an angular matching radius of 2 arcseconds and record the proper motions, parallaxes, and their uncertainties (the parallaxes are only used to determine the systematic distance bias below). We combine the proper motions with the line-of-sight velocities determined from the APOGEE spectra to obtain full 3D heliocentric velocities. We convert these to Galactocentric cylindrical coordinates using the same Galactic center and mid-plane distance as used for the positions above, and by using a solar motion of 11.1, 242., and $7.25\,\mathrm{km\,s}^{-1}$ in the radial, rotational, and vertical direction, respectively (Schönrich et al., 2010; Bovy et al., 2012c).

2.3 Distance biases

While the astroNN neural-network to determine distances returns precise distances to stars at large distances from the Sun, even a small relative bias will cause significant distortion of the observed structure at the large distance of the bar from the Sun. We determine whether there is a systematic bias in the distances determined using the neural network in the following manner. We select stars within 300 pc from the Galactic mid-plane based on their astroNN distances. To avoid outliers in the Gaia parallax measurements, we include stars with negative parallaxes only if they are larger than minus three times their uncertainty (but note that we do include negative parallaxes). We then assume that if astroNN determines that stars are at the same distance, that they truly are at the same distance, which is a reasonable assumption because a neural-network regression is very good at finding objects that have similar values of the output quantity (distance in this case) even if the output is biased. We then select stars at the same astroNN distance and determine the true distance $D$ of these stars by a maximum-likelihood fit to the parallax data, minimizing a robust objective function $\sum_{i}|\varpi_{i}-1/D|/\sigma_{\varpi,i}$, where $\varpi_{i}$ and $\sigma_{\varpi,i}$ are the parallax and its uncertainty of star $i$. We determine the uncertainty on $D$ using bootstrap resampling. In doing this fit, we add $52\,\mu\mathrm{as}$ to all Gaia parallaxes, the best-fit global parallax zero-point that we determined in a joint fit of the luminosity neural-network model and the zero-point (Leung & Bovy, 2019b).

The result of applying this procedure to stars in distance bins with a width of 200 pc is displayed in Figure 1. We see that out to $4\,\mathrm{kpc}$, the astroNN distances have a small bias of 1 % in the sense of the astroNN distances being too large. But beyond $4\,\mathrm{kpc}$, the astroNN distances are too small by a fraction that grows to 10 % at $6\,\mathrm{kpc}$ and that stays constant at larger distances. The bias determined using stars towards the direction of the bulge alone agrees with the measurement using all stars. The source of this bias is likely related to the paucity of training data at the low-temperature, high-luminosity end that makes up the large-distance sample in the Galactic mid-plane: in regions of parameter space with few training data, regularized neural networks such as astroNN (which uses drop-out regularization) tend to regress to the mean, which in our case means that predicted luminosities are too low for these stars, thus underestimating their true distance.

We approximate the behavior shown in Figure 1 using the function

[TABLE]

where $\tilde{D}$ is the raw astroNN distance. This function is shown as the smooth curve in Figure 1. We correct all distances that we use in this paper using this function.

3 The kinematics in the inner Galaxy and the bar’s pattern speed

In Figure 2, we display the median cylindrical rotational and radial velocities of 73,189 stars observed by APOGEE and Gaia within 300 pc of the Galactic mid-plane as well as their rotational frequencies. These maps clearly show the signature of the Galactic bar: the rotational velocity map within 5 kpc from the center clearly demonstrates the expected bar-shape, with the major axis at $25^{\circ}$ clockwise from the Sun–Galactic center line, the radial velocity shows the tell-tale quadrupolar inwards and outwards movement on either side of the major axis, and the rotational frequency map displays the expected dip along the major axis. In Figure 3, we show these same maps for a barred, Milky-Way-like $N$-body simulation (Kawata et al., 2017), which illustrates the striking agreement between the observations and the expected bar kinematic pattern.

The kinematic maps from Figure 2 allow for a purely-kinematic, model-free determination of the bar’s pattern speed, surface density profile, and shape through a novel application of the continuity equation. Our method is similar to the classical method of Tremaine & Weinberg (1984), but because we employ full phase-space information, we do not have to assume that the pattern speed is constant over the bar region. Assuming that the surface density $\Sigma(X,Y,t)$ is constant in a frame rotating with the bar’s pattern speed $\Omega_{b}$, $\Sigma(X,Y,t)\equiv\tilde{\Sigma}(R,\phi-\Omega_{b}t)$, the two-dimensional continuity equation becomes

[TABLE]

where $(R,\phi)$ are polar coordinates and all velocities are mean velocities in a small patch ($400\,\mathrm{pc}^{2}$ in Figure 2). Near the major axis, $\partial\tilde{\Sigma}/\partial\phi\approx 0$, which using Equation (2) allows us to kinematically determine the radial surface-density profile along the major axis as

[TABLE]

In practice, we assume that $\partial\ln\tilde{\Sigma}/\partial R$ is constant and find the best value by minimizing the median absolute difference between the left- and right-hand side of this equation. That we are able to determine the spatial profile of the inner Milky Way purely kinematically is important, because directly measuring the spatial profile by mapping the relative density of stars in our sample is complicated by the large and varying amount of interstellar extinction towards the inner Galaxy and the complex APOGEE selection function (e.g., Bovy et al. 2016a).

If we then assume that $\tilde{\Sigma}$ is stratified on similar, concentric ellipses with axis ratio $q$, all of the derivatives of $\tilde{\Sigma}$ at any $(R,\phi)$ in Equation (2) can be written in terms of $\partial\ln\tilde{\Sigma}/\partial R$ along the major axis and $q$ and we can determine $\Omega_{b}$ simply by combining the surface-density profile and a value of $q$ with the observed mean $v_{R}$ and $v_{T}$ velocities and their gradients. We re-write Equation (2) such that $\tilde{\Sigma}$ only appears as its logarithmic derivative along the major axis, $\partial\ln\tilde{\Sigma}/\partial R_{a}$, which we kinematically determine using Equation (3). We use $R_{a}$ to denote the location along the major axis for a given point $(R,\phi)$. Assuming that the surface density is constant along concentric, similar ellipses with an axis ratio $q$ and a major axis along $\phi=0$, the function $R_{a}(R,\phi)=R\,\sqrt{\cos^{2}\phi+\sin^{2}\phi/q^{2}}$. Therefore,

[TABLE]

Working out the continuity equation in terms of these derivatives then gives the pattern-speed estimator

[TABLE]

The derivative of the continuity equation similarly allows us to determine the axis ratio $q$ itself kinematically. We take the $\phi$ derivative of the continuity equation with the substitutions from Equations (4) and (5), eliminate $\Omega_{b}$ using the continuity equation, and set $\partial^{2}\ln\tilde{\Sigma}/\partial R_{a}^{2}=0$, that is, we assume that the bar’s surface density profile is exponential. This is a good assumption, because the Milky Way bar’s outer density profile constrained using red-clump star counts is found to be exponential (Wegg et al., 2015) and we see the same in the Milky-Way-like, barred $N$-body simulations below. This procedure leads to a complex equation that relates $q$ in a non-linear manner to $\partial\ln\tilde{\Sigma}/\partial R_{a}$ and to the first and second derivatives of the observed kinematics $v_{R}$ and $v_{T}$. For each spatial bin, we can simply solve this equation by evaluating it on a grid of $q$ and finding the $q$ that best satisfies the equation.

When applying Equation (3), Equation (6), and the equation for $q$ to the data and to the simulations below, we determine the first and second derivatives of $v_{R}$ and $v_{T}$ using simple finite-differencing of the maps in Figure 2, which gives the gradients in rectangular coordinates, and convert these to the gradients in cylindrical coordinates. We then apply the estimators for $\partial\ln\tilde{\Sigma}/\partial R_{a}$, $q$, and $\Omega_{b}$ to those pixels for which we can determine all of the relevant gradients within the region of interest (within 5∘ of the major axis when determining $\partial\ln\tilde{\Sigma}/\partial R_{a}$ and beyond 5∘ of the major axis when determining $q$ and $\Omega_{b}$, at $2<R/\mathrm{kpc}<5$ in both instances). This gives an estimate for each pixel within the region of interest. We robustly combine these by taking the median and estimate the random uncertainty using the median absolute deviation, scaled by 1.4826 to convert the uncertainty to an effective Gaussian uncertainty.

Applying the procedure above to the data in Figure 2, we determine that $\partial\ln\tilde{\Sigma}/\partial R=-1/(1.5\,\mathrm{kpc})$ and $q=0.41\pm 0.07$. Using Equation (6) then gives $\Omega_{b}=41\pm 1.5\,\mathrm{km\,s}^{-1}\,\mathrm{kpc}^{-1}$. Applying the procedure to various barred, Milky-Way-like $N$-body simulations (see Appendix) recovers their pattern speeds to within this uncertainty and we can therefore be confident that the procedure is robust. Our method provides a measurement of the pattern speed at each spatial location and we can therefore investigate the spatial dependence of the measured pattern speed. Figure 4 displays the radial dependence of the pattern speed, determined by combining the measurements in $\Delta R=0.33\,\mathrm{kpc}$ bins with at least 3 spatial pixels between $2\,\mathrm{kpc}$ and $5\,\mathrm{kpc}$; the combined measurement from the entire radial range is shown as the horizontal line with the gray-shaded error region. These radial measurements of $\Omega_{b}$, within their uncertainties, are consistent with the bar pattern speed being constant over $2\,\mathrm{kpc}<R<5\,\mathrm{kpc}$.

In Figure 4, we also include the location of the corotation resonance and of the 4:1 and 2:1 outer Lindblad resonances in the models of the Milky Way’s gravitational potential from Irrgang et al. (2013; model I), Bovy (2015), and McMillan (2017) as implemented in galpy (Bovy, 2015). For a bar pattern speed of $41\,\mathrm{km\,s}^{-1}\,\mathrm{kpc}^{-1}$, the 2:1 and 4:1 outer Lindblad resonance both strongly affect the kinematics of stars near the Sun (e.g., Hunt & Bovy, 2018; Hunt et al., 2019; Monari et al., 2019; Trick et al., 2019).

The pattern speed that we determine depends systematically on the parameters of the transformation from heliocentric to Galactocentric coordinates and on the value of $\partial\ln\tilde{\Sigma}/\partial R_{a}$ and $q$ that we determine. Varying the assumed value of the Sun’s distance to the Galactic center between $8$ and $8.25\,\mathrm{kpc}$ leads to values of $\Omega_{b}$ in the range $38.5$ to $41\,\mathrm{km\,s}^{-1}\,\mathrm{kpc}^{-1}$. Varying the Sun’s rotational velocity with respect to the Galactic center—the most uncertain component and the one with the largest impact on our $\Omega_{b}$ determination—by $\pm 5\,\mathrm{km\,s}^{-1}$ changes $\Omega_{b}$ by $\pm 1.25\,\mathrm{km\,s}^{-1}\,\mathrm{kpc}^{-1}$. If rather than determining $\partial\ln\tilde{\Sigma}/\partial R_{a}$ from the data, we set it to $-1/(1\,\mathrm{kpc})$, we determine $\Omega_{b}=42.5\,\mathrm{km\,s}^{-1}\,\mathrm{kpc}^{-1}$; at the other extreme, $\partial\ln\tilde{\Sigma}/\partial R_{a}=-1/(2.5\,\mathrm{kpc})$ leads to $\Omega_{b}=36.5\,\mathrm{km\,s}^{-1}\,\mathrm{kpc}^{-1}$. Similarly, setting $q$ to $0.3$ or $0.5$ gives $\Omega_{b}=43$ and $38\,\mathrm{km\,s}^{-1}\,\mathrm{kpc}^{-1}$, respectively. Considering these systematics, the recovery of the pattern speed in simulations to within $1.5\,\mathrm{km\,s}^{-1}\,\mathrm{kpc}^{-1}$ (see Appendix) and the statistical uncertainty in our measurement, our data imply that $\Omega_{b}$ is contained in the interval $[36.5,42.5]\,\mathrm{km\,s}^{-1}\,\mathrm{kpc}^{-1}$ with high confidence. To summarize this with a single value and uncertainty, we obtain $\Omega_{b}=41\pm 3\,\mathrm{km\,s}^{-1}\,\mathrm{kpc}^{-1}$, but detailed comparisons with other determinations should take note of the systematics that lead to an asymmetric uncertainty interval .

4 The formation and chemical evolution of the bulge

In Figure 5, we present maps of the median age, iron abundance $[\mathrm{Fe/H}]$, and oxygen-to-iron ratio $[\mathrm{O/Fe}]$ of stars within 300 pc of the Milky Way’s mid-plane; these maps are made with a total of 74,041, 67,308, and 62,970 stars, respectively. As in the kinematics, the bar clearly shows up in these maps, with a population that is old, metal-poor, and enhanced in oxygen compared to stars outside of the bar and in the disk. Contours of these maps clearly reveal a bar with a major axis inclined by $\approx 25^{\circ}$ clockwise from the Sun–Galactic center line in agreement with the kinematic maps of Figure 2. For example, fitting a model for the two-dimensional $[\mathrm{Fe/H}]$ that consists of concentric ellipses with a varying axis ratio, reveals a transition between an axis ratio of 0.45 near the center (in good agreement with the kinematically-determined axis ratio above) to an axis ratio of unity at $R\gtrsim 5\,\mathrm{kpc}$, therefore the size of the bar.

In Figure 6, we display the metallicity distribution and $([\mathrm{Fe/H}],[\mathrm{O/Fe}])$ track of stars in the inner Galaxy both inside and outside the bar. We select stars inside the bar as those stars within an ellipse with an axis ratio of $q=0.4$ and a semi-major axis of $5\,\mathrm{kpc}$ oriented at 25∘ clockwise from the Sun–Galactic center line (see Figure 2); stars outside the bar are those with $R<5\,\mathrm{kpc}$ that are outside of this ellipse. It is clear that stars in the inner Galaxy follow the same track, whether or not they are inside the bar, but the bar has an excess of more metal-poor stars with $[\mathrm{Fe/H}]<0$ compared to stars outside the bar. Because the barred region contains a mix of bar members and members from the axisymmetric, inner-Galaxy component traced by the ‘outside the bar’ distribution, this shows that the bar is primarily made up of metal-poor stars. In Figure 7, we similarly display the age distribution of stars inside and outside the bar in the central regions. The relation between the age distribution and the star-formation history is complex for the population of evolved, giant stars that we use as tracers. But, because the types of giants that we employ inside and outside the bar are the same, their relative formation epochs can be robustly determined. The bar is on average older than the rest of the inner Galaxy.

To connect the bar’s evolution to that of the Galactic disk, we compare the bar’s abundance and age distributions to those of the old stellar component in the Solar neighborhood, as traced by stars currently located more than 500 pc and less than 2 kpc from the mid-plane within $5<R/\mathrm{kpc}<10$, taking care to only select giants with temperatures below $4,600\,\mathrm{K}$ to match the temperatures of the stars observed in the bulge and avoid systematic biases in comparing ages and abundances determined for different stellar types. The right panel of Figure 6 compares the $([\mathrm{Fe/H}],[\mathrm{O/Fe}])$ track of stars in the bar to that of old stars in the Solar neighborhood and we find them to be the same. In Figure 7, we compare the difference between the age distributions inside and outside the bulge to the age distribution of old Solar-neighborhood stars and find that they match remarkably well. Thus, the star-formation and chemical-evolution history of stars in the central barred region is the same as that of old stars near the Sun.

5 Conclusion

In this paper, we have presented the most complete maps of the kinematics, chemical abundances, and ages of stars in the Milky Way bulge and disk. For the first time, these maps directly reveal that the Milky Way disk–bulge system is barred not only in its density profile, but also in its kinematics, chemistry, and ages. These maps herald a new era in which we can unambiguously separate the properties of stars and their orbits in the bar from those in the inner disk, bulge, and stellar halo.

We have used the maps of the kinematics of stars in the inner 5 kpc to determine the bar’s pattern speed, radial profile, and shape in a purely kinematic manner using a simple, novel application of the continuity equation. In particular, we have determined that the bar’s pattern speed $\Omega_{b}=41\pm 3\,\mathrm{km\,s}^{-1}\,\mathrm{kpc}^{-1}$. Such a pattern speed places the bar’s corotation radius at 5.1 to $5.9\,\mathrm{kpc}$ and the bar with size of $5\,\mathrm{kpc}$ therefore extends out to $\gtrsim 90\,\%$ of its corotation radius, making the Milky Way’s bar fast from a dynamical perspective. Prior measurements of the pattern speed, all based on incomplete data and/or involving complex and uncertain data modeling, range from $\approx 30-40\,\mathrm{km\,s}^{-1}\,\mathrm{kpc}^{-1}$ (Weiner & Sellwood, 1999; Rodriguez-Fernandez & Combes, 2008) to $\approx 60\,\mathrm{km\,s}^{-1}\,\mathrm{kpc}^{-1}$ (Englmaier & Gerhard, 1999), with recent measurements (Portail et al., 2017; Sanders et al., 2019) clustering around $40\,\mathrm{km\,s}^{-1}\,\mathrm{kpc}^{-1}$, in good agreement with our simple determination here. Our measurement is the first determination of the bar’s pattern speed in the bar’s outer region, as previous measurements are largely based on data at $R\lesssim 2\,\mathrm{kpc}$, and we determine the radial dependence of the bar’s pattern speed for the first time (see Figure 4). The fact that our determination of the pattern speed at $2\,\mathrm{kpc}<R<5\,\mathrm{kpc}$ is consistent with the recent measurements at $R<2\,\mathrm{kpc}$ directly shows that the entire inner Milky Way rotates as a solid body, as expected for a rigid, stable bar.

The picture that emerges from the observations of the kinematics, chemistry, and ages of stars in the Milky Way presented in this work is that the bar formed as the early chemical evolution of the Milky Way’s disk reached approximately solar abundances after a few Gyr of evolution. Studies of the age (Mackereth et al., 2017), abundance (Bovy et al., 2012a, 2016b), and kinematics (Bovy et al., 2012b; Mackereth et al., 2019) of stars in the disk indicate that the early gas disk was thick due to turbulence, as also directly observed in more massive disk galaxies at redshift $\gtrsim 1$ (Förster Schreiber et al., 2009), stabilizing it against bar formation. Direct measurements of the structure of different stellar disk components (Bovy et al., 2016b) show that as stars reached solar abundance ratios, the disk thickness had decreased to only a few hundred parsecs, at which time a bar instability could develop and the distribution of stars formed in the first few Gyr was re-arranged into a bar. Our direct observations of the chemical evolution of the bar show that this happened quickly. Chemical evolution in the inner disk then proceeded along the same chemical track to higher metallicities, while the bar rotated steadily as evidenced by our determination of its pattern speed, which is high for its size, indicating that it has not slowed down much. Our observations reveal the formation of the bar as a crucial epoch in the evolution of the Milky Way, separating the early turbulent epoch from the later, quiescent evolution (Haywood et al., 2018). A similar transition is hinted at (Kraljic et al., 2012) in the cosmic evolution of the bar fraction (Sheth et al., 2008) and disk turbulence (Wisnioski et al., 2015), but can now be directly observed in the Milky Way and placed into the wider disk-formation context.

Acknowledgements

It is a pleasure to thank Francesca Fragkoudi, Ortwin Gerhard, David Nataf, Scott Tremaine, Benjamin Weiner, and the anonymous referee for helpful comments. J. Bovy and H. Leung received support from the Natural Sciences and Engineering Research Council of Canada (NSERC; funding reference number RGPIN-2015-05235) and from an Ontario Early Researcher Award (ER16-12-061). J. Hunt is supported by a Dunlap Fellowship at the Dunlap Institute for Astronomy & Astrophysics, funded through an endowment established by the Dunlap family and the University of Toronto. J. Hunt performed part of this research at KITP with support from the Heising-Simons Foundation and the National Science Foundation (grant No. NSF PHY-1748958). We gratefully acknowledge the support of NVIDIA Corporation with the donation of a Titan Xp GPU used in this research. J.T. Mackereth acknowledges support from the ERC Consolidator Grant funding scheme (project ASTEROCHRONOMETRY, G.A. n. 772293).

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss.org.

SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

This work presents results from the European Space Agency (ESA) space mission Gaia. Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC is provided by national institutions, in particular the institutions participating in the Gaia MultiLateral Agreement (MLA). The Gaia mission website is https://www.cosmos.esa.int/gaia. The Gaia archive website is https://archives.esac.esa.int/gaia.

Appendix A Pattern speed determination for different simulations

To test the determination of the pattern speed using the continuity equation, we apply it to three $N$-body simulations of barred disk galaxies that are set up to resemble the Milky Way. Figure 3 displays the kinematics of stars within 300 pc of the mid-plane for the simulation of Kawata et al. (2017), which clearly shows the same patterns as those observed in the data in Figure 2. The other simulations are Target II and Target IV of Hunt et al. (2013). The dark-matter halo is held fixed in all of these simulations and the bar’s pattern speed is therefore quite stable, because the bar is unable to exchange angular momentum with the halo. We compute the orientation of the major axis of each simulation by fitting the surface-density profile of stars within 6 kpc from the center with a model where the density is constant along concentric, similar ellipses with an axis ratio that is a free parameter and a density profile along the major axis that is fit as a third-order polynomial. We then rotate all simulations such that their major axis is oriented at 25∘ clockwise from the Sun–Galactic center line to match the orientation of the Milky Way’s bar. As we do for the data, we then first determine the surface-density profile $\partial\ln\tilde{\Sigma}/\partial R_{a}$ along the major axis using Equation (3), using the same cuts in angle and similar cuts in radius as for the data (we only use spatial bins along the major axis that are on the same side of the galaxy as the Sun; two of the simulations have slightly longer bars than the Milky Way’s and for those we consider bins out to $R=6\,\mathrm{kpc}$). Comparing the kinematically-determined $\partial\ln\tilde{\Sigma}/\partial R_{a}$ to the actual surface-density profile of the simulation shows excellent agreement. We then use this $\partial\ln\tilde{\Sigma}/\partial R_{a}$ to determine $q$ with the equation for $q$ and then the pattern speed using Equation (6), selecting spatial bins located more than $5^{\circ}$ away from the major axis, but additionally only using bins on the same side of the galaxy as the Sun and within $55^{\circ}$ of the major axis to approximately match the populated bins in the data. We determine the following pattern speeds: (a) $\Omega_{b}=24.5\,\mathrm{km\,s}^{-1}\,\mathrm{kpc}^{-1}$ for the Kawata et al. (2017) simulation with true $\Omega_{B}=23.5\,\mathrm{km\,s}^{-1}\,\mathrm{kpc}^{-1}$, (b) $\Omega_{b}=32.2\,\mathrm{km\,s}^{-1}\,\mathrm{kpc}^{-1}$ for the Hunt et al. (2013) Target II simulation with true $\Omega_{B}=33\,\mathrm{km\,s}^{-1}\,\mathrm{kpc}^{-1}$, and (c) $\Omega_{b}=29.0\,\mathrm{km\,s}^{-1}\,\mathrm{kpc}^{-1}$ for the Hunt et al. (2013) Target IV simulation with true $\Omega_{B}=27.5\,\mathrm{km\,s}^{-1}\,\mathrm{kpc}^{-1}$. The axis ratios that we determine kinematically agree to within a few percent with those determined near the end of the bar from a density fit to the simulated bars with varying axis ratio. It is clear that the assumption of a constant axis ratio does not limit our inference, because we are able to recover the simulations’ pattern speed to high precision even though the simulations do not have constant axis ratios (Hunt & Bovy, 2018).