3D stellar motion in the axisymmetric Galactic potential and the e–z resonances

In the past decades, there has been a growing set of evidence that the disk stars of the Milky Way galaxy exhibit density and velocity structures in several phase-space planes of the Galactic disk, in both the radial and vertical directions. The long-known stellar warp in the outer disk, a vertically asymmetric distribution of stars close to the Galactic plane, is seen as a bending of the plane upwards in the first and second Galactic quadrants (longitudes $0^{\circ}\leq l\leq 180^{\circ}$) and downwards in the third and fourth quadrants (longitudes $180^{\circ}\leq l\leq 360^{\circ}$) (Momany et al., 2006). Recently, a number of stellar density substructures, appearing as overdensities or dips, have been found in the LAMOST data, at several radial and vertical positions in the Galactic disk (Wang et al., 2018).

Large-scale vertical motions as a Galactic North–South asymmetry have been revealed by spectroscopic and astrometric surveys such as SEGUE (Widrow et al., 2012), RAVE (Williams et al., 2013), LAMOST (Carlin et al., 2013), as well as the second Gaia data release (Gaia Collaboration et al., 2018). Such bulk vertical motion of stars present a wave-like behaviour, which has been referred to as bending and breathing motions, with stars coherently moving towards or away from the Galactic mid-plane, on both sides of it (Kawata et al., 2018; Ghosh et al., 2022; Khachaturyants et al., 2022). Several dynamical processes have been evoked to explain these non-zero vertical motions: the ones from external origins such as the passing of the Saggitarius dwarf galaxy through the Milky Way disk (Gómez et al., 2013), or the excitation of the disk due to interactions with dark matter subhaloes (Feldmann & Spolyar, 2015); and the ones from non-axisymmetric internal perturbations such as the breathing motion induced by the Galactic bar (Monari et al., 2015), or by the action of the spiral density waves (Faure et al., 2014; Ghosh et al., 2022; Khachaturyants et al., 2022, among others).

The density structures present in the $R$–$V_{\varphi}$ plane (Galactocentric radius versus azimuthal velocity) of disk stars, known as diagonal ridges, clearly visible in the distribution of the mean Galactic radial velocity $\left<V_{R}\right>$, can also be seen in the vertical direction from maps of the mean absolute distance from the disk mid-plane $\left<|z|\right>$ and the mean vertical velocity $\left<V_{z}\right>$ of the stellar distribution (Khanna et al., 2019; Wang et al., 2020). This attribute may be indicative of a coupling between planar and vertical motions of the stars in the disk. Several attempts to unravel the origin of these structures have been presented in the literature, with some of them relying on simulations of external mechanisms like the Sagittarius dwarf galaxy perturbation (Antoja et al., 2018; Laporte et al., 2019; Khanna et al., 2019), and others from simulations that take into account internal dynamics such as the bar or the spiral resonances (Hunt et al., 2018; Michtchenko et al., 2018; Fragkoudi et al., 2019; Barros et al., 2020).

Another striking feature associated with the vertical motion of the stars in the Galactic disk is the phase-space spiral (or the so-called snail shell) present in the $z$–$V_{z}$ plane. Many authors interpret it as evidence of an ongoing phase mixing in the vertical direction of the disk due to an influence of an outer perturbation (Antoja et al., 2018; Binney & Schönrich, 2018; Bland-Hawthorn et al., 2019; Laporte et al., 2019), or caused by the instability of a buckling bar (Khoperskov et al., 2019), or even due to the earlier discussed vertical bending waves (Darling & Widrow, 2019) or vertical breathing motion (Hunt et al., 2022, for two-armed phase spirals) in the disk. Alternatively, Michtchenko et al. (2019) showed that the phase spiral can be well-comprehended by the dynamical effects of the stellar moving groups in the solar neighbourhood.

We suggest to analyse the above-cited Galactic structures and investigate their causes by studying the stellar dynamics in the following steps:

• •

  Modeling a 3D Galactic gravitational potential. In this work, we use the Barros et al. (2016) axisymmetric Galactic potential model, adjusted to recent observational data. Firstly, the chosen potential is suitable to the study of the vertical stellar motion and, secondly, it is simple in analytical forms that gives a good representation of the radially exponential mass distribution of the Galaxy. Moreover, any non-axisymmetric mass effects (e.g., due to the central bar and/or spiral arms) can be easily added to the model posteriorly;

