THE GENESIS OF THE MILKY WAY'S THICK DISK VIA STELLAR MIGRATION

In the years since Gilmore & Reid (1983) first proposed a two-component structure to the Milky Way (MW) disk, a large body of observational work has provided supporting evidence for contrasting thin and thick disk attributes. Structurally, the thin disk scale height is shorter than the thick disk scale height (for reviews see Majewski 1993; Buser et al. 1999; Norris 1999, and references therein), and the thick disk may have a longer scale length than the thin disk (Robin et al. 1996; Ojha 2001; Chen et al. 2001; Larsen & Humphreys 2003). Kinematically, thick disk stars have larger velocity dispersions and lag the net rotation of the disk (Nissen 1995; Chiba & Beers 2000; Gilmore et al. 2002; Soubiran et al. 2003; Parker et al. 2004; Wyse et al. 2006). Additionally, thick disk stars are older and metal-poor relative to their thin disk counterparts (e.g., Majewski 1993; Chiba & Beers 2000; Bochanski et al. 2007) and at a given iron abundance, thick disk stars are α-enhanced (Fuhrmann 1998; Prochaska et al. 2000; Tautvaišienė et al. 2001; Bensby et al. 2003, 2005; Feltzing et al. 2003; Mishenina et al. 2004; Brewer & Carney 2004). Moreover, the presence of both a thin and a thick disk is not unique to the MW but a ubiquitous feature of late-type galaxies (Burstein 1979; van der Kruit & Searle 1981; Abe et al. 1999; Neeser et al. 2002; Yoachim & Dalcanton 2005, 2006; Yoachim 2007).

Recently, several Sloan Digital Sky Survey (SDSS)-based studies have provided further strong observational constraints on the structural, kinematic, and chemical properties of stars in the solar cylinder. Jurić et al. (2008, hereafter J08) used a photometric parallax method on SDSS data to estimate distances to ∼48 million stars and studied their spatial distribution. Because SDSS provides accurate photometry, which enables reasonably robust distances (10%–15%; Sesar et al. 2008), as well as faint magnitude limits (r < 22) and a large sky coverage (6500 deg²), J08 were able to strongly constrain the parameters of a model for the global spatial distribution of stars in the MW. The J08 model is qualitatively similar to previous work (e.g., Bahcall & Soneira 1980) which identifies a clear change of slope in the counts of disk stars as a function of distance from the Galactic plane; this change in slope is usually interpreted as the transition from the thin to the thick disk (Gilmore & Reid 1983; Siegel et al. 2002).

Ivezić et al. (2008, hereafter I08) further extended this global analysis of SDSS data by developing a photometric metallicity estimator and by utilizing a large proper motion catalog based on SDSS and Palomar Observatory Sky Survey data (Munn et al. 2004). I08 studied the dependence of the metallicity, [Fe/H], and rotational velocity, V_ϕ, of disk stars on the distance from the Galactic plane and detected gradients of both quantities over the distance ranging from several hundred parsecs to several kiloparsecs. Such gradients would be expected in a thin/thick disk decomposition where the thick disk stars are a separate population defined by a bulk rotational velocity lag and a lower metallicity compared to those of the thin disk. However, such a model would also predict a strong correlation between the metallicity and the velocity lag, which is excluded (∼7σ level) by the I08 analysis⁷ (see Figure 17, I08). In this work, we turn to a more sophisticated Galactic description—an N-body model—to characterize stars within the SDSS volume.

Over the past few decades, N-body simulations have been used to provide supporting evidence for three distinct theories of thick disk formation: violent relaxation (Jones & Wyse 1983), substructure disruption (Statler 1988), and heating by satellites (Quinn et al. 1993). Several works have recently revisited these ideas. Brook et al. (2004) formed a thick disk in situ at high redshift during gas-rich mergers, where star formation is triggered by the rapid accretion of gas; this result is consistent with the thick disk forming through violent relaxation of the Galactic potential. Bournaud et al. (2009) further studied the effect of violent relaxation and found that internal disk instabilities rather than mergers generate a thick disk with proper morphology. Abadi et al. (2003) showed that by tidally stripping/accreting satellites, the majority of the oldest stars in the thick disk could have formed externally. Kazantzidis et al. (2008), Villalobos & Helmi (2008), and Villalobos et al. (2010) investigated substructure disruption by using a cosmologically derived satellite accretion history to perturb an MW-like disk; subhalo-disk encounters increased the scale height of this disk at all radii effectively forming a thick disk.

In this work, we study a different mechanism of thick disk formation: radial migration (Schönrich & Binney 2009a, 2009b; Sales et al. 2009; Loebman et al. 2008; Caruana 2009). Radial migration due to transient spirals was first described by Sellwood & Binney (2002). In this model, energy and angular momentum changes occur from interactions with transient spiral arms, which move stars at the corotation resonance inward or outward in radius while preserving their nearly circular orbits. In recent years other authors have studied radial migration in simulations, due to both corotation scattering (Lépine et al. 2003) and bar–spiral resonant interactions (Lee et al. 2011). Using a semi-analytic approach, Schönrich & Binney (2009a, 2009b) demonstrated that the MW's thick disk could result from radial migration. Roškar et al. (2008a, 2008b, hereafter R08ab) studied radial migration in N-body + smooth particle hydrodynamic (SPH) simulations of disk formation and showed that migrations are possible on short timescales. They explored the implications of radial mixing for stellar populations for a variety of stellar systems, including the solar neighborhood. Here we extend their work by highlighting the vertical evolution that occurs as a result of migration.

We note that in this paper, we are not testing the validity of the other models of formation. However, recently, Sales et al. (2009) proposed using the eccentricity of orbits of stars in the thick disk to constrain the thick disk's formation mechanism; they presented eccentricity distributions that result from four N-body simulations: Abadi et al. (2003), Villalobos & Helmi (2008), R08b, and Brook et al. (2004). They found that the distributions that result from heating, radial migration, and mergers all had a strong peak at low eccentricity ( ∼ 0.2–0.3), while the distribution that results from accretion is centered at higher orbital eccentricities (〈〉 ∼ 0.5). Building on this, Wilson et al. (2011) studied the eccentricity of orbits of stars in the thick disk observed in the Radial Velocity Experiment (RAVE; Steinmetz et al. 2006) and found these results to be inconsistent with expectations for the pure accretion simulation (but see also Di Matteo et al. 2011). Ruchti et al. (2010) also utilized α measurements from RAVE to conclude that the α-enhancement of the metal-poor thick disk implies that direct accretion of stars from dwarf galaxies did not play a major role in the formation of the thick disk. Using SDSS DR7, Dierickx et al. (2010) showed that the eccentricity of stellar orbits in the thick disk implies that the thick disk is unlikely to be fully populated by radially migrated stars. We note that we cannot exclude that some fraction of the thick disk is a fossil of a past more violent history, nor can this scenario explain thick disks in all galaxies. However, in what follows, we show that a large fraction of the stars in the thick disk could have formed within the disk at smaller radii and arrived at their present location via radial migration.

