The Radial Distribution of Mono-metallicity Populations in the Galactic Disk as Evidence for Two-phase Disk Formation

Recently, the spectra of some 70,000 red giant stars from SDSS-III/APOGEE DR12 (Majewski et al. 2015) have been obtained in the H band, where the dust effects are not important, providing the element chartography of the Milky Way (MW) over an unprecedented large volume and including the Galactic plane for the first time. These very recent data opened up the possibility of studying the processes occuring along the Milky Way (MW) assembly and evolution in more detail through the imprints they have left on the stellar space distributions.

It has long been known that the stellar populations of spiral galaxies have two main components: a dynamically hot spheroid, and a cold disk. A third component, originally detected in the Milky Way through stellar counts (Gilmore & Reid 1983), has been detected as a ubiquitous excess of red flux at large galactic latitudes in external spiral galaxies (Dalcanton & Bernstein 2002; Yoachim & Dalcanton 2006). Detailed kinematic and chemical studies (e.g., Fuhrmann 1998; Bensby et al. 2003; Soubiran et al. 2003; Bensby et al. 2005; Reddy et al. 2006) confirmed that this excess was due to a distinct component, the thick disk. Later on, statistics studies of thick disks in the local Universe suggested that thick disks are ubiquous in galaxies (Comerón et al. 2011).

Observations (Gilmore et al. 1989, 1995; Fuhrmann 1998; Bensby et al. 2003; Soubiran et al. 2003; Reddy et al. 2006; Ivezić et al. 2012) indicate that the differences between the thin and thick disks involve four parameters. (i) The shape: the vertical scale-length is smaller for the thin disk than for the thick disk (e.g., Dalcanton & Bernstein 2002; Yoachim & Dalcanton 2006). (ii) The kinematics: the thin disk is colder in all velocity components, and it is more rotationally supported than the thick disk (Soubiran et al. 2003). (iii) The age: thick-disk stars are older on average (Gilmore et al. 1995). (iv) The metallicity: thin-disk stars are more metal rich than thick-disk stars (Fuhrmann 1998). Finally, (v) the $[\alpha$/Fe]: the thick disk has enhanced $\alpha$-elements compared to thin-disk populations of similar [Fe/H] abundances (e.g., Fuhrmann 1998; Bensby et al.2003; Reddy et al. 2006), suggesting shorter star formation timescales (Fuhrmann 1998; Ruchti et al. 2010).

These findings have been confirmed by observational studies made within recent or ongoing spectroscopic surveys, such as the RAdial Velocity Experiment (RAVE; Steinmetz et al. 2006), the Sloan Extension for Galactic Understanding and Exploration (SEGUE; Yanny et al. 2009), Gaia-ESO (Gilmore et al. 2012), the Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST; Cui et al. 2012), the Apache Point Observatory Galactic Evolution Experiment (APOGEE; Majewski et al. 2015), and Gaia (Gaia Collaboration et al. 2016). The full explotation of these and other planned surveys such as the GALactic Archeology with Hermes (GALAH; Martell et al. 2017) or WEAVE make this decade a Golden Age for Astroarcheology, where a dramatic advancement in the understanding of the Galaxy is expected.

While authors agree on these distinct properties of thin and thick disks, different methods for assigning a given star to the thin or thick populations can be found in literature. A kinematic classification scheme is used extensively in hydrodynamical simulations (e.g., Abadi et al. 2003; Doménech-Moral et al. 2012; Domínguez-Tenreiro et al. 2015) and in observations (Bensby et al. 2003; Soubiran et al. 2003; Kordopatis et al. 2011; Boeche et al. 2013). Other classifications rest on chemical properties such as the $[\alpha$/Fe] versus [Fe/H] (e.g., Adibekyan et al. 2012; Bovy et al. 2012b; Ramírez et al. 2013; Bensby et al. 2014; Recio-Blanco et al. 2014; Hayden et al. 2015; Martell et al. 2017, and references therein) or stellar age (see for example Kubryk et al. 2015). Haywood et al. (2013) use a classification based on the [Si/Fe] versus age relation, where the separation comes from a knee in the two-slope behavior of this relation that these authors explain in terms of a particular shape for the star formation rate history (SFRH), see Snaith et al. (2014). The two-slope behavior for [α/Fe]—age has been extended to local early-type galaxies by Walcher et al.(2015).

Much effort has so far been devoted to understand the origin of the thick disk. Basically, two different scenarios exist to explain the thick-disk emergence (see reviews in Freeman 1987; Gilmore et al. 1989; Freeman & Bland-Hawthorn 2002; Brook et al. 2004; van der Kruit & Freeman 2011; Ivezić et al. 2012; Feltzing & Chiba 2013; Minchev 2017). The first scenario links thick disks to violent formation processes from turbulent gas before the thin disk forms. The second scenario assumes a preexisting thin disk that is dynamically heated along secular evolution, with stellar migration (first proposed by Roškar et al.2008; Schönrich & Binney 2009; Loebman et al. 2011) or stellar satellite acretion (Abadi et al. 2003) as the main paths.

The first scenarios are supported in particular from cosmological simulations by Brook et al. (2004, 2012b); Stinson et al. (2013a), and Bird et al. (2013), see also Scannapieco et al. (2011). Brook et al. (2004) identified a period of fast merging at high z where chemically classified thick-disk stars form kinematically hot, before the formation of thin-disk stars. Brook et al. (2012b) confirm the previous scenario and add radial mixing to explain radial changes in thick disks. They find that radial mixing plays a comparatively minor role in disk thickening and heating. Stinson et al. (2013a) extended these results and incorporated the metallicity disk structure in their analysis.

Scannapieco et al. (2011) find that young (old) stars define thin-disk (thick-disk) structures, with kinematic and chemical properties reminescent of those of observed thin (thick) disks.

Bird et al. (2013) analyze the radial profiles of surface mass density and kinematics of mono-age stellar populations, aided by a kinematical separation of their simulated galaxy components. They find that most of the kinematically defined thick-disk stars form within the first 4 Gyr after the Big Bang, in a dynamically violent period. Their (kinematically defined) thin-disk stars formed that early would be members of a chemically thick-disk population, recovering in this way the results of Brook et al. (2012b) that thick-disk stars are born early and hot, while thin-disk stars are not. They also agree on the role radial mixing has on disk heating and thickening.

