3D test particle simulations of the Galactic disks. The kinematical effects of the bar.

3D test particle simulations of the Galactic disks. The kinematical effects of the bar.

G. Monari Kapteyn Astronomical Institute, Rijksuniversiteit Groningen, P.O. Box 800, 9700 AV Groningen, The Netherlands
11email: monari@astro.rug.nl
   T. Antoja Kapteyn Astronomical Institute, Rijksuniversiteit Groningen, P.O. Box 800, 9700 AV Groningen, The Netherlands
11email: monari@astro.rug.nl
   A. Helmi Kapteyn Astronomical Institute, Rijksuniversiteit Groningen, P.O. Box 800, 9700 AV Groningen, The Netherlands
11email: monari@astro.rug.nl
Key Words.:


Aims:To study the imprints of a rotating bar on the kinematics of stars in the thin and thick disks throughout the Galaxy.

Methods:We perform test particle numerical simulations of the thin and thick disks in a 3D Galactic potential that includes a halo, a bulge, thin and thick disks, and a Ferrers bar. We analyze the resulting velocity distributions of populations corresponding to both disks, for different positions in the Galaxy and for different structural parameters of the bar.

Results:We find that the velocity distributions of the disks are affected by the bar, and that strong transient effects are present for approximately 10 bar rotations after this is introduced adiabatically. On long (more realistic) timescales, the effects of the bar are strong on the kinematics of thin disk stars, and weaker on those in the thick disk, but in any case significant. Furthermore, we find that it is possible to trace the imprints of the bar also vertically and at least up to for the thin disk and for the thick disk.


1 Introduction

The study of kinematic groups in the Solar Neighbourhood dates back to the work of Proctor (1869) and Kapteyn (1905), with the discovery of the Hyades and Ursa Major groups. In modern times, the analysis of Hipparcos data and other surveys (Dehnen 1998; Famaey et al. 2005; Antoja et al. 2008) revealed the presence of many other local kinematic structures with varying properties, some of which are composed by stars of a wide range of age. More recently, Antoja et al. (2012) using RAVE data demonstrated that some of the known local kinematical groups extend beyond the Solar neighbourhood at least as far as from the Sun on the plane, and below it.

These main kinematic groups are in the thin disk. Their populations spread (age and metallicity) indicates that they are unlikely to be remnants of disrupted clusters (Eggen 1996). However, they can be explained in terms of the influence of the non-axisymmetric components of the Galaxy (bar and spiral arms) on the orbits of stars (Dehnen 2000; Fux 2001; De Simone et al. 2004; Quillen & Minchev 2005; Antoja et al. 2009, 2011). For instance, the existence of the Hercules stream, a local group of stars moving outwards in the disk and lagging the circular rotation, has been related to the bar’s Outer Lindblad Resonance (Dehnen 2000).

On the other hand, much of the past work on the kinematics of the thick disk has explored the presence of phase-space substructure due to minor mergers and accretion events (Helmi et al. 2006; Villalobos & Helmi 2008, 2009). For example it has been suggested that the Arcturus stream (Eggen 1971) may have an extra-Galactic origin (Navarro et al. 2004). However, it has also been advocated that this stream could have a Galactic origin (Williams et al. 2009; Minchev et al. 2009), and be a signature of a bar’s resonance (Antoja et al. 2009; Gardner & Flynn 2010). Yet this proposal may be considered somewhat speculative, as it is based on simulations of the kinematics on the plane of the Galaxy, that is ignoring the vertical motion, and it is a priori not so clear whether a resonance could affect stars of the thick disk that spend so much time far from the Galactic plane.

We thus currently have a limited understanding as to whether the non-axisymmetric components (bar and spiral arms) can induce kinematic structure in the thick disk, and how to distinguish this from the substructures associated to accretion events. It is only recently that Solway et al. (2012), using controlled N-Body simulations to study radial migration, quantified that transient spiral arms can change the angular momentum of stars in the thick disk almost as much as they do in the thin disk.

This Introduction serves to present the main motivation of our paper, namely to study the influence of the Galactic bar on the thin and thick disk in 3D, paying special attention to the vertical dimension. Our final aim is to establish if kinematic groups or other signatures of the bar may be present far from the Galactic plane, that is for orbits with large vertical oscillations, both in the thin and thick disk. To this end, we perform controlled test particle orbit integrations in a Galactic potential that includes a bar. Our simulations present a number of improvements with respect to previous studies, mainly that the integrations are done in a potential that is 3D, and that we study both the thin and thick disks.

In Section 2 we give details of the simulation techniques and choice of initial conditions. In Section 3 we analyze the resulting velocity distributions of the thin and thick disks in localized volumes around the Solar Neighbourhood, both in the Galactic plane as well as away from it. Finally, in Section 4 we present the implications of our results and our conclusions.

2 Simulation methods

Our approach to study the characteristics of the velocity distribution across the Galaxy is based on a “forward integration technique”. First we make discrete realizations of specific distribution functions, then we integrate these initial conditions forward in time, using an adaptive sizestep Bulirsch-Stoer (Press et al. 1992). Finally we analyze the resulting coarse grained distribution inside finite volumes.

This approach may be contrasted to the “backward integration technique” (e.g., Dehnen 2000; Gardner & Flynn 2010) 111Starting from a time , a point in space and a grid of velocities , the orbits are integrated backward in time. At the end of each integration, each orbit reaches a point , associated with the value of an analytic phase-space distribution function corresponding to the unbarred potential. The collisionless Boltzmann equation states that the fine grained distribution function remains constant along the orbits. It is then possible to associate to each point of the grid of velocities a value of the phase-space density, since , and to see how the velocity space is populated at the point of the Galaxy., where the velocity distribution at a single point in configuration space is obtained. Our method instead requires a large number of particles and is computationally more expensive, but it may be considered more realistic in reproducing how the phase-space is populated and is possibly more appropriate for comparison to the observations, which always probe a finite volume.

