EDP Sciences
Free Access
Issue
A&A
Volume 622, February 2019
Article Number L6
Number of page(s) 7
Section Letters to the Editor
DOI https://doi.org/10.1051/0004-6361/201834707
Published online 30 January 2019

© ESO 2019

1. Introduction

After two successful data releases, Gaia already provided to the community an unprecedented amount of high-quality knowledge about kinematics and spatial structures of the Milky Way disk and its surroundings (Gaia Collaboration 2016, 2018a). In particular, the second data release (DR2) contains the full six-dimensional phase-space information for 7.2 million stars, which made possible several discoveries about kinematics and the spatial structure of our Galaxy (Gaia Collaboration 2018b). It has become abundantly clear that nonequilibrium processes in the Milky Way such as satellite interactions, bar dynamics, and spiral structure effects have shaped the present-day structure and kinematics of the Milky Way significantly. One of the most striking new features discovered is the spiral-like distributions in the z − νz plane, which suggest that the Milky Way disk is locally out of dynamical equilibrium. In particular, Antoja et al. (2018) found spirals in the distribution of radial and azimuthal velocities in z − νz plane for stars that lie in a small volume around the Sun (< 0.1 kpc). By combining Gaia DR2 with spectroscopic surveys such as GALAH (Bland-Hawthorn et al. 2018) and LAMOST (Tian et al. 2018), this result has been quickly confirmed to hold also in a more extended volume of the Milky Way disk.

Although we already know some signatures of the nonequilibrium state of the Milky Way disk in the solar suburb (see, e.g., Siebert et al. 2011; Widrow et al. 2012; Williams et al. 2013; Bovy et al. 2015; Carlin et al. 2013; Pearl et al. 2017; Carrillo et al. 2018), the exact mechanism driving the phase-space z − νz spirals is still poorly constrained. For instance, Widrow et al. (2012) proposed that the corrugated velocity patterns in Milky Way-type galaxies can be generated by a passing satellite or dark matter subhalos (see also Feldmann & Spolyar 2015; D’Onghia et al. 2016). In particular, a slow moving satellite induces a bending-mode perturbation. However, the oscillations decay rapidly in time because of phase mixing and Landau damping, which lead to disk heating. In such a scenario, the velocity patterns seen in the observations should indicate a very recent tidal interaction (see also Widrow et al. 2014). In such a framework Antoja et al. (2018) interpreted the phase-space spirals as a signature of continuous ongoing mixing in the vertical direction, possibly excited by a recent encounter with the Sagittarius dwarf galaxy. By making use of test particle simulations, these authors estimated that the vertical phase mixing event started about 300 − 900 Myr ago. This suggestion has been confirmed qualitatively in a low-resolution N-body simulation of the interaction of Sagittarius with the Galactic disk (Laporte et al. 2018). More detailed analysis has been carried out by Bland-Hawthorn et al. (2018) who show that to produce a realistic phase-space spiral, the mass of perturber should be higher than 1010 M. Such value is in agreement with another recent dynamical study of the spirals in the (z, νz)-plane by Binney & Schönrich (2018), who noted, however, that the main weakness of such an external perturbation scenario is in the high values of mass and passage duration, which are needed to obtain a realistic model of the phase-space spirals.

Another mechanism proposed to explain the continuous phase-space mixing is the spontaneous generation of bending waves in isolated disk galaxies (Darling & Widrow 2019). By making use of test particles simulations and N-body models, these authors claimed that self-gravity is important to coupling of the in-plane and vertical motions that lead to the formation of the spirals detected by Antoja et al. (2018). We note however that a limitation of their model is the use of an ad hoc imposed bending perturbation of velocity of about 30 km s−1. In such a picture, the properties of phase-space spirals and their duration directly depend on the strength of the imposed perturbation and hence there is no difference between the isolated model with an ad hoc vertical perturbation and the excitation of phase-space spirals due to a passing Sagittarius dwarf galaxy.