The outline of this paper is as follows: in Section 2, we present our adopted simulation, which contains substantial migration; we also present a contrasting control case, which has little migration relative to the disk scale length. When we compare these two simulations, we show that migration can build a thick disk as first defined by Gilmore & Reid (1983): a component with a scale height larger than that of the thin disk. In Section 3, we qualitatively compare the adopted simulation with SDSS observations to show that they match each other sufficiently well to pursue further comparison. In Section 4, we present a detailed comparison between the simulation and the local SDSS volume focusing on the reason for the absence of a strong correlation between V_ϕ and [Fe/H]. In Section 5, we use the simulation as a proxy for the MW to show that classifying stars as members of the thin or thick disk by either velocity or metallicity leads to an apparent separation in the other property as observed. In Section 6, we compare our results to recent theoretical work that used semi-analytics to investigate how the solar neighborhood could have been shaped by radial migration and chemical evolution effects. In Section 7, we explore the correlation between [α/Fe] and age to highlight the diagnostic power of [α/Fe] as a stand-in for age. In Section 8, we summarize our results and conclusions, and in the Appendix, we reconsider recent observational claims concerning the lack of correlation between V_ϕ and [Fe/H].

We analyze the results of an N-body + SPH simulation designed to mimic the quiescent formation and evolution of an MW–mass galaxy following the last major merger. The system is initialized as in Kaufmann et al. (2007) and R08ab and consists of a rotating, pressure-supported gas halo embedded in a Navarro–Frenk–White (Navarro et al. 1997) dark matter halo. This simulation was evolved for 10 Gyr using the parallel N-body+SPH code, GASOLINE (Wadsley et al. 2004). As the simulation proceeds, the gas cools and collapses to the center of the halo, forming a thin disk from the inside out. Gas is continually infalling from the hot halo onto the disk for the duration of the simulation. Star formation and stellar feedback are modeled with subgrid recipes as described in Stinson et al. (2006). Importantly, the stellar feedback prescriptions include Type II supernova (SN II), SN Ia, and asymptotic giant branch metal production, as well as injection of supernova energy which impacts the thermodynamic properties of the disk interstellar medium (ISM). Metal diffusion is calculated from a subgrid model of eddy turbulence based on the local smoothing length and velocity gradients (Smagorinsky 1963; Wadsley et al. 2008). The simulation we utilize is nearly identical to R08ab (see R08ab for further details), but with the addition of metal diffusion between gas particles (Shen et al. 2010). We experimented with a few values of the diffusion coefficient and settled on a diffusion parameter of D = 0.05, which we find to yield reasonable metallicity distributions in the solar neighborhood.

No a priori assumptions about the disk's structure are made—its growth and the subsequent evolution of its stellar populations are completely spontaneous and governed only by hydrodynamics/stellar feedback and gravity. Although we do not account for the full cosmological context, merging in the Λ cold dark matter (ΛCDM) paradigm is a higher order effect at the epochs in question (Brook et al. 2005). Thus, our model galaxy lacks some structural components such as a stellar halo, which in ΛCDM is built up primarily during the merging process (e.g., Bullock & Johnston 2005; Zolotov et al. 2009). Our focus here, however, is disk evolution; by simplifying our assumptions, we are able to use much higher resolution and more easily study the impact of key dynamical effects on observational properties of stellar populations within the disk.

The gas disk in our model grows spontaneously from quiescent accretion of condensing halo gas. Once the baryons are arranged in the disk, they reach densities high enough for star formation, which proceeds according to our subgrid prescriptions (see Stinson et al. 2006). During most of the evolution, the stellar disk is marginally stable, with Toomre Q values ranging from 2 to 3 in the main part of the disk (Q is higher for R < 2 kpc where velocity dispersion is high and in the outskirts where the density drops quickly). This marginal (in)stability allows for the continuous triggering of transient spirals.

To quantify the spirals, we expand the density distribution in a Fourier series and extract pattern speeds and amplitudes of the dominant patterns at different radii. We find that multiple spirals are always present in the disk, each with a different pattern speed dominating a different part of the disk. The m = 2 spirals are always the dominant component, although spirals of higher multiplicity are also present. Figure 1 shows the time evolution of the m = 2 Fourier amplitude in four different radial bins. The amplitude of m = 2 spirals varies from 0.1 to 0.4, though there are not many instances where the Fourier amplitude exceeds 0.2. Further details regarding the nature of spirals in our simulations are beyond the scope of this paper, but will be presented in a forthcoming publication.

Zoom In Zoom Out Reset image size

Based on their simulations, R08ab presented the implications of stellar radial migration resulting from the interactions of stars with transient spiral arms (Sellwood & Binney 2002) on the observable properties of disk stellar populations. Radial migration efficiently mixes stars throughout the disk into the solar neighborhood, resulting in a flattened age–metallicity relation, a broader metallicity distribution function (R08ab), and significantly altered kinematics (Sales et al. 2009). Figures 2 and 3 illustrate the basic premise of this paper—stars that presently occupy the solar neighborhood were not all born there, but have come from a wide range of radii. As the stars migrate radially outward, their vertical motions carry them to much larger heights; because the disk mid-plane potential is shallower at the new location, the amplitude of stars' vertical oscillations increases. This effect is clearly seen in Figure 3—as we consider distributions of formation radius and age at increasing heights above the plane, the stars are predominantly old and have originated in the inner disk. Previous studies have shown (Schönrich & Binney 2009a, 2009b; Sales et al. 2009; Loebman et al. 2008; Caruana 2009) that the vertical evolution resulting from radial migration can influence the characterization of the thick disk.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this paper, we use an improved version of the simulation utilized by R08ab because the simulation has roughly the same mass and size as the MW. As such we often refer to it, for brevity, as the MW simulation, to distinguish it from the second simulation we present below. We stress, however, that the MW simulation is not meant to represent the MW in detail—the simulated galaxy is slightly more massive (V_circ ∼ 250  km s⁻¹) and only has a weak oval whereas the MW is barred.