All the particles of our simulations are integrated in a potential composed by two parts: an axisymmetric part (formed by a bulge, a thin disk, a thick disk and a dark halo) and a non-axisymmetric part (the bar). The bar is introduced in the potential gradually with time and replaces the bulge. The mass of the bulge is continuously transferred to the bar, and becomes null when the bar has grown completely (see details in Section 2.2)

Here and in the rest of the paper and are the Galactocentric spherical and cylindrical radial coordinate respectively, is the polar angle measured from the long axis of the bar in the direction of Galactic rotation, and the vertical coordinate. Moreover, is the radial Galactocentric velocity (positive towards the center of the Galaxy), is the azimuthal velocity (positive in the direction of the Galactic rotation) and .

2.1 Axisymmetric part of the potential

Table 1: Parameters of the two axisymmetric models A0+B1 and A0+B2.

Below we give a detailed description of the individual components contributing to the axisymmetric part of the potential. Table 1 lists their characteristic parameters.

The bulge follows a Hernquist (1990) potential,


with . We choose two different values for the mass of the bulge (and, therefore, of the bar): (B1) and (B2, see Section 2.2).

We represent the thin and thick disk potentials with Miyamoto & Nagai (1975) models,


where “” stands for “thin” or “thick”. The mass ratio between the two disks is and this gives a thick-to-thin disk density normalization of near the Sun, similarly to what is observed in our Galaxy (Jurić et al., 2008). We have chosen the Miyamoto-Nagai functional form to represent the disks because of their mathematical simplicity which results in computational convenience. Since in reality the Galactic disks follow more closely an exponential form in the density (Binney & Tremaine 2008), we set the characteristic parameters of the Miyamoto-Nagai model to resemble the gravitational potential of two exponential disks with radial scale lengths (for both disks) and vertical scale lengths and . Although the potentials are similar, the forces differ. At the Miyamoto-Nagai thin disk overestimates the exponential disk vertical force up to at , while it underestimates the force by at . The exponential thick disk force in the Miyamoto-Nagai model is underestimated at all heights we consider. These small differences in the force field can presumabxly be regarded as being smaller than the uncertainties in the true form and values of the characteristic parameters of the Milky Way gravitational field.

The dark halo follows the NFW potential (Navarro et al. 1997)


with , and (Battaglia et al. 2005).

The axisymmetric potential consisting of halo, thin and thick disks is referred as “A0” (and is the same for all models explored), while when the bulge is added we refer to A0+B1 or A0+B2 depending on the mass of the bulge.

2.2 Bar potential

Parameter GB LB
(GB1) (LB1)
(GB2) (LB2)
(A0+GB1) (A0+LB1)
(A0+GB2) (A0+LB2)
(A0+GB1) (A0+LB1)
(A0+GB2) (A0+LB2)
(A0+GB1) (A0+LB1)
(A0+GB2) (A0+LB2)
Table 2: Parameters of the bar and location of the main resonances.

The 3D bar used in our simulations is a Ferrers (1870) bar, whose density profile is given by


where . We choose and explore two sets of axis lengths taken from Bissantz & Gerhard (2002) and López-Corredoira et al. (2007), respectively. We dub them “Galactic Bar” (GB) and “Long Bar” (LB), in the same fashion as Gardner & Flynn (2010). The bar parameters are listed in Table 2.

For the pattern speed of the bar we use the value , i.e. in the range of current determinations, which is between and (Gerhard 2011). It is important to notice that the effects on the velocity distribution in a certain volume are mainly set by the proximity to the main resonances. This implies that a faster (slower) pattern speed yields similar effects as the original bar in outer (inner) volumes.

As anticipated, the bar potential grows adiabatically in the simulations. Specifically, each simulation starts with an axisymmetric potential (either A0+B1 or A0+B2), and the mass of the bulge is transferred to the bar using the function


which runs smoothly from to (Dehnen 2000). In particular:

  • for , and ,

  • for , and .

In our models (B1) and (B2), which is consistent with the mass estimates by Dwek et al. (1995) and Zhao & Mao (1996).

In all our simulations the bar grows following Eq. (5), up to time . The choice of is not particularly relevant, as it was shown in Minchev et al. (2010) that a longer growth time of the bar only linearly delays the effects of the bar (in this case, the formation of the kinematic structures). Throughout this paper we always use as in Fux (2001), where .

Figure 1: Circular velocity curves of the components of A0 (halo, thin and thick disk), bulge B2 and the Fourier component of the bar GB2. The solid black line is the resulting circular velocity for model A0+GB2.

In Fig. 1 we show the circular velocity curves of the components of A0 and of the bulge B2. In the same figure we plot the circular velocity of the term of the Fourier decomposition of the bar potential GB2 (GB with ). Note that, although the total mass of the bulge and the bar are the same, inside their velocity curves (blue and purple lines, respectively) differ because the enclosed masses are different. While the Hernquist bulge extends to infinity, all the bar’s mass is confined inside a radius equal to its semi-major axis (see Eq. 4). As a consequence, when the bar growth has been completed, the resulting circular velocity curve has changed. For this reason, for our analysis we use the circular velocity curve obtained by adding the contributions of the dark halo, thin and thick disk and component of the bar potential (A0+GB2). This curve is represented by the black line in Fig. 1. We see how the total rotation curve is mostly influenced by the bar for approximately , by the thin disk for and by the dark halo for . The value of the circular velocity at Solar radius, and for is for A0+GB1 and for A0+GB2.

2.3 Initial conditions

We use two sets of initial conditions for our simulations: ICTHIN and ICTHICK. ICTHIN mimics the typical kinematics and density distribution of the intermediate age population of the thin disk, while ICTHICK represents the thick disk. We generate low (LR) and high resolution (HR) realizations, that have and particles respectively. We perform HR simulations for our standard Galactic model (A0+GB2, see below), to help to distinguish the resonant features that are in the wings of the distribution and often hidden in the noise.