Internal processes, external tides, and local galaxy disk perturbations are the main physical factors that could be responsible for the observed phase-space mixing process in the solar neighborhood, but we still do not know the dominant cause of the observed deviations from equilibrium state. In this Letter, we want to study and discuss an alternative scenario for the formation of the phase-space spirals found in Gaia DR2, related to internal mechanisms, and not to external perturbations, by showing in particular that these spirals may be the last vestiges of the bar buckling instability. This study is performed using an extreme high-resolution three-dimensional N-body simulation. The outline of the Letter is as follows: in Sect. 2, we describe our model; in Sect. 3, we examine the phase-space structure of the simulated Milky Way-type galaxy, which develops a bar, that then buckles and forms a boxy/peanut bulge; and finally, in Sect. 4, we summarize the conclusions that can be drawn from our study.

2. Model

In order to investigate the phase-space mixing in the Milky Way type galaxy, we performed a single, extremely high-resolution, N-body simulation of disk galaxy evolution embedded into a live dark matter halo. We adopted the initial conditions from a recent study of the bar/bulge regions of the Milky Way by Fragkoudi et al. (2017) and details are presented in Appendix A.

For a deeper understanding of the evolution in the phase-space distribution of disk stars, we require simulations that are capable of resolving the dynamics down to few tens of parsecs scale; this resolution is much higher than most of the current Milky Way-like galaxies simulations. In order to get inside the rich substructures in the phase space, we used 108 particles to model the stellar disk component and 41 451 200 particles to model the live dark matter halo. To our knowledge, this is one of the highest resolution, N-body simulations of a disk galaxy.

3. Results

The global evolution of our simulated galaxy is in good agreement with many previous N-body models in which bars are usually weakened by the buckling instability (see, e.g., Raha et al. 1991; Combes et al. 1990; Debattista et al. 2004; Athanassoula 2005). In our simulation, the initially axisymmetric disk develops a prominent bar in about two-three rotation periods. Between about 0.5 and 1 Gyr, the bar becomes vertically unstable and buckles, breaking the symmetry with respect to the equatorial disk midplane (see Fig. 1). After ≈1.2 Gyr, the bar profile across the disk plane becomes again symmetric, but the inner disk region now has a larger vertical thickness, and its isodensity-contours show a clear boxy morphology.

thumbnail Fig. 1.

Face-on and edge-on projected density evolution (in log-scale) of the galactic disk. Buckling phase of the bar is depicted in second and third frames. We note that beyond the bar radius the disk density distribution stays axisymmetric during the entire simulation.

Open with DEXTER

We may state that, so far, in relation to the study of the bar buckling process, attention has mostly been drawn to collective phenomena associated with the formation of bulges, while its possible impact on the outer disk regions – that is essentially outside the bar corotation – has not received particular attention. The bending instability has been associated with nonlinear processes during the bar evolution. The growth of the bar involves stars with low-velocity dispersion from the outer to inner regions, where the disk is already significantly hotter. Such a kinematical pressure anisotropy appears to drive the fire-horse instability across the bar (Raha et al. 1991), which in turn leads to the thickening of the disk, as a way to remove the source of instability in the disk. In addition to this, we find that the process of the vertical bar instability (see Fig. 1) is also reflected in the kinematical properties of stars in the outer (relative to the bar) disk.

As follows from our simulation, indeed the development of the boxy/peanut structure also affects the disk areas outside the bar region. To demonstrate this in Fig. 2 we present the evolution of the bar strength (Fourier harmonics A2), median vertical offset ⟨z⟩ and median vertical velocity ⟨νz⟩ measured for stars at different galactocentric radii from 1 to 10 kpc. During the early linear phase of the bar growth (< 0.5 Gyr) our model is stable against the bending instability at all radii. However, once the disk starts to buckle and the bar strength decreases, bending waves appear initially in the central part of the disk, the bar region, forming bell-like structures, which later propagate farther out toward the galactic outskirts. After a relatively short vertical impulse in which the edge of the buckling bar acts as a source of bending waves, quasi-stationary ring-like structures are established over the entire outer galaxy. In particular, a symmetric m = 0 perturbation moves outward, and its amplitude weakly increases in the outer parts of the disk. It takes roughly 0.5 Gyr for perturbations to reach the outer disk part (> 10 kpc). Such a bending wave propagation process is very much similar to those found in a number of various theoretical and numerical works about the excitation of internally driven bending waves in the galactic disks (see, e.g., Hunter & Toomre 1969; Kulsrud et al. 1971; Merritt & Sellwood 1994; Khoperskov et al. 2010; Rodionov & Sotnikova 2013; Khoperskov & Bertin 2017; Chequers & Widrow 2017).