In order to further illustrate the importance of radial migration, we repeated much of our analysis on a control case. The control simulation is a system with the same initial conditions as the MW simulation except for having a higher angular momentum content with a dimensionless spin parameter λ = 0.1 (Bullock et al. 2001). This results in a more extended disk (final disk scale length =3.71 kpc, versus 2.54 kpc for the MW simulation), possibly similar to a low surface brightness galaxy; we therefore refer to this simulation as the low surface brightness (LSB) simulation. Due to its lower surface density, the disk forms weaker spirals and as a result the stellar populations at all radii are less affected by radial mixing. When we compare migration as a function of scale lengths, we find that there is significantly less migration in the LSB simulation than in the MW simulation. This can be seen in Figure 2: the dotted lines show the distribution of formation radii for stars in the LSB simulation (versus the MW simulation shown in solid lines); the LSB distribution shows a much smaller offset from the stars' present radius, highlighted in gray, than the MW simulation's. Thus, the weaker spirals of the LSB simulation cause significantly less radial migration than in the MW simulation.

We explicitly checked that the weaker spiral structure in the LSB case is not a result of numerical suppression due to limited spatial resolution. Calculating the ratio of the softening length to λ_crit = 4π²GΣ_s/κ², where G is the gravitational constant, Σ_s is the stellar surface density and κ is the epicyclic frequency, we find that for the LSB case it is ≪1 at all radii of interest. For the fiducial simulation, the ratio only approaches unity beyond 15 kpc. Because λ_crit is indicative of the spatial scale of expected spirals, this test argues that we are adequately resolving the structure in our simulated disks.

The distribution of stellar mass away from the mid-plane is strongly affected by radial migration; this can be seen in Figure 4, which contrasts the MW simulation against the LSB case. Here the normalized mass density distribution within four analogous cylindrical volumes drawn from a variety of radii is presented. At larger radii, the steepness of the profiles is quite different; the LSB simulation has a constant slope while the MW simulation shows a transition from a steep to a shallow density distribution. Thus the MW simulation cannot be characterized by a single-exponential or sech² component in the vertical direction, as we show explicitly in the following section and in Figure 6. It is this double-component nature which first led to the identification of the thick disk (Gilmore & Reid 1983); we have thus shown that internal evolution may also give rise to this important feature.

Zoom In Zoom Out Reset image size

The two models not only yield very different vertical density profiles, but also the resulting stellar population properties are drastically different in the model with little migration. Comparing Figure 5 with Figure 3, it is clear that unlike in the MW simulation, where virtually all stars found at significant heights above the plane come from the inner disk, most stars in the LSB disk simply stay put over time. The properties of stars away from the mid-plane in the MW simulation discussed in the forthcoming sections are therefore a unique consequence of the radial migration process.

Zoom In Zoom Out Reset image size

In the following section, we compare SDSS observations with the MW simulation to demonstrate its usefulness as a model for understanding the MW thick disk. Here we study the stellar mass distribution, rotational velocity, and metallicity as functions of distance from the Galactic plane, |z|, and galactocentric cylindrical radius, R. We draw qualitative comparisons between the data sets by examining their mass-weighted metallicity and kinematic distributions in this R–|z| space.

The observed MW disk is best fit by a two-component model that is exponential both in the R and z directions (see Table 10, bias-corrected results, J08). The top panel of Figure 7 shows the mass-weighted density distribution of the entire MW simulation at its final time step. This distribution is in qualitative agreement in both the R and z directions with J08 for up to ∼2.5 kpc above the disk's plane and ∼15 kpc from the galactic center.

Figure 6 shows our best fits to the vertical density profiles for the MW simulation for radial bins R = 2–4 kpc, R = 4–6 kpc, R = 7–9  kpc, and R = 9–12 kpc. As in J08, the innermost radial bin is best fit by a double-exponential function; however, all other bins are better fit by the sum of two double sech² profiles, which is in agreement with the theoretical work by Spitzer (1942) and observational results for other galaxies by Yoachim & Dalcanton (2006). Despite the small discrepancy between sech² versus exponential fits, the simulation is in good qualitative agreement with the SDSS-based analysis of the MW (J08).

Zoom In Zoom Out Reset image size

In addition, Figure 6 shows the vertical density profiles as divided into two populations: stars that have migrated greatly (|R − R_form| > 2 kpc, blue dots) and stars that have not migrated much (|R − R_form| ⩽ 2 kpc, red dots). Note that we have plotted error bars reflecting Poisson noise for each bin. We note that the second-component fit for radial bins R = 7–9 kpc and R = 9–12 kpc is entirely dominated by star particles that have migrated. We expect that stars that moved from elsewhere combined with stars that were born locally should not naturally conspire to produce a single sech² profile.

For the solar cylinder, R = 7–9 kpc, we find that the model distribution of stars as a function of |z| resembles a double sech² profile, with the "transition" height of |z|  ∼ 0.75 kpc, comparable to ∼1 kpc found by J08 for the Galactic disk. We have found the best-fit scale heights to be 381 pc and 913 pc, in qualitative agreement with the best-fit scale heights of 270 pc and 1200 pc for the SDSS data (J08). The scale height ratio suggested by the simulation is slightly low—∼2.4, instead of ∼3 from the data. Moreover, the simulation's thick disk to thin disk normalization is slightly discrepant—∼0.14, rather than ∼0.12 suggested by the SDSS data. Overall, we find that the thick disk in the MW simulation, formed through the process of radial migration, is qualitatively very similar to the observed Galactic thick disk. In contrast, as Figure 4 shows, this conclusion does not apply to the LSB simulation.

The median metallicity of the MW disk exhibits a clear vertical gradient (see bottom two panels of Figure 9 in I08). Notably, the spatial variation of the median metallicity does not follow the distribution of the stellar number density (I08). The middle panel of Figure 7 shows that the MW simulation reproduces a qualitatively similar metallicity distribution; we note that a constant additive offset of 0.1 dex has been applied to [Fe/H] throughout the simulation so that the median metallicity in the plane of the disk in the solar cylinder matches observations (see the metallicity distribution function in R08b). As expected, at low Galactic latitudes and small radii, the volume is dominated by high (near solar) metallicities. At higher latitudes the volume is increasingly metal-poor.