The positions of the particles in both disks are distributed following Miyamoto-Nagai densities, corresponding to the potentials described in Section 2.1. The density distribution of each disk, , is derived from Eq. (2) through Poisson’s equation (for the complete expression see Miyamoto & Nagai 1975). The positions of the particles are randomly drawn from these profiles with a method based on the Von Neumann’s rejection technique (Press et al. 1992) and is generated uniformly between and .

Once the positions are generated, the velocities are assigned in the following way. We describe the radial velocity dispersion , as


(Hernquist 1993), which implies that, far enough from the center, . In the center the smoothing parameter reduces the velocity dispersion. Hernquist (1993) uses . We prefer , which makes the smoothing much more gradual. We relate the tangential velocity dispersion to the radial through the epicyclic approximation, i.e.,


where and are the circular and epicyclic frequencies (Binney & Tremaine 2008).

Finally, the vertical velocity dispersion is obtained by solving the following Jeans equation for an axisymmetric system


where is the total potential of the Galaxy. Assuming a velocity ellipsoid aligned with the and axes, i.e., , the first term in the left-hand side of Eq. (8) vanishes, and


In practice, we compute Eq. (9) on a discrete grid of equispaced points on the meridional plane and, to obtain in a generic point, we linearly interpolate between the nearest grid points.

Each velocity component of a particle is assigned by drawing a random number from a Gaussian distribution with the respective dispersion. The tangential component of the velocity, , is corrected for the asymmetric drift ,


(Binney & Tremaine 2008) where, again, we set .

The initial conditions generated in this way are not fully consistent with the potential. In particular we have assumed the velocity distribution to be locally Gaussian. Therefore we expect these to change in time, while reaching equilibrium with the potential. These transient effects are described in Minchev et al. (2009) and produce kinematic arcs in the Solar Neighbourhood. To avoid confusion between these effects and those induced by the bar we first let our initial conditions evolve in the axisymmetric potential for a time . We find that suffices to reach a stationary state around and inside the Solar radius. We use these evolved initial conditions as effective initial conditions for our simulations with the bar.

Figure 2: Velocity dispersions before and after the evolution of the initial conditions in the A0+B2 axisymmetric potential for and for stars with final . The final velocity dispersions at are also indicated.

The evolution of the velocity dispersions for ICTHIN and ICTHICK is shown in Fig. 2 for the axisymmetric potential A0+B2 (the results for A0+B1 are similar). We see that the radial and tangential velocity dispersions increase in the outer parts of the Galaxy and decrease in the inner regions, both for ICTHIN and ICTHICK. In contrast, the vertical velocity dispersion does not evolve for ICTHIN and only slightly increases for ICTHICK. This is likely because the vertical velocity dispersion is obtained directly by solving the Jeans equation, while the radial dispersion and tangential profiles are only imposed later using Eq. (7) and related through the epicyclic approximation .

An intermediate age population in the thin disk () has local velocity dispersions (Robin et al., 2003; Holmberg et al., 2007). After evolution, our set of initial conditions ICTHIN has comparable characteristics with slightly larger and smaller . The measured local dispersions of the stars in the thick disc are according to Robin et al. (2003) or smaller in the and directions as reported in Bensby et al. (2003). Our initial conditions ICTHICK after the evolution in the imposed potential are consistent with these measurements.

2.4 Orbits and resonances

Let us, for a moment, consider orbits in an axisymmetric potential. In the epicyclic approximation (that is a good description only of orbits that are almost circular and with small vertical motion amplitudes) the azimuthal velocity of a star, located at in the Galaxy, is


where is the radius of the guiding center of its orbit. When higher order terms are considered, depends also on the value of (e.g., Dehnen 1999). Of particular interest are the orbits with such that


where and are two integer numbers. For these orbits the epicyclic frequency resonates with the frequency of the circular orbit at , in the reference frame rotating with pattern speed . In other words, orbits well described by the epicyclic approximation (i.e., not too eccentric) that satisfy this condition are closed in this reference frame. Important for the rest of the paper are, especially, the Outer Lindblad Resonance (OLR, ) where and , the Resonance (), where and , and the Corotation Resonance (), where (i.e., ). In Fig. 3 we show the positions of the main resonances for the A0 +GB2 potential but considering only the Fourier component of GB2, for (Table 2).

Figure 3: Positions of the main resonances for the model A0+GB2 ( Fourier component) with (horizontal dashed line). The position is given by the at which the curves cross the horizontal line.

Now, let us assume a small non-axisymmetric perturbation (the bar) to the axisymmetric potential, rotating with pattern speed . If the perturbation is small, the motion of most stars will still be approximately described by the epicyclic theory. However, some stars will be more strongly affected by the perturbation, namely those that satisfy Eq. (12) in the axisymmetric potential. As stated by Binney & Tremaine (2008), at the resonances the “perturbation is acting with one sign for a long time. If the effects of a perturbation can accumulate for long enough, they can become important, even if the perturbation is weak”.

2.5 Vertical motion

The simple dynamical picture that we described in the previous section is complicated in two ways in the Milky Way: Galactic disks are 3D structures, i.e. stars have vertical motions, and especially the thick disk population is kinematically hot.