The ability of radially migrated inner stars to cause thick dyanmically hot extended disks is carefully analyzed in Minchev et al. (2012) using cosmological and pre-prepared simulations. In both cases, the authors also conclude that this process does little for disk heating and thickening, but can change the spatial chemical structure of the disks, as shown by Kubryk et al. (2013). These changes, however, do not wash out the chemical imprints left by disk formation according to the pre-prepared simulation results of Curir et al. (2014). Similar results on the effects of radial migration along secular evolution have been presented by Vera-Ciro et al. (2014).

Some observational results suggest or are consistent with the previous scenario. For example, Haywood et al. (2013), Snaith et al. (2014) and Haywood et al. (2015) conclude from an analysis of the chemical properties of a set of solar-vicinity FKG stars (Adibekyan et al. 2012), combined with a careful determination of their age, that the thick-disk stars have formed at early times out of a turbulent gas. This sets the chemical conditions for a latter thin-disk formation in a more quiescent situation.

Different results from different surveys point in this direction. For example, Kordopatis et al. (2015) found in their analyses of data from the Gaia-ESO survey that the mixing of metals in the young Galaxy (e.g., turbulent gas in disks or radial stellar migration) was more efficient at early times before the (current) thin disk started forming (see also Mikolaitis et al. 2014). Other authors also prefer this scenario on the basis of their results on radial and/or vertical abundance gradientes from different surveys (see, for example Cheng et al. 2012; Xiang et al. 2015, from RAVE and LAMOST data, respectively) or consider this possibility (Boeche et al. 2014). Comerón et al. (2015) present the first Integral Field Unit spectroscopy of an edge-on galaxy, ESO 533-4, with enough depth and quality to study the thick disk. Even if not conclusive, their results suggest that the thick disk of ESO 533-4 formed in a relatively short event before the thin disk.

As for the second scenario, authors highlight different processes acting on a preexisting thin disk as the origin of the thick disk: radial migration (Roškar et al. 2008; Schönrich & Binney 2009; Loebman et al. 2011), stars accreted through satellites and dragged into the plane of a preexisting disk (Abadi et al. 2003; Meza et al. 2005), heating of a preexisting disk by satellite accretion (Quinn et al. 1993), which was later on confirmed by other authors (see, e.g., Hayashi & Chiba 2006; Kazantzidis et al. 2008; Villalobos & Helmi 2008; Bekki & Tsujimoto 2011; Qu et al. 2011), more recently by Moetazedian & Just (2016) and Ruiz-Lara et al. (2016). The secular thickening of a disk embedded in a fluctuating potential has also been investigated (Fouvry et al. 2016a, 2016b), as well as the effects of early massive clump formation in unstable gas-rich disks (e.g., Noguchi 1999; Bournaud et al. 2007; Agertz et al. 2009; Ceverino et al. 2010), see however Buck et al. (2017) or the popping of stellar clusters (Kroupa 2002; Assmann et al. 2011) on disk thickening. These two last scenarios link thickening to formation processes rather than to the heating of a preexisting thin disk.

We see that there is currently a lively debate on the origin of the thick disk. The metallicity cartography of the Galactic disk provided by the APOGEE project contributes new very interesting possibilities to this debate.

Of particular interest for studies of the thin and thick disk are mono-metallicity populations (MMPs, i.e., stars in narrow bins in [Fe/H] within wider $[\alpha$/Fe] ranges). Analyzing the red giant sample from SDSS-III/APOGEE DR12, Bovy et al. (2016) present a careful determination of their radial structure, and show that radial MMPs distributions show a bimodal behavior according to their $[\alpha$/Fe] content: (i) MMPs with enhanced $[\alpha$/Fe] show disk radial surface densities, ${\rm{\Sigma }}$(${R}_{\mathrm{cyl}}$), that are well described by single exponential distributions, with a unique scale length regardless of the MMP, confirming with a much better statistics and radial coverage previuos results by Bovy et al. (2012c) in the SDSS/SEGUE survey; see also Kordopatis et al. (2012). (ii) MMPs with low-[$\alpha$/Fe] show a more complex behavior, with their respective MMP radial surface densities exhibiting continuously varying shapes, more centrally concentrated as [Fe/H] increases, whose overall envelope is an exponential disk.

Taking enhanced and low-[α/Fe] stellar populations as thick- and thin-disk populations, respectively, Bovy et al. (2016) maintain that the results above cast doubts on the classical thin-disk thick-disk dichotomy because of the complexity of thin disk splitting into MMPs with different surface densities. This issue has not yet received a complete answer so far (but see Haywood et al. 2016; Minchev et al. 2017), in particular when a cosmological context is envisaged. The aim of this paper is to analyze this question within the wider context of disk formation inside non-isolated halos, more particularly whithin the two-phase assembly scenario of halo formation. Analytical models as well as N-body simulations show that two different phases can be distinguished along the halo mass assembly, as first proposed by Wechsler et al. (2002), Zhao et al. (2003), Brook et al. (2005), Salvador-Solé et al. (2005), see also Griffen et al. (2016) for recent results. These are (i) first a violent rapid phase with high mass aggregation rates, resulting from collapse-like events in the Cosmic Web environment, implying high merger rates, and (ii) later on a slow phase with lower mass aggregation rates. Small-box hydrodynamical simulations (Domínguez-Tenreiro et al. 2006) as well as larger box simulations (Oser et al. 2010; Domínguez-Tenreiro et al. 2011) confirmed this scenario and its implications on a possible scenario for thick-disk formation (Brook et al. 2004), elliptical properties at low redshift (Cook et al. 2009), and on classical bulges (Obreja et al. 2013).

In this paper we employ hydrodynamical simulations in a cosmological context in order to address the question to which extent the two-phase halo assembly scenario can determine the detailed MMP disk chartography that is now available in terms of differentiating the thin versus the thick disk. To answer to this question, we need to determine the physical basis of the thin- versus thick-disk differentiation within the two-phase scenario. This issue is addressed in depth in this work.

