© P. Di Matteo et al. 2019

1. Introduction

The Milky Way (hereafter MW) bulge is a boxy or peanut-shaped bulge (hereafter B/P; see Okuda et al. 1977; Maihara et al. 1978; Weiland et al. 1994; Dwek et al. 1995; Binney et al. 1997; Babusiaux & Gilmore 2005; López-Corredoira et al. 2005; Rattenbury et al. 2007; Cao et al. 2013; Wegg & Gerhard 2013; Portail et al. 2015, 2017a,b; Ciambur et al. 2017). In its vast majority, it consists of stars that originated in its disk. Several works indeed now agree that the contribution of a classical bulge needs to be limited to a low percentage of its total mass (Shen et al. 2010; Kunder et al. 2012, 2016; Di Matteo et al. 2014; Debattista et al. 2017; Gómez et al. 2018). The contribution of the stellar halo, whose density is expected to peak in the inner few kiloparsecs of the Galaxy, also seems to be marginal and to not exceed a few percent of the bulge total mass (Ness et al. 2013a). While there is growing consensus on these general results, the details remain more controversial. One of the points that is still debated concerns the origin of the metal-poor (−1 dex  ≤  [Fe/H]  ≤  0 dex), α-enhanced population observed in the bulge. Compared to the metal-rich population ([Fe/H]> 0), metal-poor stars appear kinematically hotter at latitudes b ≲ −4 deg, and colder closer to the midplane (Ness et al. 2013b; Rojas-Arriagada et al. 2014; Zoccali et al. 2017). They are not strongly peanut-shaped (Rojas-Arriagada et al. 2014, but show an evident difference between stars with metallicities below and above ∼ − 0.5 dex, see Ness et al. 2013a), and dominate the outer bulge region. Their proportion, relative to the total population, decreases with decreasing height from the Galactic plane until at least latitudes b ≃ −5 deg (Ness et al. 2013a), and then it seems to show an inversion, becoming dominant again at very low latitudes (Zoccali et al. 2017).

Di Matteo et al. (2014, 2015) proposed that the metal-poor population observed in the bulge is the Galactic thick disk, mapped in the B/P-shaped bulge differently from the metal-rich stars because of its hotter kinematics. In particular, we have suggested (see also Ness et al. 2013a, for this same interpretation) that stars with metallicity between −1 ≤ [Fe/H]≤ − 0.5 dex can be associated with the old Galactic thick disk (ages greater than about 10 Gyr and metallicities below −0.5 dex, see Haywood et al. 2013), while bulge stars with metallicities in the range −0.5 ≤ [Fe/H]≤0 would correspond to stars of the young Galactic thick disk (ages between 8 and 10 Gyr and metallicities between −0.5 dex and solar, see Haywood et al. 2013). The reasons behind this suggestion have been reviewed in Di Matteo (2016), and we refer to that paper for a discussion of the disk-bulge connection in this scenario. In particular, Di Matteo (2016) and Fragkoudi et al. (2017a) have shown that a composite disk galaxy made initially of a kinematically cold, intermediate, and hot disk, which mimick the thin disk, the young disk, and the old thick disk of the Galaxy, respectively, can reproduce the morphological and populations trends observed in the MW bulge. These are that the colder the population, the stronger the bar, and at all heights where a B/P morphology is observed, the hotter the population, the weaker the B/P structure. The chemical patterns and metallicity distribution function of bulge stars can also be understood as a consequence of the differential mapping of the composite Galactic (thin+thick) disk into the B/P bulge, modulo their scale lengths and heights (see Fragkoudi et al. 2017b, 2018; Haywood et al. 2018). In this context, the agreement of the metallicity map that this scenario predicts with that of the MW bulge, as derived from APOGEE data, is remarkable (Fragkoudi et al. 2018).

Our findings, based on models tailored to match the MW bulge and inner-disk populations, agree with the general result that in a composite disk, populations with different kinematics are mapped differently into a bar and B/P bulge. This process is referred to as kinematic fractionation by Debattista et al. (2017). This result, first shown by Bekki & Tsujimoto (2011), has since then indeed been confirmed and extended in a number of other works (Di Matteo 2016; Fragkoudi et al. 2017a; Athanassoula et al. 2017; Buck et al. 2018). However, Debattista et al. (2017) suggested that the main driver of this different response to the bar perturbation is the radial, in-plane random motion of stars, not their vertical motions. In their view, a composite disk made of stellar populations with the same initial thickness but different in-plane random motions constitutes a sufficient condition to reproduce the trends observed in the MW bulge, and they conclude that as a consequence, a thick disk, that is, a disk both radially and vertically hot, is not necessary to explain the main characteristics of the MW bulge.

In this paper we revisit this issue about the driver of the trends observed in the MW bulge. First, we aim at understanding whether the in-plane, radial random motions of disk stars constitute the main drivers of the trends observed in the MW bulge, as suggested, and the role played by vertical motions. Second, by comparing some of the current properties of the MW bulge with the model predictions, we aim at finding signatures and suggesting criteria to constrain the kinematic conditions that were present in the early MW disk.

For this study, we analyzed two dissipationless N-body simulations of composite-disk galaxies. These two models represent two extreme cases: one where all three components of the disk have a fixed vertical velocity dispersion and different planar/radial (in-plane) dispersion, as in Debattista et al. (2017), and another where they all have a fixed radial dispersion and different vertical random motions (thickness). Our analysis reveals that disk populations with different initial vertical random motions and the same in-plane random motions are at first order subject to a similar mapping into a B/P-shaped bulge as is experienced by stellar populations with equal initial thickness and different in-plane motions. More specifically, both initial conditions lead to composite B/P bulges, where former disk stars with initially the highest velocity dispersions redistribute in a thicker B/P structure and weaker bar than disk stars that initially have a cold kinematics. As a consequence of this similar mapping, the B/P bulges in both models show vertical metallicity and age gradients, and a pinching of the metallicity maps along the bulge minor axis, as also observed Gonzalez et al. (2017). Despite these same trends, we note some significant differences in the morphology of the resulting B/P structure in the two cases, however. In particular, when the disk is made of stellar populations with the same initial thickness, which are different only in their in-plane motions, as in Debattista et al. (2017), the peanut-shaped structure appears at the same height above the plane for all populations, independent of their kinematics, and its strength and shape is remarkably similar. This is at odds with what is observed in the MW bulge. We discuss the reasons behind this different behavior and how the coupling of radial and vertical motions acts in the two cases.