thumbnail Fig. 2.

Evolution of the bar strength A2 (top panel), median vertical coordinate of particles (middle panel), and mean vertical velocity (bottom panel). All frames show the evolution of parameters at different radii in a circular annulus of 0.5 kpc width: from blue in the center to brown lines at the outskirts. To improve the visibility of the ⟨z⟩,⟨νz⟩ time-dependent variations we shift each curve relative to the previous by a constant value of δvi = 0.2 and δri = 0.05 for vertical velocity and vertical coordinate respectively. The end of the bar buckling phase corresponds to ≈1 Gyr.

Open with DEXTER

The consequences of the bar buckling can persist over a significant time – more than 10 rotational periods at the solar radius – and it can have observational counterparts in the (z, νz) phase-space distribution. Indeed, independent of the mechanism, the bending waves propagating across the disk suppress the mixing process in the phase space. It is thus natural to investigate the kind of features left by this mechanism in the z − νz plane. In order to compare our simulation with the results by Antoja et al. (2018), in Fig. 3 we present the evolution of the phase-space mixing patterns in the disk at R = 8 kpc from the center in a small arc of 15° of 1 kpc radial width oriented with 30° relative to the bar. This region contains of about 547 856 to 624 231 particles at different times. At early times, before the perturbation from the buckling bar reaches the solar radius, (z, νz) phase-space structure is highly symmetric in the innermost part (low |⟨z⟩| and |⟨νz⟩|), and only a weak quadrupole pattern dominates in the outer region of the distribution. Such a picture is known as a tilt of the velocity ellipsoid (see, e.g., Sect. 5.1 in Bland-Hawthorn et al. 2018, and references therein), which is the manifestation of a fully mixed galactic disk. Once the vertical perturbations arrive at the outer disk regions (T ≈ 0.5 Gyr), in the νr(z, νz) distribution, we start to observe an asymmetric two-arm spiral structure in (z, νz)-plane, which corresponds to a transition phase between a tilted ellipse and a single arm spiral. Figure 3 demonstrates a clear similarity to the Gaia-like phase-space spiral at T = 1.3 Gyr. The spiral structure is more clearly seen in the distribution of radial velocity (z, νz) space; this structure is less evident in azimuthal velocity, but still distinguishable in density distribution. Thus we find that internally driven bending waves trigger a clear correlation between in-plane and vertical motions that lead to spiral structures that are in good agreement to the phase-space spirals observed at the solar neighborhood.

thumbnail Fig. 3.

Evolution of phase-space spiral in z − νz plane for a solar neighborhood-like region in the disk, at R = 8 kpc. At all times (except the first frame) the bar is rotated by 30° relative to the direction on the galactic center. Top row: stellar density distribution log(ρ(z, νz)/max(ρ(z, νz))), middle row: mean radial velocity, and bottom row: mean residual azimuthal velocity νφ(z, νz)−⟨νφ(z, νz)⟩.

Open with DEXTER

Surprisingly, the phase-space spiral pattern is still visible in a solar neighborhood-like region at later times (T >  1.5 Gyr) after the end of the boxy/peanut bulge formation. The reason for this is that in our model the bar-driven bending waves do not dissipate even after the initial source (buckling of the bar) disappears and the weave-like perturbations persist in the disk for tens of rotational periods (up to 4 Gyr). The bending waves (see Fig. 2) are not kinematical because they do not dissipate on the dynamical timescale. This has been predicted by Darling & Widrow (2019), who showed that self-gravity is crucial for the persistence of long-lived spontaneous bending waves in disks. Hence disk self-gravity is essential for the support of bending waves triggered by the bar buckling. That is why, even a very long time (≈3 Gyr) after the buckling of the bar – z − νz spiral patterns are clearly seen. However, the morphology of spirals is very much time-dependent, and the outer part of the distribution always shows the existence of a quadrupole structure, which is visible as a top right to bottom left reddish asymmetry in the radial velocity distribution νr(z, νz) (see Fig. 3).