The epicyclic approximation assumes that the orbits of the stars are not too eccentric and that the amplitudes of the vertical motions are so small that horizontal and vertical motions are nearly decoupled. In an axisymmetric potential the motion is exactly decoupled when the radial and vertical force ( and ) are respectively functions of and only. This is only true very near to the Galactic plane. In Fig. 4 we have plotted the radial (top) and vertical (bottom) forces for the total potential (colour) and the contribution of the bar ( terms of the Fourier decomposition, black solid lines) for model A0+GB2. The bottom panel shows, as just discussed, that is independent of only at small (where it runs parallel to the cylindrical radius axis). However, already beyond , it varies significantly with , over the range explored by many orbits that visit the Solar Neighbourhood. On the other hand, the top panel shows that the isocontours are less curved and that this force decreases more slowly with (and this is even more true if we consider the bar’s radial force only as indicated by the black curves). For this reason the radial motion may be expected to be more similar with height, than the vertical motion with radius. Furthermore, it is important also to note that orbits with angular momentum and vertical oscillations larger than live in a range of that is hundred of parsecs different from that of a planar orbit with the same angular momentum (see Binney & McMillan 2011).

Also interesting is that for all our models the bar induces horizontal forces that, at , are almost independent on (top panel of Fig. 4, solid black curves). In fact, the relative importance of the perturbation due to the bar slightly increases with , as a consequence of the decreasing radial force of the axisymmetric part of the potential. Therefore, we expect stars far from the Galactic plane to experience a slightly stronger relative effect of the bar compared to those on the plane.

The exact details of how the vertical motions, particularly for the populations with larger kinematic temperature, affect the simple dynamical picture previously sketched, require however, as especially the bottom panel of Fig. 4 shows, suitable and fully 3D orbital integrations.

Figure 4: Isocontours of the (top panel) and (bottom panel) forces for the the potential A0+GB2 (colors) and terms of the Fourier decomposition of GB2 (black solid lines). The meridional plane is shown at (along the long axis of the bar).

3 Results

Motivated by the above discussion, we turn to numerical simulations to study how the bar impacts the kinematics of stars near the Sun and neighbouring volumes.

We conventionally place the Sun on the Galactic plane, at and we consider cylindrical volumes throughout the Galaxy of radius and height , centered at some . For each volume we study the vs. velocity distribution.

3.1 General effects of the bar and time evolution

Figure 5: Evolution in time of the kinematics of the simulated LR thin disk in the Solar Neighbourhood volume for the potential A0+GB2. The dashed horizontal lines correspond to (yellow), (red), (purple) and (blue). The density field is obtained with the modified Breiman estimator, described in Ferdosi et al. (2011), with optimal pilot window . The contours enclose of the probability, where . In the lower left part of the panel we indicate the number of particles in the volume.

In Fig. 5 we show the evolution in time of the kinematics in the Solar Neighbourhood volume for the model A0+GB2 and the LR thin disk initial conditions (LR ICTHIN). Specifically, to describe the kinematics we use a probability density function of , determined with the modified Breiman estimator, described in Ferdosi et al. (2011), with optimal pilot window . We show density contours that enclose (from inside to outside) of the probability, where . The -th contour from inside out corresponds, therefore, to the contour. The evolution seen in Fig. 5 is qualitatively similar for all the other potentials (and respective values of the parameters) explored. As we can see, the form of the distribution changes with time (see Minchev et al. 2010 for a colder disk). Initially, the distribution of stars is featureless, symmetric in , and with a tail of stars at low due to the asymmetric drift. At the distribution has changed completely: it is split up in two regions and presents other deformations that make it asymmetric, both in and .

Before proceeding with the analysis of the evolution in this plane, it will be useful for the rest of the paper to introduce a terminology describing the features seen in velocity space. In Fig. 5 we have indicated the of four important families of resonant orbits, using Eq. (11) for : (yellow dashed line), (red dashed line), (purple dashed line) and (blue dashed line). As we can see, the diagonal valley in the velocity distribution corresponds roughly to , confirming the large effect of the Outer Lindblad Resonance on the kinematics of the stars. Following Dehnen (2000), we call the peak in the density of stars above the valley the “LSR Mode”, and the one below the “OLR Mode”. Furthermore, in the right part of the LSR Mode, at , there is an elongation in , which we call “Horn”.

Previous studies have used default integration times of (Dehnen 2000). We notice here that between and there is still significant evolution of the distribution, which only reaches a stationary configuration later. For instance, the valley between the LSR Mode and the OLR is progressively filled with stars. The OLR and LSR Modes change shape in time.

In general, it is clear from Fig. 5 that all the features associated to the bar are most prominent for short integration times. The features that we find in the velocity distribution for these times are similar to those observed by Dehnen (2000) in 2D integrations with a quadrupole bar. This is the first time that they are seen in 3D simulations. On the other hand, the time evolution confirms, again in 3D, the results obtained by Fux (2001) in his study of the phase-mixing of orbits in 2D and with the quadrupole bar.

As our thin disk corresponds to an intermediate-old population of stars, it is unlikely that they experienced the effects of the bar only for (). On the other hand, Cole & Weinberg (2002), using properties of infrared carbon stars from the Two Micron All Sky survey, state that the Milky Way Bar is likely to be younger than , fixing an upper limit for its formation to ago (however studies of the bulge indicate that it is dominated by a stellar population older than , e.g., Minniti & Zoccali 2008). Given this, we decide to choose a default integration time of , when the kinematics have roughly reached stationarity. Note that we are effectively assuming that the properties of the bar, such as its pattern speed did not change much during this time.

3.2 Thin disk

Figure 6: The effect of different bar potentials on the kinematics of the simulated thin disk in 3 volumes centered at an angle and for .