• •

  Studying the Galactic model in the whole phase space, in order to obtain all the possible regimes of stellar motion that could provide reasonable explanations to the observed Galactic phase-space structures. This is, in part, the main objective of the present paper;

• •

  Performing numerical simulations through numerical integrations of stellar orbits, to verify whether our achievements are satisfactory. This is a topic for a future work.

In the present paper, we focus on the study of the stellar orbits in the Galactic disk assuming a zeroth-order approximation of a time-independent and axisymmetric Galactic gravitational potential. The whole phase space around the Sun is investigated through techniques widely used in the Celestial Mechanics. Resonances between the radial and vertical independent modes of the stellar motion are characterized in terms of resonant zones, formed by islands of stability that capture and trap stars inside them, enhancing the stellar density, and separatrices that deplete the objects close to these regions. Our goal, for a subsequent work, is to investigate possible associations between the observed Galactic phase-space structures and the resonant zones originated from the commensurabilities between the radial and vertical frequencies of the stellar motion.

This paper is organized as follows. In Section 2, we describe the 3D axisymmetric Galactic gravitational potential model used for the construction of the Hamiltonian function and the calculus of the stellar orbits. In Section 3, we study the 3D stellar dynamics on the representative plane $R$–$V_{\varphi}$ in the radial and vertical directions. In Section 4, regimes of stellar motion are analysed in terms of the beating between the planar and vertical frequencies of oscillation, and, in Section 5, we extend the analysis for a wide range of initial vertical velocities $V_{z}$. Concluding remarks are drawn in the closing Section 6.

Here, we briefly describe the 3D model for the axisymmetric Galactic potential and refer the reader to Barros et al. (2016) for more details. The model considers the contributions of three Galactic components into the potential; they are the axisymmetric disk, the central spheroidal bulge and the spherical halo of dark matter.

The axisymmetric disk is composed of the stellar (thin and thick disks) and gaseous (H I and $\mathrm{H}_{2}$ disks) components. Each of these components is a superposition of three Miyamoto–Nagai disks (MN, Miyamoto & Nagai, 1975), which reproduces a radially exponential mass distribution (Smith et al., 2015). We adopt the 1st–, 2nd– and 3rd–order expressions for the potential of the MN–disks at the position $R$ and $z$ (cylindrical coordinates), respectively:

where $\zeta=\sqrt{z^{2}+b^{2}}$, with $b$ being the vertical scale length, while $M$ and $a$ are the mass and the radial scale length of each disk component, respectively.

The potential of the thin disk is then written as a combination of three third-order MN–disks

while the potential of the thick disk is a combination of three first-order MN–disks

The contributions of the gaseous disks components H I and $\mathrm{H}_{2}$ into the total axisymmetric potential are, respectively,

The physical and structural parameters of the components of the Galactic axisymmetric potential are obtained by applying a fitting procedure, which adjusts the analytical rotation curve Eq. (10) to the observed one. For the observation-based rotation curve, we use data of H I-line tangential directions from Burton & Gordon (1978) and Fich et al. (1989), CO-line tangential directions from Clemens (1985), and maser sources data associated with high-mass star-forming regions from Reid et al. (2019) and Rastorguev et al. (2017). From these data, the Galactic radii and rotation velocities were calculated taking the distance of the Sun from the Galactic center as $R_{0}=8.122$ kpc (GRAVITY Collaboration et al., 2018) and the Local Standard of Rest (LSR) velocity at the Sun $V_{0}=233.4$ km s^-1 (Drimmel & Poggio, 2018). As additional constraints to the fitting procedure, we used the value of the local angular velocity $\Omega_{0}$, the local disk surface density $\Sigma_{0}$, and the local surface density within $|z|\leq 1.1$ kpc. For details of the data used and the fitting procedure, see Barros et al. (2016, 2021).

The parameters obtained are shown in Table 1 and Table 2. The contributions of each component of the axisymmetric potential to the Galactic rotation curve are plotted in Fig. 1 (top panel), together with the calculated rotation curve (red line) and the data points. Figure 1 (bottom panel) also shows the radial profile of the $K_{z}$ vertical force resulting from the Galactic model calculated at $z=1.1$ kpc. We note a good agreement between the modeled force and the measured $K_{z}$–force, based on measurements by Bovy & Rix (2013) (black dots with error bars).