These two models constitute a first step to explore to what extent in-plane and vertical random motions reshape the properties of B/P bulges and their surrounding disks, but none of them evolves into properties that are compatible with those of MW disk. As we discuss, it is indeed necessary that disk populations both radially and vertically warm, such as those currently found in the MW thick disk, were present in the Galaxy at the beginning of its secular evolution, in order to reproduce some of the trends observed in the Galaxy.

The paper is organized as follows: in Sect. 2 we describe the two models we analyze in this paper, their properties, and the adopted numerical methods; in Sect. 3 we present the results, by first discussing some main general trends (morphology, metallicity, and ages) of their B/P bulges and their surrounding disks (Sect. 3.1), then the strength of the B/P shape and its dependence on the initial disk kinematics (Sect. 3.2), and finally, we outline our conclusions in Sect. 4.

2. Models

In this paper, we analyze two dissipationless simulations of galaxies that consist of the superposition of disks with different initial kinematics. In Model 1, the galaxy consists of three stellar disk components of equal initial thickness (and hence initial vertical velocity dispersion), but different initial in-plane motions (i.e., different radial and tangential velocity dispersions). In Model 2, the galaxy also consists of three stellar disks, but with equal initial in-plane motions and different initial thickness. These two models do not aim at representing a realistic model for the MW at the time of the formation of the bar, but are intended to constitute two significantly different initial conditions, which are required to explore the drivers of the main characteristics observed in the MW bulge. In particular, Model 1 has a similar setup to the models analyzed in Debattista et al. (2017) (i.e., equal thickness and different in-plane motions) and therefore allows a comparison with this work.

In each model, the three disks are represented by Miyamoto-Nagai density profiles (Miyamoto & Nagai 1975), with characteristic heights and scale lengths as given in Table 1. In each of the two models, the cold, intermediate, and hot disks have masses that represent 50%, 30%, and 20% of the total disk mass, respectively. In this context, the choice of the mass of each of these disks is arbitrary, and the conclusions of this work do not depend on it. However, we use these percentages because they are reminiscent of those of the thin, young, and old thick disks, respectively (for a derivation of the Galactic mass growth and relation with the thin and thick disks, we refer to Snaith et al. 2014, 2015).

These composite disks are then embedded in a dark matter halo, which is modeled as a Plummer sphere (Plummer 1911), and whose parameters are reported in Table 1. The resulting rotation curves are shown in Fig. 1. Face-on and edge-on maps of stars in the bulge region at the initial time are given in Appendix A, where we also discuss further details about the adopted parameters and their effect on the final bulge morphology.

Fig. 1. Initial conditions of our composite stellar disks. First row: Circular velocity (black curve) of the model, together with the rotation curves of the cold (blue), intermediate (green), and hot (red) disks. Second row: In-plane radial velocity dispersion of the cold (blue), intermediate (green), and hot (red) disks as a function of the distance R from the galaxy center. Third row: In-plane tangential velocity dispersion of the cold (blue), intermediate (green), and hot (red) disks as a function of the distance R from the galaxy center. Bottom row: Vertical velocity dispersion of the cold (blue), intermediate (green), and hot (red) disks as a function of the distance R from the galaxy center. Model 1 is shown in the left column, and Model 2 in the right column, as indicated.

To generate initial conditions, we employed the iterative method described in Rodionov et al. (2009). This method allows generating initial conditions at equilibrium with the required density profiles and/or kinematics constraint, avoiding the relaxation processes (and departure from initial conditions) that are often observed in N-body models of disk galaxies. This is particularly useful for building composite disks, each with a specific kinematic and/or density profile, allowing a full control on their initial state. For the models presented here, the initial velocity dispersion profiles of the cold, intermediate, and hot populations are reported in Fig. 1. In Model 1 the initial thickness of the three disks is the same, as is the corresponding vertical velocity dispersion profile; the in-plane motions (radial and tangential) are different, however. In contrast, in Model 2, the vertical velocity dispersions are different, while the in-plane motions are initially the same. We investigate in the following how these two models evolve secularly, and how these (initially different) composite disks are mapped into the B/P bulge. We refer to Appendix B for a discussion of the velocity dispersion structure that is imposed for the initial disks and of the final convergence of the method.

The two simulations were run with a recently developed parallel MPI tree-code that takes into account the adaptive spatial decomposition of particle space between nodes. The multi-node tree-code is based on the 256 bit AVX instructions, which significantly speed up the floating point vector operations and sorting algorithms (Khoperskov et al., in prep.). Fifteen million particles were employed for each model, 10 million in the disk components, and 5 million in the dark matter halo. A time step Δt = 2 × 10⁵ yr was adopted, together with a gravitational smoothing length of 50 pc. Both simulations were run over a timescale of 5 Gyr. In less than 1 Gyr from the beginning of the simulations, a bar and a B/P bulge form in both models. We then analyzed the final configuration, at t = 5 Gyr, to quantify the properties of the B/P bulge and surrounding disk.

3. Results

3.1. Multiple populations with different in-plane or vertical kinematics: relative contribution to the B/P bulge, metallicity, and age trends

The face-on and edge-on projection of stars in the bulge region for Model 1 are shown in Fig. 2. We recall that in this model, the simulated galaxy consists of three disks with different in-plane kinematics (the radial velocity dispersion increases moving from the cold disk to the hot disk, see Fig. 1), but same vertical velocity dispersions. This setup is thus similar to those analyzed in Debattista et al. (2017), and the morphology of the bulge in Fig. 2 indeed agrees well qualitatively with the idealized model presented in their paper. We note, however, that the two models remain different in the details and thus illustrate the possible variety of models and morphologies present in each class.

Fig. 2. Model 1, morphology of the bar-B/P region at t = 5 Gyr. First and third column, from top to bottom: Face-on (Col. 1) and edge-on (Col. 3) absolute stellar densities of all stars (top panel), and stars initially in the cold (second panel), intermediate (third panel), and hot (bottom panel) disks. Second and fourth column, from top to bottom: Fraction of cold- (second panel), intermediate- (third panel), and hot-disk (bottom panel) stars in the B/P bulge, seen face-on (Col. 2) and edge-on (Col. 4). The fraction is defined as the ratio of cold-, intermediate-, or hot-disk stars to the total (i.e., cold+intermediate+hot) disk stars in the B/P bulge. In the edge-on maps, only stars in the bar (|x| ≤ 5.5 kpc and |y| ≤ 3 kpc) have been selected.