The paper is organized as follows. In Section 2 we describe the codes and the simulated galaxies. The occurrence of a two-phase mass assembly for these galaxies is examined in Section 3, where the separating event between the phases is identified as halo virialization, defining a timescale ${t}_{\mathrm{vir}}$. In this section we propose a classification scheme for stars into thick and thin populations based on ${t}_{\mathrm{vir}}$, and in Section 4 we study its relationship with a chemical classification based on age. In Section 5 we confirm that this classificaction scheme leads indeed to thin- and thick-disk populations whose properties differ in agreement with observational data. The radial MMP distributions for the disks of simulated galaxies are analyzed in Section 6, where we show that we recover the results of Bovy et al. (2016). To decipher the physical processes underlying this behavior, the stellar birth places are linked to their birth times and metallicity in Section 7, and they are analyzed in terms of the gas metallicity structure before and after ${t}_{\mathrm{vir}}$. With these findings in mind, we return to the two-phase galaxy assembly in Section 8 and analyze the effects each phase has on the physical conditions for star formation, gas metallicity gradient maintenance and removal, and the matching of the MMP distributions at stellar birth to those observed at z = 0. The summary, discussion, and conclusions are presented in Section 9.

Accurate conservation of angular momentum, a detailed implementation of chemical evolution, and the effects of discrete energy injection by stellar physics are fundamental issues in cosmological hydrodynamical simulations. Comparisons of results using different codes is advisable because we search for effects coming from a generic and fundamental level of physical description. In view of these considerations, in this paper we present results of simulations that were run with two different SPH codes: P-DEVA, and GASOLINE.

P-DEVA (Martínez-Serrano et al. 2008) is an entropy-conserving AP3M-SPH code in which angular momentum conservation was the main concern in its design. Chemical evolution implementation makes use of the ${Q}_{\mathrm{ij}}$ formalism (Talbot & Arnett 1973), which relates each nucleosynthetic product to all its different sources. Ejecta from supernovae (SNe) as well as from low- and intermediate-mass stars have been taken into account. The stellar evolution data described in Gavilán et al. (2005) have been used for low- and intermediate-mass stars, and those in Woosley & Weaver (1995) for high-mass stars. Supernova Ia (SNIa) rates were computed according to Ruiz-Lapuente et al. (2000) and their element production according to Iwamoto et al. (1999). We have considered the evolution of the following elements: H, He⁴, C¹², C¹³, N¹⁴, O¹⁶, Ne²⁰, Mg²⁴, Si²⁸, S³², Ca⁴⁰, and Fe⁵⁶. Stellar feedback is implicitly implemented through (inefficient) star formation parameters, as discussed in Agertz et al. (2011).

In GASOLINE (Wadsley et al. 2004), SNe feedback is implemented using the blastwave formalism (Stinson et al. 2006). In the MaGGIC runs (Brook et al. 2012a) feedback from massive stars is also taken into account (Stinson et al. 2013b). The chemical evolution implementation follows Raiteri et al. (1996), and we track nine elements: H, He⁴, C¹², N¹⁴, O¹⁶, Ne²⁰, Mg²⁴, Si²⁸, and Fe⁵⁶. Metals are ejected from SNII, SNIa, and the stellar winds driven from asymptotic giant branch (AGB) stars, using literature yields for SNIa (Nomoto et al. 1997). Each SN releases ${10}^{51}$ erg, and the cooling is delayed for 4 Myr for the neighboring gas particles. Ejected mass, metals, and thermal energy from SNe are distributed to the nearest-neighbor gas particles using the smoothing kernel (Stinson et al. 2006).

In both codes the star formation recipe follows a Kennincutt-Schmidt-like law with a given density threshold, ${\rho }_{* }$, and star formation efficiency ${c}_{* }$. In both codes it is assumed that stars more massive than 8 ${M}_{\odot }$ produce type II supernovae (SNII), and the yields of Woosley & Weaver (1995), with the Fe ejecta divided by two according to Timmes et al. (1995), have been adopted. In addition, stars take and retain their progenitor gas particle abundances. The abundance increments synthetyzed by a given stellar particle within an integration timestep are transferred to its nearest gas particle, allowing stars to share their yields with the surrounding interstellar medium.

In both codes, abundance diffusion within the gas is included based on (unresolved) turbulent mixing (see, for example, Scalo & Elmegreen 2004), ensuring that gas abundance increment reaches distances farther away from the stellar position and explosion sites. This metal mixing tends to equilibrate the element content of spatially close gas particles.

A zoom-in simulation technique has been used, with initial baryon particle masses $\delta {M}_{\mathrm{bar}}$, a minimum smoothing length ${h}_{\mathrm{soft}}$, and a periodic box length ${L}_{\mathrm{box}}$ as given in Table 1, where we also list the star formation efficiency c_* and density threshold ${\rho }_{* }$ parameters. We explore two P-DEVA galaxies (HD-5004A and LD-5101A) and one GASOLINE—MaGICC galaxy (g1536-L*). Their stellar and gaseous masses are given in Table 1. The main difference between the two P-DEVA galaxies is that by construction, HD-5004A forms and evolves in a dense environment, while the environment of LD-5101A has a lower density. For the g1536-${{\rm{L}}}^{* }$ galaxy, the halo to be zoomed-in has been chosen to have low merging activity.

Table 1.  The Initial Mass of Gas Particles ($\delta {M}_{\mathrm{bar}}$, it is Fixed for P-DEVA Runs), Minimum SPH Smoothing Length (${h}_{\mathrm{soft}}$), Density Threshold (${\rho }_{* }$), Star Formation Efficiency (${c}_{* }$), Initial Mass Function (IMF), Stellar (${M}_{* }$) and Gas (${M}_{\mathrm{gas}}$) Mass of the Simulated Galaxies at z = 0, Universe Age used to Separate the Thick from the Thin Stellar Populations (${t}_{\mathrm{chem}}$), and Periodic Box Length (${L}_{\mathrm{box}}$)