The dynamical model defined by the Hamiltonian in Eq. (12) is of two-degree-of-freedom. The resulting stellar motion can be formally represented by two coupled oscillations: one, in the radial (equatorial) direction, is described by the pair of the variables $R$–$p_{r}$, and other, in the vertical direction, described by the pair $z$–$V_{z}$. The phase space of the system is four-dimensional, that makes it difficult to visualize the stellar orbits. Yet, in this paper, we introduce a representative plane, which allows us to present the main features of the 3D stellar dynamics. This plane is the $R$–$V_{\varphi}$ plane, where $V_{\varphi}$ is the tangential velocity defined as $V_{\varphi}=L_{z}/R$; we show this plane in Fig. 2.

To study the stellar motions on the $R$–$V_{\varphi}$ plane, we fix, in this paper, the initial values of the planar linear momentum at $p_{r}=0$ and the vertical height at $z=0$. It is worth noting that this choice preserves the generality of the presentation, since all closed stellar motions under the potential given by Eq. (11) pass through these conditions. The initial value of the vertical velocity is chosen as $V_{z}=20$ km s^-1, except when the dependence of the motion on the initial vertical configuration is analysed (Sect. 5). The chosen velocity value is close to the value of the standard deviation of the $V_{z}$ distribution of the stars from the Gaia DR3 (Gaia Collaboration et al., 2016, 2023), with $|z|<1$ pc (we obtained $\sigma_{V_{z}}=18$ km s^-1).

First, we plot on the $R$–$V_{\varphi}$ plane in Fig. 2 the main characteristics of the Hamiltonian dynamics: the orbital energy and the angular momentum, which are both constants of motion in our model. The energy levels and $L_{z}$–levels are shown by dashed and continuous lines, respectively. The rotation curve obtained using Eq. (10), is shown by the red curve.

By definition, the rotation curve is a place of the circular orbits, which are orbits of the minimal energy, for a given $L_{z}$-value. The geometrical interpretation of this condition is that, given the location of a circular orbit on the rotation curve, the corresponding $L_{z}$-level is tangent to the minimal energy level on the $R$–$V_{\varphi}$ plane. This fact can be observed in Fig. 2, considering that the energy is continuously increasing with the increasing radial distances $R$, as shown on the top panel in Fig. 2. The same $L_{z}$-level intersects the levels of higher energy always at two points, which are turning points of the radial oscillations at the condition $p_{r}=0$.

We plot in Fig. 2 the direct projections of three 3D stellar orbits: the one is the Sun’s orbits (cyan dots) calculated with initial values $R=8.122$ kpc, $p_{r}=-12.9$ km s^-1, $V_{\varphi}=245.6$ km s^-1, $z=0.015$ kpc and $V_{z}=7.78$ km s^-1 (Drimmel & Poggio, 2018). The orbits of other two fictitious stars (black dots) were starting at points 1 and 2, respectively, with $p_{r}=0$, $z=0$ and $V_{z}=20$ km s^-1. We can verify that all orbits are closed, that is, each one evolves between two turning points, which belong to the same energy levels, along the corresponding $L_{z}$–levels (constant of motion), around the circular orbits, which lie at the intersection of the corresponding energy and $L_{z}$–levels.

The detailed analysis of the main features of the stellar motions on the representative $R$–$V_{\varphi}$ plane can be done in two steps: first, focusing on the radial oscillations on the Galactic equatorial plane, and, then, on the vertical oscillations around the equatorial plane. By understanding the behaviour of each component of the stellar motion in the axisymmetric potential, we can consider new effects produced by their interaction. We will do it in the next sections.

To understand the features of the resonant motion, we start plotting the amplitude of the vertical oscillation ($z_{\mathrm{max}}$) as a function of the initial radial distance $R$ on the top panel in Fig. 6. In contrast to one another shown on the top panel in Fig. 4, the amplitude $z_{\mathrm{max}}$ was now calculated with initial conditions chosen along the $L_{z}$–level corresponding to the Sun’s orbit. In this case, all orbits are eccentric, with exception of one located at $R_{c}=8.522$ kpc, at the intersection of the $L_{z}$–level with the rotation curve in Fig. 4.

For the orbits starting along the same $L_{z}$–level at the same initial configurations in the vertical direction ($z=0$ and $V_{z}=20$ km s^-1), the minimal value of $z_{\mathrm{max}}$ is associated with the circular orbit, at $R_{c}=8.522$ kpc. The vertical amplitude is generally increasing when the orbital eccentricity grows with both increasing and decreasing distances from the circular orbit. However, this evolution of $z_{\mathrm{max}}$ is not monotonous, but shows sudden increase/decrease at some radial distances. The nature of this behaviour can be analyzed using the dynamical power spectrum method, which allows us to detect the change of a regime of motion through the analysis of the evolution of the main frequencies of the dynamical system (detailed description of the method and its applications to different systems are found in e.g. Michtchenko et al., 2002, 2017).