Fig. 6 shows the kinematics of the thin disk in the Solar Neighbourhood and in two nearby volumes aligned with the Solar Neighbourhood () but at different radii ( and ) for various potentials and LR simulations. The first and second column show the distributions with the bar potential with and , respectively. Not surprisingly, the effect of the least massive bar is weaker compared to a bar twice as massive. The density contrast between the OLR Mode and the LSR Mode is smaller, and the Horn is slightly less prominent for the least massive bar. The overall distribution has also a more symmetric appearance. These small differences are explainable by simply noting that the force of the perturbation only slightly smaller in the case of the least massive bar. The non-axisymmetric part of the force (i.e. excluding the monopole term associated to the bar) differs only by at between the GB1 and GB2. Some of the differences between the distributions may also be ascribed to the difference between the circular velocity curves of A0+B1 and A0+B2, which have the resonances at slightly different positions. For instance, for , for A0+GB1 and for A0+GB2 (see Table 2). Since the resonances are mostly responsible for the effects of the bar on the kinematics of the stars, we expect some differences at fixed 222These differences may also reflect a dissimilar evolution of the distributions. As explained in Minchev et al. (2010), the typical time of libration of the stars trapped to the resonances depends on the bar strength; since the shape of the distributions is defined mostly by the resonances, we may expect, at the same snapshot, different evolutionary stages for different bar strengths..

The second and third columns of Fig. 6 allow us to compare more massive bar models () with different geometrical properties (GB and LB). The Long Bar (third column) has slightly different effects on the thin disk velocity distributions than the less elongated Galactic Bar. It produces sharper features such as a Horn with a more defined lower edge, especially for and . This was already noticed by Gardner & Flynn (2010) when they compared the responses to these two Ferrers bar models, but we confirm here their results in our 3D simulations. Nevertheless, we stress that, for this long integration time, the differences between the effect of the different bars on the kinematics inside the various volumes explored is much less significant than for shorter integration times. Furthermore, for short integration times the structures are much more evident in all our bar models.

Given the uncertainties in having a satisfactory model for the bar (recent studies even suggest that it could be a superposition of something similar to our Galactic Bar and Long Bar, see Robin et al. 2012), we use from now on as a default model the Galactic Bar and the potential A0+GB2, confident that the effects will be similar for all the other models, especially for long integration times, as illustrated in Fig. 6.

Figure 7: Kinematics in the Solar Neighbourhood at for two subsets of the thin disk HR simulation with A0+GB2 potential. Left panel: the most planar particles in the simulation. Right panel: a random subsample of the simulation with the same number of stars.

Since most works so far have used 2D simulations, it is important to understand if significant differences exist for the kinematics near the plane in the 2D vs 3D simulations. To this end we consider a subset of the HR 3D thin disk simulation in the potential A0+GB2, constituted by those particles that, at , have and . These are the particles whose vertical oscillations have the smallest amplitude in the simulation: at their dispersion is still less than and their dispersion less than . It is then reasonable to expect they behave in a similar fashion to particles in a 2D simulation, distributed on the Galactic plane with similar density and kinematics at , and subject to the same potential A0+GB2.

Fig. 7 shows the kinematics of the subset of most “planar” particles (left panel) and the kinematics of a random subsample of the whole 3D simulation with the same number of stars (right panel), in the Solar Neighbourhood volume at the default integration time. We note that there are no large differences between the velocity distributions for these two data sets. The two sets do have slightly differing velocity dispersions: the planar particles are slightly colder (, while for the 3D subsample), but this is a reflection of the initial conditions. This agreement stems likely from the fact that the planar and vertical motions are decoupled for most thin disk particles. Since our 3D thin disk is particularly cold in ( at ), it is also possible that this is more true in our simulations than in reality since the observed local (Robin et al., 2003; Holmberg et al., 2007).

Figure 8: Kinematics of the simulated HR thin disk in the potential A0+GB2 in different volumes at . Different columns correspond to volumes centered at different angles from the bar in the direction of antirotation, from (leftmost column) and with step . Different rows correspond to volumes centered at different Galactocentric radii, from (topmost row) and with step . The dashed lines correspond to (yellow), (red), (purple) and (blue). Notice how for the volumes in the bottom row and are superposed.

In Fig. 8 we show the kinematics of the thin disk in the vicinity of the Solar Neighbourhood at and in HR. From left to right, the columns correspond to volumes centered at angles , and (aligned with the bar). From bottom to top, the rows correspond to different Galactocentric radii from to in steps of . The yellow, red, purple and blue dashed lines represent of the main resonances (as in Fig. 5), obtained by using Eq. (11) with the of the corresponding volume. This figure shows that the main kinematic features associated to the bar vary as a function of position in the Galaxy, similarly to the effects already reported in the literature for 2D simulations. First, we note that the main features shift in as varies. At (first row), is in the low tail of the distribution (), the OLR Mode has a small extension, and the Horn is placed at the lower right edge of the contour (-th contour). For and , and in consequence the separation between LSR Mode and OLR Mode, lies at the center of the distribution (), and the LSR Mode and OLR Mode become clearly separated. For (bottom row), the LSR mode can hardly be seen and the Horn is now at the upper right edge of the distribution. Since these features are mostly associated with the resonant orbits, their shift can be explained in light of Eq. (11), for each of the resonant families ( constant).

Figure 9: Average of thin disk particles in volumes centered at as a function of for the HR simulation A0+GB2 , at (solid red). For comparison we show with the solid black line the average over all angles for the axisymmetric initial conditions ICTHIN, and with dashed lines twice the standard deviation around this average. The top and bottom panels are for and respectively. The blue line in the top panel represents the average of thin disk particles in volumes centered at as a function of for the LR simulation A0+GB2, at (solid blue).

We can also observe a change in the mean the of the OLR Mode as the orientation of the volume with respect to the bar () is varied. While the OLR Mode is almost completely symmetric with respect to the axis for , it moves towards for more negative . At the OLR Mode is completely in the half plane at almost all radii. We quantify this change in Fig. 9. The red solid line shows the mean for particles with , in volumes at as function of (plotted in steps ). The black solid line represents the same quantity, but for the axisymmetric initial conditions and averaged between all the angles. The dashed lines represents twice the standard deviation of the set of means obtained for all the angles. This figure shows that the OLR Mode varies periodically in as a function of , and that the amplitude of this oscillation is significantly above the noise level. Notice however how the curve is not perfectly symmetric: there is a wiggle for , and the amplitude is different between positive and negative . Because of the symmetry of the potential, we would expect instead , for each . This happens because the disk has not reached a fully stationary equilibrium configuration at at this time (Fux 2001; Mühlbauer & Dehnen 2003). A much longer (and unrealistic) integration time of () is required to reach complete symmetry (Fig. 9, blue line). For volumes approximately aligned or perpendicular to the bar the OLR Mode has a null mean, that is the OLR Mode is centered. For these cases, however, the behaviour is still significantly different from the axisymmetric case, as one can see directly in the right column of Fig. 8, at , where the bimodality, despite being centered with respect to , is still present, and as shown in Fig. 10. In this figure we restrict the attention to those stars with , where the differences from the axisymmetric initial conditions are enhanced.