Zoom In Zoom Out Reset image size

As with metallicity, previous studies have found a gradient in the median V_ϕ with respect to z in the MW (see Figures 5, 8, and 9 in Bond et al. 2010 (hereafter B10), and references therein). These authors concluded that V_ϕ is also well characterized by a non-Gaussian distribution (see Binney 2010 for the implied distribution function). The bottom panel of Figure 7 shows that the MW simulation reproduces qualitatively similar V_ϕ properties, including a |z| gradient.

Thus on a gross scale the MW simulation qualitatively matches patterns observed in mass density, metallicity, and rotational velocity in the MW disk. We now look at the solar cylinder in greater detail.

4.1. Vertical Gradients

Here we study in detail the distributions of age, stellar density, rotational velocity, and metallicity as functions of |z| for the solar cylinder in the MW simulation. To be consistent with the analysis of high Galactic latitude SDSS data by J08 and I08, we select model particles from an annulus with 7 kpc ⩽R ⩽ 9 kpc. This radial cut spans 2.8–3.6 scale lengths from the center of the model (scale length 2.5 kpc) and covers the Sun's location of ∼3 disk scale lengths from the MW center (scale length ∼2.6 kpc; see bias-corrected value; Table 10, J08).

The behavior of the MW simulation in this particular volume is illustrated in Figure 8. The top panel shows the age distribution as a function of |z|. As expected, near the mid-plane, the population is dominated by young stars, and the mean age monotonically increases with increasing distance from the mid-plane. Only very old stars are found at large |z|: by |z| ∼ 0.5 kpc the mean age is already ∼5 Gyr. Recall that Figure 3 shows that with increasing distance from the mid-plane, the stellar population becomes dominated by older stars that formed closer to the center of the disk. These two figures taken together give a coherent picture of the net dynamic effect on the system: on average stars move both radially outward and away from the mid-plane over time.

Zoom In Zoom Out Reset image size

Radial migration is able to change the extent of a star's vertical oscillation, as well as its rotational velocity (middle panel, Figure 8). In the MW simulation, we find that the rotational velocity slows by 18  km s⁻¹ for every kpc (roughly comparable to the gradient of −30  km s⁻¹ kpc⁻¹ found by I08). Note that stars which significantly lag the rotation of the disk in the mid-plane have traditionally been regarded as members of the thick disk.

Another observed trend that is physically well motivated by the net outward and upward movement of stars over time is the decline of metallicity with increased height (bottom panel, Figure 8). In the MW simulation, the metallicity distribution changes with |z| with a best-fit gradient of ∼0.18  dex  kpc⁻¹, again in qualitative agreement with the measured value for the MW of ∼0.30  dex  kpc⁻¹ (I08).

To understand how this trend arises, recall that the top panel of Figure 8 shows that the stellar population away from the plane is dominated by old stars. As Figure 3 shows, the oldest stars at large |z| mostly formed in the inner 2–4 kpc. Moreover, at early times, the radial metallicity gradient was steep (see Figure 2 in R08b); stars that are now in the thick disk in the solar cylinder in the MW simulation were once at the outer edge of the forming thin disk and hence formed at a low metallicity. As a result, stars that formed at a radius of 2–4 kpc early on in the galactic history are necessarily metal-poor. These stars have been subsequently moved out and up over time. At later times, the metallicity gradient flattened out as the disk grew. Again, looking at Figure 2 (top right panel), one can see a significant fraction of young stars (i.e., less than 4 Gyr old) formed at R  ∼ 6 kpc where there was a relatively more metal-rich ISM. These stars had less time to migrate and as a result remain closer to the mid-plane of the disk. This complex co-dependence of radial migration, birth location, and metallicity gradient evolution then gives rise to a vertical metallicity gradient in the solar neighborhood.

4.2. Thin/Thick Disk Transition

We now turn our attention to a particular region within the solar cylinder of the MW simulation: |z| = 0.5–1.0 kpc, R = 7–9 kpc. This region in the simulation is analogous to the thin–thick disk "transition zone" considered by I08: within this volume a roughly equal number of thin and thick disk stars is expected (see the bottom left panel, Figure 6). The transition region within the MW simulation occurs in approximately the same place as in the data (0.75 kpc versus ∼1 kpc) which allows us to draw a direct comparison to the slice analyzed in I08. We show here that all the trends observed in the SDSS data can be explained by a continuous distribution rather than two distinct populations.

In Figure 9, we show the distribution of observable properties as a function of age within this slice of the MW simulation. The top left panel of Figure 9 shows the distribution of ages in this thin volume slice: the region is well populated (∼30,000 star particles) and predominantly old (36% older than 7 Gyr, 63% older than 5 Gyr). The remaining panels in Figure 9 illustrate probability densities of formation radius, metallicity, and rotational velocity versus age. In all cases, the distributions do not suggest distinct populations. For the given |z| slice, stars older than about 4 Gyr are both numerous and formed significantly interior to their present location (bottom left panel). These stars also show a rotational velocity lag of ∼20 km s⁻¹, as shown in the top right panel. The oldest stars, those with ages ≳ 8 Gyr, have significantly lower metallicities than younger stars (bottom right panel). As noted previously, the oldest stars formed at a range of interior radii, when the ISM metallicity was low and had a steep radial gradient. These effects drive the behavior of various observables, as discussed below.

Zoom In Zoom Out Reset image size

4.3. Rotational Velocity versus Metallicity

We have shown in Section 4.1 that the MW simulation yields vertical trends in metallicity and rotational velocity similar to those found in SDSS observations presented in I08 and B10. SDSS data also revealed a surprising lack of strong correlation between V_ϕ and [Fe/H], contrary to the expectations based on a traditional two-disk model (see the Appendix for a discussion of an observational claim that a weak correlation does exist for −1 <[Fe/H] < − 0.5). In the top left panel of Figure 10, we show V_ϕ versus [Fe/H] for the vertical slice considered above (|z| = 0.5–1.0 kpc) corresponding to the transition between the thick and thin disk in the simulation. The MW simulation also yields a lack of strong correlation between these quantities: although both rotational velocity and metallicity show robust vertical gradients, when stars are selected from a thin z slice, the correlation between velocity and metallicity is much weaker.

Zoom In Zoom Out Reset image size