Object	$\delta {M}_{\mathrm{bar}}$	${h}_{\mathrm{soft}}$	${\rho }_{* }$	${c}_{* }$	IMF	${M}_{* }$	${M}_{\mathrm{gas}}$	${t}_{\mathrm{chem}}$	${L}_{\mathrm{box}}$
	${10}^{5}$ ${M}_{\odot }$	${{\rm{h}}}^{-1}$ kpc	${\mathrm{cm}}^{-3}$	%		${10}^{10}$ ${M}_{\odot }$	${10}^{10}$ ${M}_{\odot }$	Gyr	Mpc
g1536-L	1.90	0.15	9.4	3.3	Chab03	2.32	1.97	8.0	34
HD-5004A	3.94	0.20	6.0	1.0	Salp55	3.26	0.67	7.7	10
HD-5101A	3.79	0.20	12.0	0.8	Salp55	1.29	0.33	7.2	10

Download table as:  ASCIITypeset image

These P-DEVA galaxies have been studied by Doménech-Moral et al. (2012) at z = 0, who analyzed some of their thin- versus thick-disk population properties, GASOLINE galaxies by Brook et al. (2012b), Stinson et al. (2013a), and Obreja et al. (2014). In addition, Obreja et al. (2013), who addressed bulge formation, Domínguez-Tenreiro et al. (2014), who analyzed their HI and H2 content and panchromatic spectral energy distributions (SEDs), and Miranda et al. (2016), who studied abundance gradients, included galaxies run using both codes. In all these cases, the consistency with observational data is very satisfactory, including the differences between the thin and thick disk described in Section 1 that are analyzed in these papers.

A key issue in this paper is to elucidate if the halos of the analyzed systems assembled their mass through a two-phase process, as argued in Section 1.

The answers lie in Figure 1, where the evolution of the halo mass assembly is illustrated through the plot of the halo virial mass, ${M}_{\mathrm{halo}}$(${t}_{{\rm{U}}}$), as a function of time (black curves). In this figure mergers appear as discontinuities, ${\rm{\Delta }}$ ${M}_{\mathrm{halo}}$(${t}_{{\rm{U}}}$), spanning their respective merger time interval. As expected (see Section 1), we see that two phases clearly stand out, a fast phase where the mass-assembly rate (i.e., the curve slope at given times) is very fast (high slopes), and a slow phase in which this rate is low or very low (low slopes), see, i.e., Griffen et al. (2016) for a recent similar result. Note also that at high redshifts, mergers have their respective merger time intervals superimposed in some cases (resulting in a rather continuous slope), with ${\rm{\Delta }}$ ${M}_{\mathrm{halo}}$(${t}_{{\rm{U}}}$) showing a wide range of values. In the slow phase the discontinuities are less frequent and tend to be smaller, with no major mergers in any of the three galaxies. An analysis of these curves in detail, taking into account the cosmological context, indicates that the frequent mergers at early times are but the effects of the collapse of the region surrounding the halo.⁷ This explains the violence of the mass-building events at high redshift, and it is the basis for the so-called two-phase formation scenario.

Zoom In Zoom Out Reset image size

As a second proof of the two-phase behavior shown by haloes, in each of the three panels of Figure 1 we plot the radii, ${r}_{\mathrm{halo}}$(${t}_{{\rm{U}}}$), that enclose the particles forming the respective z = 0 haloes at different ${t}_{{\rm{U}}}$ (blue curves). The radii have been normalized to their respective maximum values. We see that these radii first increase, until they reach a maximum (turn-around), and then they decrease, until they come to a quasi-equilibrium value that is retained up to z = 0. This ${r}_{\mathrm{halo}}$(${t}_{{\rm{U}}}$) behavior indicates that at high redshift, after expansion and turn-around, we witness halo collapse and virialization in consistency with the predictions of the spherical collapse model, see Padmanabhan (1993), followed by a slower evolution period. In each panel of this figure the short time interval separating the decreasing from the quasi-equilibrium behavior of ${r}_{\mathrm{halo}}$(${t}_{{\rm{U}}}$) curve (marked with a pink band) is close to, and in some cases rather coinciding with, the short time interval separating the halo mass-assembly rates from fast to slow. This is expected because both dynamical timescales mark halo virialization. Therefore, only a dynamical timescale ${t}_{\mathrm{vir}}$, based on the pink bands, is considered in this paper.

This two-phase halo mass assembly translates into a two-phase baryon mass assembly at the galactic scale. This is illustrated in Figure 2, where we plot again ${M}_{\mathrm{halo}}$(${t}_{{\rm{U}}}$) (black curves) and ${M}_{\mathrm{bar}}$(${t}_{{\rm{U}}}$), the baryon mass in the galactic object along the evolution (red curves). We see that ${M}_{\mathrm{bar}}$(${t}_{{\rm{U}}}$) roughly follows ${M}_{\mathrm{halo}}$(${t}_{{\rm{U}}}$), and it also shows a two-phase behavior, with no major mergers along the slow phase in any of the three panels. For more details, see for example Figure 3 in Obreja et al. (2013), where the mass aggregation tracks along the main branch of the merger tree (MATs) are given at different fixed radii for its cold baryons and stars. The ${M}_{\mathrm{halo}}$(${t}_{{\rm{U}}}$) curves are also plotted, and a correspondence between the mass increments (i.e., merger events) at both scales, halo and galaxy, clearly stands out.⁸

Zoom In Zoom Out Reset image size

Therefore we can conclude that not only do our simulated galaxies show a two-phase halo mass assembly, but the baryon mass assembly of the whole disk galaxy shows a two-phase behavior as well, which is closely linked to that of its respective halo.

Let us stress that the physical conditions for star formation are very different before and after ${t}_{\mathrm{vir}}$. Before virialization, stellar activity is very high, either as star formation (Figure 2), winds, or as SNe explosions. The many fast mergers inject mechanical energy into the early galactic system, and together with the energy from discrete sources (i.e., energy feedback), they increase the gas turbulence. At the same time, the systems experience high rates of infall of low-metallicity gas. The mixing of infalling with in situ gas, aided by turbulence, is an effective way to remove metallicity gradients. In addition, in this violent phase the gravitational potential is time dependent, and therefore scattering of the already-formed stellar populations can also be expected. Thus, the spatial distribution of stars can change considerably along the fast phase relative to the distribution at their birth time.