Figure 10: The distribution of thin disk stars in for the volume centered at and for stars with . The red line represents the simulation A0+GB2, at . The solid black line represents the axisymmetric initial conditions.
Figure 11: Vertex deviation in volumes centered at as a function of . The description of the lines is the same of Fig. 9.

The vertex deviation quantifies the orientation of the velocity ellipsoid on the plane. Following Binney & Merrifield (1998) this is:


The corresponding values are quoted in all panels of Fig. 8. In the Solar Neighbourhood vicinity is generally small. This means that the velocity ellipsoid is only slightly misaligned with the respect to the and axis. The vertex deviation behaviour is fully quantified in Fig. 11. For , we see that is positive for , null for (volumes aligned with the bar), and negative for , until it changes again sign at . This is consistent with the results of Mühlbauer & Dehnen (2003), who showed that the bar can cause significant periodic perturbations in the velocity moments throughout the Galaxy for the 2D case.

Figure 12: Kinematics of the thin disk HR simulation in the potential A0+GB2 in volumes centered at and and and (first, second and third column respectively).

Our 3D simulations offer us the unique opportunity to trace the disk kinematics far from the Galactic plane. In Fig. 12 we show the velocity distribution of the thin disk for the default bar model, for volumes at , at angles (top row) and (bottom row), but centered at different . For both angles it is possible to trace the characteristic features (LSR and OLR Mode, Horn) present at , at least up to . This is confirmed by the periodic oscillation of the average of the OLR Mode which is significant, even at this height. Periodicity is present also at (see lower panel of Fig. 9).

We mentioned in earlier sections how the relative strength of the bar slightly increases with . It might be then not surprising to see effects of the bar at these heights. However the volumes at these heights are populated very differently than the ones located on the plane: the orbits are more eccentric, they have larger vertical amplitudes and they rotate more slowly in average (). It seems, however, that the differences in the velocity distributions at various heights are not much driven by the different types of orbits but mostly by the much smaller numbers of stars (and related Poisson noise) far from the plane.

3.3 Thick disk

Figure 13: Kinematics of the simulated thick disk in the potential A0+GB2 in different volumes at . As in Fig. 8 the columns correspond to volumes centered at different angles, and the rows to volumes centered at different Galactocentric radii. The dashed lines indicate to (yellow), (red), (purple) and (blue).

Even in a hot and vertically extended population like the thick disk, it is still possible to recognize the dynamical imprints of the bar on its kinematics. Fig. 13 (equivalent for the thick disk to Fig. 8) shows that, for , we can find the Horn in all volumes (see the contour) at larger for decreasing . At the same time, we can recognize the OLR Mode as a distinct feature with the respect to the LSR mode in several volumes. Specifically, for and the fourth contour shows a lack of stars at positive in the range .

Fig. 13 also confirms that the bar effects are systematic and vary with the location of the volumes in the Galaxy. We see that the main structures associated with the OLR Resonance (Horn, separation between LSR and OLR Mode), are still present as deformations in the 3rd, 4th and 5th contours of the distributions. Fig. 14 (equivalent to Fig. 9 for the thin disk) indicates that, also in this case, the of the OLR Mode varies significantly, as a function of . However, in the cases of volumes aligned with the bar ( and ), where as in the axisymmetric case, the imprints of the bar on the distribution are not so clear as in the thin disk case, as Fig. 15 shows.

Figure 14: Average of thick disk particles in volumes centered at as a function of , for simulation A0+GB2, at (solid red). As in Fig. 9 we show for comparison the average (solid black) and twice the standard deviation (dashed black) for the axisymmetric initial conditions ICTHICK. The top, middle and bottom panels correspond to , and respectively.
Figure 15: The -distribution of thick disk stars for the volume centered at and for stars with . The red line represents the simulation A0+GB2, at . The solid black line represents the axisymmetric initial conditions.
Figure 16: Kinematics of the thick disk simulation in the potential A0+GB2 in volumes centered at and and and (first, second and third column respectively). Notice that for , we can double the resolution using the volumes centered at .

The vertical coherence of the kinematic structures in the thick disk is depicted in Fig. 16. The vertical extent of the thick disk allows us to have enough stars, even at . The distributions are similar on the plane and for and . The Horn and the OLR Mode are present at the two angles, visible especially in the 2nd, 3rd and 4th contours. This is confirmed by the analysis of the mean of the OLR Mode (shown in the lower panels of Fig. 14). Even for , the periodicity of the mean of the OLR Mode is still noticeable.

We have repeated these analyses on the OLR Mode, both for the bar models GB1 and LB2, and found that the same periodicity can be observed, with a larger amplitude of the oscillations in the LB2 case.

3.4 Combination of thin and thick disk

Figure 17: Kinematics of the thin and thick disk simulations together in the potential A0+GB2, in volumes centered at and and and (first, second and third column respectively). This is obtained combining thin disk particles with a a random subsample of thick disk particles, so that the thick-to-thin density normalization at is .

