Constraining the Milky Way assembly history with Galactic Archaeology. Ludwig Biermann Award Lecture 2015

TL;DR

Galactic Archaeology aims to reconstruct the Milky Way's history by combining stellar chemical compositions, ages, and dynamical models to understand its formation and evolution.

Contribution

This paper reviews methods to constrain the Milky Way's assembly history using chemo-dynamical modeling and phase-space analysis.

Findings

Constraints on Galactic bar, spiral structure, and satellites from phase-space data

Impact of multiple patterns and radial migration on disc dynamics

Role of mergers in thick disc formation

Abstract

The aim of Galactic Archaeology is to recover the evolutionary history of the Milky Way from its present day kinematical and chemical state. Because stars move away from their birth sites, the current dynamical information alone is not sufficient for this task. The chemical composition of stellar atmospheres, on the other hand, is largely preserved over the stellar lifetime and, together with accurate ages, can be used to recover the birthplaces of stars currently found at the same Galactic radius. In addition to the availability of large stellar samples with accurate 6D kinematics and chemical abundance measurements, this requires detailed modeling with both dynamical and chemical evolution taken into account. An important first step is to understand the variety of dynamical processes that can take place in the Milky Way, including the perturbative effects of both internal (bar and spiral structure) and external (infalling satellites) agents. We discuss here (1) how to constrain the Galactic bar, spiral structure, and merging satellites by their effect on the local and global disc phase-space, (2) the effect of multiple patterns on the disc dynamics, and (3) the importance of radial migration and merger perturbations for the formation of the Galactic thick disc. Finally, we discuss the construction of Milky Way chemo-dynamical models and relate to observations.

Equations

$${\Phi }_0(r)=v_0^2\log (r),$$