The large number of particles in our simulation enable us to derive a global view of the phase-space disk structure. This allows us to understand the uniqueness of the phase-space structure recently discovered in the extended solar vicinity. To do this, in Figs. B.1, B.2, we plot the distributions of stellar density, radial, and azimuthal velocities in the z − νz plane for different parts of the disk at different times: 0.76 Gyr (during the bar buckling) and at 2 Gyr (when the boxy-bulge is formed). Figure B.2 shows the same distributions but at very late times 4 Gyr. These plots show a clear evidence of ongoing phase-space mixing occurring across the entire galaxy. In particular, at early times, when the bending waves are just triggered by the bar, the inner disk (R ≤ 5) shows the presence of one-arm spiral structures. While the outer disk (R >  5) is not vertically perturbed yet and various symmetric quadrupole-like structures for radial velocity distribution in z − νz plane persist there. Thanks to the homogeneous outer disk (beyond the bar size, see Fig. 1), we obtain clear periodic variations of the phase-space patterns depending on the bar position angle (see black contours in Figs. B.1, B.2). Thus, we find that the phase-space structure of disk stars in an isolated, barred Milky Way-type galaxy varies across the disk and the spiral pattern in z − νz plane is nonsteady, but long-lived. This is because the disk self-gravity is important to support and enhance the bending waves beyond the bar radius, even few gigayear after the end of the bar buckling phase.

4. Conclusions

In this paper, which is part of a series exploring various chemo-kinematical features found by Gaia DR2 and large-scale spectroscopic surveys (GALAH, LAMOST), we examine the rich phase-space structure of stellar populations in the Milky Way-type disk. To this aim, we present a model of the natural formation of the spiral patterns in z − νz- plane in an isolated galaxy with a boxy/peanut bulge, such as that of the Milky Way, where the mechanism driving the vertical oscillations is the buckling instability of the bar. The existence of long-lived vertical oscillations in stellar disk has been explored by numerous N-body simulations, but in this work we demonstrate a clear impact of the bar buckling phase on the generation of bending waves propagating across the disk. As a consequence, coherent vertical and in-plane motions are both responsible for generating various phase-space patterns in the galactic disk.

Our model is the first self-consistent, isolated galaxy model that is successful in reproducing the spiral structures in the phase-space discovered by Antoja et al. (2018) and it reveals a rich phase-space structure over the entire Milky Way-type galaxy disk. We have been able to show that the phase-space spirals are related to a nonequilibrium state of the Milky Way disk, and by means of comprehensive analysis of a high-resolution N-body simulation, we demonstrate clearly that it is hard to complete the phase-space relaxation even in an isolated galaxy. We confirm the conclusions of the studies by Widrow et al. (2014) and Darling & Widrow (2019) that self-gravitating disks are able to support long-lived vertical oscillations even when the source of bending instability is gone. That is why in our pure N-body simulation the bending waves do not dissipate quickly after the end of the bar buckling.