In practice, surveys generally have a mix of both the thin and thick disk populations. This is why in Fig. 17 we show the velocity distribution of the thin and the thick disk together. Since the thick disk has of the mass of the thin disk (Section 2.1), we combine the particles of the thin disk for the A0+GB2 potential, with a random subsample of thick disk particles run in the same potential. As a result, in the volumes on the Galactic plane the number of stars of the thick disk is that of the thin disk. Fig. 17 shows how, in a sample of stars selected independently of the parent disk population, the imprints of the dynamics of the bar are easily detectable. In particular, we see how for the distribution completely resembles that of the thin disk, while the thick disk, dominates the distribution at . All the main features, described in the previous figures (Horn, OLR Mode, LSR Mode) are easily detectable at every height.

4 Summary and discussion

In this work we have studied the 3D response of test particles, representing the thin and thick disk component of the Milky Way, to the Galactic gravitational field including a rotating central bar. Our goal was to understand if this non-axisymmetric component could induce moving groups, especially in the case when the finite vertical thickness of the Galaxy is taken into account.

We have found that the imprints of the bar are apparent both in a population resembling the intermediate/old population of the thin disk of the Milky Way and in a hot and vertically extended population resembling the thick disk. Thanks to our 3D simulations, this is the first time that such imprints have been traced far from the Galactic plane, up to in the thin disk and up to in the thick disk. These results are not strongly dependent on the bar models used, as all the simulations explored with different structural parameters (semi-major axes, vertical thickness and masses) yield similar results.

The effects of the bar that are seen in our simulations are clearly related to the resonant interaction between the rotation of the bar and the orbits of the stars in the disk. Many of the stellar particles in the vicinity of the Sun have orbits strongly affected by the bar’s Outer Lindblad Resonance of the bar. This OLR is apparent in a splitting of the simulated velocity distributions in two main groups. Our simulations also show that the impact of the Outer Lindblad Resonance varies with position in the Galaxy, depending on Galactocentric radius and angle from the major axis of the bar. On a larger scale, the characteristics of the velocity distributions are periodic with respect to the orientation angle of the bar, tracing the symmetries of the bar’s potential. This is manifested for instance, on the vertex deviation around the Solar radius and in the mean radial velocity of the Outer Lindblad Resonance groups.

In this sense our work agrees with previous results on the effects of the Outer Lindblad Resonance of the bar on the kinematics of stars in the Solar Neighbourhood (Dehnen 2000; Fux 2001; Mühlbauer & Dehnen 2003; Gardner & Flynn 2010). The main difference is that our simulations are more realistic as they are 3D and incorporate both a thin and a thick disc. For the analyses we have considered volumes of finite spatial extent (as in the observations) in contrast to the previous studies, which studied velocity distributions in a single point in the configuration space. Due to this difference as well as our longer integration times, we find more diffuse kinematic structures.

The presence of structures at large heights above the plane may also be understood from the fact that near the Sun the forces along the polar radius direction do not vary strongly with height, and certainly much less than the vertical forces. This implies that resonant features in the velocity plane, even at large heights from the plane, are expected. Our 3D simulations have shown that even in these regions of the Galaxy where a large fraction of stars have horizontal and vertical motions that are not decoupled, and the orbital eccentricity is large (especially for the thick disk), the bar significantly transform the velocity distribution, in a similar way as on the plane.

We detect strong transient effects for approximately 10 bar rotation periods after the bar is introduced. Short integration times, like in e.g., Dehnen (2000), produce much clearer bar signatures. However, we believe that it is unlikely that an old population would have experienced the bar for such a short time scale.

Our simulations offer an explanation for the discovery that several of the observed kinematic groups in the Solar Neighbourhood are also present far below the plane of Galaxy (, Antoja et al. 2012). If some of these kinematic groups are direct consequence of the bar gravitational force (like, e.g., the Hercules group seems to be), our results show that it is not surprising to be able to recognize them far from the Galactic plane. Moreover, our simulations predict that, depending on the location of the Sun in the Galaxy and with the respect of the bar, some imprints of the bar should be recognizable even beyond these heights.

Although we find kinematic structures in the thick disk, the Arcturus stream does not appear in our simulations for our default (long) integration times. This moving group is expected at in the Solar Neighbourhood of our simulations ( in Williams et al. 2009). However, for short integration times and the long bar model, we can identify the same structures and distortions in the low tail of the kinematic distribution in our thin disk, already associated to Arcturus by Gardner & Flynn (2010). Our simulations seem to suggest that, if the Arcturus stream is an imprint of the bar, the bar should be very young. However, the lack of overdensity in the Arcturus region could also be a reflection of the choice of the initial conditions of our simulations if we are to associate it to a non-axisymmetry of the potential rather than an accreted origin.

One important assumption of this work is that the the pattern speed (and other properties) of the bar has not changed in the last . A change in the pattern speed should result in a change in the resonances, and in which stars are affected. But whether this evolution would produce more or less bar signatures on the disc phase-space remains to be seen.

Finally, as shown in Solway et al. (2012), the spiral arms induce angular momentum changes that affect not only the thin but also the thick disc. This suggests that spiral arms may also produce kinematic signatures on the thick disc like those induced by the bar and studied here. It would be interesting to quantify the relative importance of the effects of spiral arm resonances with respect to the bar’s.

In our upcoming paper we will present a more exhaustive quantification of the substructures and the bar effects and a dynamical interpretation of the results based on the orbital frequency analysis.

The results obtained in this study unveil a more complex and structure rich thick disk, which has been likely affected not only by external accretion events but also by secular evolution induced by the disk non-axisymmetries such as the bar and the spiral arms.