$$\Phi_{\rm b}=A_{\rm b}(\epsilon_{\rm b})\cos[2(\phi-\Omega_{\rm b}t)]\times\left\{\begin{array}[]{cclcr}\left(r_{\rm b}\over r\right)^{3}&,&r&\geq&r_{\rm b}\\ 2-\left(r\over r_{\rm b}\right)^{3}&,&r&\leq&r_{\rm b}\end{array}\right.$$

$$\frac{{\overline{v}}_d}{\overline{d}}$$

$${\overline{\mu }}_l$$

$$2C$$

$$2K$$

$${\Phi }_s(r,\varphi ,t)={ϵ}_s\cos [\alpha \ln \frac{r}{r_0}-m(\varphi -{\Omega }_{s}t)],$$

$${ϵ}_s\approx -2\pi G{\Sigma }_s{r}_0/(\alpha v_c^2),$$

$$\begin{gathered} {ϵ}_s\approx -\frac{{\Sigma }_s}{{\Sigma }_0}\frac{r_0}{r_d}\frac{e^{\frac{-r_0}{r_d}} \\ }{0.39\alpha }. \end{gathered}$$

$$\Delta {\sigma }_{\mathrm{z}}=({\sigma }_{\mathrm{z},\mathrm{all}}-{\sigma }_{\mathrm{z},non\_mig})/{\sigma }_{\mathrm{z},\mathrm{all}},$$

Full text

1 Introduction

The goal of Galactic Archeology (Freeman & Bland-Hawthorn, 2002) is to dissect the Milky Way into its various components (discs, bulge, bar and halo) and thus to disentangle the various processes that contributed to their formation and evolution. Galactic Archaeology relies on the assumptions that (i) the dynamics of formation is locked in the phase-space structure of stellar populations and that (ii) stellar atmospheres preserve the chemical imprint of their birth cloud for most of their lifetime. Chemical elements synthesized inside stars are later injected into the interstellar medium (ISM) and incorporated into the next generations of stars. Because different elements are released into the ISM by stars of different masses and on different timescales, stellar abundance ratios are thus directly related to the star formation and gas accretion history. Because stars move away from their birthplaces (a process known as radial migration), chemical information is crucial for understanding the Galactic formation history.

The importance of this topic is manifested in the number of Galactic surveys dedicated to obtaining spectroscopic information for a large number of stars, e.g., RAVE (Steinmetz et al., 2006), SEGUE (Yanny et al., 2009), APOGEE (Majewski et al., 2010), HERMES (Freeman, 2010), Gaia-ESO (Gilmore et al., 2012), and LAMOST (Zhao et al., 2006). This effort will soon be complemented by more than a billion stars observed by the Gaia space mission (Perryman et al., 2001). Millions of these will have accurate proper motions and parallaxes, which together with existing spectroscopic data, and especially with the advent of the dedicated Gaia follow-up ground-based surveys WEAVE (Dalton et al., 2012) and 4MOST (de Jong et al., 2012), will enable Galactic Archaeology as never before.

Before we can hope to understand the past Milky Way history we need a good understanding of its present state, in particular of its disc, where the majority of baryons are concentrated. This is already not a trivial task, due to the Sun’s position close to the Galactic disc midplane - we cannot simply observe the disc morphology as we do in external face-on galaxies. Therefore, mostly indirect methods have been used to constrain the Galactic bar and spiral structure.

While in axisymmetric discs energy and angular momentum are conserved quantities and are, thus, integrals of motion (Binney & Tremaine, 2008), this is not true for the more realistic case of potentials including perturbations from a central bar and/or spiral arms. In the case of one periodic perturbation there is still a conserved quantity in the reference-frame rotating with the pattern – the Jacobi integral $J=E-L\Omega_{p}$, where $E$ is the energy of the particle, $L$ is its angular momentum, and $\Omega_{p}$ is the pattern angular velocity. This is no longer the case, however, when a second perturbation with a different patterns speed is included.

It has now been well established that the Milky Way disc contains both a bar (as in more than 50% of external disc galaxies) and spiral structure moving at different pattern speeds, making it difficult to solve such a dynamical system analytically. Instead, different types of numerical methods are usually employed, from simple test-particle integrations, to preassembled N-body and SPH systems, to unconstrained, fully cosmological simulations of galaxy formation. All of these techniques have their strengths and weaknesses. Test particles are computationally cheap, allow for full control over the simulation parameters (such as spiral and bar amplitude, shape, orientation and pattern speed) but lack self-gravity. N-body simulations offer self-consistency but bar and especially spiral structure parameters are not easy to derive and not well controlled. Finally, in addition to being very computationally intensive, the outcomes of hydrodynamical cosmological simulations are even less predictable, with merging satellites and infalling gas making it yet harder to disentangle the disc dynamics; these are, however, much closer to reality in their complexity and a necessary ultimate step in the interpretation of observational data.

Before we consider more complex systems, we first present simple test-particle models that illustrate the effect of bar and spiral density waves.

1.1 Resonances in galactic discs

Galactic discs rotate differentially with nearly flat rotation curves, i.e., constant circular velocity as a function of galactic radius. In contrast, density waves, such as a central bar and spiral structure, rotate as solid bodies. Therefore stars at different radii would experience different forcing due to the non-axisymmetric structure. Of particular interest are locations in the disc where the stars are in resonance with the perturber. The corotation resonance (CR), where stars move with the pattern, occurs when the angular rotation rate of stars equals that of the perturber. The Lindblad resonances (LRs) occur when the frequency at which a star feels the force due to a perturber coincides with the star’s epicyclic frequency, $\kappa$. As one moves inward or outward from the CR circle, the relative frequency at which a star encounters the perturber increases. There are two values of $r$ for which this frequency is the same as the radial epicyclic frequency. This is where the inner and outer Lindblad resonances (ILR and OLR) are located. Quantitatively, LRs occur when the pattern speed $\Omega_{p}=\Omega\pm\kappa/m$, where $m$ is the multiplicity of the pattern111$m=2$ for a bar or a two-armed spiral structure and $m=4$ for a four-armed spiral.. The negative sign corresponds to the ILR and the positive to the OLR. While Bertil Lindblad defined these for the case of an $m=2$ pattern (thus strictly speaking the ILR/OLR are the 2:1 resonances), for an $m=4$ pattern the ILR/OLR must be the 4:1 resonances.

Since second order resonances, i.e., 4:1 for a two-armed spiral or bar, or 2:1 for a four-armed spiral, can also be quite important (as will be shown in §3.1), we need a convenient way to refer to them. It is somewhat confusing and unclear how the 4:1 resonances are referred to in the literature. The inner 4:1 resonance for an $m=2$ pattern is known as the Ultra-harmonic resonance (UHR). Some also describe the inner and outer 4:1 resonances as the IUHR and the OUHR, others as the inner and outer $m=4$ resonance. If the pattern multiplicity is $m=4$, then these become the ILR and OLR. To our knowledge, there is no terminology for the 2:1 resonances of an m=4 pattern. Minchev & Famaey (2010) proposed to generalize the standard notation of Lindblad resonances by allowing to refer to both 2:1 and 4:1 resonances, regardless of the multiplicity of the pattern. By adopting this nomenclature, we will refer to the 2:1 ILR/OLR and the 4:1 ILR/OLR for both two-armed (or bar) and four-armed spiral structure. Naturally, other resonances can also be described in this manner, e.g., 3:1, 5:1, 6:1 ILR/OLR.

Fig. 1 illustrates the relationship between the pattern angular velocity and the radii at which resonances occur for a flat rotation curve and nearly circular orbits. The red horizontal line indicates a bar pattern speed $\Omega_{b}=55.5$ km/s/kpc (likely for the Milky Way). Inner and outer spiral structures moving with different pattern speeds are also shown by the blue and green lines, respectively. The actual extent of the patterns is indicated by the solid horizontal lines. This is a typical situation seen in N-body simulations for multiplicity of $m=2$ (e.g., Quillen et al. 2011) and a possible configuration for the Milky Way. The vertical red, blue, and green lines give the radial positions of each resonance for the bar, inner, and outer spiral structure, respectively (solid lines: 2:1; dashed lines: 4:1, dotted lines: CR). Note that for two (or more) non-axisymmetric patterns moving at different angular velocities there will always exist regions in the disc where resonances overlap.

1.2 The Milky Way bar and spirals

Due to our position in the Galactic disc, the properties of the Milky Way bar are hard to observe directly. Hence its parameters, such as orientation and pattern speed, have been inferred indirectly from observations of the inner Galaxy (e.g., Blitz & Spergel 1991; Weinberg 1992). The bar has also been found to affect the local velocity distribution of stars. The way this works is as follows. If the Sun happened to lie close to a Lindblad resonance, then the local stellar velocity distribution would exhibit clumps belonging to two different orbital families, which are on nearly closed orbits in the reference frame moving with the pattern (bar or spiral structure). In the case of an OLR, a family of nearly closed orbits supporting the spiral/bar orientation exists, while a second family is misaligned with the structure outside the OLR. These are the $x_{1}(1)$ and $x_{1}(2)$ orbits, respectively (e.g., Dehnen 2000; Fux 2001; Minchev et al. 2010; see Fig. 2, left). Near the ILR the situation is similar but reversed, where the inner orbital family supports the structure but the outer one is misaligned with it.

Hipparcos (Perryman et al., 1997) + GCS (Nordström et al., 2004) data revealed more clearly a stream of old disc stars with an asymmetric drift of about 45 km/s and a radial velocity $u<0$, with $u$ and $v$ positive toward the Galactic center and in the direction of Galactic rotation, respectively. This concentration of stars in the ”$u-v$ plane” is known as the Hercules stream. The numerical work of Dehnen (1999, 2000); Fux (2001); Minchev et al. (2010), and more recently Antoja et al. (2012), has shown that this stream can be explained as the effect of the Milky Way central bar if the Sun is placed just outside the 2:1 Outer Lindblad Resonance (OLR). Due to the inhomogeneity in age and metallicity of Hercules stars, a number of works (e.g., Famaey et al. 2005, 2007; Bensby et al. 2007) have concluded that a dynamical effect, such as the influence of the bar, is a more likely explanation than a dispersed cluster. Most estimates agree on a bar orientation222See Fig. 2, second panel, for the definition of the bar angle. of $\phi_{b}=30^{\circ}\pm 10$ and a pattern speed of $\Omega_{b}/\Omega_{0}=1.9\pm 0.1$, where $\Omega_{0}$ is the local standard of rest (LSR) rotation rate. A drastically different bar pattern speed has recently been suggested by the longer bar half-length measured by Wegg et al. (2015), compared to previous works – $r_{b}=5.0\pm 0.2$ kpc, which places the CR at 5-7 kpc. In such a case the bar 2:1 OLR would lie in the range 8.5-12 kpc and, thus, would not work as an explanation for the Hercules stream.

The Milky Way spiral structure is more poorly known than the bar. Cepheid, HI, CO and far-infrared tracers suggest that the Milky Way disc contains a four-armed tightly wound structure (see also Vallée 2016), whereas Drimmel & Spergel (2001) have shown that the near-infrared observations are consistent with a dominant two-armed structure. A dominant two-armed and a weaker four-armed structure has been proposed by Amaral & Lepine (1997).

Similarly to the effect of the bar, the spirals can be linked to clumps in the $u-v$ plane, as first shown by Quillen & Minchev (2005). Using an orbital weighting function technique, this work showed that a two-armed spiral density wave with pattern speed placing the Sun near the 4:1 inner Lindblad resonance can account for two major clumps in the solar neighborhood’s velocity distribution: the Pleiades/Hyades moving group corresponds to the one family of orbits, and the Coma Berenices moving group corresponds to another family. Similar patter speed estimate was obtained by Pompéia et al. (2011) and Siebert et al. (2012).

1.3 Multiple patterns in galactic discs

Multiple patterns in N-body simulations have been known to exist since the work of Sellwood (1985) and Sellwood & Sparke (1988), who found that a bar can coexist with a spiral pattern moving at a much lower angular velocity. Tagger et al. (1987) and Sygnet et al. (1988) explained this as the non-linear mode coupling between the bar and the spiral wave. These findings were later confirmed by the numerical studies of Masset & Tagger (1997) and Rautiainen & Salo (1999). According to the theoretical work by Tagger et al. (1987) and Sygnet et al. (1988), two patterns can couple non-linearly as they overlap over a radial range, which coincides both with the CR of the inner one and the ILR of the outer one. This coincidence of resonances results in efficient exchange of energy and angular momentum between the two patterns. The coupling between the two patterns generates beat waves (as we describe below), also found to have LRs at the interaction radii, resulting in a strong non-linear effect even at relatively small amplitudes. Rautiainen & Salo (1999) showed that coupling between a CR and 4:1 ILR, as well as ILR and OLR is also possible in N-body simulations. Waves couple with a selection of frequencies which optimizes the coupling efficiency. Strong exchange of energy and angular momentum is then possible among the coupled waves.

How do multiple patterns affect the dynamics of galactic discs? Quillen (2003) considered the dynamics of stars that are affected by perturbations from both spiral arms and a central bar by constructing a one-dimensional Hamiltonian model for the strongest resonances in the epicyclic action-angle variables. Quillen pointed out that when two perturbers with different pattern speeds are present in the disc, the stellar dynamics can be stochastic, particularly near resonances associated with one of the patterns. Similar findings were presented more recently by Jalali (2008). All these results are not surprising since it has already been shown by Chirikov (1979) that in the case of resonance overlap the last KAM surface between the two resonances is destroyed, resulting in stochastic behavior. It is therefore expected that resonance overlap could give rise to both velocity dispersion increase and radial migration. Both of these possibilities are explored in §4.1 and §4.2, respectively.

A comprehensive discussion of the Galactic bar and spiral parameters can be found in the recent review by Bland-Hawthorn & Gerhard (2016).

1.4 Radial migration

The power of Galactic Archaeology has been threatened both by observational and theoretical results, showing that stars most probably move away from their birthplaces, i.e, migrate radially. Observational signatures of radial migration (or mixing) have been reported in the literature since the 1970’s, with the pioneering works by Grenon (1972, 1989). Grenon identified an old population of super-metal-rich stars (hereafter SMR), presently at the Solar vicinity, but with kinematics and abundance properties indicative of an origin in the inner Galactic disc. SMR stars show metallicities which exceed the present day ISM and those of young stars at the solar vicinity. The metallicity of the solar vicinity, however, is not expected to increase much since the Sun’s formation, or in the last $\sim$4 Gyr, due to the rather inefficient star formation rate (SFR) at the solar radius during this period (e.g., Chiappini et al. 2003; Asplund et al. 2009). Hence, pure chemical evolution models for the Milky Way thin disc cannot explain stars more metal rich than $\sim$0.2 dex and the effect of radial migration needs to be considered.

N-body simulations have also long shown that radial migration is unavoidable. Raboud et al. (1998) studied numerical simulations aimed at explaining the results by Grenon (1989) of a mean positive $u$-motion, which the authors interpreted as metal-rich stars from the inner galaxy, wandering in the solar neighborhood. However, Raboud et al. (1998) interpreted their findings as stars on hot bar orbits, not recognizing that permanent changes to the stellar angular momenta are possible. It was not until the work by Sellwood & Binney (2002) that radial migration was established as an important process affecting the entire disc, where stars shift guiding radii due to interaction with transient spiral structure. A number of works using N-body and N-body/SPH simulations (e.g., Roškar et al. 2008; Grand et al. 2012; Minchev et al. 2011) have confirm that migration always takes place in numerical galactic discs.

A different radial migration mechanism was proposed by Minchev & Famaey (2010), who considered the simultaneous propagation of a bar and a long-lived spiral density wave. In such a setup angular momentum redistribution arises from the overlap of resonances associated with different modes or from the stellar mass transiency as perturbers with different pattern speeds interfere constructively (see Comparetta & Quillen 2012), but not from the growth and decay of transient modes. This work, along with studies of diffusion coefficients in barred discs (Brunetti et al., 2011; Shevchenko, 2011), predicts a variation in migration efficiency with time and disc radius, establishing that the dynamical influence of the bar plays an integral part of Milky Way disc modeling. Aside from internal structure, perturbations due to minor mergers have also been shown to be effective at mixing the outer discs (Quillen et al., 2009; Bird et al., 2012), but also can, indirectly, affect the entire disc by inducing (reinforcing) spiral and bar instabilities (e.g., Purcell et al. 2011). Considering the established presence of a central bar, spiral structure and evidence for merger activity in the Milky Way, it is clear that all of the above mentioned radial migration mechanisms would have an effect on the Galactic disc.

Stars found today at a given small Galactic disc region (e.g., the solar neighborhood) can have a range of birth radii but are mostly indistinguishable in their kinematics from locally born stars. Therefore, chemical abundances need to be invoked in identifying migrators. Probably one of the best ways of quantifying radial migration in the Milky Way is using the technique of chemical tagging (Bland-Hawthorn et al., 2010), where stars born in the same cluster (now dispersed) are expected to appear as a clump in the multi-dimensional chemical space. This is one of the main objectives of the ongoing GALAH survey (Freeman, 2010).

1.5 The Galactic thick disc

The formation of galactic thick discs has been an important topic ever since their discovery in external galaxies (Burstein, 1979; Tsikoudi, 1979) and in the Milky Way (Gilmore & Reid, 1983). The large uncertainties in important observational constraints in the Milky Way, such as the age-velocity-metallicity relation, abundance gradients and their evolution, have led to different scenarios to be proposed for the formation of the Galactic thick disc.

One possibility is that the stars comprising thick discs are born thick at high redshift from internal gravitational instabilities in gas-rich, turbulent, clumpy discs (Bournaud et al., 2009; Forbes et al., 2012) or in the turbulent phase associated with numerous gas-rich mergers (Brook et al., 2004, 2005). They could also have been created through accretion of galaxy satellites (Meza et al., 2005; Abadi et al., 2003), where thick disc stars then have an extragalactic origin.

Another possibility is that thick discs are created through the heating of preexisting thin discs by minor mergers (Villalobos & Helmi, 2008; Di Matteo et al., 2011). Evidence for merger encounters can be found in the phase-space structure of Milky Way disc stars (e.g., Minchev et al. 2009; Gómez et al. 2012b, 2013).

Several works (e.g., Schönrich & Binney 2009b; Loebman et al. 2011) have also proposed that radial migration can give rise to thick disc formation by bringing out high-velocity-dispersion stellar populations from the inner disc and the bulge. More detailed dynamical studies (Minchev et al., 2012a; Martig et al., 2014b; Vera-Ciro et al., 2014) have more recently shown that migration does not contribute to any significant level to disc thickening, but on the opposite, it suppresses flaring when external perturbations are included (Minchev et al., 2014a; Grand et al., 2016). This is discussed further in §6.1.

Finally, a recently proposed thick disc formation mechanism is the superposition of coeval flaring disc subpopulation (Minchev et al., 2015). Disc flaring can result from a number of different sources, the most effective most likely being the perturbative effect of mergers on the host disc. In an inside-out formation scenario, the old populations dominate in the inner disc and younger in the outer disc, thus the flaring disappears when the total stellar population is taken into account. By making the distinction between a thick disc defined as the [$\alpha$/Fe]-high, old stellar population in the Milky Way and a geometrically defined thick discs in observations of edge-on galaxies, this new view of the formation of thick discs resolves a number of apparent contradictions. More on this can be found in §7.

2 Constraining the Galactic bar

2.1 The effect of a recently formed bar

While earlier work has shown the bar to affect the $u-v$ plane at higher velocities, some of the low-velocity moving groups in the solar vicinity have been explained only with the effect of spiral density waves, such as the splitting of the Pleiades/Hyades and Coma Berenices moving groups (e.g., Quillen & Minchev 2005). The reason the bar’s OLR should not affect stars with kinematics cooler than the Hercules stream is that the solar radius lies about 0.8 kpc outside the OLR (dashed circle), thus, stars on near circular orbits at the OLR cannot reach the solar neighborhood. This situation can be seen in the leftmost panel of Fig. LABEL:bar1, where the Hipparcos volume is shown by the yellow-filled small circle. Hotter stars near the OLR, however, can appear close to the Sun on their apocenters.

Studying the effect of the bar, we explored a different scenario in Minchev et al. (2010), where the time evolution of the $u-v$ plane was examined just after the bar formation. These test-particle simulations were performed in a Milky Way-like potential to which a perturbation due to a central bars was added, similar to Dehnen (2000). The background axisymmetric potential due to the disc and halo has the form

[TABLE]

corresponding to a flat rotation curve. We model the nonaxisymmetric potential perturbation due to the Galactic bar as a pure quadrupole

[TABLE]

Here $A_{\rm b}(\epsilon_{\rm b})$ is the bar’s gravitational potential amplitude, identical to the same name parameter used by Dehnen (2000); the strength is specified by $\epsilon_{\rm b}=-\alpha$ from the same paper. The bar length is $r_{\rm b}=0.8r_{\rm cr}$ with $r_{\rm cr}$ the bar corotation radius. The pattern speed, $\Omega_{\rm b}$ is kept constant. The bar amplitude $\epsilon$ is initially zero, grows linearly with time at $0<t<t_{1}$ and transitions smoothly to a constant value after $t=t_{1}=4$ bar rotations. This insures a smooth transition from the axisymmetric to the perturbed state.

Previous work has used as a measure of bar strength the parameter $Q_{T}$ (Combes & Sanders, 1981). This is the ratio of the maximum tangential force to the azimuthally averaged radial force at a given radius. From eq. 2 this definition yields $Q_{T}=2A_{b}/v_{c}^{2}$. We examine bar amplitudes in the range $0.1<Q_{T}<0.4$ as expected from observations of various galaxies and from N-body simulations (Combes & Sanders, 1981). This corresponds to $0.013<|\epsilon_{b}|<0.05$.

In our units the solar neighborhood radius is $r_{0}=1$; the circular speed is $v_{0}=1$ everywhere since the rotational curve is flat. To convert to real units we use a LSR tangential velocity of 240 km/s, and Galactocentric distance of 8 kpc. The 2:1 OLR with the bar is achieved when $\Omega_{\rm b}/\Omega_{0}=1+\kappa/2\approx 1.7$, where $\kappa$ is the epicyclic frequency. For a flat rotation curve $\kappa=\sqrt{2}\Omega_{0}$, where $\Omega_{0}$ is the angular velocity of the LSR.

Interestingly, we found that a steady state bar induced transient features at low velocities in the simulated solar neighborhood velocity distribution due to the initial response of the disc to the bar formation. We associate these velocity streams with two quasi-periodic orbital families librating around the stable $x_{1}(1)$ and $x_{1}(2)$ orbits near the bar’s OLR (Fig. 2, left two panels). In a reference frame moving with the bar these, otherwise stationary, orbits precess on a timescale dependent on the strength of the bar. The effect of this precession can be seen in Fig. 2, where the left panel illustrates the $x_{1}(1)$ and $x_{1}(2)$ orbits and their orientation with respect to the bar. The second panel shows the results of a test-particle simulation from Minchev et al. (2010) at time $t=7.5$ (in units of solar rotations). Only stars initially on circular orbits close to the OLR are shown, with those inside/outside the OLR colored in blue/red. Note that, because of the precession, at this time the $x_{1}(1)$ and $x_{1}(2)$ orbits orientation is not as in the left panel but offset by about $30^{\circ}$ in the clockwise direction. This behavior allows the two (kinematically cold) orbital families to reach the solar neighborhood and manifest themselves as clumps in the $u-v$ plane moving away from ($x_{1}(2)$), and toward ($x_{1}(1)$) the Galactic center.

The nine panels on the left in Fig. 3 show the simulated $u-v$ plane for nine neighborhoods at different Galactic disc positions, as indicated by the yellow circles in the x-y plot on the right. The solar neighborhood volume is shown as a black circle. For this particular simulation, only test-particles on nearly circular orbits were subjected to the potential in order to avoid the sea of hot orbits, which we already know give rise to Hercules stream-like feature (see Fig. 2 third panel).

For a given time, the tangential velocity $v$ of resonant features in the $u-v$ plane is set by the bar pattern speed and the radial velocity $u$, is set by the bar’s orientation. These variations are illustrated by the blue lines in the figure. Thus, assuming some features at low velocities in the Hipparcos velocity distribution are of resonant origin, we can match stream positions in the $u-v$ plane and estimate the bar pattern speed $\Omega_{b}$ and the orientation of its major axis with respect to the solar Galactocentric line, $\phi_{0}$. In addition, structure varies with time due to the libration of the quasi-periodic orbits around the fixed points, which allows us to constrain the bar formation time.

Depending on the bar parameters and time since its formation, this model is consistent with the Pleiades and Coma Berenices, or Pleiades and Sirius moving groups seen in the Hipparcos stellar velocity distribution, when the Milky Way bar angle is $30^{\circ}\la\phi_{0}\la 45^{\circ}$ and its pattern speed is $\Omega_{b}/\Omega_{0}=1.82\pm 0.07$. Since the process is recurrent, a good match could be achieved about every six LSR rotations. However, to be consistent with the fraction of stars in the Pleiades, we estimated that the Milky Way bar formed $\sim 2$ Gyr ago. This model argues against a common dynamical origin for the Hyades and Pleiades moving groups.

We the advent of Gaia accurate proper motions and parallaxes, combined with high precision line-of-sight velocities from APOGEE and RAVE (already after Gaia DR1 and DR2) and ultimately by the 4MOST survey (scheduled to start operation in 2021), we will be able to test these predictions by trying to reproduce the variation of structure in the $u-v$ plane at different positions in the disc (as in the left half of Fig. 3).

2.2 Relation to the Oort C constant

As discussed above, previous work has related the Galactic bar to structure in the local stellar velocity distribution.

Another indirect way to constrain the bar was presented by Minchev et al. (2007), where we showed that the bar also influences the spatial gradients of the velocity vector via the Oort constants.

We can linearize the local velocity field about the LSR and write the mean radial velocity $\overline{v}_{d}$ and longitudinal proper motion $\overline{\mu}_{l}$ as functions of the Galactic longitude $l$ as

[TABLE]

where $\overline{d}$ is the average heliocentric distance of stars, $A$ and $B$ are the usual Oort constants, and $C$ and $K$ are given by

[TABLE]

Here $r$ and $\phi$ are the usual polar coordinates and $v_{\phi}=v_{0}+v$, where $v_{0}$ is the circular velocity at the Solar radius, $r_{0}$. Considering a flat rotation curve, the derivatives of $v_{\phi}$ in the above equations are identical to the derivatives of $v$. $C$ describes the radial shear of the velocity field and $K$ its divergence. For an axisymmetric Galaxy we expect vanishing values for both $C$ and $K$333Note, however, that $C$ and $K$ would also be zero in the presence of nonaxisymmetric structure if the Sun happened to be located on a symmetry axis.. Whereas $C$ could be derived from both radial velocities and proper motions, $K$ can only be measured from radial velocities, in which case accurate distances are also needed.

The study of Olling & Dehnen (2003) not only measured a non-zero $C$, implying the presence of non-circular motion in the local disc, but also found that $C$ is more negative for older and redder stars with a larger velocity dispersion, while it is roughly zero for the cold (thus young) population. More recent determination of the $C$ constant using RAVE data confirmed these results (Siebert et al., 2011). This variation of $C$ with velocity dispersion/color/age is surprising as a hotter stellar population is expected to have averaged properties more nearly axisymmetric, and hence, a reduced value of $|C|$ (e.g., Minchev & Quillen 2007).

Assuming the Galactic bar affects the shape of the distribution function of the old stellar population in the solar neighborhood, an additional constraint on the bar can be provided by requiring that a model reproduces the observed value of the Oort constant $C$. In other words, in addition to relating the dynamical influence of the bar to the local velocity field (see § 2.1 above), $C$ provides a link to the gradients of the velocities as well.

In an effort to reproduce these observational results, using test-particles integrations we simulated measurements of the Oort $C$ value in a gravitational potential including the Galactic bar, as in § 2.1 above. For all other parameters fixed, we defined a cold and a hot sample resulting from an initial radial velocity dispersion of 10 km/s and 40 km/s, respectively.

In Fig. 4 we present our results for $C$ as a function of the bar angle, $\phi_{0}$ (the angle by which the Sun’s azimuth lags the bar’s major axis). Each panel shows a simulation with a different pattern speed and sample depth of $\overline{d}=200$ pc. Solid and dotted lines represent the results for cold and hot discs, respectively. $C$ is presented in units of $\Omega_{0}=v_{0}/r_{0}$. To make the discussion less cumbersome, we write $C_{\rm h}$ and $C_{\rm c}$ to refer to the values for $C$ as estimated from the hot and cold discs, respectively.

$C_{\rm h}$ (dotted lines in Fig. 4) varies with galactic azimuth as $C_{\rm h}(\phi_{0})\sim\sin(2\phi_{0})$ for all of the $\Omega_{\rm b}$ values considered. On the other hand, the cold disc values (solid lines) exhibit different variations, depending on the bar pattern speed or equivalently, on the ratio $r_{0}/r_{OLR}$. Closer to the OLR (left panels of Fig. 4), $C_{\rm c}(\phi_{0})$ approaches the functional behavior of $C_{\rm h}(\phi_{0})$. Away from the OLR (right panels), $C_{\rm c}(\phi_{0})$ is shifted by $90^{\circ}$ compared to $C_{\rm h}(\phi_{0})$, i.e., $C_{\rm c}(\phi_{0})\sim-C_{\rm h}(\phi_{0})$. Moreover, the hot disc yields an increase in the amplitude of $C_{\rm h}(\phi_{0})$ as the pattern speed nears the OLR. Whereas the same trend is apparent for the cold disc, the gradient of the amplitude of $C_{\rm c}(\phi_{0})$ is much larger. This is consistent with our expectation that the cold disc is affected more by the bar, especially near the OLR. While close to the OLR $|C_{\rm h}(\phi_{0})|<|C_{\rm c}(\phi_{0})|$, we observe the opposite behavior away from it. This could be explained by the results of Mühlbauer & Dehnen (2003), where it was found that high velocity dispersion stars tend to shift the “effective resonance” radially outwards.

The red vertical lines in Fig. 4 indicate possible solar azimuths where the observed Oort $C$ trends are well matched. By comparing measurements of $C$ with our simulations we constrained the pattern speed as $\Omega_{\rm b}/\Omega_{0}=1.87\pm 0.04$, where $\Omega_{0}$ is the local circular frequency, and found the bar angle to lie in the range $20^{\circ}\leqslant\phi_{0}\leqslant 45^{\circ}$.

3 Constraining the Milky Way spiral structure

3.1 Constraints using the U-V plane

Using a spiral pattern speed of $\sim 20$ km/s/kpc (assuming LSR rotation of 28 km/s/kpc), Quillen & Minchev (2005) used an orbital weighting function technique to show that a two-armed spiral density wave can split the solar neighborhood’s velocity distribution into two major clumps. A different approach but with a similar result was presented by Pompéia et al. (2011) and discussed below. In that work a test-particle simulation by Minchev & Quillen (2008) was employed, in which initially a disc was populated by stars on circular orbits and then numerically integrated in the axisymmetric potential from eq. 1, but including a spiral, instead of a bar perturbation term.

The spiral potential is given by

[TABLE]

where $\epsilon_{s}$ is the spiral strength, related to the amplitude of the mass surface density of spirals, $\Sigma_{s}$, as

[TABLE]

as shown in Binney & Tremaine (2008). The parameter $\alpha$ is related to the pitch angle of the spirals, $p$, by $\alpha=m\cot(p)$. The azimuthal wavenumber $m$ is an integer corresponding to the number of arms. We consider both two-armed and four-armed spiral structure with $\alpha=-4$ and $-8$, respectively, where the negative sign corresponds to trailing spirals. Elmegreen (1998) found that grand-design spirals have arm-interarm contrasts of 1.5-6, corresponding to a fractional amplitude of $0.2<\Sigma_{s}/\Sigma<0.7$, which is in agreement with Rix & Zaritsky (1995) who estimated $0.15<(\Sigma_{s}/\Sigma)<0.6$.

For a maximum exponential disc the peak circular speed in the disc has the value $v_{c}\simeq 0.622\sqrt{GM_{d}/r_{d}}$ at $r\simeq 2.15r_{d}$, where $M_{d}$ is the disc mass inclosed by $r_{d}$. The surface density of the disc at radius $r$ is $\Sigma(r)=({M_{d}/2\pi r_{d}^{2}})e^{-r/r_{d}}$. Eliminating $M_{d}$ from these expressions, we find at $r_{0}$, $v_{c}^{2}\simeq 0.39\times 2\pi Gr_{d}\Sigma_{0}e^{r_{0}/r_{d}}$. Substituting this expression for $v_{c}$ in eq. 7, the relation between the relative potential and the relative overdensity becomes

[TABLE]

The perturbation is grown from zero to its maximum strength in four rotation periods at $r_{0}$. In order to improve statistics, positions and velocities are time averaged for 10 spiral periods. We distribute particles (stars) between in inner and outer galactic radii $(r_{in},r_{out})=(0.3r_{0},2.0r_{0})$. New particles are added until the final number of outputs is $2.5\times 10^{6}$. In addition, the two-fold symmetry of our model galaxy is used to double this number. More details about these simulations can be found in Minchev & Quillen (2008).

The top panel of Fig. 5 shows the effect of a two-armed spiral structure on orbits near the 4:1 ILR444Also known as the inner ultra-harmonic resonance or IUHR.. This 4:1 resonance gives rise to square orbits, in the frame moving with the spiral pattern, even though the imposed spiral wave is two-armed. Similarly to the case of the bar (see §2), two families of closed orbits are excited by the resonance, where one supports the spiral structure (blue particles) and the other one is misaligned with it (red particles). $\phi_{0}$ indicates the orientation of spiral structure with respect to the Sun’s azimuth, given by the angle between the dotted Galactocentic arrows. For a Sun orientation at $20^{\circ}$ with respect to a concave arm, both orbital families can enter the solar neighborhood (yellow circle). The Galactocentric axes are in units of $r_{0}$, the Galactocentric radius of the Sun.

The bottom panel of Fig. 5 shows the effect on the $u-v$ plane for the configuration shown in the top panel. To match spatially the Hipparcos stellar sample, only particles in a 200 pc circle around the fiducial Sun are selected. Each orbital family gives rise to a stream (or moving group) in velocity space. The dense clump at $(u,v)\approx(-35,-17)$ km/s can be associated with the Hyades and the shallow one at $(u,v)\approx(10,0)$ km/s with the Sirius moving groups. The contour levels correspond to 0.2, 0.31, 0.43, 0.55, 0.67 and 0.8 of the maximum value at the centre of the Hyades clump.

This model constrains the spiral pattern speed to $\Omega_{s}/\Omega_{0}=0.65$ and orientation with respect to the Sun, $\phi_{0}=20^{\circ}$ to an uncertainty of less than 5%.

3.2 The effect of spiral structure on the global disc phase space

Another way to obtain constraints on the Milky Way spiral structure is through stellar samples covering large parts of the Galactic disc. To facilitate the interpretation of the huge amounts of data expected from Gaia and spectroscopic follow-up surveys, studies of how spiral structure affects the global disc phase space are needed. One such attempt was made by Minchev & Quillen (2008), where test-particle simulations of a galactic disc perturbed by a steady-state spiral density wave used to relate structure in velocity and morphology to the spiral parameters.

In the top row of Fig 6 we present stellar number density contour plots for two such simulations with different spiral pattern speeds, as indicated. The background axisymmetric disc is subtracted to emphasize the spiral structure. The quantity plotted is $\Delta\Sigma=(\Sigma-\Sigma_{axi})/\Sigma_{axi}$, where $\Sigma$ and $\Sigma_{axi}$ are the perturbed and axisymmetric stellar number densities. The resonances get closer together for the faster spiral on the right. Distances are in units of the solar radius, $r_{0}$. Darker colors correspond to lower density. The inner $0.3r_{0}$ disc is not plotted since we do not model the Galactic center. Note the disruption of the spirals near the 2:1 LRs (dashed circles).

The second and third rows of Fig 6 show maps of the residual mean line-of-sight velocity, $\Delta V_{d}$, and the residual stellar density, $\Delta\Sigma$, for the disc in the top-left panel. Here, however, these quantities are plotted versus Galactic longitude $l$ (x-axis), and heliocentric distance $d/r_{0}$ (y-axis). The orange curve is the projection of the solar circle. The solar phase angle is $\phi_{0}=20^{\circ}$ in the second row and $\phi_{0}=40^{\circ}$ in the third row. The minimum (black) and maximum (white) contour values are displayed on top of each panel. Well defined structure is present for both $\Delta V_{d}$ and $\Delta\Sigma$ and that is seen to change with the change in solar orientation with respect to the spirals.

The fourth and fifth rows of Fig 6 are similar to the two above, but show the results for a four-armed structure with a phase angle of $20^{\circ}$ and a pattern speed $\Omega_{s}=0.7\Omega_{0}$ (fourth row) and $\Omega_{s}=0.9\Omega_{0}$ (fifth row). Thus, the variation in structure seen here is due to the different pattern speeds used.

Fig 6 indicates that (1) the solar orientation with respect to the spiral, (2) the spiral pattern speed and (3) the number of spiral arms can be associated with structure in these observationally motivated maps.

Minchev & Quillen (2008) found that the axisymmetric potential needs to be known to $\sim 10\%$, line-of-sight velocities to $\sim 20$ km/s, and distance uncertainties need to be less than $\sim 30\%$, in order to be able to constrain spiral structure. The mean line-of-sight velocity and the velocity dispersion are affected by up to $\sim 35$ km/s which is well within the detectable limit of even low-resolution spectroscopic surveys.

A survey of stars close to the disc midplane and covering large area is need to apply this method of constraining the spirals. Such large-scale Galactic surveys have now become a reality, most notably, the APOGEE-1 near-infrared SDSS project, which covers Galactocentric azimuth of about $60^{\circ}$ in the radial range $3<r<14$ kpc for some $1.5\times 10^{5}$ stars. The surveyed disc area will soon increase with the forthcoming APOGEE-2 (north and south) and the addition of $\sim 3\times 10^{5}$ more stars. In the very near future Gaia and follow-up surveys 4MOST and WEAVE will increase this number to millions of stars with accurate proper motions, parallaxes, radial velocities, photometry, and chemical composition.

4 The importance of multiple patterns in galactic discs

So far we have only considered simple test-particle models, but which allow a high degree of control and the means of cheaply sweeping parameter space. Only the effect of a bar or that of a spiral perturbation was used for the models discussed in § 2 and § 3. Using test-particle simulations, we could isolate, understand, and quantify better the effect of an individual perturber, i.e., the bar or a spiral wave. It is well known, however, that the Milky Way disc, as well as more than 50% of disc galaxies in general, harbor both types of non-axisymmetric components. Moreover, both analyses of N-body simulations (e.g., Tagger et al. 1987; Sygnet et al. 1988) and expanding galaxy images in Fourier components (e.g, Elmegreen et al. 1992; Rix & Rieke 1993) have shown the existence of multiple spiral density waves, which can propagate simultaneously in galaxy discs.

Complicating the disc dynamics by the inclusion of multiple patterns and studying their combined effect on the kinematical heating of, and redistribution of angular momentum (radial migration) in, galactic discs is the topic of the following §4.1 and §4.2, respectively. In §4.3 we show that, indeed, N-body SPH simulations support the idea that spiral patterns can be long-lived features, justifying the assumptions we have been making so far.

4.1 Stochastic heating from multiple spiral waves

A new mechanism for increasing the stellar velocity dispersion with time in galactic discs (known as the age-velocity relation) was described by Minchev & Quillen (2006), where we studied the effect of two spiral density waves propagating with different pattern speeds. The simulations performed were similar to the ones described in §3, however, two perturbation terms in the form of eq. 6 were added to the axisymmetric potential (eq. 1) and evolved in time. While the two spiral perturbations were described the same way, their number of arms, strengths (or amplitudes) $\epsilon_{s,1}$, $\epsilon_{s,2}$ and pattern speeds $\Omega_{s,1}$, $\Omega_{s,2}$ were allowed to differ. Both spiral amplitudes were grown in four solar rotations, as in §3, after which they were kept constant.

Minchev & Quillen (2006) adopted a configuration of a primary two-armed and a weaker four-armed spiral wave perturbations, as done previously (e.g., Amaral & Lepine 1997) . The novelty in this model was the assumption that there was a non-zero relative angular velocity between the two spiral patterns. This introduces an additional parameter - the pattern speed of the secondary spiral wave.

Fig. 7 illustrates how two long-lived spiral density waves propagating with different pattern speeds affect initially circular test-particle orbits. The first column shows a face-on view of particles starting initially on a ring just inside (red) and a ring just outside (green) the CR (dashed black circle) of a two-armed spiral wave of an intermediate strength. Top and bottom plots are identical for the first column. Note the unbroken banana, or “horseshoe”, shape of the orbits.

The second column of Fig. 7 presents the same initial conditions, but near the 4:1 OLR (top) or the 4:1 ILR (bottom) of a secondary, four-armed wave. This spiral has an amplitude only 1/3 that of the primary, two-armed wave. Only small distortions in the particles’ initially circular orbits are apparent. For all three single-spiral simulations no significant time variation was observed.

In the third column of Fig. 7 the two rings of particles are perturbed by both spiral waves shown in each row, shown at t=10 rotation periods, which is about 2.5 Gyr of evolution. Both top and bottom plots are in the reference frame of the corotating, two-armed spiral (not shown). Remarkably, the secondary wave, which has an amplitude of only 1/3 that of the first one, is enough to disrupt the banana orbits near the CR. Moreover, the distances travelled by the particles are now significantly larger, especially in the top panel. The irregular shapes of these orbits suggests the presence of stochasticity, which causes an increase in the particle velocity dispersion, i.e., the disc is heated kinematically with time.

The fourth column of Fig. 7 shows the same setup as in the third column but at t=40 rotations, amounting to about 10 Gyr of evolution. At this time the particles appear well mixed, especially at the CR+OLR region (top panel). Concentration of particles at the CR Lagrange points is seen in the bottom panel.

Previous studies of heating mechanisms have concentrated on heating from transient spiral density waves (e.g., Carlberg et al. 1985; De Simone et al. 2004; Jenkins 1992). However in that case changes in the stellar velocity dispersion only occur during the spirals’ growth and decay. Minchev & Quillen (2006) observed an increase in the velocity dispersion even when the spiral density waves have ceased to grow. Hence, this is a different heating mechanism compared to those explored by these previous investigations.

We note that the above described process causes not only heating but also radial migration, i.e., exchange of angular momentum across the resonance overlap region (discussed in §4.2). In other words, some of the particles away from their initial radii in Fig. 7 are there on their apo- or pericenters, while others have shifted their guiding radii altogether.

Minchev & Quillen (2006) examined how different combinations of pattern speeds, at a given radius, affect the stellar velocity dispersion, finding that this heating mechanism is strongly dependent on the relative velocity between stars and the spiral patterns and thus on the Galactocentric distance. This expected radial variation may be erased when radial mixing of stars is taken into account.

4.2 Migration due to the overlap of long-lived patterns

As shown by Minchev & Quillen (2006), when two perturbers moving at different pattern speeds are imposed on a galactic disc, an increase in the random motions of stars is expected, i.e., the disc heats. In addition, it was found that in regions of resonance overlap stars drift radially with time from their birth radii (see Fig. 7). This suggests that such a mechanism could be responsible for radial migration in the disc. Minchev & Famaey (2010) quantified this idea more clearly for the case of a bar and one set of spiral structure, by examining the effect of the simultaneous propagation of the two perturbers. Test-particle simulations were used with a setup similar to that described in §4.1 above, but where the axisymmetric disc potential from eq. 1 is perturbed by a central bar (eq. 2) and a spiral density wave (eq. 6), instead of two spiral waves.

In the first row of Fig. 8 we first show how much migration is induced in a galactic disc by a single perturber – either a bar or a spiral density wave. Plotted are contours of stellar density for the changes in guiding radius, $\Delta r_{g}$, as a function of the initial guiding radius, $r_{g,0}$. Bar and spiral amplitudes are indicated in each panel. The first two panels form left to right present simulations with a bar of an intermediate strength ($Q_{T}=0.25$) and a strong one ($Q_{T}=0.4$), respectively. The third and fourth panels show a four-armed spiral structure with relative overdensity $\Sigma_{s}/\Sigma_{0}=0.2,0.35$. The rightmost panel presents a simulation with a two-armed spiral structure with an amplitude giving rise to $\Sigma_{s}/\Sigma_{0}=0.25$. The different vertical lines depict the 2:1 ILR, 4:1 ILR, CR, 4:1 OLR and 2:1 OLR, as indicated. Depending on the pattern speed, some resonances might not be present in the disc. $|\Delta r_{g}|$ increases significantly only near the corotation of each perturber.

In the second row of Fig. 8 we show the time development of the changes of guiding radius in a stellar disc perturbed by both a central bar and spiral structure. Panels from left to right show the temporal evolution of the system for 0.3, 0.6, 1.2, 1.8 and 2.5 Gyr. Both the bar and the four-armed spiral structure have intermediate strengths: $\epsilon_{b}=-0.03$ ($Q_{T}=0.25$) and $\epsilon_{s}=-0.02$ ($\Sigma_{s}/\Sigma_{0}=0.35$ for m=4). At the beginning of the simulation $|\Delta r_{g}|$ increases mainly at the CR of each perturber (dotted lines), which would be the case of linearly adding the individual effects of the bar and spirals seen in the top row. At later times, however, large changes in angular momentum occur, suggesting that non-linear effects become important. The changes in guiding radii (angular momentum) increase throughout the simulation timespan manifesting the effect of multiple overlapping patters.

What this experiment demonstrated for the first time was that transient patterns are not a necessity for efficient migration (e.g., Sellwood & Binney 2002), but also long-lived, multiple perturbers can be very effective. Because of the hard task of inferring the longevity of spirals in N-body systems, our controlled simulations were the ideal tool for this proof of concept, by allowing us to keep the patterns unevolving.

4.3 Spiral structure longevity in N-body simulations

At present the nature of galactic disc spiral structure is not well understood. Though it is generally accepted that spirals are density waves there exist two competing theories: (i) transient/recurrent spirals, and (ii) long-lived, steady-state spirals.

Recurrent spiral instabilities have been reported by Sellwood & Carlberg (1984) and Sellwood & Lin (1989) in their simulations of self gravitating discs. It was argued by Toomre & Kalnajs (1991) that these transient spiral density waves are due to the swing-amplification mechanism as first formulated by Toomre (1981).

Alternatively, the concept of quasi-stationary density waves was developed (mostly analitically) by Lin et al. (1969) and culminated in the work by Bertin et al. (1989a, b) and Lowe et al. (1994). While thought to always produce short-lived spirals, N-body simulations have been constructed to yield long-lived spiral density waves lasting for over five rotation periods, by introducing an inner Q-barrier to shield the 2:1 ILR (Thomasson et al., 1990; Elmegreen & Thomasson, 1993; Donner & Thomasson, 1994; Zhang, 1996).

Minchev et al. (2012b) studied the longevity of spiral structure in N-body Tree SPH simulations from the GalMer database (Di Matteo et al., 2007). These simulations develop strong bar and spiral structure. The full description and simulation setup can be found in Minchev et al. (2012b).

Minchev et al. (2012b)’s novel way of estimating the longevity of spirals in N-body discs consisted of following the time evolution of Fourier power spectrograms, which allowed to assessing the development of radial extent, amplitude, and patterns speed. Unlike in most previous works, where power spectrograms are typically computed over time periods of 0.5-1 Gyr, in our analyses we decreased the time window to 300 Myr and computed spectrograms every 150 Myr from $t=150$ to $t=900$ Myr (using outputs every 10 Myr).

We present the results for the GalMer gSa galaxy model in Fig. 9. The horizontal and vertical axes show galactic radius, $r$, and angular frequency, $\omega$, respectively; to convert between $\omega$ and the patter speed we use $\Omega=\omega/m$, where $m$ is the multiplicity. The series of panels show the time evolution of the $m=2$ component. At $t=300$ Myr the central bar is seen to extend from the galactic center to $r\approx 5.5$ kpc and has a frequency $\omega\approx 0.12~{}{\rm Myr^{-1}}$. A spiral wave of frequency $\omega\sim 0.05~{}{\rm Myr^{-1}}$ is also seen, extending between the bar’s CR and the disc break at that time. At later times this wave is always present, increasing in length up to $r\sim 18$ kpc at $t=750$ Myr. When animated, this two-armed feature is aways present in the spectrograms, smoothly changing from one time output in the figure to the next. We, therefore, conclude that this pattern has a lifetime $\gtrsim 600$ Myr. The smooth decrease in strength, as well as the radial extent with time, is also in agreement with the conclusion that we see the same pattern in all snapshots. In contrast, an inner $m=2$ spiral, seen mostly in the middle panel, exhibits strong variations with time. When animated, this feature appears to be driven by the interaction of the bar with the outer $m=2$ pattern, bouncing back and forth between the two on a timescale consistent with their beat frequency. In other words, every time the (faster) bar encounters the (slower) outer spiral wave, this inner wave is regenerated, speeding up to catch up with the bar ($t=450$ and $t=600$ Myr) and later on slowing down to reconnect with the outer spiral ($t=750$ Myr). Such an inner structure, connecting the bar with the dominant spiral has been reported before and has been proposed to provide an explanation for the nature of the “long” bar in our Galaxy (Athanassoula, 2005; Martinez-Valpuesta & Gerhard, 2011; Romero-Gómez et al., 2011; Athanassoula, 2012).

We estimated above that the lifetime of the two-armed outer spiral wave in our gSa model is $\gtrsim 600$ Myr. At its average rotational frequency of $\omega\sim 0.04~{}{\rm Myr^{-1}}$, this corresponds to $\gtrsim 4$ rotations, which is relatively long-lived and comparable to estimates found in previous work (e.g., Thomasson et al. 1990; Elmegreen & Thomasson 1993; Donner & Thomasson 1994; Zhang 1996). It should be noted that the model we consider here has a substantial bar, which may be related to the longevity of spirals, especially if mode coupling is present. This spiral lifetime is most likely a low limit because of the quick disc heating taking place in these simulations, due to the lack of gas infall, which, if present, would rejuvenate the spiral structure (Binney & Tremaine, 2008).

5 Constraining recent merger events

In addition to the moving groups seen in the Hipparcos velocity distribution discussed earlier, high-velocity groups of stars (i.e., $|v|\gtrsim 60$ kpc, with $v$ the residual tangential component) have also been identified in various local surveys. Because the bar and spirals have been found to only affect the $u-v$ plane at $|u|,|v|\lesssim 50$ km/s. These moving groups are thus usually attributed to accreted populations. For instance, the Arcturus stream (Williams et al., 2008) at $v=-100$ km/s, has been interpreted as originating from the debris of a disrupted satellite (Navarro et al., 2004; Helmi et al., 2006). Two other groups of stars moving at $v\sim 80$ km/s (Arifyanto & Fuchs, 2006) and $v\sim-160$ km/s (Klement et al., 2008) were assigned similar origin, based on their kinematics.

Another interpretation for these high-velocity streams was proposed by Minchev et al. (2009), showing that these can be associated with the dynamical effect on the disk by a recent merger event. The distinction here is that these stars are not accreted but belong to the disk and are grouped together in velocity because of the external perturbation.

How this works is illustrated in Fig. 10, which shows the time development of the simulated solar neighborhood velocity distribution for an axisymmetric disk with initial conditions emulating an impact following a minor merger perturbation. Contours show the $u-v$ plane and the solid line is the tangential velocity distribution.

It is evident that as time increases features get closer together because phase wrapping takes place. There is a particular time when the separation between clumps is similar to that observed, or about 20 km/s. This is shown in Fig. 11 where the model at $\sim 1.9$ Gyr after the perturbation is compared to a combination of the observational samples by Nordström et al. (2004) and Schuster et al. (2006).

This model provides an explanation for the following high-velocity streams seen in the solar neighborhood: $v\approx-60$ (HR 1614, De Silva et al. 2007)555Although HR 1614 is established to be a dissolving cluster (De Silva et al., 2007), it can nevertheless be explained by this model, as long as it is older than the time of merger event. , $v\approx-80$ (Arifyanto & Fuchs, 2006), $v\approx-100$ (Arcturus), and $v\approx-160$ (Klement et al., 2008). In addition, it predicts four new features at $v\approx-140,-120,40$ and 60 km/s.

The separation of the induced velocity streams due to this ”disk ringing” could be related to the time of the last merger impact, suggesting an event $\sim 1.9$ Gyr ago. The validity of this diagnostic was confirmed with N-body simulations (Gómez et al., 2012b) and found to be consistent with structure in the SEGUE G-dwarf sample (Gómez et al., 2012a). Gómez et al. (2013) showed that just such a disturbance from the Sagittarius dwarf galaxy can explain the recent discovery of vertical waves near the Sun in SEGUE (Widrow et al., 2012) and RAVE (Williams et al., 2013).

6 Effects of radial migration

6.1 Migration and disc thickening

Several works have previously suggested that radial migration can give rise to thick disc formation by bringing out high-velocity-dispersion stellar populations from the inner disc and the bulge. Such a scenario was used, for example, in the analytical model of Schönrich & Binney (2009b), where the authors claimed to explain the Milky Way thick- and thin-disc characteristics (both chemical and kinematical) without the need of mergers or any discrete heating processes. Similarly, the increase of disc thickness with time found in the simulation by Roškar et al. (2008) has been attributed to migration in the work by Loebman et al. (2011).

6.1.1 Migration induces disc flaring in isolated discs

A first effort to demonstrate how exactly radial migration affects disc thickening in dynamical models was done by Minchev et al. (2012a). It was shown that stellar samples arriving from the inner disc have slightly higher velocity dispersions, which will result in them being deposited at higher distances above the galactic midplane. However, the opposite effect arrises from samples arriving from the outer disc (with lower velocity dispersions). Therefore, the overall migration effect on the disc thickening is minimal throughout most of the disc extend, except in the very inner/outer parts of the disc, where only inward/outward migrators are deposited. This naturally results in disc flaring, as shown in Fig. 7 by Minchev et al. (2012a). We explained this as the conservation of vertical action as opposed to conservation of the vertical energy assumed before. Several independent groups, using different simulation techniques and setups, have confirmed that migration does not thicken the disc (Martig et al. 2014b; Vera-Ciro et al. 2014; Grand et al. 2016; Aumer et al. 2016).

This effect is illustrated in the left column of Fig. 12 for an isolated N-body Tree-SPH simulation (the GalMer gSb model) of a barred disc. The top left panel of Fig. 12 shows the changes of guiding radius, $\Delta r_{g}$, vs the initial guiding radius, $r_{g,0}$, during a time period of 1 Gyr; $r_{g}$ for each star is estimated using the values of the specific angular momentum and rotation curve as $r_{g}=L/v_{\phi}$. The percentage of stars in each contour level is given by the color bar on the right. The migration seen in the top panel is solely due to internal (secular) evolution, i.e., the effect of the central bar and spiral arms.

Next we separated migrating from non-migrating stars (in the considered time period) by applying the technique described by Minchev et al. (2012a). This consists of separating stars in a given radial bin into migrators and non-migrators as follows. Non-migrators are those particles found in the selected radial bin at both the initial and final times, while migrators are those that were not present in the bin initially but are there at the final time. We distinguish between outward and inward migrators – those initially found at radii smaller or larger than the annulus considered, respectively. This is done for radial bins over sampling the entire radial extent of the disc. See Minchev et al. (2012a) for more details.

In the second top-to-bottom left panel of Fig. 12 we plot the vertical velocity dispersion profiles of inward migrators, outward migrators and the non-migrating population, as indicated. We find that stars arriving from the inner disc are slightly hotter than the non-migrating populations but the ones arriving from the outer disc are cooler. In the third top-to-bottom left panel it can be seen that the overall effect of all migrators is negligible to the overall vertical velocity dispersion and, thus, to the disc thickening. Some cooling inside $r=3$ kpc is seen, which is related to the accumulation of inward migrating stars with small vertical actions in that region.

Finally, to quantify the changes to the disc vertical velocity dispersion resulting from migration in the given period of time, we plot the fractional changes in the bottom left panel of Fig. 12. As in Minchev et al. (2012a), we estimate these as

[TABLE]

where $\sigma_{\rm z,all}$ and $\sigma_{\rm z,non\_mig}$ are the vertical velocity dispersions for the total population and the non-migrators, respectively. Fluctuations of less than 10% are seen around $\Delta\sigma_{z}=0$. The positive slope of $\Delta\sigma_{z}$ with radius indicates an overall cooling/heating in the galactic center/disc outskirts. This is what gives rise to the flaring described by Minchev et al. (2012a). The degree of flaring cause by migration is much less than that by infalling satellites (40% vs 800%, respectively, over $\sim 4$ scale-lengths; compare Minchev et al. 2012a to Bournaud et al. 2009).

6.1.2 Migration suppresses disc flaring when infalling satellites are present

It is well known from both observations and cosmological simulations that minor mergers take place in the formation of galactic discs. The intensity of these interactions with the host disc decrease with redshift but can persist until today, as evident in the Milky Way (e.g., the Sagittarius dwarf galaxy, Ibata et al. 1994, 1995). Such interactions will have the effect of heating more the disc outskirts at any given time of the disc growth (more so at high redshift), because of the low mass density there.

Interestingly, and to complicate matters, when this more realistic scenario is considered, migration has the opposite effect on the disc vertical profile compared to the effect of an isolated galaxy – the role of outward and inward migrators is reversed in that they now cool and heat the disc, respectively. We next demonstrate this in the right column of Fig. 12, which presents a simulation by Martig et al. (2012) used for the chemo-dynamical model presented by (Minchev et al. 2013, 2014a; see §9). We consider all stars born right before the last massive merger encounters the disc at $t=1.4$ Gyr. The strong redistribution of angular momentum seen in the $r_{g,0}-\Delta r_{g}$ plane is caused both by the tidal effect of the satellite, which plunges through the galactic center, and the strong spiral structure induced in the gaseous component.

Examining the right second top-to-bottom panel, it is remarkable that inward and outward migrators during the merger have reversed roles compared to the case of an isolated disc (left column), where stars migrating inwards have positive contribution to $\sigma_{z}$ and those migrating outward cool the disc.

The net effect of migrators during the merger can be seen in the third top-to-bottom right panel of Fig. 12. The overall contribution to the vertical velocity dispersion from the migrating stars during the merger is negative, in the sense that it is lower than that of the stars which did not migrate. We emphasize that we only considered stars born before the merger took pace, therefore, the effect seen is not related to the accreted population.

In the bottom right panel of Fig. 12 we can see a negative slope for $\Delta\sigma_{z}$ vs $r$, in contrast to the left panel, indicating that the disc flaring induced by the merging satellite is being counteracted by outward migrating stars. A decrease in $\sigma_{z}$ of up to 30% in the disc outskirts is seen.

Fig. 12 showed that during a massive merger sinking deep into the disc center, migrators cool the outer disc, thus working against disc flaring. This is related to the stronger effect of mergers on the outer disc, owing to the exponential decrease in the disc surface density. Stars arriving from the inner parts during a merger, would therefore be cooler than the rest of the population. Note that flaring for a given mono-age population, nevertheless, still results from infalling satellites and that gives rise to the formation of a thick disc, as discussed next.

7 Formation of galactic thick discs by the flaring of mono-age populations

Stellar disc density decomposition into thinner and thicker components in external edge-on galaxies find that thicker disc components have larger scale-lengths than the thin discs (e.g., Yoachim & Dalcanton 2006; Pohlen et al. 2007; Comerón et al. 2012). While this is consistent with results for the Milky Way when similar morphologically (or structural) definition for the thick disc is used (e.g., Robin et al. 1996; Ojha 2001; Jurić et al. 2008), it is in contradiction with the more centrally concentrated older or [$\alpha$/Fe]-enhanced stellar populations (e.g., Bensby et al. 2011; Cheng et al. 2012; Bovy et al. 2012). This apparent discrepancy may be related to the different definition of thick discs - morphological decomposition or separation in chemistry.

Additionally, while no flaring is observed in external edge-on discs (van der Kruit & Searle, 1982; de Grijs, 1998; Comerón et al., 2011), numerical simulations suggest that flaring cannot be avoided due to a range of different dynamical effects. The largest source is most likely satellite-disc interactions (e.g., Bournaud et al. 2009; Villalobos & Helmi 2008), which have been found to increase an initially constant scale-height by up to a factor of $\sim 10$ in 3-4 disc scale-lengths. We showed in the left row of Fig. 12 that purely secular evolution (in the absence of external perturbations) also causes flared discs, due to the redistribution of disc angular momentum. Other sources of disc flaring include misaligned gas infall (Scannapieco et al., 2009) and reorientation of the disc rotation axis (Aumer & White, 2013).

Minchev et al. (2015) studied the formation of thick discs using two suites of simulations of galactic disc formation, one of which we present here. This is a full cosmological zoom-in hydro simulation, using initial conditions from one of the Aquarius Project haloes (Springel et al., 2008; Scannapieco et al., 2009). The technique used here is Tree-PM SPH with 300 pc spatial and $4.4\times 10^{5}$ M*⊙* mass resolution. This is a non-barred galaxy with stellar mass $5.5\times 10^{10}M_{\odot}$ and a disc scale-length $h_{d}=4$ kpc. Further details about this simulation can be found in Aumer et al. (2013), their model Aq-D-5.

The simulation forms an initial central component during an early epoch of violent merger activity. Gas-rich mergers supply the initial reservoir of gas at high redshift and merger activity decreases with redshift, similarly to what is expected for the Milky Way. This inside-out disc formation results in a centrally concentrated older stellar population. The general formation and evolutionary behavior of this model is similar to many recent simulations in the cosmological context (e.g., Brook et al. 2012; Stinson et al. 2013).

We decomposed the stellar discs into mono-age populations, i.e., narrow bins of age, where we used $\Delta$age=1 Gyr. It was found that single exponentials provided good fits to the column density in the vertical direction for all age bins and at all radii, in agreement with Martig et al. (2014a). In contrast, to properly fit the vertical density of the total stellar population required the sum of two exponentials.

In the top panel of Fig. 13 we plot the scale-height variation with galactocentric radius, $r$, in the region 1-5 disc scale-lengths, $h_{d}$. Both the radius and scale-height, $h_{z}$, are in units of $h_{d}$. It can be seen that significant flaring is present, which increases for older coeval populations.

In contrast to the flaring found for all but the youngest mono-age populations, the thin and thick disc decomposition of the total stellar population including all ages, results in no apparent flaring. This is shown by the triangle and square symbols overlaid in the top panel of Fig. 13.

What is the reason for the flaring of mono-age discs? In numerical simulations flaring is expected to result from a number of mechanisms related to galactic evolution in a cosmological context (e.g., Bournaud et al. 2009; Kazantzidis et al. 2008; Villalobos & Helmi 2008; Aumer & White 2013). Even in the absence of environmental effects, flaring is unavoidable due to secular evolution alone (radial migration caused by spirals and/or a central bar, Minchev et al. 2012a). It should be stressed here that, while migration flares discs in the lack of external perturbations, during satellite-disc interactions it works against disc flaring (Minchev et al. 2014a, see §6.1). Yet, this is not sufficient to completely suppress the flaring induced by orbiting satellites, as evident from the top row of Fig. 13. This suggests that external effects are much more important for the disc flaring in this simulation. Because the mass and intensity of orbiting satellites generally decreases with decreasing redshift, so does the flaring induced. It can be expected that at a certain time secular evolution takes over the effect of external perturbations.666Minchev et al. (2014b) suggested that the time at which internal evolution takes over can also be inferred from the shape of the [$\alpha$/Fe]-velocity dispersion relation of narrow metallicity samples.

What is the reason for the lack of flaring in the total disc population? In an inside-out formation scenario, the outer disc edge, where flaring is induced, moves progressively from smaller to larger radii because of the continuous formation of new stars in disc subpopulations of increasing scale-length. At the same time the frequency and masses of perturbing satellites decreases. Because of the inside-out disc growth, which results in more centrally concentrated older samples (see Fig. 13, middle panel), the younger the stellar population, the further out it dominates in terms of stellar mass. The geometrically defined thick disc, therefore, results from the imbedded flares of different coeval populations, as seen in the top panel of Fig. 13.

The bottom panel of Fig. 13 shows that a geometrical thick disc is expected to have a negative age gradient. For this particular simulation the mean age decreases from $\sim 10.5$ to $\sim 6$ Gyr in four disc scale-lengths. Such an age drop of mean stellar age at high distances, $|z|$, from the disc midplane explains the inversion in $[\alpha$/Fe] gradients with increasing mean $z$ found by Anders et al. (2014) in APOGEE data.

8 Chemo-dynamical modeling of the Milky Way

So far we have focused on the dynamics of the Milky Way disc, which tells us mostly about its current state. To be able to go back in time and infer the Milky Way evolutionary history, however, we need to include both stellar chemical and age information.

8.1 Summary of Milky Way chemo-dynamical evolution modeling techniques

A major consideration in a disc chemo-dynamical model is taking into account the effect of radial migration, i.e., the fact that stars end up away from their birth places. Below we briefly summarize models which include radial migration.

$\bullet$ Semi-analytical models tuned to fit the local metallicity distribution, velocity dispersion, and chemical gradients, etc., today (e.g., Schönrich & Binney 2009a; Kubryk et al. 2015) or Extended distribution functions (Sanders & Binney, 2015):

$-$ Easy to vary parameters

$-$ Provide good description of the disc chemo-kinematic state today

$-$ Typically not concerned with the Milky Way past history

$-$ Time and spatial variations of migration efficiency due to dynamics resulting from non-axisymmetric disc structure is hard to take into account.

$\bullet$ Fully self-consistent cosmological simulations (e.g., Kawata and Gibson 2003; Scannapieco et al. 2005; Kobayashi and Nakasato 2011; Brook et al. 2012):

$-$ Dynamics self-consistent in a cosmological context

$-$ Can learn about disc formation and evolution

$-$ Not much control over final chemo-kinematic state

$-$ Problems with SFH and chemical enrichment due to unknown subgrid physics

$-$ Much larger computational times needed if chemical enrichment included.

$\bullet$ Hybrid technique using simulation in a cosmological context + a classical (semi-analytical) chemical evolution model (Minchev et al. 2013, hereafter MCM13):

$-$ Avoids problems with SFH and chemical enrichment in fully self-consistent models

$-$ Can learn about disc formation and evolution

$-$ Not easy to get Milky Way-like final states.

Below we focus on the results of the latter model.

9 The MCM13 hybrid chemo-dynamical model

To properly model the Milky Way it is crucial to be consistent with observational constraints at redshift $z=0$, for example, a flat rotation curve, a small bulge, a central bar of an intermediate size, gas to total disc mass ratio of $\sim 0.14$ at the solar vicinity, and local disc velocity dispersions close to the observed ones.

It is clear that cosmological simulations would be the natural framework for a state-of-the-art chemo-dynamical study of the Milky Way. Unfortunately, as discussed by MCM13, a number of star formation and chemical enrichment problems still exist in fully self-consistent simulations. We have, therefore, resorted to the next best thing – a high-resolution simulation in the cosmological context coupled with a pure chemical evolution model.

The simulation used is part of a suite of numerical experiments presented by Martig et al. (2012), where the authors studied the evolution of 33 simulated galaxies from $z=5$ to $z=0$ using the zoom-in technique described by Martig et al. (2009). This technique consists of extracting merger and accretion histories for a given halo in a $\Lambda$-CDM cosmological simulation and then re-simulating at much higher resolution (150 pc spatial, and 104-5 M*⊙* mass resolution).

Originally, our galaxy has a rotational velocity at the solar radius of 210 km/s and a scale-length of $\sim 5$ kpc. To match the Milky Way in terms of dynamics, at the end of the simulation we downscale the disc radius by a factor of 1.67 and adjust the rotational velocity at the solar radius to be 220 km/s, which affects the mass of each particle according to the relation $GM\sim v^{2}r$, where $G$ is the gravitational constant. This places the bar’s CR and 2:1 OLR at $\sim 4.7$ and $\sim 7.5$ kpc, respectively, consistent with a number of studies (e.g., Dehnen 2000; Minchev et al. 2007, 2010; see §2). At the same time the disc scale-length, measured from particles of all ages in the range $3<r<15$ kpc, becomes $\sim 3$ kpc, in close agreement with expectations in the Milky Way.

9.1 Birth radius distribution of solar vicinity stars

To illustrate the migration efficiency in the MCM13 model, in the left panel of Fig. 15 we show the birth radii of stars ending up in a solar neighborhood-like location ($7<r<9$ kpc, green shaded strip, and $|z|<3$ kpc) after 11.2 Gyr of evolution. The solid black line shows the total population, which peaks close to $r_{0}=6$ kpc due to radial migration. The entire sample is also divided into six age-groups, shown by the curves of different colors and line-styles. The strongest effect from radial migration is found for the oldest stars (red curve), whose distribution has a maximum at $r\approx 3$ kpc, or inside the bar’s CR (at $\sim 4.7$ kpc). Note that locally born stars of all ages can be found in the solar neighborhood. A relatively smooth transition of the peak, from older to younger groups of stars, is observed; this is expected, since even for a constant migration efficiency older stars would be exposed longer to perturbations. While a wide range of birth radii is seen for all age groups, the majority of youngest stars were born at, or close to the solar neighborhood bin.

9.2 The age-metallicity relation

In the middle panel of Fig. 15 shows stellar density contours of the resulting age-metallicity relation for the “solar” cylinder shown in the left panel. Overlaid on top of the contours is our input chemical evolution model for different initial radii, as indicated, giving an insight into the origin of stars found in this localized region. The excess of stars above the local curve (cyan) is due to migrators coming from the inner disc. Similarly, contours below the local curve result from stars originating in the outer disc.

The mean metallicity, shown by the pink dashed curve, is found to follow closely the in-situ born population (cyan curve). Some flattening is observed, mostly for ages $\gtrsim$9 Gyr, however the final distribution is by no means flat. The reason for this minor effect on the local metallicity gradient, despite the strong migration, is the fact that at the Sun’s intermediate distance from the Galactic center the change in metallicity arising from stars migrating from the inner regions is mostly compensated for by stars arriving from the outer disc.

9.3 The birth place of the Sun

With the care taken into defining a proper solar radius in a simulation with Milky Way characteristics we can made an estimate for the solar birth radius. The right panel of Fig. 15 displays density contours of the $r_{0}$-[Fe/H] plane for all local stars. The cyan curve shows the input solar-age metallicity gradient. This is taken to be 4.6 Gyr look-back time, consistent with, e.g., Bonanno et al. 2002; Christensen-Dalsgaard 2009; Houdek & Gough 2011. Assuming an error of $\pm 0.06$ dex in [Fe/H], we find a possible Sun birth radius in the range $5.0<r_{\odot,birth}<6.8$ kpc (where the horizontal green, transparent strip meets the cyan curve)777This estimate is also dependent on the solar metallicity error assumed, e.g., MCM13 used $\rm[Fe/H]_{err}\pm 0.1$ dex obtaining $4.4<r_{\odot,birth}<7.7$ kpc.. This result is dependent on the migration efficiency in the simulation and is in good agreement with the estimate of $6.6\pm 0.9$ kpc by Wielen et al. (1996).

9.4 Constraining the entire Milky Way merger history

Using data from the RAVE and SEGUE surveys, Minchev et al. (2014b) presented a previously unknown chemo-kinematic relation, where the variation of stellar velocity dispersion with [Mg/Fe] ratios is reversed beyond 0.4 dex (Fig. 15, panel d). This is unexpected because stars with higher [Mg/Fe] ratios, at a given radius, should be older and, thus, should possess larger random energies. An excellent match by the MCM13 chemo-dynamical model (Fig. 15, panel e) revealed that this new relation offers a unique way to recover the Milky Way merger history, where the peak of each metallicity subpopulation gives away the time and mass of the perturber and the decline on the right is related to older cooler stars migrating from the inner disk. This observational relation has been also found in Gaia-ESO data (Guiglion et al., 2015) and is currently being investigated with APOGEE data.

By using our range of models, we will quantify global structures in the Milky Way disk caused by mergers as a function of satellite mass, inclination, and compactness. We may discover, for example, that bar formation is related to a previous merger, tidal deformation, or gas accretion event. We may also find evidence for kinematic and metallicity variations in the disk directly related to episodes of gas accretion tied to merger events.

9.5 Comparison with observations

The MCM13 model has been found to match simultaneously the following observational constraints, some of which are reproduced in Fig. 16:

$\bullet$ The disc morphology, scale-length and rotation curve today (see MCM13)

$\bullet$ The local age-velocity dispersion relation (Sharma et al. 2014)

$\bullet$ The more centrally concentrated [$\alpha$/Fe]-enhanced (old) disc (Bensby et al. 2011; Bovy et al. 2012; see Minchev, Chiappini, and Martig 2014b, Fig.11)

$\bullet$ The distribution of scale-heights for mono-abundance subpopulations found in SEGUE G-dwarfs (MCM13, Fig.13)

$\bullet$ The reversal of the radial [$\alpha$/Fe] and metallicity gradients (e.g., APOGEE - Anders et al. 2014), when sample distance from the disc mid-plane is increased (Minchev et al. 2014b)

$\bullet$ The MDF for stellar samples at different distance from the disc midplane (Fig. 16, a, b)

$\bullet$ The metallicity variation with vertical distance from the plane (Fig. 16, c)

$\bullet$ The inversion in velocity dispersion relation in RAVE and SEGUE (Fig. 16, d, e; Minchev et al. 2014; Guiglion et al. 2015)

$\bullet$ The [Fe/H]-[$\alpha$/Fe] plane (Fig. 16, f)

$\bullet$ The age-[$\alpha$/Fe] relation (Fig. 2 by Minchev et al. (2016))

$\bullet$ The age-[Fe/H] relation (Fig. 2 by Minchev et al. (2016))

$\bullet$ The flaring of mono-abundance populations (assuming similarity to mono-age populations) found by Bovy et al. (2015) and predicted earlier by Minchev et al. (2015), where Model 1 in the latter paper presents the same galaxy as the one used for the MCM13 chemodynamical model.

It is impressive that this simple match between a simulation in the cosmological context with Milky Way characteristics and a classical Milky Way chemical evolution model, where for both the SFH is such that the disc grows inside-out, was able to account for so many chemo-kinematic relations. This model needs to be improved to account also for the gap in the [$\alpha$/Fe]-[Fe/H] relation seen for a long time in high-resolution spectroscopic data in the solar vicinity (e.g., Fuhrmann 2004; Bensby et al. 2013; Adibekyan et al. 2013) and more recently found in a range of Galactic radii using APOGEE data (Anders et al., 2014; Nidever et al., 2014). This separation into low- and high-[$\alpha$/Fe] sequences (sometime referred to as the thin and thick discs)888Note that this is the chemical definition of thin and thick disks, which is different than a definition using morphology (or geometry). See discussion in §7. is most likely the result of a gap in the SFH at high redshift, as first suggested in the Two-Infall Model by Chiappini et al. (1997).

10 Discussion and conclusions

In this work we reviewed some methods used to constrain the dynamics of the Milky Way disc, in particular, its central bar, spiral structure, and the formation of the thick disc. While kinematics alone can be used to constrain the current dynamical state of the disc, chemical information and accurate age estimates are needed to understanding the Milky Way formation, with the help of detailed chemo-dynamical evolution modeling.

We emphasize the importance of taking into account today’s Galactic disc morphology in the construction of chemo-dynamical models, which ensures that resonances are placed at the appropriate Galactic radii. Although there is no guarantee that a model displaying similar final state to that of the Milky Way had the same evolutionary history, taking as a constraint the only phase-space snapshot we have of the Galaxy is a must.

The field of Galactic Archaeology will soon be transformed by the Gaia mission, which will provide precise positions, proper motions, radial velocities, and stellar parameters for up to 1 billion stars of the Milky Way and will map distances out to 15 kpc from the Sun. This can be contrasted to a total of a few million Milky Way stars observed to date, only a fraction of which lie farther than  2-3 kpc from the Sun. Stars with accurate 6D kinematics are currently available only within  0.1 kpc from the Sun (yellow dot in Fig. 17) provided by Hipparcos and GCS. This will soon change already with the first and second Gala data releases expected in 2016 and 2017, respectively, which will complement existing spectroscopic data (e.g., RAVE-Gaia and APOGEE-Gaia synergies led by the Leibniz-Institut für Astrophysik Potsdam AIP) and will dramatically increase the disc coverage of stars with precise kinematics (see Fig. 17).

Availability of accurate ages is also very important to make progress in the field of Galactic Archaeology by breaking degeneracies in chemo-dynamical models. This has recently become evident with the unexpected results of Chiappini et al. (2015) and Martig et al. (2015), who used CoRoT (Baglin et al., 2006) and Kepler (Gilliland et al., 2010) asteroseismic ages, respectively, combined with APOGEE chemical information, to show the existence of significantly young high-[$\alpha$/Fe] stars (but see also Yong et al. 2016). This is a largely unexpected result for chemical evolution modeling, because [$\alpha$/Fe] has been thought to always be a good proxy for age. Stellar ages for much larger samples and broader disc coverage are expected in the very near future from Kepler-2 and the Gaia mission. These will help break degeneracies and refine chemo-dynamical models, thus bringing us a step closer to understanding the formation of our Galaxy.