In the understanding of the origin of the phase-space spirals, we are thus left with two possible mechanisms: an external, recent perturbation by a massive satellite, as suggested by Antoja et al. (2018), Binney & Schönrich (2018), and Bland-Hawthorn et al. (2018), and an internal mechanism, related to the secular evolution of the disk, as shown in this paper. The main limitations of the accretion scenario are currently related to the high mass of the perturber and the recent onset of the perturbation. As discussed by Binney & Schönrich (2018), indeed, the perturbing satellite would need to be massive and the duration of the effect brief in time. A relatively high mass for the perturber is also found by Bland-Hawthorn et al. (2018), at the level of 3 × 1010M. Perhaps the mass argument is less important though, since there have been suggestions that the progenitor of Sagittarius dwarf galaxy could have been more massive than initially envisaged (see Gibbons et al. 2017, 1011M), but still this is the original mass of the satellite and not the mass after several Gyrs of stripping. The internal scenario does not suffer of any of these two limitations. In our model, we have not fine-tuned either the strengths of the bar and of the boxy/peanut bulge or the timescale of the buckling instability. Interestingly, we find that the timing argument can be relaxed, that is we do not need the perturbation to be as recent as predicted in the case of the satellite accretion. However, it will be necessary to address in future works the upper limit of this timescale in the case in which the phase-space spirals originate from the buckling instability. Could the buckling have occurred about 8 Gyr ago at the time when bars are expected to appear in galaxies of the size of the Milky Way (see Sheth et al. 2008; Kraljic et al. 2012), or do we need a more recent or maybe recurrent (see Martinez-Valpuesta et al. 2006) buckling to explain the observed phase-space patterns? Also, the role of a dissipative component in counteracting the secular heating of the disk, the role of spiral asymmetries, and the Gaia selection function will need to be taken into account in future models. Despite these limitations, our work shows once more that the role of the bar, in our Galaxy, may have been fundamental also in generating of long-lasting phase-space patterns and incomplete mixing.

Acknowledgments

We thank the anonymous referee for her/his valuable comments that helped to improve this paper. This work was granted access to the HPC resources of CINES under the allocation 2017-040507 (PI : P. Di Matteo) made by GENCI. This work has been supported by the ANR (Agence Nationale de la Recherche) through the MOD4Gaia project (ANR-15-CE31-0007, PI: P. Di Matteo). Numerical simulations are carried out partially with support by the Russian Foundation for Basic Research (16-32-60043, 16-02-00649) and using the equipment of the shared research facilities of HPC computing resources at Lomonosov Moscow State University supported by the project RFMEFI62117X0011. PB acknowledge the support by Chinese Academy of Sciences through the Silk Road Project at NAOC, through the “Qianren” special foreign experts program and also under the President’s International Fellowship for Visiting Scientists program of CAS. PB acknowledges the partial financial support of this work by the Deutsche Forschungsgemeinschaft through Collaborative Research Center (SFB 881) “The Milky Way System” (subproject Z2) at the Ruprecht-Karls-Universitat Heidelberg. PB acknowledge the support of the Volkswagen Foundation under the Trilateral Partnerships grant No. 90411 and the special support by the NASU under the Main Astronomical Observatory GRID/GPU computing cluster project.

References

Appendix A: Model description

In this work we employ a multicomponent model for a galaxy consisting of three cospatial disk populations with different parameters. Each disk component was represented by a Miyamoto-Nagai density profile (Miyamoto & Nagai 1975) that has a characteristic scale length of 2, 2, and 4.8 kpc, vertical thicknesses of 0.15, 0.3, and 0.6 kpc and masses of 1.86, 2.57, and 4.21 × 1010 M, respectively. Our simulation also includes a live dark matter halo whose density distribution follows a Plummer sphere (Plummer 1911), which has a total mass of 3.81 × 1011 M and a radius of 21 kpc. The choice of parameters leads to a galaxy mass model with a circular velocity of ≈220 km s−1 and Toomre stability parameter QT ≈ 1.5 at 8 kpc. The initial setup has been generated using the iterative method by Rodionov et al. (2009).

To avoid small-scale bends due to initial artificial nonaxisymmetries, we forced a proper symmetry for the initial positions and velocities of all particles in the simulation (e.g., similar to Debattista 2014). As a consequence, the stellar disk is highly symmetric about the midplane. This model, where the thick disk is massive and centrally concentrated, can reproduce the morphology of the metal-rich and metal-poor stellar populations in the Milky Way bulge, as well as the mean metallicity and [α/Fe] maps as obtained from the APOGEE data (Fragkoudi et al. 2018). Hence we emphasize that our model is meant to be representative of the Milky Way disk, its stellar populations, and global kinematics.