Figure 6 bottom shows two proper frequencies, $f_{r}$ and $f_{z}$, which are the frequencies of the radial and vertical oscillations, respectively, as functions of the radial distances $R$. The evolution of the $f_{r}$–frequency (red dots) and its harmonics is monotonous over the whole $R$–range; it reaches the maximal value at $R_{c}=8.522$ kpc, that is expected since the circular orbit is an equilibrium solution of the effective potential $\Phi_{\mathrm{eff}}=\frac{L_{z}^{2}}{2R^{2}}+\Phi_{0}(R,z)$, with the constant $L_{z}$. On the other hand, the behaviour of the $f_{z}$–frequency (black dots) on the dynamical spectrum is highly non-harmonic and we can observe the appearance of bifurcations, stable islands and an erratic scattering of the dots, which is characteristic of the chaotic motion.

Analyzing the evolution of the proper frequencies in Fig. 6 bottom, we verify that bifurcations occur when the beats between the two independent frequencies occur. This condition is known as resonances between distinct modes of motion, which lead to changes in the regime of motion. One beat occurs around $R=14$ kpc, where $f_{R}\cong 2\,f_{z}$, giving origin to the island of the 2/1 resonant motion. The periodic orbit of the 2/1 resonance is shown in the subspaces $R$–$p_{R}$ (top) and $z$–$V_{z}$ (bottom) in Fig. 7 (column a); the projections of this orbit in the Poincaré maps are shown by red dots in Figs. 3 and 5.

The beating frequencies are also observed at $R=17.8$ kpc in the dynamical spectrum in Fig. 6 bottom; this event is associated to the 4/1 resonance. The radial and vertical oscillation modes of the corresponding periodic orbit are shown in Fig. 7 (column b); the projection of this orbit in the Poincaré maps is shown by cyan dots in Figs. 3 and 5. In Fig. 5, we also show, by green dots, the projection of the 6/1 resonant orbit (see Fig. 7 c); the location of this resonance in the dynamical spectrum is shown at $R=19.2$ kpc.

We denote these resonances as the $e$–$z$ resonances, where $e$ and $z$ denote planar eccentricities and vertical heights of stellar orbits, respectively. The most important from the $e$–$z$ resonances is the 1/1 resonance, whose domain is large and separated by the layers of chaotic motion (separatrices). For a given value of $L_{z}$, it appears at very large Galactic distances; in the dynamical map in Fig. 6, its domain is extended from 20 kpc to 35 kpc. The motion inside the 1/1 resonance is stable, with the amplitude of the vertical oscillation reaching its minimal value at the center (locus) of the resonance, at $\sim$27 kpc; this locus appears as a fixed point (black dot) in the Poincaré map in Fig. 5. The locus is a periodic orbit of the 1/1 resonance shown in Fig. 7 (column d). The neighborhood of the 1/1 resonance is generally a domain of highly unstable motion. This is due to the overlap of the high-order resonances of the type $f_{r}/f_{z}\cong m/n$ ($m$ and $n$ are integers).

Finally, it is worth emphasizing that, comparing the evolution of the two frequencies in Fig. 6 bottom, we note the qualitative difference in the behaviour of two independent modes: while the radial motion seems to be unaffected by the passages through the resonances, their impact on the vertical motion is significant. This could be explained by the peculiar characteristics of the massive disk potential.

In this section, we investigate the resonant behaviour as a function of the initial vertical velocity, $V_{z}$. For this, we introduce the representative plane $R$–$V_{z}$ of the initial conditions and fix the rest of the variables at $p_{r}=0$, $z=0$ and $L_{z}=1995$ kpc km s^-1, this last corresponding to the Sun’s angular momentum.

Figure 8 presents the amplitude of the vertical oscillation in the form of the $z_{\textrm{max}}$–levels on the representative $R$–$V_{z}$ plane. The red curve shows the positions of the circular orbits of the constant $L_{z}$, which are slightly dislocated in the direction of the larger Galactocentric distances, when the initial vertical velocities increase. The panel beside the $R$–$V_{z}$ plane allows us to quantify the $z_{\textrm{max}}$–amplitude, showing its values calculated along the loci of the circular orbits (red curve). The smoothed evolution of this quantity with the increasing initial $V_{z}$ can be observed.