In contrast with this situation, during the slow-assembly phase, the potential at the halo scale is nearly spheroidal and axial at galactic scales. Axial symmetry can occasionally be broken, for example owing to disk disruption caused by satellites or when a bar develops. On the other hand, gas turbulence is expected to be low during this phase, because energy injection events are now scarce at any scale.

The expected imprints of these different physical conditions on the properties of the stellar populations born before and after virialization are in line with the properties of the thick- and thin-disk stars, respectively, properties highlighted on the basis of an empirical classification. In Section 5 we show that this is indeed the case. However, the physical timescale t_vir is not an observable. In order to enable a comparison to observations, we first need to determine an observable timescale that is linked to ${t}_{\mathrm{vir}}$ to classify disk stars. In the next section we show that there is a chemical timescale that matches this requirement.

To identify three components in the simulations, we proceed in two steps. First, stars are split into their spheroid (bulge and stellar halo) and disk components using the kk-means method in the space of kinematical variables (see Dhillon et al. 2004; Doménech-Moral et al. 2012). The following kinematical variables are used: energy E, eccentricity $\tfrac{{j}_{z}}{{J}_{c}(E)}$, and $\tfrac{{j}_{p}}{{J}_{c}(E)}$, where j_z and j_p are the projections of a given particle angular momentum on the disk axis and plane, respectively, and J_c(E) is the angular momentum of the circular orbit with E energy. With this method, the spheroid and disk star properties found in observations are recovered (Doménech-Moral et al. 2012; Obreja et al. 2013; Domínguez-Tenreiro et al. 2015). A similar kinematic selection is made for the GASOLINE galaxy.

In a second step, disk stars have to be assigned to either the thin or the thick disk based on an empirical timescale linked to ${t}_{\mathrm{vir}}$. As explained in Section 1, the empirical classification methods resting on a timescale take advantage of the two-slope behavior of the [α/Fe]–age relation, based on solar-neighborhhod stellar data and recently extended to early-type galaxies.

The SFRH for the whole galaxies⁹ are given in Figure 2. We see that early SFRs are high (i.e., they show a low characteristic timescale for SF), and then they decrease to low values at later times. For the g1536-${{\rm{L}}}^{* }$ galaxy, the SF timescale contrast between early and late SFRs is less marked than for the other two galaxies.

The chemical [O/Fe]-age relation for disks is drawn in Figure 1, where we can see that they show a two-slope behavior (rather steep at high z and much flatter at low z), as their zeroth-order shape.¹⁰ The slope contrast is less marked for the g1536-${{\rm{L}}}^{* }$ galaxy than for the other two (because the SF timescale contrast between early and late SFRs is less important here), and in addition, it shows a lurch around ${t}_{{\rm{U}}}\sim 7$ Gyr, i.e., when its SFRH (Figure 2) shows an excess after a dip (see Figure 10 in Snaith et al. 2015, for an explanation).

This global two-slope behavior can be used to define a chemical scale separating the thin-disk from the thick-disk populations. The simplest option is to make a two-slope fit to the [O/Fe]–age plot. Among the different possible solutions, we have chosen the solution ${t}_{\mathrm{chem}}$ that is the closest to ${t}_{\mathrm{vir}}$. The ${t}_{\mathrm{chem}}$ values are observational proxies for ${t}_{\mathrm{vir}}$. Their values for the simulated galaxies are given in Table 1, and they are marked as vertical magenta lines in Figures 1 and 2.

The next task is to ensure that the stellar populations classified as belonging to the thick and thin disks (i.e., older and younger than ${t}_{\mathrm{chem}}$, respectively) do show different elemental abundance distributions as well as different kinematics, and indeed different trends in kinematics with age, with metallicity, and with $[\alpha$/Fe], similar to those found in observational data, see Section 1.

A detailed analysis of the different properties shown by stars in the thick- and thin-disk components in a sample of P-DEVA galaxies can be found in Doménech-Moral et al. (2012). These authors have analyzed the respective sizes and shapes, as well as the age and $[\alpha$/Fe] and [Fe/H] distributions of the different galaxy components according to different classification schemes. Their conclusions are that the properties shown by thin- and thick-disk stars in this simulated sample of galaxies are nicely consistent with observations, regardless of the classification scheme. Similar analyses have been performed for the GASOLINE g1536-L* galaxy (see, for example Stinson et al. 2013a).

In what follows we consider these different properties more closely by focusing on the different trends in three-dimensional kinematics with age, metallicity, and $[\alpha$/Fe] shown by the disks, and by comparing them to observational data. It is worth noting that detailed data for stellar populations are currently only available for the MW, while our simulations have been run from random, i.e., non-constrained, initial conditions. Therefore, only qualitative consistency between the former and the latter can so far be required to pass the validation test, with quantitative agreements adding further strenghths, but they are not strictly necessary.

The particle velocity field is usually expressed in cylindrical coordinates (the symbols ${V}_{\phi }$, ${V}_{{\rm{z}}}$, and ${V}_{{\rm{R}}}$ are used hereafter for the tangential, vertical, and radial velocities, respectively, and ${\sigma }_{\phi }$, ${\sigma }_{{\rm{z}}}$, and ${\sigma }_{{\rm{R}}}$ for their corresponding dispersions). The velocities are calculated relative to the galaxy center of mass, and the system is oriented such that the z-axis at ${t}_{{\rm{U}}}$ is parallel to the angular momentum of the gaseous disk at the same epoch. All the results shown in this section refer to the 3 kpc < ${R}_{\mathrm{cyl}}$ < 20 kpc disk shells, where ${R}_{\mathrm{cyl}}$ is the radial distance in the cylindrical coordinate system. Results for the 6 kpc < ${R}_{\mathrm{cyl}}$ < 10 kpc galaxy shells, somewhat similar to the MW solar neighborhood in its geometric characterization, show no remarkable differences, but have higher statistical noise. For this reason, we present only the results in 3 to 20 kpc cylindrical shells.