We can understand why there is no strong correlation between V_ϕ and [Fe/H] if we consider the thin slice as occupied by an age ensemble that was brought together by radial migration. In this light, we can decompose the sample by age to see how the behavior of each population is modified with time (see the top left panel of Figure 9 for the histogram of the age subsamples considered).

The top right panel of Figure 10 shows that a strong correlation between rotational velocity and metallicity exists for young stars; the gradient is −71  km s⁻¹ dex⁻¹ (more metal-rich young stars rotate more slowly). The young stars constitute a small fraction of the overall population in the selected volume (10%) and, as the top right and bottom right panels of Figure 9 illustrate, young stars have correlated metallicity and formation radius. Thus if a young star has a metallicity different than the surrounding ISM, it must originate from somewhere other than the solar cylinder and be at or near perigalacticon (if R_g ⩾ 9 kpc, where R_g is the radius of the stellar guiding center) or at or near apogalacticon (if R_g ⩽ 7 kpc); consequently, it will either lead or lag the local standard of rest (LSR).

For increasingly older age bins, the gradient diminishes and ultimately fades away. The bottom left and right panels of Figure 10 show V_ϕ versus [Fe/H] for intermediate (4 Gyr ⩽ Age ⩽ 6 Gyr) and old (8 Gyr ⩽ Age ⩽ 10 Gyr) populations, respectively. The intermediate age stars constitute 24% of the stellar population in the thin slice and retain a slight gradient =−19  km s⁻¹ dex⁻¹. The oldest stars make up 21% of the stellar population and show a nearly flat/slightly positive slope of 8  km s⁻¹ dex⁻¹.

The peak in the distribution of intermediate and oldest age star particles is offset to progressively lower metallicity and rotational velocity. The metallicity of intermediate age stars is in the range −0.5 < [Fe/H] <0.5 while that of old stars is −1.0 < [Fe/H] <0.0. We can think of these three age bins as dominating different portions of the overall V_ϕ versus [Fe/H] space, with the peak of intermediate and old stars being perpendicular to the correlation for young stars.

The peak mean value of the oldest stars is [Fe/H]  =   − 0.5, V_ϕ = −210 km s⁻¹; this is significantly lower than the peak mean value of the full age sample [Fe/H]  =   − 0.1, V_ϕ = −225 km s⁻¹. The full age sample's peak mean value matches the intermediate age values, which is in turn lower than the peak mean value of the youngest stars [Fe/H]  =   − 0.05, V_ϕ = −245 km s⁻¹. When the entire population is considered as a whole, the correlation becomes much weaker.

Why does the V_ϕ–[Fe/H] gradient diminish with increasing age? To understand this, we return to Figure 9 and recall that most of the stars within this thin cut are intermediate to old age, and they did not form within 7 kpc ⩽ R ⩽ 9 kpc. A range of formation radii corresponds to a range of formation environments; this maps to a range of metallicities for stars of a given old age bin. At the same time, depending on a star's formation location and subsequent migration and scattering off the disk substructure, it can end up with a range of possible rotational velocities (see the bottom left panel, Figure 9). Hence, even a single old age bin samples a wide range of formation environments and metallicities, in addition to a wide range of rotational velocities resulting from unique dynamical histories, which are not directly correlated.

We stress the significant difference between these results and the traditional two-disk interpretation: in the latter case a correlation between V_ϕ and [Fe/H] is expected in an older, metal-poor, thick disk population. In contrast, here we find that a trend is present in the young stars but absent in the older stars. The trend in the young stars can be easily understood as arising from epicyclic motions of stars with their birth radii imprinted into their metallicities. As migration moves the guiding centers of stars, this metallicity encoding is erased, and the correlation between V_ϕ and [Fe/H] disappears.

4.4. Geneva–Copenhagen Survey

As we demonstrated above, we can understand why there is no strong correlation between V_ϕ and [Fe/H] when we decompose the sample by age. Older populations have had more time to radially mix; thus in the older age bins the expected trend is much weaker. That is, the evolution of the trends shown in Figure 10 is a unique signature of radial mixing taking place in the disk.

We can verify whether this is the case in the solar neighborhood by utilizing observational data taken from the Geneva–Copenhagen Survey (GCS; Holmberg et al. 2009). The GCS samples a wide range of stellar ages and is reasonably well populated with old stars to make qualitative comparisons to the MW simulation in the plane of the disk. Because stellar age estimates are fairly uncertain, we compare only the lowest and the highest thirds of the GCS distribution. While what we have shown in Figure 10 is the relationship at 0.5–1 kpc above the plane, these trends evolve in the same way as a function of age in the mid-plane.

We have selected the non-binary stars from the GCS and repeated the analysis we presented above. Figure 11 shows the results: the top panel shows that there is no net trend when the sample is considered as a whole. Splitting the sample into two broad age bins (0–3 Gyr and 6–15 Gyr) yields a trend only in the young stars—the correlation is absent for the older stars. The evolution of the observed V_ϕ versus [Fe/H] trend in the GCS sample therefore matches our expectations based on the simulation and suggests that the stellar populations in the solar neighborhood have been influenced by radial mixing.

Zoom In Zoom Out Reset image size

We now turn our attention to assigning thin or thick disk membership in a manner analogous to observational studies, i.e., either based on kinematics or metallicity. We find that classifying stars as members of the thin or thick disk by either velocity or metallicity leads to an apparent separation in the other property, as observed.

5.1. Membership Based on Kinematic Criteria

One of the observational differences between kinematically selected thin and thick disk stars is that the latter have higher abundance of α elements at a given [Fe/H] (Bensby et al. 2005; Feltzing 2006, and references therein). This difference is often interpreted as evidence for different formation histories. Our simulations include a calculation of the oxygen abundance following the prescription by Raiteri et al. (1999); given that SNe II mostly yield oxygen (Hoffman et al. 1999), we use oxygen as a proxy for all α elements. Because metal yields are not precisely known, the absolute values for [α/Fe] within the MW simulation do not match the MW. However, we believe the relative differences in [α/Fe] to be qualitatively accurate; thus the MW simulation gives a good qualitative perspective on distributions in [α/Fe] space within the solar cylinder.