Our analysis for this model confirms the findings of Debattista et al. (2017; see also Combes et al. 1990; Athanassoula 2013, for earlier works on the strength of bars in single-component disks and their dependence on the initial disk kinematics). Stars with different in-plane kinematics are mapped differently in the B/P bulge. When viewed face-on, the colder the initial disk, the stronger the bar. This can be appreciated by looking at the isodensity contours, which are more elongated for the initially cold disk, but it is also quantified by means of the m = 2 coefficient, A₂, of the Fourier analysis of the face-on density maps (Fig. 3, left column), which shows that the colder the disk initially, the higher the corresponding A₂ value. In the edge-on projection, stars with initially the coldest kinematics show the thinnest distribution, as we verified by estimating ⟨h_z⟩, the median of the absolute value of the z-component of their positions, as a function of the distance R from the galaxy center (Fig. 2). Because of this different response to the bar perturbation, stars that initially belong to disks with different in-plane kinematics contribute differentially to the B/P bulge in both its face-on and edge-on projection. In agreement with Debattista et al. (2017), the colder the disk, the higher its relative contribution to the B/P structure. Stars in the initially kinematically hottest disk dominate at larger heights above the disk. As a consequence, their relative contribution along the B/P minor axis increases with height above the plane (see Fig. 2, fourth column. This trend, that is, the larger thickness of the initially in-plane hot disk with respect to the colder ones, arises because in the inner regions of the simulated disk the azimuthal and vertical frequencies, Ω and ν, respectively, are similar (see Appendix C), and thus in-plane and vertical motions are not decoupled, as is usually assumed in the case of very thin or infinitely thin disks (see Binney & Tremaine 1987). As a result, a disk that initially has only a hot in-plane kinematics will rapidly relax into a disk that is also vertically hot, because of this coupling.

Fig. 3. First row: Bar strength A2 as a function of time for all stars (dashed black curve), and for stars initially in the cold (blue curve), intermediate (green curve), and hot (red curve) disks. Second row: Height profiles in the B/P region of all stars, and stars initially in the cold, intermediate, and hot disks. Models 1 and 2 are shown in the left and right columns, respectively. In the panels in the first row, the gray area indicates the time of disk bending. At the end of this instability phase, a B/P-shaped bulge is formed in both models, and its morphology does not change significantly after that time.

The analysis of Model 2 interestingly reveals that all the trends observed for Model 1 are satisfied in this case as well (Fig. 4). When viewed face-on, the lower the initial vertical velocity dispersion of the disk, the stronger the bar, as is measured by means of the A₂ coefficient as a function of time (Fig. 3, right column). When viewed edge-on, the lower the initial vertical velocity dispersion of the disk, the thinner the B/P structure. We note, however, that the thickness of cold and hot populations in the two cases is different: outside the B/P region, the hot population in Model 2 has a higher thickness than that measured for the cold population, while the difference is milder in Model 1. For example, at R = 6 kpc, the thickness of the hot population is about three times higher than that measured for the cold population, while in Model 1 the hot population is only 1.3 times thicker than the cold one at the same radius. Finally, because the three disks in Model 2 have, by construction, a different initial thickness, their relative contribution to the B/P bulge depends on their initial kinematics: in particular, the weight of stars with initially the hottest kinematics increases with the height above the plane.

Fig. 4. Same as Fig. 2, but for Model 2.

The similarity in the trends of cold and hot populations with height above the plane that we found for the two models also has some important implications on the mean metallicity and age of stars throughout the B/P bulge, as we investigate in the following. To this aim, we assign to stars in each disk a metallicity and age as follows: the cold, intermediate, and hot disks have Gaussian metallicity distributions with means at [Fe/H] = 0.3, −0.2, and −0.6 dex, respectively, and dispersions σ_[Fe/H] = 0.1, 0.3, and 0.3; for the ages, we assume for the cold, intermediate, and hot disks uniform distributions in the intervals [0, 8], [8, 10], and [10, 13] Gyr, respectively. The metallicity mean values and dispersions are similar to those found by Ness et al. (2013a) for components A, B, and C in the Galactic bulge. Since, as recalled in the introduction, in our scenario we associate these three bulge components. By “components” here we mean stars that are associated with different metallicity intervals, as defined by Ness et al. (2013a) for A, B, and C. with the Galactic disk(s) (the thin, young thick, and old thick disks, respectively), which, as robustly established, show a decreasing metallicity with increasing (in-plane and vertical) velocity dispersions, the choice to assign these metallicity values to our three disks is natural. We emphasize, however, that other choices would have been possible, and that the trends discussed in the following are not affected by the specific metallicity values we employed, because the three modeled disks have different chemistry and ages (i.e., a trend between metallicity/age and velocity dispersion must exist). Since age distributions in the Galactic bulge as well as in most of the Galactic disk are still lacking, we decided to assign them in the simplest way by relating our modeled disks to the thin, young, and old thick disks, respectively, and making use of the age estimates of (thin and thick) disk stars at solar vicinity, as given by Haywood et al. (2013). With these choices of metallicities and ages, half of the disk stellar mass has metallicities below solar and ages younger than 8 Gyr, in agreement with the results presented in Snaith et al. (2014) and Haywood et al. (2016) for the inner Galactic disk (inside 6–7 kpc from the Galactic center).

Before the formation of the bar and the B/P structure, the edge-on metallicity and age maps of Models 1 and 2 are shown in Fig. 5. Because in Model 1 all disks initially have the same thickness, no trend is visible when the modeled galaxy is seen edge-on, and the vertical gradient along the bar major axis is null. For Model 2, however, because the disks have different initial vertical velocity dispersions, a vertical gradient is in place ab initio in both the metallicity and age maps, with mean metallicities(/ages) decreasing(/increasing) with height above the plane. As an example, for Model 2 along the bulge minor axis, the initial metallicity gradient is −0.20 dex kpc⁻¹, and the age gradient is 1.72 Gyr kpc⁻¹. The two gradients were estimated by comparing the metallicity values (/age) at z = 0 and z = 2 kpc from the galaxy plane. We note, however, that in Model 2 the metallicity (/age) gradient is mostly due to the presence of very metal-rich (/young) stars close to the galaxy midplane (at heights below 500 pc), while for larger heights the gradient is nearly flat in this case as well. Once the B/P bulge is formed, the metallicity and age maps of Models 1 and 2 appear rather different from their initial state. In both models, X-shaped metallicity and ages maps are clearly visible (as found also in the models by Athanassoula et al. 2017; Debattista et al. 2017), with metal-rich (young stars) associated with the initially cold disk mostly redistributed in the peanut-shaped configuration. Along the bulge minor axis, vertically away from the plane, the metallicity and age maps become rapidly dominated by the metal-poor, old intermediate and hot disks, giving rise to a vertical negative (/positive) metallicity (/age) gradient. It is striking to note that not only are the global metallicity and age trends similar in the two models (Fig. 5), but the absolute values and strengths of the gradients are similar as well. Except for the outermost bulge region outside the B/P lobes, which we discuss in the following section, at all x ≲ |4| kpc, the vertical metallicity and age profiles appear very similar, with comparable slopes. As an example, along the bulge minor axis, between z = 0 and z = 2 kpc from the galaxy midplane, the metallicity gradient is equal to −0.089 dex kpc⁻¹ in Model 1 and to −0.091 dex kpc⁻¹ in Model 2. The age gradients are 0.721 Gyr/kpc for Model 1 and 0.782 Gyr kpc⁻¹ for Model 2. The age and metallicity maps generated in these two models are thus remarkably similar, despite the significantly different initial conditions that we adopted, and this is valid for all choices of the initial metallicity and age distributions in the disk populations, since here metallicity and age are simply tags of the stellar particles that trace their spatial distribution.

Fig. 5. Edge-on metallicity and age maps for stars in Model 1 (first and third row) and for Model 2 (second and fourth row) at time t = 0 (first column) and at time t = 5 Gyr (second column). Only stars in the bar (|x| ≤ 5.5 kpc and |y| ≤ 3 kpc) have been selected. In all these maps, the bar when present is oriented side-on. Isodensity contours are shown in black.

As a final remark, we emphasize that both models produce B/P bulges whose metallicity/age maps appear more pinched or peanut-shaped than the stellar density distribution itself (cf. metallicity/age maps with isodensity contours in Fig. 5). This feature has been observed in external bulges (Gonzalez et al. 2017), and has also been reported in N-body simulations (Debattista et al. 2017; Aumer & Binney 2017), and here we show that it is not necessarily the signature of a pre-existing disk with stellar populations differentiated in their in-plane random motions, as previously suggested.

From this analysis we thus conclude that initial vertical random motions are as important as in-plane random motions in determining the overall bulge population, metallicity, and age trends, and that previous statements emphasizing the dominant role of in-plane motions in determining these trends are not confirmed by our analysis.

Finally, we would like to briefly comment on the metallicity/age trends that these two models generate outside the B/P bulge in the surrounding disk. In Fig. 6, we show the large-scale (i.e., over the whole extent of the galaxy) edge-on metallicity and age maps for both models. As shown in Appendix D, during the whole evolution in Model 1, the populations outside the bulge region are never differentiated in their vertical velocity dispersion and thus scale heights. This implies that the vertical metallicity and age profiles of the disk do not exhibit any vertical gradient as a result of combining different populations. The consequence of the absence of any vertical differentiation in kinematics present in Model 1 is that the mean metallicity and age of the disk do not vary with z: this feature is already evident at the borders of the B/P bulge and persists throughout the disk. This means that while in the B/P bulge the two models predict similar trends of metallicity/age with height above the plane, significant differences are expected outside the B/P region and in the surrounding disk. In particular, the absence of any dependence of metallicity or age on the height above the plane, as observed in Model 1, is in sharp contrast to what we know about the MW disk: both in the solar vicinity and on a scale of several kiloparsec, metallicity and ages vary with height above the plane and with the vertical velocity dispersion of the corresponding populations (Strömberg 1946; Spitzer & Schwarzschild 1951; Nordström et al. 2004; Seabroke & Gilmore 2007; Holmberg et al. 2007, 2009; Bovy et al. 2012, 2016; Haywood et al. 2013; Sharma et al. 2014; Martig et al. 2016; Ness et al. 2016; Mackereth et al. 2017). To generate vertical metallicity/age gradients in the disk and in the region outside the lobes of the B/P bulge, it is necessary that populations of different metallicities/ages are differentiated in their vertical velocity dispersion ab initio before they are trapped in the bar instability, otherwise, this relation cannot be put in place afterward by secular evolution processes alone.

Fig. 6. Edge-on metallicity and age maps at the final time of the simulation for Model 1 (first and second row), and Model 2 (third and fourth row). The B/P region and the surrounding disk are shown in the maps. Isodensity contours are shown in black.

3.2. Multiple populations with different in-plane or vertical kinematics: on the shape and strength of the B/P bulge

Although Figs. 2 and 4 show similar trends with height from the plane, we note that the initial disk kinematics imprints some significant differences in the intrinsic morphology of the final B/P bulges. A first way to quantify these differences is shown in Fig. 7, where we estimate the strength of the B/P structure by means of the sixth Fourier component (B₆) of edge-on density isophotes in the bulge region (see Ciambur 2015; Ciambur & Graham 2016, for further details¹). We observe that while in Model 2 the amplitude of the B₆ coefficient depends on the initial kinematics of the population (the hotter the kinematics, the lower B₆, and this throughout the whole extent of the bulge), in Model 1 the B₆ amplitude appears to be the same for all populations, independently of their initial kinematics, both as a function of the projected radial length of the peanut shape, R_Π, and of its projected height, z_Π (for a definition of R_Π and z_Π, see Fig. 1 in Ciambur & Graham 2016). In the regions where all three populations coexist in the bulge (about z_Π≤1 kpc), the trends of B₆ with R_Π and z_Π are indeed indistinguishable for Model 1.

Fig. 7. First row: B/P strength, quantified by means of the B6 parameter, as a function of the projected radial length of the peanut shape, RΠ, for Model 1 (left panel) and Model 2 (right panel). Second row: B/P strength as a function of the projected height of the peanut shape, zΠ, for Model 1 (left panel) and Model 2 (right panel). In all plots, the strength of the peanut shape is shown at t = 5 Gyr for all stars (black curves) and for stars initially belonging to the cold (blues curves), intermediate (green curves), and hot (red curves) disks. In all panels, error bars are represented by shaded areas and indicate the intensity rms scatter along each isophote after subtracting the best-fit pure ellipse with harmonic terms up to and including the sixth order.

A second way to appreciate the different B/P shapes in the cold and hot populations of the two models is illustrated in Fig. 8, where we show isodensity contour maps of the B/P bulge in Models 1 and 2. For both models and for all populations, we normalized the isodensity contours to their central density value. The normalized density levels shown by the contours were chosen in order to probe similar regions of the B/P bulge for hot and cold populations. The difference between the two models appears very clearly: while in Model 1 the isodensity shapes inside the B/P lobes, that is, x ≤ 2 kpc, are the same at all heights for both cold and hot populations, in Model 2 the shape of the isodensity contours depends on the initial kinematics of the population: the hotter it is initially, the rounder and less peanut-shaped the contour is at all z. We refer to Appendix A for a more extensive discussion on the populations that dominate the innermost bulge regions of Model 2. Figure 8 thus shows that in Model 1 the thickness of the B/P structure changes moderately with the initial kinematics of the population, but not its intrinsic shape. The trends observed for Model 1 in Fig. 2 (the decreasing contribution of the cold-disk population with height above the plane when compared to the hot population) and in Fig. 5 (the X-shaped metallicity/age maps) are rather due to final different vertical density gradients of the cold and hot populations, and thus to different thicknesses, rather than to an intrinsic difference in their isodensity shape. This result is interesting because it shows that X-shaped metallicity maps are not only the result of the spatial superposition of populations with different metallicities and different intrinsic peanut or X shapes, but can also be produced when multiple populations with different metallicities, but similar X shapes, coexist in a B/P bulge. In this latter case, the X-shape metallicity and age maps result simply because of a different vertical density decay of these populations.

Fig. 8. Isodensity contours of the B/P bulge in Model 1 (top panel) and Model 2 (bottom panel) for the initially kinematically cold (blue contours) and hot (red contours) disk populations at t = 5 Gyr. The contours of the intermediate disk are not shown, but they are bracketed by those of the cold and hot disk populations. The numbers at the top left of each panel indicate the values of the isodensity contours for cold and hot populations normalized to their central density. Only stellar particles in the bar region (|x| ≤ 5.5 kpc and |y| ≤ 3 kpc) have been selected for these plots.

Another consequence of the similarity in the isodensity contours found for Model 1 is that the peanut shape appears at the same height above the plane, at about 600 pc, for all populations, independent of their initial in-plane kinematics (see also Fig. 9). This trend seems at odds with what is found in the MW bulge, where the strength and appearance of the B/P shape depends on the metallicity of bulge stars (Ness et al. 2013a; Rojas-Arriagada et al. 2014). According to our experiment, this observational result favors a scenario where metal-rich and metal-poor populations in the MW bulge initially had different initial vertical random motions, that is, before they were trapped by the B/P bar. If the metal-poor and metal-rich populations were different only in their in-plane random motions, but not in the vertical ones, the bimodality in the stellar density distribution would show a weak or even null dependence on the height above the plane, as we find for Model 1.

Fig. 9. Distribution of stars along the bar major axis at different heights above the plane, as indicated at the top of each panel. The distributions of cold, intermediate, and hot stars in the B/P structure are indicated with blue, green, and red histograms, respectively. The distribution appears bimodal for all populations at the same height above the plane in Model 1 (left panel), while in Model 2 (right panel), the bimodality appears at larger heights for the hot population than for the colder ones.

The question remains why the B/P shape of hot and cold populations so similar in a disk whose populations are only differentiated by their in-plane motions, as in Model 1, and how general this result is. While it is premature to generalize, and a larger statistics of Model 1-like simulations is required before deriving final conclusions, we can give some first explanations here of the observed similarities in the B/P morphology of cold and hot populations in Model 1. In this model, all populations initially have the same scale lengths and thus similar guiding radii, and they differ only in their initial radial velocity dispersions. In the epicyclic approximation, orbits in an axisymmetric potential, with different radial velocity dispersions σ_R but the same guiding radius R_g, have the same azimuthal frequency Ω, where Ω = V_c/R_g, V_c being the circular velocity at R = R_g. Furthermore, when the epicyclic approximation is not valid because the departure from the orbit circularity is significant, kinematically hot and cold populations with similar guiding radii still have similar azimuthal frequencies (with differences typically smaller than 10%), as we show in Fig. C.1 in Appendix C. The distribution of Ω being similar for hot and cold populations, and because in Model 1 the distribution of vertical frequencies ν is the same by construction (all populations have the same thickness), this implies that the fraction of stars that satisfies the inequality ν ≥ 2(Ω − Ω_p) (Merritt & Sellwood 1994), that is, which can respond to the bending/buckling instability that initiates the B/P formation, is similar in all populations for any given value of the bar pattern speed Ω_p. Moreover, because of the initial similarity of Ω and ν at any given radius, the percentage of resonant material with the vertical inner Lindblad resonance (Combes et al. 1990) is similar for all populations. Finally, in Model 1 we observe that the bar formation rapidly increases the radial velocity dispersion of the cold disk to values similar to those of the hot population (see Fig. D.1 in Appendix D), thus erasing their initial kinematic differences. The similarity of Ω and ν and the convergence of the in-plane velocity dispersions to similar values implies similar B/P shapes in all populations, independent of their initial in-plane kinematics.

4. Conclusions

We have made use of dissipationless N-body simulations of boxy/peanut-shaped bulges formed from composite stellar disks that consist of kinematically cold and hot stellar populations to discuss the main drivers of the trends observed in a B/P bulge. To this aim, we have analyzed two extreme models, the first consisting of disk populations with the same initial vertical random motions, σ_z, but different in-plane random motions, σ_R, the second consisting of disk populations with different initial thickness, but equal in-plane random motions. We chose such extreme initial conditions to separately study the role of in-plane versus vertical kinematics in determining the final properties of a B/P bulge.

While a larger sample of N-body simulations is needed to understand how representative our findings are of the two classes of models we studied, we can derive some first conclusions. At first order, disk populations with different initial σ_z and same σ_R are subject to a similar differential mapping into a B/P-shaped bulge as is experienced by stellar populations with equal initial σ_z and different σ_R. In both cases,

• when viewed face-on, the colder the disk population initially, the stronger the bar; when viewed edge-on, the hotter the disk population initially, the thicker the B/P structure.

• at a given height above the plane, the relative contribution of stars in the B/P-shaped bulge depends on their initial kinematics, and the relative weight of stars initially in the cold disk population decreases with height (and vice versa, that of the hot disk population increases with height).

• because of this differential distribution of initially cold- and hot-disk populations in the B/P-shaped bulge, vertical metallicity(/age) gradients are naturally formed along the bulge minor axis, and X-shape metallicity/age distributions as well, if a velocity dispersion-metallicity(/age) relation initially exists in the disk.

As a consequence, the vertical as well as the radial velocity dispersions are both responsible for generating similar trends in a B/P-shaped bulge like that of the Milky Way, and we cannot conclude that in-plane random motions dominate more strongly than vertical random motions in determining these trends, as has recently been suggested.

However, some important differences exist in the characteristics of the B/P bulges and their surrounding disks, as generated with these two models.

• We find that the model where disk populations initially have only different in-plane random motions, but similar thickness, results in a B/P-shaped bulge where all populations have a similar peanut shape, independent of their initial kinematics or metallicity. As a consequence, the peanut shape appears at the same height above the plane for all populations. This is at odds with the trends observed in the Milky Way bulge.

• Finally, a model where the disk stellar populations are initially differentiated only in their in-plane kinematics generates disks with no relation between metallicity/age and vertical velocity dispersion and where all populations converge to the same final radial velocity dispersion. This model is thus not applicable to the Milky Way.

On the basis of these two models, we conclude that a metal-poor, kinematically (radial and vertical) hot component, that is, a thick disk, is necessary in the Milky Way before bar formation.

Acknowledgments

The authors are grateful to the referee for their comments that helped improve the clarity of the manuscript and the presentation of the results. This work has been supported by the ANR (Agence Nationale de la Recherche) through the MOD4Gaia project (ANR-15-CE31-0007, P.I.: P. Di Matteo). FF has been supported by a postdoctoral grant from the center National d’Etudes Spatiales (CNES). SK is supported by a postoctoral grant by ANR. This work was granted access to the HPC resources of CINES under the allocation A0020410154 made by GENCI. PDM and MH acknowledge the generous hospitality of Valerie de Lapparent and the “Galaxies” team at the Institut d’Astrophysique de Paris, where this work has been partially developed.

References

Appendix A: Initial conditions and their effect on the final bulge morphology

In Fig. A.1 the face-on and edge-on stellar density maps of stars in the inner regions are shown for Models 1 and 2. To generate the initial positions of the particles in the two models, we randomly extracted a number N of particles for each disk component according to the corresponding Miyamoto-Nagai density distribution of characteristic radius a and characteristic height h (see Table 1 for the number of particles and characteristic spatial scales adopted for the disks in the two models). For this random realization, a cut in the in-plane, R, and vertical, z, distances has been adopted, so that only particles with positions inside (R_max, z_max) = (10 × a_hot disk, 10 × h_hot disk) are retained. This choice corresponds to a density cut at approximately 10⁻³ times the central stellar density for all disks in the two models (see Fig. 1 in Miyamoto & Nagai 1975). The adoption of these initial cuts explains why in Model 1 the initial vertical extension of all disks is smaller than that of the corresponding disks of Model 2, since in Model 1 all disks have the same characteristic height. The limiting density induced by this cut is low enough to guarantee that no significant artificial effect in the further evolution of the system is observed. As an example of this assertion, the cold disks in Models 1 and 2 initially have both the same mass and same characteristic radii, but differ only in the vertical cut adopted for the generation of the initial conditions, z_max being equal to 2.5 kpc and 9 kpc for Models 1 and 2, respectively. Despite this different initial vertical extension, the final B/P bulge morphology for the two cold populations is remarkably similar in the two models at the final time of the simulations (cf. Figs. 2 and 4 ).

Fig. A.1. Morphology of the bulge region at time t = 0 for Models 1 and 2. First and third column, from top to bottom: Face-on absolute stellar densities of all stars (top panel), and stars initially in the cold (second panel), intermediate (third panel), and hot (bottom panel) disks for Model 1 (Col. 1) and Model 2 (Col. 3). Second and fourth column, from top to bottom: Edge-on absolute stellar densities of all stars (top panel), and stars initially in the cold (second panel), intermediate (third panel), and hot (bottom panel) disks for Model 1 (Col. 2) and Model 2 (Col. 4). In the edge-on maps, only stars with |x| ≤ 5.5 kpc and |y| ≤ 3 kpc have been selected.

As pointed out by Miyamoto & Nagai (1975; see, again, their Fig. 1), and as is also visible in the edge-on maps shown in our Fig. A.1, for high h/a ratios (b/a ratios, in Miyamoto & Nagai (1975) nomenclature), a central bulge-like part appears in the initial edge-on density maps. In the following, we estimate the contribution of this central bulge-like structure to the final B/P morphology of the bulge for the hot-disk component of Model 2. Because of its h/a ratio is the highest, this is indeed the component with the most pronounced central bulge-like part in our models. To estimate the contribution of the central bulge-like component to the final B/P morphology, we selected all particles in the hot disk of Model 2 that initially have |x| ≤ 5.5 kpc and |y| ≤ 3 kpc, corresponding to the final B/P bulge region, and which, in particular, lie inside the bulge-like part of the hot disk, which we define as the region delimited by a density level equal to 0.1 times the central density of the hot disk (see Fig. A.2, left panel). At time t = 0, the fraction of particles inside this bulge-like region constitutes on average about 80% of the whole hot disk inside the final B/P region, but only 15% of the whole stellar distribution. The edge-on distribution of the bulge-like particles at the final time is shown in the second panel of Fig. A.2. It appears that despite their initial bulge-like morphology, these particles do also respond to the B/P instability and show a (weak) peanut-like morphology at the end of the simulation. It is also interesting to note that these particles tend to preferentially populate the innermost regions of the final B/P bulge described by the hot-disk population. This is similar to the finding of Di Matteo et al. (2014) and Gómez et al. (2016), who showed that stars born in the innermost regions of a disk tend to populate the innermost parts of a B/P bulge, while stars born at larger distances tend to be preferentially redistributed in the outskirts of the B/P structure, preserving their Jacoby energy, as has been shown by Martinez-Valpuesta & Gerhard (2013). Overall, the bulge-like component represents a significant fraction of the hot disk in the innermost B/P region, especially for |z| ≤ 1 kpc (see Fig. A.2, third panel). However, it does represent only a marginal fraction of the whole B/P bulge populations, which also includes particles formerly in the cold and intermediate disks (see Fig. A.2, fourth panel). Thus, while this central bulge-like component certainly dominates the morphology of the kinematically hot part of the B/P bulge that consists of hot-disk stars only, it only marginally affects the final whole B/P morphology and its spatial extension.

Fig. A.2. From left to right, first panel: Edge-on stellar density map of particles in the hot disk of Model 2 at t = 0. The density contours trace the whole density distribution of the hot disk in the central regions, the colored map corresponds to the bulge-like part of this distribution (see text for details). Second panel: Edge-on stellar density map at the final time of the simulation of hot-disk particles initially inside the bulge-like region. Third panel: Relative contribution of hot-disk particles initially inside the bulge-like region to the final B/P bulge, as defined by hot-disk particles only. Fourth panel: Relative contribution of hot-disk particles initially inside the bulge-like region to the final B/P bulge as defined by all disk particles, i.e., cold and intermediate disks also included.

Appendix B: Convergence of the iterative method used to generate the initial conditions

To generate initial conditions at equilibrium for the two models, we employed the iterative method described in Rodionov et al. (2009). All disks of Models 1 and 2 have Myamoto-Nagai density profiles, with characteristic scale lengths and heights given in Table 1. The dark matter halos follow a Plummer distribution, with characteristic parameters also given in Table 1. We did not impose any restriction on the kinematics of the system, as was done in some of the examples reported in Rodionov et al. (2009). That is, none of the galactic components was required to converge to any specific velocity dispersion profile. To initialize the procedure, we assigned the initial azimuthal velocities of each disk component in such a way as to have an initial constant (i.e., independent of the radial distance from the center) rotational lag² equal to 0, 30, and 60 km s⁻¹ for the cold, intermediate, and hot disk of Model 1, respectively. This guarantees the different in-plane kinematics of the three disks in Model 1 at the end of the iterative procedure. For Model 2, no initial rotational lag was imposed to any of the disk components, and their different vertical kinematics is the consequence of their different characteristic heights. Twenty iterations were sufficient to obtain initial conditions at equilibrium in the galactic potentials. As an example of this convergence, in Figs. B.1 and B.2 we show the radial, azimuthal, and vertical velocity dispersion profiles of the cold, intermediate, and hot disks in Model 1 and 2, respectively, at each step of the iteration from 1 to 20. In both figures, the insets show the evolution of these profiles for 20 supplementary iterations, from step = 21 to step = 40, to demonstrate that all our models converged after 20 iterations and that no further significant change in the velocity dispersion profiles is observed.

Fig. B.1. From top to bottom: Evolution of the radial, azimuthal, and vertical velocity dispersion profiles as a function of radius R during the iterative procedure adopted to generate the initial conditions of Model 1. For each plot and each inset in the plot, different colors represent different steps of the iteration, as indicated by the the color bar. In all plots and insets, the black dotted curve indicates the velocity dispersion after 20 iterations, when we consider the system converged to an equilibrium solution. Cold, intermediate, and hot disks are shown in the left, middle, and right columns, as indicated.

Fig. B.2. Same as Fig. B.1, but for Model 2.

Appendix C: Initial distribution of azimuthal and vertical frequencies

In Fig. C.1 we show the distributions of guiding radii, azimuthal, Ω, and vertical frequencies, ν, for the two simulations at t = 0, when the disk is still axisymmetric. While Model 2 shows a significant difference in the distribution of the vertical frequencies for the three populations, in Model 1 all distributions appear remarkably similar. As a consequence, in Model 1 the distribution of Ω − ν/2 is also similar for the three disks, implying that the fraction of stars that can respond to the buckling instability, that is, the fraction of stars that satisfies the inequality Ω_p ≥ Ω − ν/2, is similar in the three disks. In the previous inequality, Ω_p indicates the bar pattern speed.

Fig. C.1. From left to right: Distribution of guiding radii, azimuthal frequencies Ω, vertical frequencies ν, and of the differences Ω − ν/2 for Model 1 (top panels) and Model 2 (bottom panels) at t = 0. In each plot, the blue, green, and red curves correspond to the cold-, intermediate-, and hot-disk components, as indicated. In each plot, we also give the 25th, 50th, and 75th percentile of the distributions, with colors corresponding to each of these components.

Appendix D: Evolution of the in-plane and vertical velocity dispersions with time

In this section we show the temporal evolution of the radial and vertical velocity dispersion of the initially kinematically cold and hot disks for Models 1 and 2 in Figs. D.1 and D.2, respectively. In both plots, these velocity dispersions are shown for eight different radial annuli that cover the whole extent of the simulated galaxy, from its inner regions to its outskirts. For clarity, we do not show the intermediate disk in these plots, but we checked that its temporal behavior is always bracketed between that of the initially hot and cold populations. We refer to Sect. 3.1 for a discussion of the implications of these plots on the vertical metallicity gradients outside the B/P bulges of Model 1 and 2, and to the comprehensive works of Aumer et al. (2016), Aumer & Binney (2017) and Aumer et al. (2017) for the evolution of vertical and in-plane random motions in composite thin and thick stellar disks.

Fig. D.1. Top and middle rows: Radial (left panels) and vertical (right panels) velocity dispersions as a function of time for stars initially in the cold (thin curves) and hot (thick curves) disks of Model 1. Velocity dispersions are shown for eight different radial annuli in the modeled disk, as indicated in the top left panel, which are also represented by different colors. In the middle panels, velocity dispersions are normalized to their corresponding value at the initial time of the simulation. Bottom row: Ratio of the radial (left panel) and vertical (right panel) velocity dispersions of stars in the initially hot and cold disks for the eight different radial annuli as a function of time. The dashed line indicates the ratio equal to unity.

Fig. D.2. Same as Fig. D.1, but for Model 2.

All Tables

All Figures

Fig. 1. Initial conditions of our composite stellar disks. First row: Circular velocity (black curve) of the model, together with the rotation curves of the cold (blue), intermediate (green), and hot (red) disks. Second row: In-plane radial velocity dispersion of the cold (blue), intermediate (green), and hot (red) disks as a function of the distance R from the galaxy center. Third row: In-plane tangential velocity dispersion of the cold (blue), intermediate (green), and hot (red) disks as a function of the distance R from the galaxy center. Bottom row: Vertical velocity dispersion of the cold (blue), intermediate (green), and hot (red) disks as a function of the distance R from the galaxy center. Model 1 is shown in the left column, and Model 2 in the right column, as indicated.

In the text

Fig. 2. Model 1, morphology of the bar-B/P region at t = 5 Gyr. First and third column, from top to bottom: Face-on (Col. 1) and edge-on (Col. 3) absolute stellar densities of all stars (top panel), and stars initially in the cold (second panel), intermediate (third panel), and hot (bottom panel) disks. Second and fourth column, from top to bottom: Fraction of cold- (second panel), intermediate- (third panel), and hot-disk (bottom panel) stars in the B/P bulge, seen face-on (Col. 2) and edge-on (Col. 4). The fraction is defined as the ratio of cold-, intermediate-, or hot-disk stars to the total (i.e., cold+intermediate+hot) disk stars in the B/P bulge. In the edge-on maps, only stars in the bar (|x| ≤ 5.5 kpc and |y| ≤ 3 kpc) have been selected.

In the text

Fig. 3. First row: Bar strength A2 as a function of time for all stars (dashed black curve), and for stars initially in the cold (blue curve), intermediate (green curve), and hot (red curve) disks. Second row: Height profiles in the B/P region of all stars, and stars initially in the cold, intermediate, and hot disks. Models 1 and 2 are shown in the left and right columns, respectively. In the panels in the first row, the gray area indicates the time of disk bending. At the end of this instability phase, a B/P-shaped bulge is formed in both models, and its morphology does not change significantly after that time.

In the text

	Fig. 4. Same as Fig. 2, but for Model 2.
In the text

	Fig. 5. Edge-on metallicity and age maps for stars in Model 1 (first and third row) and for Model 2 (second and fourth row) at time t = 0 (first column) and at time t = 5 Gyr (second column). Only stars in the bar (|x| ≤ 5.5 kpc and |y| ≤ 3 kpc) have been selected. In all these maps, the bar when present is oriented side-on. Isodensity contours are shown in black.
In the text

	Fig. 6. Edge-on metallicity and age maps at the final time of the simulation for Model 1 (first and second row), and Model 2 (third and fourth row). The B/P region and the surrounding disk are shown in the maps. Isodensity contours are shown in black.
In the text

Fig. 7. First row: B/P strength, quantified by means of the B6 parameter, as a function of the projected radial length of the peanut shape, RΠ, for Model 1 (left panel) and Model 2 (right panel). Second row: B/P strength as a function of the projected height of the peanut shape, zΠ, for Model 1 (left panel) and Model 2 (right panel). In all plots, the strength of the peanut shape is shown at t = 5 Gyr for all stars (black curves) and for stars initially belonging to the cold (blues curves), intermediate (green curves), and hot (red curves) disks. In all panels, error bars are represented by shaded areas and indicate the intensity rms scatter along each isophote after subtracting the best-fit pure ellipse with harmonic terms up to and including the sixth order.

In the text

Fig. 8. Isodensity contours of the B/P bulge in Model 1 (top panel) and Model 2 (bottom panel) for the initially kinematically cold (blue contours) and hot (red contours) disk populations at t = 5 Gyr. The contours of the intermediate disk are not shown, but they are bracketed by those of the cold and hot disk populations. The numbers at the top left of each panel indicate the values of the isodensity contours for cold and hot populations normalized to their central density. Only stellar particles in the bar region (|x| ≤ 5.5 kpc and |y| ≤ 3 kpc) have been selected for these plots.

In the text

Fig. 9. Distribution of stars along the bar major axis at different heights above the plane, as indicated at the top of each panel. The distributions of cold, intermediate, and hot stars in the B/P structure are indicated with blue, green, and red histograms, respectively. The distribution appears bimodal for all populations at the same height above the plane in Model 1 (left panel), while in Model 2 (right panel), the bimodality appears at larger heights for the hot population than for the colder ones.

In the text

Fig. A.1. Morphology of the bulge region at time t = 0 for Models 1 and 2. First and third column, from top to bottom: Face-on absolute stellar densities of all stars (top panel), and stars initially in the cold (second panel), intermediate (third panel), and hot (bottom panel) disks for Model 1 (Col. 1) and Model 2 (Col. 3). Second and fourth column, from top to bottom: Edge-on absolute stellar densities of all stars (top panel), and stars initially in the cold (second panel), intermediate (third panel), and hot (bottom panel) disks for Model 1 (Col. 2) and Model 2 (Col. 4). In the edge-on maps, only stars with |x| ≤ 5.5 kpc and |y| ≤ 3 kpc have been selected.

In the text

Fig. A.2. From left to right, first panel: Edge-on stellar density map of particles in the hot disk of Model 2 at t = 0. The density contours trace the whole density distribution of the hot disk in the central regions, the colored map corresponds to the bulge-like part of this distribution (see text for details). Second panel: Edge-on stellar density map at the final time of the simulation of hot-disk particles initially inside the bulge-like region. Third panel: Relative contribution of hot-disk particles initially inside the bulge-like region to the final B/P bulge, as defined by hot-disk particles only. Fourth panel: Relative contribution of hot-disk particles initially inside the bulge-like region to the final B/P bulge as defined by all disk particles, i.e., cold and intermediate disks also included.

In the text

Fig. B.1. From top to bottom: Evolution of the radial, azimuthal, and vertical velocity dispersion profiles as a function of radius R during the iterative procedure adopted to generate the initial conditions of Model 1. For each plot and each inset in the plot, different colors represent different steps of the iteration, as indicated by the the color bar. In all plots and insets, the black dotted curve indicates the velocity dispersion after 20 iterations, when we consider the system converged to an equilibrium solution. Cold, intermediate, and hot disks are shown in the left, middle, and right columns, as indicated.

In the text

	Fig. B.2. Same as Fig. B.1, but for Model 2.
In the text

Fig. C.1. From left to right: Distribution of guiding radii, azimuthal frequencies Ω, vertical frequencies ν, and of the differences Ω − ν/2 for Model 1 (top panels) and Model 2 (bottom panels) at t = 0. In each plot, the blue, green, and red curves correspond to the cold-, intermediate-, and hot-disk components, as indicated. In each plot, we also give the 25th, 50th, and 75th percentile of the distributions, with colors corresponding to each of these components.

In the text

Fig. D.1. Top and middle rows: Radial (left panels) and vertical (right panels) velocity dispersions as a function of time for stars initially in the cold (thin curves) and hot (thick curves) disks of Model 1. Velocity dispersions are shown for eight different radial annuli in the modeled disk, as indicated in the top left panel, which are also represented by different colors. In the middle panels, velocity dispersions are normalized to their corresponding value at the initial time of the simulation. Bottom row: Ratio of the radial (left panel) and vertical (right panel) velocity dispersions of stars in the initially hot and cold disks for the eight different radial annuli as a function of time. The dashed line indicates the ratio equal to unity.

In the text

	Fig. D.2. Same as Fig. D.1, but for Model 2.
In the text