The authors gratefully acknowledge support from the European Research Council under ERC Starting Grant GALACTICA-240271.


  • Antoja et al. (2008) Antoja, T., Figueras, F., Fernández, D., & Torra, J. 2008, A&A, 490, 135
  • Antoja et al. (2011) Antoja, T., Figueras, F., Romero-Gómez, M., et al. 2011, MNRAS, 418, 1423
  • Antoja et al. (2012) Antoja, T., Helmi, A., Bienayme, O., et al. 2012, MNRAS, 426, L1
  • Antoja et al. (2009) Antoja, T., Valenzuela, O., Pichardo, B., et al. 2009, ApJ, 700, L78
  • Battaglia et al. (2005) Battaglia, G., Helmi, A., Morrison, H., et al. 2005, MNRAS, 364, 433
  • Bensby et al. (2003) Bensby, T., Feltzing, S., & Lundström, I. 2003, A&A, 410, 527
  • Binney & McMillan (2011) Binney, J. & McMillan, P. 2011, MNRAS, 413, 1889
  • Binney & Merrifield (1998) Binney, J. & Merrifield, M. 1998, Galactic Astronomy
  • Binney & Tremaine (2008) Binney, J. & Tremaine, S. 2008, Galactic Dynamics: Second Edition, ed. Binney, J. & Tremaine, S. (Princeton University Press)
  • Bissantz & Gerhard (2002) Bissantz, N. & Gerhard, O. 2002, MNRAS, 330, 591
  • Cole & Weinberg (2002) Cole, A. A. & Weinberg, M. D. 2002, ApJ, 574, L43
  • De Simone et al. (2004) De Simone, R., Wu, X., & Tremaine, S. 2004, MNRAS, 350, 627
  • Dehnen (1998) Dehnen, W. 1998, AJ, 115, 2384
  • Dehnen (1999) Dehnen, W. 1999, AJ, 118, 1190
  • Dehnen (2000) Dehnen, W. 2000, AJ, 119, 800
  • Dwek et al. (1995) Dwek, E., Arendt, R. G., Hauser, M. G., et al. 1995, ApJ, 445, 716
  • Eggen (1971) Eggen, O. J. 1971, PASP, 83, 271
  • Eggen (1996) Eggen, O. J. 1996, AJ, 112, 1595
  • Famaey et al. (2005) Famaey, B., Jorissen, A., Luri, X., et al. 2005, A&A, 430, 165
  • Ferdosi et al. (2011) Ferdosi, B. J., Buddelmeijer, H., Trager, S. C., Wilkinson, M. H. F., & Roerdink, J. B. T. M. 2011, A&A, 531, A114
  • Ferrers (1870) Ferrers, N. M. 1870, Royal Society of London Philosophical Transactions Series I, 160, 1
  • Fux (2001) Fux, R. 2001, A&A, 373, 511
  • Gardner & Flynn (2010) Gardner, E. & Flynn, C. 2010, MNRAS, 405, 545
  • Gerhard (2011) Gerhard, O. 2011, Memorie della Societa Astronomica Italiana Supplementi, 18, 185
  • Helmi et al. (2006) Helmi, A., Navarro, J. F., Nordström, B., et al. 2006, MNRAS, 365, 1309
  • Hernquist (1990) Hernquist, L. 1990, ApJ, 356, 359
  • Hernquist (1993) Hernquist, L. 1993, ApJS, 86, 389
  • Holmberg et al. (2007) Holmberg, J., Nordström, B., & Andersen, J. 2007, A&A, 475, 519
  • Jurić et al. (2008) Jurić, M., Ivezić, Ž., Brooks, A., et al. 2008, ApJ, 673, 864
  • Kapteyn (1905) Kapteyn, J. C. 1905, Reports of the British Association for the Advancement of Science, 13, 264
  • López-Corredoira et al. (2007) López-Corredoira, M., Cabrera-Lavers, A., Mahoney, T. J., et al. 2007, AJ, 133, 154
  • Minchev et al. (2010) Minchev, I., Boily, C., Siebert, A., & Bienayme, O. 2010, MNRAS, 407, 2122
  • Minchev et al. (2009) Minchev, I., Quillen, A. C., Williams, M., et al. 2009, MNRAS, 396, L56
  • Minniti & Zoccali (2008) Minniti, D. & Zoccali, M. 2008, in IAU Symposium, Vol. 245, IAU Symposium, ed. M. Bureau, E. Athanassoula, & B. Barbuy, 323–332
  • Miyamoto & Nagai (1975) Miyamoto, M. & Nagai, R. 1975, PASJ, 27, 533
  • Mühlbauer & Dehnen (2003) Mühlbauer, G. & Dehnen, W. 2003, A&A, 401, 975
  • Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493
  • Navarro et al. (2004) Navarro, J. F., Helmi, A., & Freeman, K. C. 2004, ApJ, 601, L43
  • Press et al. (1992) Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical recipes in FORTRAN. The art of scientific computing
  • Proctor (1869) Proctor, R. A. 1869, Royal Society of London Proceedings Series I, 18, 169
  • Quillen & Minchev (2005) Quillen, A. C. & Minchev, I. 2005, AJ, 130, 576
  • Robin et al. (2012) Robin, A. C., Marshall, D. J., Schultheis, M., & Reylé, C. 2012, A&A, 538, A106
  • Robin et al. (2003) Robin, A. C., Reylé, C., Derrière, S., & Picaud, S. 2003, A&A, 409, 523
  • Solway et al. (2012) Solway, M., Sellwood, J. A., & Schönrich, R. 2012, MNRAS, 422, 1363
  • Villalobos & Helmi (2008) Villalobos, Á. & Helmi, A. 2008, MNRAS, 391, 1806
  • Villalobos & Helmi (2009) Villalobos, Á. & Helmi, A. 2009, MNRAS, 399, 166
  • Williams et al. (2009) Williams, M. E. K., Freeman, K. C., Helmi, A., & RAVE Collaboration. 2009, in IAU Symposium, Vol. 254, IAU Symposium, ed. J. Andersen, Nordströara, B. m, & J. Bland-Hawthorn, 139–144
  • Zhao & Mao (1996) Zhao, H. & Mao, S. 1996, MNRAS, 283, 1197
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description