The top left panel of Figure 12 shows the overall dependence of [α/Fe] on [Fe/H] for R = 7–9 kpc, |z| = 0.0–0.3 kpc within the MW simulation (i.e., the "solar neighborhood").⁸ A similar behavior is seen for the stars at |z| ∼ 1 kpc in agreement with Bensby et al. (2005). When the same sample is separated by age, distinct portions of the parameter space are covered. In particular, at low [Fe/H], old stars show an enhancement of [α/Fe] relative to young stars.

Zoom In Zoom Out Reset image size

Locally, observed thick disk stars are selected kinematically, rather than by age (Prochaska et al. 2000; Reddy et al. 2003, 2006; Allende Prieto et al. 2004). We reproduce the qualitative behavior of observations by following similar steps with the MW simulation. The top right panel in Figure 12 shows a Toomre diagram for simulated stars, with selection cuts motivated by Nissen & Schuster (2009). The bottom left panel shows that these kinematically selected thin and thick disk stars from the solar neighborhood show similar bifurcation of [α/Fe] versus [Fe/H] behavior as young and old subsamples shown in top left panel. In the bottom right panel, stars from just the shaded [Fe/H] cut in the bottom left plot are examined.

Stars falling within the kinematically selected thick disk region have a higher fraction of old stars relative to the overall population. These old stars are α-enhanced as they formed in the interior of the disk and radially migrated to their present location. Thus kinematically dividing the stars locally yields a sample that is older and α-enhanced. Although the simulation and the data are not in detailed quantitative agreement, these qualitative results imply that the differences in [α/Fe] versus [Fe/H] for kinematically selected thin and thick disk stars may be another consequence of mixing effects that result from radial migration.

5.2. Membership Based on [α/Fe] Criteria

Although not as commonly done, it is equally plausible to take the converse approach; we can assign membership to the thin and thick disk based on an [α/Fe] cut and then study the kinematic and [Fe/H] distributions that result (Navarro et al. 2011; Lee et al. 2010). Here we follow the technique outlined by Lee et al. (2010), who chemically divided the Galactic disks using SDSS SEGUE-1 data; we select star particles within R = 7–11 kpc, |z| = 0.3–2.0 kpc and split them so that stars with [α/Fe] ⩾−0.1 are assigned thick disk membership.

Assigning membership based on [α/Fe] effectively divides the disk into two populations: old stars and young-intermediate age stars, as can be seen clearly in the top left panel of Figure 13 (see also Schönrich 2010; Haywood 2008). Why is this the case? Star particles are born α-enhanced if they form in a region with a high local star formation rate (SFR) and little SN Ia pollution. Since most old stars originated near the center of the disk, and that region is where the SFR was high, these stars are naturally α-enhanced.

Zoom In Zoom Out Reset image size

We compare our results to three trends discussed in Lee et al. (2010): the radial metallicity gradient, distribution in V_ϕ and distribution in [Fe/H], which we show clockwise from the top right in Figure 13. The top right panel illustrates the best-fit radial metallicity gradients as derived from the mass-weighted mean value of [Fe/H]. Cutting by [α/Fe] results in no trend (slope ∼ 0 kpc dex⁻¹) in the thick disk and a negative trend in the thin disk (slope ∼ −0.2 kpc dex⁻¹); this is similar to the thick disk slope =0 kpc dex⁻¹ and thin disk slope =−0.3 kpc dex⁻¹ observed by Lee et al. (2010). The bottom right panel illustrates the cumulative V_ϕ distributions that result for the MW simulation when stars are separated by [α/Fe]: the thin and thick disk trends are offset by ∼22 km s⁻¹, which is qualitatively similar to, if quantitatively smaller than, the ∼50 km s⁻¹ offset found by Lee et al. (2010). Finally, the bottom left panel shows the cumulative [Fe/H] distributions; we find the thin and thick disk trends offset by ∼0.35 dex, not dissimilar from the observed ∼0.4 dex offset (Lee et al. 2010).

Therefore, by adopting a kinematic or [α/Fe] selection criteria used by observers, we are able to reproduce an apparent separation in the other property, despite the fact that there are no distinct populations in the model galaxy.

Recently, Schönrich & Binney (2009a, 2009b) investigated how radial migration and chemical evolution shape the solar neighborhood, by incorporating for the first time a prescription for radial migration in a semi-analytic model of Galactic chemical evolution. Their model represented a disk in which star formation commenced at all radii simultaneously (i.e.,without inside-out growth), with radially varying star formation rates set to yield a disk with an appropriate scale length. The guiding center radii of the stars in their model changed (i.e., stars migrated radially) according to a parameterized probabilistic prescription whose normalization was left as a free parameter in the model. The vertical structure of their disk was determined based on the assumption that locally coeval stars comprise an isothermal population with a velocity dispersion given by observational constraints. Because stars migrating from the inner disk retain their velocity dispersions but encounter a lower restoring potential in the outer disk, they populate the disk away from the plane. Hence, Schönrich & Binney (2009b) showed that a thickened component is a natural consequence of radial migration, and by fitting the model they also reproduced many of the canonical features of the thick disk (i.e., enhanced [α/Fe] ratios, older ages, lower metallicities, rotational lag), similar to what we have shown in the previous section.

It is therefore reassuring that our results presented here agree qualitatively with those of Schönrich & Binney (2009b), given that our modeling methods are entirely different, and that our model was not specifically tuned to the Galaxy. However, we also find subtle yet crucial differences. Figure 4 of Schönrich & Binney (2009b) shows that due to radial mixing, a population that shares the same average rotational velocity can show immense variations in its chemical composition. The top left panel of Figure 14 shows that our model also yields stellar populations with the similar average rotational velocities but very different chemistry. However, we note the contours of mean V_ϕ in this plane have a different orientation in our simulation compared to those shown in Figure 4 of Schönrich & Binney (2009b), even though the age structure is very similar; at fixed [Fe/H] age increases monotonically with rising [α/Fe] (middle panel of Figure 14 and Figure 5 of Schönrich & Binney 2009b).

Zoom In Zoom Out Reset image size

A hint of a reason for this discrepancy is provided by scrutinizing the metal-rich end of these figures. In our model, the [α/Fe]-deficient metal-rich population originated in the interior of the disk (see the bottom left panel of Figure 14) and has migrated to the present radius without very much heating. This is apparent from the fact that the mean V_ϕ for this population is only very slightly lagging the LSR (∼−240 km s⁻¹). On the other hand, in the model of Schönrich & Binney (2009b), that same population shows considerable lag from the LSR. This discrepancy implies that there are qualitative differences in the treatment of radial mixing between their prescription and our simulation.