5.1. Kinematics Versus Stellar Age

Figure 3 shows the tangential, ${V}_{\phi }$(age), vertical, ${V}_{{\rm{z}}}$(age), and radial, ${V}_{{\rm{R}}}$(age), velocities (dotted lines) and their respective velocity dispersions (solid lines) for the three simulated galaxies. Curves corresponding to the thick-disk (thin-disk) components are plotted using thick (thin) line types. The horizontal axis bin sizes are 0.66 Gyr. No qualitative changes result when half this value is used.

Zoom In Zoom Out Reset image size

We see that in all the cases, as expected, the thin disk appears to be more rotationally supported than the thick disk. The differences at equal ${R}_{\mathrm{cyl}}$ (not drawn) are in the range ∼35–15 km ${{\rm{s}}}^{-1}$, consistent with Veltz et al. (2008), who find an asymmetric drift of ${V}_{\mathrm{lag}}=33\pm 2\,\mathrm{km}\,{{\rm{s}}}^{-1}$ for the Galactic thick disk. Similar values were found for example by Chiba & Beers (2000) and Dambis (2009). The ${V}_{\phi }$(age) curves show a decreasing behavior with increasing age, consistent with results of Haywood et al. (2013), see their Figure 16.¹¹ The g1536-${{\rm{L}}}^{* }$ galaxy shows a remarkable decrement due to major mergers around ${t}_{{\rm{U}}}\sim 7$ Gyr, which changes the ordered radial velocity into velocity dispersion. This shows that our simulated galaxies recover the qualitative behavior of local stars regarding ${V}_{\phi }$(age).

The deviations of the ${V}_{{\rm{z}}}$ and ${V}_{{\rm{R}}}$ curves from zero are a measure of the misalignments of the stellar disk plane from the cold gaseous plane at each ${t}_{{\rm{U}}}$. No remarkable differences are detected, except for weak fluctuations that in some cases become more important at the old age end, where disks are not yet well defined. It is worth noting that these curves inform about occurrences of stellar disks as early as 12 Gyr ago.

Concerning the velocity dispersions, the ${\sigma }_{\phi }$(age), ${\sigma }_{{\rm{z}}}$(age), and ${\sigma }_{{\rm{R}}}$(age) curves show that the thick disk has always higher velocity dispersion values than the thin disk. More specifically, the ${\sigma }_{\phi }$(age) and ${\sigma }_{R}$(age) curves (blue and green solid lines) show common shape patterns as we go from young to old stellar populations (patterns not found in LD-5101A). After a slow increase with age for young stars, a marked dispersion increase within a short age interval is apparent, again followed by a low slope or even flat behavior. Even if the issue is beyond the scope of this paper, we note that these g1536-L* and HD-5004A abrupt velocity dispersion increments occur at ages when the disks suffer a particularly violent and/or destabilizing dynamical event. These are major mergers around ${t}_{{\rm{U}}}\sim 7$ Gyr (age ∼ 6.5 Gyr) for g1536-L*, see Figure 1, or a loss of axial symmetry around ${t}_{{\rm{U}}}\sim 9.5$ Gyr (age ∼ 4.0 Gyr) for galaxy HD-5004A, see Figure 9. Galaxy LD-5101A is not involved in any such events, and it consequently does not show such remarkable slope changes either.

The vertical velocity dispersion (solid red curves) also tends to increase with age. However, the changes in slope are less important and less abrupt than the in-plane dispersions.

These patterns are consistent with the results of Haywood et al. (2013), see the velocity dispersion curves shown in their Figures 11 and 12, where an almost flat behavior at the young age end, followed by a fast increase, stand out.

Comparing Figure 3 with the MW values of these dispersions listed in Table 2, we can see that not only the trends but also the numerical values of the dispersions are satisfactorily consistent with observational data. To assess the degree of consistency, we recall that the LD-5101A galaxy is less massive that the MW and that the g1536-${{\rm{L}}}^{* }$ galaxy has suffered a major merger near ${t}_{\mathrm{vir}}$. Moreover, at the old age end (say, age older than 10 Gyr) the stellar disk of g1536-${{\rm{L}}}^{* }$ galaxy is poorly populated.

Table 2.  Milky Way Velocity Dispersion Values from Observational Data

	${\sigma }_{{\rm{R}}}$ (km ${{\rm{s}}}^{-1}$)	${\sigma }_{\phi }$ (km ${{\rm{s}}}^{-1}$)	${\sigma }_{{\rm{z}}}$ (km ${{\rm{s}}}^{-1}$)	References
Low age end	25	⋯	10	Haywood et al. (2013)
	25 ± 2	20 ± 4	10 ± 1	Vallenari et al. (2006)
Old age end	70	⋯	45	Haywood et al. (2013)
	74 ± 11	50 ± 7	38 ± 7	Vallenari et al. (2006)
	63 ± 6	39 ± 4	39 ± 4	Soubiran et al. (2003)

Download table as:  ASCIITypeset image

5.2. Kinematics Versus [$\alpha /{Fe}$]

The kinematics versus [O/Fe] relations are drawn in Figure 4 for the three simulated galaxies, where the horizontal bin sizes are 0.05 dex. No significant changes occur when we use a bin size of 0.025 dex. We see that the disk tangential velocities, ${V}_{\phi }$, (dotted blue curves) for thick disks are always lower than those corresponding to thin disks. They tend to decrease with increasing $\alpha$ enrichment, in consistency with Recio-Blanco et al. (2014) for chemically classified thin-disk and thick-disk stars from GAIA-Giraffe data, see their Figure 18. It is worth noting that a gap exists between thin- and thick-disk population rotations ${V}_{\phi }$([$\alpha$/Fe]), even in the $[\alpha$/Fe] ranges where both populations have similar values, recovering observational trends found by Haywood et al. (2013) (their Figure 14).

Zoom In Zoom Out Reset image size