The bifurcation structures, which are characteristic of the resonances, appear outside the loci of the circular orbits in Fig. 8 and are strengthened with the increasing planar eccentricities of the orbits. The black curve on the beside panel shows the evolution of the $z_{\textrm{max}}$–amplitude calculated along the constant initial $R=13$ kpc. We can observe the behaviour which is similar to that shown on the top panel in Fig. 6, indicating the passages through some $e$–$z$ resonances. The most prominent passage is through the 1/1 resonance, whose loci are shown by the blue dots in Fig. 8. The projection of the Sun’s trajectory on the $R$–$V_{z}$ plane is shown by the cyan dots. Its proximity to the rotation curve avoids the capture in one of the resonances.

In this paper, we report the existence of the resonant motion of the kind $e$–$z$ ($e$ and $z$ being the planar eccentricity and the vertical altitude of stellar orbits, respectively), produced by the axisymmetric Galactic 3D potential.

We investigate the spacial motion of the stars applying the elaborated model for the axisymmetric potential of the Galactic disk (Barros et al., 2016). The model accounts for the combined gravitational effects due to the two stellar disks (thin and thick disks) and two gaseous disks (H I and $\mathrm{H}_{2}$ disks). Moreover, to account for a radially exponential Galactic mass distribution, each of the disks is approximated by the composition of the three Miyamoto-Nagai disk profiles (Smith et al., 2015). The physical and structural parameters of the modeled components were chosen to correspond to the observable rotation curve of the Galaxy.

The motion of stars in the axisymmetric Galactic potential is investigated in the whole phase space applying several techniques from the Celestial Mechanics. The one of those is the dynamical mapping of the representative plane chosen here as the $R$–$V_{\varphi}$ plane. Based on the conservation laws, we identify two independent modes of the stellar motion, radial and vertical ones. The radial motion is an oscillation around the circular solution, which is a circular orbit belonging to the rotation curve and characterized by the same angular momentum of the radial mode. The vertical motion is an oscillation around the equatorial Galactic plane.

For nearly circular orbits, with small planar eccentricities, two oscillations are weakly coupled. It is worth noting that the circular orbits survive even in the case when vertical velocities are very high. However, when the planar eccentricities increase, the coupling between two modes of motion also increase. Their interaction is better visualised in the dynamical power spectra, which clearly show the regions in the phase space where two proper frequencies are commensurable. This beating frequencies indicate the resonant zones, which we denote as $e$–$z$ resonances.

In the Hamiltonian theories, resonances are fundamental properties of dynamical systems with two or more degrees of freedom. They occur in domains of the phase space where the frequencies of the independent modes of motion are commensurable. The locations and the sizes of the resonant domains are strongly dependent on the physical parameters of the model adopted to describe the system under study. Thus, associating some observable structures of the objects to the resonances, we can assess the reliable ranges of the parameters of the applied model. One example of this approach is shown in Michtchenko et al. (2018), where we relate the moving groups in the solar neighbourhood to the Lindblad resonances produced by the spirals perturbations, and estimate the spiral strength and pattern speed.

The main features of the $e$–$z$ resonant motion are sudden increase/decrease in the amplitude of the vertical oscillation, capture and trap of the stars inside the stable resonant zones, enhancing the density, and deplete of the regions close to the separatrices of the resonances. This behaviour is characteristic of the vertical oscillations, while the planar radial motion is only slightly affected by the $e$–$z$ resonances; this is probably due to peculiar properties of a disk potential at all.

It is worth emphasizing that, despite the fact that the resonances occur at very large Galactic distances, the high eccentricities of the resonant orbits allow us their detection in the solar neighbourhood, according to Fig. 7. Indeed, analyzing the behaviour of the Gaia-eDR3 sample from the solar vicinity ($R_{\odot}\pm 0.2$ kpc and $|p_{R}|\geq 100$ km s^-1), we have detected 10% (from one thousand randomly chosen stars) of the objects showing the resonant oscillations, mainly, inside the 2/1 resonance. Of course, whether these stars are really evolving inside the $e$–$z$ resonances, it depends strongly on the parameters adopted in our model. Thus, we should be able to observe the described manifestations of the $e$–$z$ resonances analysing the distributions of the observable proper elements of the stars. The comparison to theoretical predictions could allow us to improve unknown values of the parameters, which describe the observable Galactic mass distribution.