The differences in the velocity structure in this plane may also be a result of the assumption in the Schönrich & Binney model that the entire disk begins forming stars at once with a peak in star formation occurring ∼10 Gyr ago everywhere. In our model, it is impossible to have an old metal-poor star which formed at the solar radius—all of these stars must have migrated to their present position because the solar neighborhood in our simulation does not exist ∼9 Gyr ago.

In an observational study, Haywood (2008) argued that the existing solar neighborhood samples (e.g., the GCS) show signatures of radial mixing as proposed by Sellwood & Binney (2002). In our sample shown in Figure 14, the metal-poor, low [α/Fe] stars have high velocities and young ages and are at or near perigalactic passage from the outer disk into the solar neighborhood sample. This is consistent with the arguments put forth by Haywood (2008) as observational evidence of radial mixing. However, the large velocities of this tail in the distribution signify that these stars are on fairly eccentric orbits and therefore do not necessarily represent a population that migrated via corotation resonance.

We have demonstrated that we can understand the observed V_ϕ, [Fe/H], and [α/Fe] distributions as a consequence of radial migration effects; the significance of these effects became clear when we decomposed the MW simulation within the solar cylinder by age. Unfortunately, age is not easily accessible observationally. However, α is a directly measurable quantity and, as we have shown, for the oldest stars, α-enhancement is a reasonable proxy for age (see the top left panel of Figure 13, as well as Schönrich 2010; Haywood 2008). Thus we now reassess the relationship between V_ϕ and [Fe/H] by decomposing Figure 10 in cuts of [α/Fe] so that we present testable model results for upcoming observational surveys.

Figure 15 shows V_ϕ versus [Fe/H] for the "transition zone" considered in Section 4 (|z| = 0.5–1.0 kpc, R = 7–9 kpc), as split by two broad bins in [α/Fe]. Below each V_ϕ versus [Fe/H] plot is a histogram of the ages represented in the given bin. In the right panel, the high [α/Fe] sample contains exclusively old stars (99% older than 7 Gyr); the corresponding V_ϕ versus [Fe/H] figure shows a weak trend between the quantities (best-fit slope ∼ −14 km s⁻¹ dex⁻¹ for −0.8 < [Fe/H]  <  0; when the fit is restricted to [Fe/H]  <  − 0.5, the slope is steeper and more consistent with the results from Spagna et al. 2010). In contrast to this behavior, the low [α/Fe] content bin, shown in the left panel, has few old stars (5% older than 7 Gyr) while sampling young to intermediate aged stars relatively equally (22% between 0 and 2 Gyr old, and 31% between 4 and 6 Gyr old). Here there is a strong trend between V_ϕ and [Fe/H], with an overall best-fit value of ∼25 km s⁻¹ dex⁻¹. We note that this [α/Fe] decomposition works equally well in the mid-plane as in the "transition zone." It is worthnoting that the slopes we get in V_ϕ versus [Fe/H] are somewhat dependent on the [α/Fe] cut we use to define the thin and thick disk; as we show here, the high [α/Fe] stars show a weak correlation, and this correlation has a slope in the opposite sense as the low [α/Fe] stars.

Zoom In Zoom Out Reset image size

For the low [α/Fe] cut, we have fit the mass-weighted mean values with a single linear fit; however, it is equally plausible to fit two lines here, one at −0.6 dex ⩽ [Fe/H] ⩽−0.2 dex and the other at −0.2 dex ⩽ [Fe/H] ⩽0.4 dex. In that case, the linear fit to low [Fe/H] shows no trend while the linear fit to high [Fe/H] shows a strong trend. Notably, the portion spanning low [Fe/H] is dominated by the intermediate aged stars while the portion spanning high [Fe/H] is dominated by a younger population. This "knee" is a persistent feature; when the upper limit on this [α/Fe] cut is lowered, fewer intermediate age stars are sampled, and the knee in the trend shifts to lower [Fe/H]. Note that the trend is evident here because on average the stars in this [Fe/H] space have experienced less radial mixing than older stars within the same spatial volume. Thus, it is possible to recover this signature of radial mixing even in the absence of age estimates, but with knowledge of [α/Fe], [Fe/H], and V_ϕ for an unbiased population of stars located out of the mid-plane.

We also reconsider here the gradient in V_ϕ as a function of distance from the mid-plane, as shown in Figure 16. As with Figure 15, we decompose V_ϕ by cuts in [α/Fe]; we find that each subsample displays a similar slope (9 versus 11  km s⁻¹ kpc⁻¹). When the full sample is considered, the slope is larger (18  km s⁻¹ kpc⁻¹) because, in addition to the intrinsic slope, there is also transition from thin disk to thick disk, with the latter having a velocity lag of 19  km s⁻¹. Thus, we suggest that if radial migration is significant within the solar cylinder, observational data divided by cuts in [α/Fe] should show similar slopes in V_ϕ versus |z|, but the mean values in the mid-plane should be shifted relative to each other.

Zoom In Zoom Out Reset image size

To date, [α/Fe] measurements have only been accessible for small targeted samples; our initial comparison with a compilation of all currently available data is particularly encouraging (Navarro et al. 2011). We note current work, like that of the SDSS SEGUE collaboration, aims to obtain a large, well-calibrated [α/Fe] data set (Lee et al. 2011). We eagerly anticipate the application of cuts on [α/Fe] to an unbiased population of stars that fall within the region considered by I08 or indeed any large sample within the solar cylinder; such an analysis would further elucidate whether radial mixing has played an important role in shaping the distribution of the MW stars over time.

We have used an N-body model, designed to mimic the quiescent formation and evolution of an MW-type galactic disk (Roškar et al. 2008a, 2008b) to interpret recent SDSS-based observational constraints on the structure of the MW disk. While not a detailed model of the MW, it produces good qualitative agreement with the data and sheds new light on the origin and evolution of the observed disk structure. Of particular importance is the role of radial migration in mixing stars born throughout the disk into the solar neighborhood.

The properties of the Roškar et al. model for the overall spatial, metallicity, and kinematic distributions of the MW stars are in qualitative agreement with the SDSS data (Jurić et al. 2008; Ivezić et al. 2008; Bond et al. 2010). While there are quantitative differences in the spatial gradients of these distributions, as well as in their detailed behavior, even this qualitative agreement is remarkable because the simulation was not fine tuned to match the MW. The model not only reproduces the observed change of slope in the counts of disk stars as a function of distance from the Galactic plane (the original motivation for introducing a separate thick disk component), but also predicts gradients in metallicity and rotational velocity.