The thick-disk vertical and radial dispersion curves, ${\sigma }_{{\rm{z}}}$ and ${\sigma }_{{\rm{R}}}$, respectively, are flat or increase very slowly in Figure 4. The thin-disk population curves increase for [O/Fe] higher than ∼0.0 dex in two galaxies. Thick-disk populations are always dynamically hotter than the corresponding thin-disk populations.

Lee et al. (2011) have analyzed the kinematics–$[\alpha$/Fe] relations in the SEGUE G-dwarf sample. By plotting dispersions versus $[\alpha$/Fe] with no thin-disk versus thick-disk splitting, they showed that the relation is flat at lower [$\alpha /\mathrm{Fe}$], and then it increases. The same behavior is shown in particular by the HD-5004A galaxy when no classification of stars is taken into account. The authors also found that the dispersions at the low $[\alpha$/Fe] end (mostly young stars, that is, stars not affected by specific dynamical events) increase from ${\sigma }_{{\rm{z}}}$([$\alpha$/Fe]) to ${\sigma }_{\phi }$([$\alpha$/Fe]), and further to ${\sigma }_{{\rm{R}}}$([$\alpha$/Fe]), in consistency with our findings. The numerical values are also consistent for HD-5004A. Similar trends have been found in RAVE data (Steinmetz et al. 2006) by Minchev et al. (2014).

5.3. Kinematics Versus [Fe/H]

In Figure 5 we plot the kinematic curves versus [Fe/H]. We used a horizontal bin size of 0.15 dex, and the results are stable when we halve this value, for instance. In this figure, the tangential velocity curves, ${V}_{\phi }$([Fe/H]), for the thick disks of HD-5004A and g1536-${{\rm{L}}}^{* }$ (thick dotted blue lines) increase with [Fe/H]. The tangential velocity thin-disk curve clearly decreases for g1536-${{\rm{L}}}^{* }$. In any case, the thin-disk ${V}_{\phi }$([Fe/H]) curves show a gap relative to the thick-disk curves in the [Fe/H] intervals that are populated by both thin- and thick-disk stars. These patterns of the tangential velocity curves are consistent with the findings of Lee et al. (2011), their Figure 3; Haywood et al. (2013), their Figure 14; Recio-Blanco et al. (2014), their Figure 17; and Allende Prieto et al. (2016) from GAIA DR1 data combined with APOGEE, their Figure 3.

Zoom In Zoom Out Reset image size

As for the dispersions in our simulated galaxies, their shapes for the thick disks are rather flat, except for g1536-${{\rm{L}}}^{* }$, whose ${\sigma }_{\phi }$([Fe/H]) (thick solid blue lines) and ${\sigma }_{{\rm{R}}}$([Fe/H]) (thick solid green lines) decrease slightly for [Fe/H] ≳ −0.5 dex. The thin-disk dispersions decrease with [Fe/H] for HD-5004A and are rather flat in the other two cases. A gap in the velocity dispersion between thin- and thick-disk populations is found for the three galaxies, with those corresponding to the thick disk always higher than those corresponding to the thin disk. The highest dispersions are always in the radial direction. These behaviors are consistent with the results of Recio-Blanco et al. (2014), see their Figure 20.

To summarize, the well-established fact that thin- and thick-disk stars show different kinematics and indeed different trends in kinematics with stellar age, metallicity, and $[\alpha$/Fe], is nicely recovered by the galaxies we simulated. This completes previous results on thick- versus thin-disk star properties. We therefore have thin and thick disks in our simulated galaxies, with the correct behavior as expected from observations. Thus, our simulations can be used as a testbed to analyze the origin of the thin- versus thick-disk differentiation, as well as the origin of the MMP radial distributions found by Bovy et al. (2016). This is the aim of this paper, and the next sections are devoted to study this issue.

The star surface number densities ${\rm{\Sigma }}$(${R}_{\mathrm{cyl}}$) have been calculated for MMPs, and they are depicted in the first row of Figure 6 for the thin-disk (left) and thick-disk (right) stars, respectively. Colors stand for different [Fe/H] bins, which for thin disks are navy blue [Fe/H] < −0.4 and brick red for [Fe/H] > 0.2 ([Fe/H] > 0.3), for HD-5004A and LD-5101A (g1536-L*), the other colors stand for five (six) intermediate 1 dex wide [Fe/H] bins. For the thick disks, color bins are pale green [Fe/H] < −0.8, sky blue for [Fe/H] > −0.25, and the other colors stand for three intermediate bins ∼0.2 dex wide.

Zoom In Zoom Out Reset image size

The MMP surface number densities show a dichotomy: (i) the ${\rm{\Sigma }}$(${R}_{\mathrm{cyl}}$) for thick disks have similar scale lengths regardless of the MMP (except for the most metal-rich one in g1536-L*), and no breaks (see Stinson et al. 2013a, for similar results), (ii) for the thin disks, ${\rm{\Sigma }}$(${R}_{\mathrm{cyl}}$) have different shapes depending on the MMP: those corresponding to high [Fe/H] bins peak at low ${R}_{\mathrm{cyl}}$, while those corresponding to low [Fe/H] values are rather flat. Moreover, some slope changes appear that tend to be placed at increasing ${R}_{\mathrm{cyl}}$ for lower [Fe/H]. Qualitatively, we recover the APOGEE results for the Galactic MMP radial structure for high and low $[\alpha$/Fe] populations (a proxy for old and young populations, respectively).

To decipher the physical processes underlying this behavior, the birth places for each star belonging to either the thin or the thick disks have been determined. In the bottom panels of Figure 6 we plot the surface number densities that the thin (left) and thick (right) disks would have if stars had not moved away from their birth places. The [Fe/H] bins are the same as those corresponding to their respective upper panels. These plots indicate that different MMPs are preferentially born at different locations, characterized by an increasing ${R}_{\mathrm{cyl}}$ as the metallicity decreases, the effect being much more marked for the thin- than for the thick-disk populations. We can also see that the stars in thick disks diffused outward between their birth time and z = 0.

Different MMPs are also preferentially born at different epochs, ${t}_{{\rm{U}}}$, as Figure 7 shows for the disk of g1536-${{\rm{L}}}^{* }$, see the color code in the caption. Note that they are not mono-age populations, see Minchev et al. (2017) for similar results.¹² At fixed times, i.e., horizontal cuts in this figure, stars reflect the gas metallicity structure at their birth.

