3D test particle simulations of the Galactic disks. The kinematical effects of the bar.
Key Words.:
Abstract
Context:
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.
Conclusions:
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 nonaxisymmetric 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 phasespace 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 extraGalactic 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 nonaxisymmetric 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 NBody 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 BulirschStoer (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) ^{1}^{1}1Starting 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 phasespace 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 phasespace 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 phasespace 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 nonaxisymmetric 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
Parameter  

(B1)  
(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,
(1) 
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,
(2) 
where “” stands for “thin” or “thick”. The mass ratio between the two disks is and this gives a thicktothin disk density normalization of near the Sun, similarly to what is observed in our Galaxy (Jurić et al., 2008). We have chosen the MiyamotoNagai 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 MiyamotoNagai 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 MiyamotoNagai 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 MiyamotoNagai 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)
(3) 
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) 
The 3D bar used in our simulations is a Ferrers (1870) bar, whose density profile is given by
(4) 
where . We choose and explore two sets of axis lengths taken from Bissantz & Gerhard (2002) and LópezCorredoira 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
(5) 
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 .
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 semimajor 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 MiyamotoNagai 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
(6) 
(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.,
(7) 
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
(8) 
where is the total potential of the Galaxy. Assuming a velocity ellipsoid aligned with the and axes, i.e., , the first term in the lefthand side of Eq. (8) vanishes, and
(9) 
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 ,
(10) 
(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.
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
(11) 
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
(12) 
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).
Now, let us assume a small nonaxisymmetric 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.
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
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 phasemixing of orbits in 2D and with the quadrupole bar.
As our thin disk corresponds to an intermediateold 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
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 nonaxisymmetric 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 ^{2}^{2}2These 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.
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).
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).
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.
The vertex deviation quantifies the orientation of the velocity ellipsoid on the plane. Following Binney & Merrifield (1998) this is:
(13) 
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.
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
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.
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
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 nonaxisymmetric 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 (semimajor 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 nonaxisymmetry 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 phasespace 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 nonaxisymmetries such as the bar and the spiral arms.
Acknowledgements.
The authors gratefully acknowledge support from the European Research Council under ERC Starting Grant GALACTICA240271.References
 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., RomeroGó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ópezCorredoira et al. (2007) LópezCorredoira, M., CabreraLavers, 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. BlandHawthorn, 139–144
 Zhao & Mao (1996) Zhao, H. & Mao, S. 1996, MNRAS, 283, 1197