For the N-body system integration, we used our parallel version of the TREE-GRAPE code (Fukushige et al. 2005, Sect. 4.1) with multithread usage under the SSE and AVX instructions. We also ported the original GRAPE-based (GRAvity PipEline) tree code to few different hardware platforms, including the CPU and the recent IA Graphics Processing Unit platform using the Compute Unified Device Architecture. In recent years we already used and extensively tested our hardware-accelerator-based gravity calculation routine in several galaxy dynamics studies where we obtained accurate results with a good performance (Shumakova & Berczik 2005; Bien et al. 2008; Pasetto et al. 2010, 2011). In the simulation we adopted the standard opening angle θ = 0.5 and a gravitational softening parameter equal to 10 pc. For the time integration, we used a leapfrog integrator with a fixed step size of 0.1 Myr (Khoperskov et al. 2014).

thumbnail Fig. A.1.

Circular velocity (blue curve) of the model, together with the rotation curve (red).

Open with DEXTER

Appendix B: Additional plots

thumbnail Fig. B.1.

Spatial variations of the phase-space patterns in simulated galaxy at 0.76 Gyr (top row) and at 2 Gyr (bottom row). Each frame corresponds to number density (left panels), mean azimuthal (middle panels), and radial velocities (right panels) in z − νz plane for stars selected in small arcs |R − Ri|< 0.5 kpc and |ϕ − ϕi|< 7.5°, where cylindrical coordinates (Ri, ϕi) are indicated on global axes. Each frame axes size is [ − 1; 1] kpc × [−50; 50] km s−1. Total galaxy stellar surface density contours are over-plotted by black lines, which clearly show the position of the bar in galactic cylindrical coordinates. Number of particles varies from ≈4 × 106 in the innermost frames to ≈0.3 × 106 in the outermost frames.

Open with DEXTER

thumbnail Fig. B.2.

Same as in Fig. B.1, but for latter time of 4 Gyr.

Open with DEXTER

All Figures

thumbnail Fig. 1.

Face-on and edge-on projected density evolution (in log-scale) of the galactic disk. Buckling phase of the bar is depicted in second and third frames. We note that beyond the bar radius the disk density distribution stays axisymmetric during the entire simulation.

Open with DEXTER
In the text
thumbnail Fig. 2.

Evolution of the bar strength A2 (top panel), median vertical coordinate of particles (middle panel), and mean vertical velocity (bottom panel). All frames show the evolution of parameters at different radii in a circular annulus of 0.5 kpc width: from blue in the center to brown lines at the outskirts. To improve the visibility of the ⟨z⟩,⟨νz⟩ time-dependent variations we shift each curve relative to the previous by a constant value of δvi = 0.2 and δri = 0.05 for vertical velocity and vertical coordinate respectively. The end of the bar buckling phase corresponds to ≈1 Gyr.

Open with DEXTER
In the text
thumbnail Fig. 3.

Evolution of phase-space spiral in z − νz plane for a solar neighborhood-like region in the disk, at R = 8 kpc. At all times (except the first frame) the bar is rotated by 30° relative to the direction on the galactic center. Top row: stellar density distribution log(ρ(z, νz)/max(ρ(z, νz))), middle row: mean radial velocity, and bottom row: mean residual azimuthal velocity νφ(z, νz)−⟨νφ(z, νz)⟩.

Open with DEXTER
In the text
thumbnail Fig. A.1.

Circular velocity (blue curve) of the model, together with the rotation curve (red).

Open with DEXTER
In the text
thumbnail Fig. B.1.

Spatial variations of the phase-space patterns in simulated galaxy at 0.76 Gyr (top row) and at 2 Gyr (bottom row). Each frame corresponds to number density (left panels), mean azimuthal (middle panels), and radial velocities (right panels) in z − νz plane for stars selected in small arcs |R − Ri|< 0.5 kpc and |ϕ − ϕi|< 7.5°, where cylindrical coordinates (Ri, ϕi) are indicated on global axes. Each frame axes size is [ − 1; 1] kpc × [−50; 50] km s−1. Total galaxy stellar surface density contours are over-plotted by black lines, which clearly show the position of the bar in galactic cylindrical coordinates. Number of particles varies from ≈4 × 106 in the innermost frames to ≈0.3 × 106 in the outermost frames.

Open with DEXTER
In the text
thumbnail Fig. B.2.

Same as in Fig. B.1, but for latter time of 4 Gyr.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.