The robust qualitative agreements between the data and model predictions motivate the use of model quantities inaccessible to observations, such as the stellar age and the ISM metallicity at the time and position of stellar birth, to interpret recent SDSS results. In particular, the relationship between the metallicity and rotational velocity at ∼1 kpc from the Galactic plane can be understood as due to complex interplay between the ISM metallicity at the time and position of stellar birth, and the subsequent secular evolution largely driven by spiral arms.

No a priori assumptions about the disk's structure are incorporated in the model and the disk evolves entirely in isolation—yet it reproduces the main observational results which for years have motivated decomposition of the disk into two presumably distinct components. The absence of mergers in this model implies that they are not required to explain the overall disk structure. And while merger remnants are detected within the MW disk (e.g., J08), their influence is apparently well localized. Moreover, even if a primordial thick disk is present, having formed via accretion/external heating, it is likely to be substantially polluted by migrating disk stars.

The same mixing effects are likely responsible for the observed differences in α element abundance between kinematically selected thin and thick disk stars. By adopting kinematic selection criteria used by observers, we are able to reproduce distinctive [α/Fe] versus [Fe/H] trends similar to those seen in the data, despite the fact that there are no distinct populations in the model galaxy. We find that cutting on [α/Fe] is to a very good approximation to an age cut, particularly for the α-enhanced population; moreover, separating stars into very old and not-so-old components gives rise to observed correlations.

We look forward to the improved data derived from the emerging generation of surveys such as SEGUE (Rockosi et al. 2009) and APOGEE (Majewski et al. 2010). Key to garnering a deeper understanding of the importance of radial migration in the MW evolution is gathering both precise age determinations and detailed chemical compositions. We are optimistic that this study will lead to further observational and theoretical work.

We thank our numerous SDSS (www.sdss.org) collaborators for their valuable contributions and helpful discussions. We thank the referee, Jerry Sellwood, for a thorough report which greatly improved the original manuscript. S. Loebman and Ž. Ivezić acknowledge support by NSF grants AST-0707901 and AST-1008784 to the University of Washington, and by NSF grant AST-0551161 to LSST for design and development activity. Ž. Ivezić acknowledges support by the Croatian National Science Foundation grant O-1548-2009. This research was supported in part by the NSF through TeraGrid resources provided by TACC and PSC. R. Roškar and T. R. Quinn were supported by the NSF ITR grant PHY-0205413 at the University of Washington.

Recently, Spagna et al. (2010) found that the rotational velocity for disk stars is correlated with metallicity for −1 < [Fe/H] <−0.5 at 1 kpc < z < 3 kpc (and no correlation for [Fe/H] > − 0.5). Notably, they found a gradient within 40–50  km s⁻¹ dex⁻¹, such that more metal-poor stars rotate more slowly. Their claim is in conflict with several other observational studies (Carollo et al. 2010; I08; B10).

We reconsider the Spagna et al. (2010) findings using SDSS DR7 data to try to understand the differences between these observational findings. We note that a direct comparison between I08 and Spagna et al. (2010) is non-trivial as Spagna et al. (2010) use their own proper motion measurements based on the GSC-II catalog and do not provide a comparison to the Munn et al. (2004) proper motions on a star-by-star basis. Additionally, their color selection is more generous (0.0 < g − r < 0.9) than previous studies (0.2 < g − r < 0.6; I08). This more generous selection lets in blue horizontal branch stars, as well as red stars where [Fe/H] reliability decreases. Despite differences in selection criteria, we can reproduce Spagna et al. (2010) results when using the SDSS spectroscopic sample (see the left panel, Figure 17).

Zoom In Zoom Out Reset image size

While the left panel of Figure 17 only considers a narrow bin of z = 1.0–1.5 kpc, we find that the best-fit lines are reproducible in other bins as well. Notably, the median V_ϕ as a function of [Fe/H] closely follows the Spagna et al. (2010) result. We also reproduce a bimodal distribution of stars in the metallicity direction, with modes at [Fe/H] ∼ − 0.65 and [Fe/H] ∼ − 0.4 (see Section 3.2, Spagna et al. 2010).

However, we detect no V_ϕ–[Fe/H] correlation when we consider a complete sample selected in the meridional plane (l ∼ 0 ° or l ∼ 180 °), where proper motion alone suffices to measure rotational velocity (see the right panel, Figure 17). This sample is essentially complete in selected color–distance limits and thus not subject to strong selection effects present in the SDSS spectroscopic sample. This sample shows negligible dependence of the median V_ϕ on [Fe/H] for [Fe/H]  >   − 1.0 (one can also see the bias due to halo stars for [Fe/H]<−1.0). We note that the slope of the V_ϕ versus [Fe/H] relation for the spectroscopic sample decreases when using photometric [Fe/H] instead of spectroscopic [Fe/H] (as well as corresponding distance and V_ϕ) to about 20 km s⁻¹ dex⁻¹(from 40 to 50  km s⁻¹ dex⁻¹). This change in the V_ϕ versus [Fe/H] slope indicates that it might be possible that the relation does have some slope which is masked by systematic errors in the photometric metallicity estimator; however, even if this is so, such a slope must be smaller than that implied by the spectroscopic subsample.

It is possible that the Spagna et al. (2010) results are caused by selection biases in the SDSS spectroscopic sample. Notably, when we fit two Gaussians (one for the disk and halo) to the V_ϕ and [Fe/H] distributions, about 45% of the spectroscopic sample are halo stars (consistent with Spagna et al. (2010) within ∼5% errors), while halo stars make up only 8% of the complete sample. When we compare the two panels in Figure 17 we are using exactly the same volume and exactly the same measurements: the only difference between the samples is that the spectroscopic sample includes only ∼2% of all the stars, with the selection probability about 10 times higher for halo stars than for disk stars. On the other hand, ongoing work with SDSS G dwarf sample suggests that photometric metallicity errors may cause flattening of the V_ϕ versus [Fe/H] (Y. S. Lee 2011, private communication). Irrespective of the resolution of this puzzle, we note that even if the correlation is as strong as suggested by Spagna et al. (2010), the drop in rotational velocity is only ∼20  km s⁻¹ between [Fe/H] = −1 and [Fe/H] = −0.5.