Zoom In Zoom Out Reset image size

The gas distributions at different times, ${t}_{{\rm{U}}}$, are shown in Figure 8 for HD-5004A, within the range of thick- and thin-disk formation ages (top and bottom rows, respectively). Colors stand for different [Fe/H] bins, with the same color code as in Figure 7. In this figure we can see a segregation in galactocentric distance, according to the gas metallicity in the thin disk. This segregation is not removed by metal enrichment due to stellar evolution. This segregation and its preservation across time implies that thin-disk stars belonging to different MMPs are born at preferential times and disk locations.

Zoom In Zoom Out Reset image size

In the early period of the thick disk, when its stars are born, the gas ${R}_{\mathrm{cyl}}$ segregation due to metallicity is not as important as later on. In fact, this segregation does not appear until z ∼ 1.5, and once it appears, it is almost removed each time a dynamically violent event occurs (because they favor metal mixing, see Section 3). These sequences of order-disorder can be appreciated in the top rows of Figure 8. Similar situations arise for the other two simulated galaxies. As a consequence, the stellar MMPs differentiation at birth for thick-disk stars are not as strongly marked as in the case of the thin disk.

Violent dynamical events occur much more often at early times than later on during the thin-disk formation period because of the two-phase characteristic of mass assembly. This has been proven and discussed in Section 3. We now focus on how the dichotomy of fast/slow rates of mass assembly affects the characteristics of thin-disk and thick-disk stellar populations and, particularly, those of their MMPs radial distributions.

As stated in Sections 3 and 4, the thin disk begins to form after the galaxy halo reaches its virial equilibrium, becoming, to a first approximation, a quasi-stationary system with almost spherical symmetry at the halo scale, and an axial symmetry at the galaxy scale. We have found that in this situation, gaseous disks form inside out and chemical evolution proceeds within them, building up element gradients as a consequence of stellar evolution. These gradients are detected as the gas ${R}_{\mathrm{cyl}}$ segregation we see in the second row of Figure 8.

In contrast, the thick disk forms its stars while the galaxy halo is still in its fast-assembly phase, when the halo has a time-dependent potential. In this period, many fast mass increments occur, carrying with them matter (dark matter, gas, and some stars) that is either diffuse or clumpy, and energy. In this unstable situation, the chemical gradient is built and removed, as we see in the sequence plotted in the top panels of Figure 8. Gradient removal occurs not only as a result of the infalling gas mixing with the local disk gas, but presumably is also helped by an increase in turbulent diffusion caused by the injection of dynamical and stellar feedback energy into the system.

Finally, we explore why the ${R}_{\mathrm{cyl}}$ distributions of stellar disk birth places do not match the current distributions at z = 0, for either the thin or the thick disks. To see this, in Figure 6 we compare the first with the second rows on the left (right) panels for the thin (thick) disks. These differences must be due to stellar radial mixing (blurring, radial migration, or scattering,¹³ see, e.g., Sellwood & Binney 2002 and Halle et al. 2015).

Even if the precise distinction of these mechanisms is beyond the scope of this paper, the importance of radial mixing and its correlation with destabilizing dynamical events is not difficult to asses.

For the thin disk, Figure 6 indicates that the radial mixing decreases from HD-5004A to LD-5101A to g1536-${{\rm{L}}}^{* }$. An analysis of the HD-5004A evolution during the slow phase, shown in Figures 1 and 2, indicates that a burst of star formation occurs around ${t}_{{\rm{U}}}\sim 9.5$ Gyr. Within the same time interval, this galaxy develops a non-axial configuration, which clearly stands out at z = 0.33 in Figure 9, where we can also see that the destabilizing agent is a double minor merger. In this way, thin-disk stars suffer an important radial migration.

Zoom In Zoom Out Reset image size

No important breaking of axial symmetry appears during the slow phase of LD-5101A. This means that radial migration and its effects are milder than in the previuos case. The effects of radial migration are even milder in the case of galaxy g1536-${{\rm{L}}}^{* }$.

The different unstability patterns shown by the three simulated galaxies during their slow phase can be understood in terms of their respective environments. Indeed, while HD-5004A lives in a dense environment, the environment of LD-5101A has a lower density. In addition, the host halo of g1536-${{\rm{L}}}^{* }$ has low merging activity at late times.

These correlations strongly suggest that radial mixing in the slow phase is an effect of disk instabilities, caused by dynamical activity that temporarily breaks the axial symmetry of the system. Hence, we can say that the final MMP ${R}_{\mathrm{cyl}}$ distributions of the thin disk is shaped by disk instabilities. The formation of thick-disk stars in a time-varying potential implies a high degree of dynamical scattering, which explains the uniformity of their ${R}_{\mathrm{cyl}}$ distributions, regardless of the MMPs.

It this paper we addressed the issue of the physical processes underlying the origin of the thick disk relative to the thin disk. Recent results on radial distributions of mono-metallicity populations (MMPs, i.e., stars in narrow bins in [Fe/H] within wider [α/Fe] ranges) in the Galactic disk by Bovy et al. (2016) have cast doubt on the classical bimodality of thin versus thick disks. Our work aimed at explaining such MMPs radial distributions in terms of these physical processes. We showed that these distributions are straightforward consequences of the two-phase mass-assembly scenario for disks.

It is important to note that in this work similar results have been obtained from the analyses of galaxies simulated using two very different codes, placing our results on the level of fundamental physical processes beyond particular code implementations.

9.1. Summary of Results

9.2. Discussion

9.3. Conclusions

This work was partially supported through MINECO/FEDER AYA2012-31101 and AYA2015-63810-P grants (Spain). C.B. thanks the Ramon y Cajal program. A.O. has been funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—MO 2979/1-1, and thanks the MICINN and MINECO (Spain) for financial support through a FPI fellowship. We also acknowledge the Centro de Computación Científica (Universidad Autónoma de Madrid, Red Española de Supercomputación) for computational support.