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

The full phase space information on the kinematics of a huge number of stars provided by the Gaia third Data Release raises the demand for a better understanding of the 3D stellar dynamics. In this paper, we investigate the possible regimes of motion of stars in the axisymmetric approximation of a Galactic potential model. The model consists of three components: the axisymmetric disk, the central spheroidal bulge and the spherical halo of dark matter. The axisymmetric disk is divided into stellar and gaseous disk subcomponents, each one modeled by three Miyamoto-Nagai profiles. The physical and structural parameters of the Galaxy components are adjusted by observational kinematic constraints. The phase space of the two-degrees-of-freedom model is studied by means of the Poincar\'e and dynamical mapping, the dynamical spectrum method and the direct numerical integrations of the Hamiltonian equations of motion. For the chosen physical parameters, the nearly-circular and low-altitude stellar behaviour is composed of two weakly coupled simple oscillations, radial and vertical motions. The amplitudes of the vertical oscillations of these orbits are gradually increasing with the growing Galactocentric distances, in concordance with the exponential mass decay assumed. However, for increasing planar eccentricities and the altitudes over the equatorial disk, new regimes of stellar motion emerge as a result of the beating between the radial and vertical oscillation frequencies, which we refer to as e-z resonances. The corresponding resonant motion produces characteristic sudden increase or decrease of the amplitude of the vertical oscillation, bifurcations in the dynamical spectra and the chains of islands of stable motion in the phase space. The results obtained can be useful in the understanding and interpretation of the features observed in the stellar 3D distribution around the Sun.


Introduction
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 • ≤ l ≤ 180 • ) and downwards in the third and fourth quadrants (longitudes 180 • ≤ l ≤ 360 • ) (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 nonaxisymmetric 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 φ plane (Galactocentric radius versus azimuthal velocity) of disk stars, known as diagonal ridges, clearly visible in the distribution of the mean Galactic radial velocity ⟨V R ⟩, can also be seen in the vertical direction from maps of the mean absolute distance from the disk mid-plane ⟨|z|⟩ and the mean vertical velocity ⟨V z ⟩ 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 φ 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.  (Reid et al. 2019;Rastorguev et al. 2017), and H I and CO tangent-point data (Burton & Gordon 1978;Clemens 1985;Fich et al. 1989).The red curve shows the analytical rotation curve expressed by Eq. ( 10).The contributions of the modelled disk (continuous curve), bulge (short-dashed line) and dark halo (longdashed line) to the axisymmetric Galactic potential are also shown.Bottom: Vertical force K z at z = 1.1 kpc as a function of the Galactic radius.The black dots with error bars are from observation measurements by Bovy & Rix (2013), while the grey solid curve shows the K z force radial profile at z = 1.1 kpc resulting from our Galactic model.

3D Galactic model
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 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 ζ = √ 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 H 2 into the total axisymmetric potential are, respectively, The potential of the Galactic bulge is derived by a Hernquist density distribution profile, in the form (Hernquist 1990): where two free parameters, M b and a b , are the total mass and the scale radius of the bulge, respectively.Finally, we consider a spherical dark halo, whose potential is modelled with a logarithmic potential in the form (e.g.Binney & Tremaine 2008): where r h is the core radius and v h is the circular velocity at large R (i.e., relative to the core radius), respectively.The analytical rotation curve is then calculated using the expression.
where the axisymmetric potential Φ 0 is just the sum of the potentials of the Galactic components described above, that is: 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.4km s −1 (Drimmel & Poggio 2018).As additional constraints to the fitting procedure, we used the value of the local angular velocity Ω 0 , the local disk surface density Σ 0 , and the local surface density within |z| ≤ 1.1 kpc.For details of the data used and the fitting procedure, see Barros et al. (2016Barros et al. ( , 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 stellar orbits are calculated through the numerical integrations of the equations of motions defined by the Hamiltonian function given as where p r and V z are the linear radial momentum and vertical velocity, respectively, while L z is the angular momentum, which is a constant in the axisymmetric approximation (all momenta are given per mass unit)

3D stellar dynamics on the representative plane
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 φ plane, where V φ is the tangential velocity defined as V φ = L z /R; we show this plane in Fig. 2.
To study the stellar motions on the R-V φ 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 σ V z = 18 km s −1 ).
First, we plot on the R-V φ 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 φ 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 -Fig.3. Poincaré map of the stellar 3D orbits calculated with the initial conditions chosen along the Sun's L z -level.The initial vertical conditions were chosen at z = 0 and V z = 20 km s −1 .The projections are calculated at instants when the star crosses the equatorial plane with V z > 0. The orbits are concentric, with the circular orbit at the center (black dot).This orbit lies at the intersection of the rotation curve and the Sun's L z -level on the R-V φ plane, at R c = 8.522 kpc and V φ = 234 km s −1 .Red and cyan dots show the projections of two periodic orbits on the Poincaré map corresponding to the 2/1 and 4/1 resonances, respectively (see later Fig. 7).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.9km s −1 , V φ = 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 φ 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.

Radial oscillations
The R-V φ plane shown in Fig. 2 is commonly used for understanding the radial (or planar) motion of stars.Indeed, due to the conservation of the angular momentum during the motion, the projections of the 3D orbits are aligned along the corresponding L z -levels (continuous lines), as observed in the case of our three examples.Since the orbital energy is also conserved in our problem, the orbit's projections must match the corresponding energy level at conditions p r = 0, which was chosen in the construction of the plane.By definition, this condition defines two turning points of a closed orbit, which delimit the path of the radial oscillation on the R-V φ plane and, consequently, the maximal (R max ) and minimal (R min ) values of the radial orbital variation (and the tangential velocity V φ , for a given L z -value).
Any L z -level crosses the rotation curve (red curve) at one point at R c on the R-V φ plane in Fig. 2, which corresponds to the circular orbit of the minimal energy for the corresponding L z .At this condition, two turning points merge at one point, characterizing the circular orbit with the zero-amplitude R-oscillation, that is, the periodic orbit of the two-degree-of-freedom problem.All other orbits calculated along the same L z -level, are quasiperiodic orbits oscillating with the non-zero R-amplitudes, in such a way that larger the deviations from the rotation curve, larger are the amplitudes of the radial variations.
The amplitude of the radial oscillation is related to the planar eccentricity of the orbit as which is uniquely defined by the values of the orbital energy and angular momentum of the stellar orbit.In the conservative problem considered in this paper, it is a constant of motion.The gain/loss of the orbital energy due to some external physical processes, without changing the angular momentum of the system, will increase/decrease the amplitude of R-oscillation and, consequently, the planar eccentricity.This effect, known as 'heating' in the stellar dynamics, is analogous to the tidal gravitational interactions in the planet-satellite systems.
The equatorial motions of the stars are shown in Fig. 3, where we plot the orbits in the subspace R-p R .The orbits were calculated through numerical integrations of the equations of motion defined by the Hamiltonian (Eq.12), with initial conditions chosen along the level L z = 1995 kpc km s −1 , which corresponds to the Sun's orbit (see Fig. 2).To obtain the projection of the 3D motion on the Galactic plane, we gather the orbital coordinates at the instants when the orbits cross the equatorial plane (the condition z = 0), in the direction of the positive z-values (the condition V z > 0).This approach is known as the Poincaré surfaces of section method, which allows us to plot the radial motion separately from the vertical one.
Figure 3 shows that all orbits in the R-p R subspace oscillate around a circular periodic orbit (black dot), which lies at the intersection of the corresponding L z -level and the rotation curve, at R c = 8.522 kpc and V φ = 234 km s −1 on the R-V φ plane in Fig. 2.This periodic orbit appears as a fixed point on the R-p R plane, indicating that the amplitude of the radial oscillation is zero (or e = 0).At the same time, the vertical motion started with the initial V z = 20 km s −1 has non-zero amplitude; it is a simple oscillation around the equatorial plane (z = 0).
In the vicinity of the fixed point (black dot), between 7.0 kpc and 10.5 kpc, the concentric low-eccentricity orbits, with e < 0.2, evolve in good agreement with the epicyclic approximation (e.g.Binney & Tremaine 2008).The orbit of the Sun, with eccentricity of 0.061, is located in this region of the phase space.The values of the radial and vertical frequencies for the Sun's L z level, derived under the epicyclic approximation, are equal to 0.00542 and 0.01387 (in units of 1/Myr), respectively, that corresponds to the radial and vertical periods of 184.5 Myr and 72.07 Myr.These periods derived from the power spectrum of the Sun's orbit, calculated numerically through our model, are of 161.3 Myr and 75.2 Myr, respectively.
The orbits with the increasing eccentricity, however, deviate from the epicyclic approximation and, when the eccentricity approaches 0.5, a notable feature appears in Fig. 3. Indeed, at the radial distances beyond 12 kpc, we can clearly observe a lack of circulating orbits, but there are some islands on the Poincaré   3, except the projections were calculated at conditions p R = 0 and ṗR > 0 on the z-V z plane.Red, cyan and green dots show the projections of three periodic orbits on the Poincaré map corresponding to the 2/1, 4/1 and 6/1 resonances, respectively.The large island of the orbits oscillating around the center at z = 0 kpc and V z = 103.0km s −1 (black dot) is related to the strong 1/1 resonance (see later Fig. 7).
map.As shown below, the orbits in this region are different structurally from the low-eccentricity orbits and are related to the new kind of motion mode.The planar axisymmetric model is not able to explain this feature, thus, we proceed investigating the vertical motion and its interaction with the planar component.

The vertical oscillations
The same representative R-V φ plane shown in Fig. 2 can be used to study the stellar motion in the vertical direction.For this, we plot in Fig. 4, by continuous lines, the levels of the maximal altitude over the equatorial plane that a star reaches during its vertical oscillation (z max ).To obtain z max , we integrate numerically, over several billions years, the equations of motion defined by the Hamiltonian (Eq.12), over the 200×200 grid of the initial conditions covering the R-V φ plane.All orbits started with p R = 0 on the galactic equatorial plane (at z = 0) and with the same initial vertical velocity V z = 20 km s −1 .The panel on the top of Fig. 4 shows z max obtained in this way for the circular orbits located along the rotation curve (red curve).We plot also the projection of the Sun's orbit (cyan dots) along the L z -level (dashed line) on the R-V φ plane.
To understand the z max -evolution on the Galactic R-V φplane, we consider first the nearly circular orbits, distributed along the rotation curve.For these low-eccentricity orbits, the amplitude of the vertical oscillation is smoothly increasing the growing Galactocentric distance, as shown on the top panel in Fig. 4.This result is consistent with the radially decaying exponential mass distribution model applied and is in agreements with the observed increase of the scale height with the increase of the Galactic radius of the disk stars, which is the well-known flaring of the Galactic disk (López-Corredoira et al. 2002;Amôres et al. 2017;Robin et al. 2022).
However, the continuous evolution of the z max -levels is interrupted in the regions of the increasing planar eccentricities when we move away from the rotation curve on the R-V φ plane in Fig. 4, in the direction of both the higher and lower tangential velocities.To understand this behaviour, it is worth noting that, in domains of regular oscillations, small changes in the initial conditions lead to small changes of the elements of the orbits, in particular, the amplitude z max .Consequently, the levels suffer slight displacements on the map when the initial conditions are gradually changed.However, in the vicinity of the resonances, small changes in the initial configurations produce large changes of the orbital elements, forming singular structures, such as resonant islands and stochastic layers.On the dynamical map in Fig. 4, these structures appear in the form of stalactites of different widths.
On the Poincaré map in Fig. 5, the same structures appear as chains of islands contrasting with the orbits, which regularly circulate around the origin.To construct this map, we pick up the orbital coordinates z and V z at the instants when the star is in the turning point of its orbit (p r = 0) of the minimal radial distance ( ṗr > 0).The initial configurations of the orbits were chosen along the Sun's L z -level and z and V z fixed at 0 and 20 km s −1 , respectively, as in Fig. 3.However, in this case, the initial Rvalues were extended to the very large radial distances, up to 27 kpc.
Comparing this map to the map of the equatorial motion shown in Fig. 3, we verify that the behaviour in the vertical direction is more complicated.We can observe the projections of the qualitatively different types: there are circulating orbits of regular motion at smaller stellar altitudes.With the growing height, the chains of islands of different width appear; they are separated by the stochastic layers.A large island dominates the region of high vertical velocities in the positive half-plane of the map and is surrounded by the sea of chaotic motion.To identify the cause of this behaviour, we apply a special spectral analysis method and describe the results obtained in the next section.

The e-z resonances
To understand the features of the resonant motion, we start plotting the amplitude of the vertical oscillation (z 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 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 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 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. 2002Michtchenko et al. , 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 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.3 and 5).Column (b): the 4/1 resonant orbit with e = 0.72 (cyan dots in Figs. 3 and 5).Column (c): the 6/1 resonant orbit with e = 1.01 (green dots in Fig. 5).Column (d): the 1/1 resonant orbit with e = 0.802 (black dot in Fig. 5).
Analyzing the evolution of the proper frequencies in Fig. 6 bottom, we verify that bifurcations occur when the beats between the two 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 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.2kpc.
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 ∼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 highorder resonances of the type f r / f z 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.

The dependence on the initial vertical velocity V z
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 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 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 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.8.The R-V z plane of initial conditions chosen along the Sun's L zlevel, with p R = 0 and z = 0.The equidistant levels of the maximal vertical deviation from the equatorial plane, z max , are shown by continuous curves; its value is mainly increasing with the increasing V z (see beside panel to associate the absolute z max -values).The red curve shows the location of the circular orbits with L z = 1995 kpc km s −1 as a function of the initial V z .The structures, which are characteristic of the e-z resonances, appears on both sides of the red curve, expanding their domains with increasing eccentricities.The blue curve shows loci of the very strong 1/1 resonance.The projection of the Sun's orbit is shown by cyan dots; its maximal z-amplitude is around 0.85 kpc.Beside panel: z max -values as a function of the initial V z , for circular orbits (red curve) and orbits with the initial R = 13 kpc (black curve).

Conclusions
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 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 φ 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 ⊙ ± 0.2 kpc and |p R | ≥ 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.

Fig. 1 .
Fig. 1.Top: Rotation curve of the Galaxy.The observed rotation curve is represented by the points with error bars, which indicate masers data from high-mass star-forming regions(Reid et al. 2019;Rastorguev et al. 2017), and H I and CO tangent-point data(Burton & Gordon 1978;Clemens 1985;Fich et al. 1989).The red curve shows the analytical rotation curve expressed by Eq. (10).The contributions of the modelled disk (continuous curve), bulge (short-dashed line) and dark halo (longdashed line) to the axisymmetric Galactic potential are also shown.Bottom: Vertical force K z at z = 1.1 kpc as a function of the Galactic radius.The black dots with error bars are from observation measurements byBovy & Rix (2013), while the grey solid curve shows the K z force radial profile at z = 1.1 kpc resulting from our Galactic model.

Fig. 2 .
Fig.2.Topology of the Hamiltonian (Eq.12) on the representative R-V φ plane.Continuous curves are the levels of the angular momentum L z = R × V φ , while dashed curves are the energy levels, calculated with the fixed initial values p R = 0, z = 0 and V z = 20 km s −1 .Rotation curve is shown by red dots.The projections of the three 3D orbits on the plane are: the Sun's orbit (cyan dots), one orbit starting at the configuration 1 and other at the configuration 2 (black dots).Note that the motion occurs always along a fixed L z -value, due to conservation of the angular momentum L z .The conservation of the energy during the motion defines two turning points of each orbit, both lying at the intersections of the corresponding L z -and energy levels.Top panel: the total energy (see Eq. 12), in arbitrary units, calculated along the rotation curve, as a function of R.

Fig. 4 .
Fig. 4. Representative R-V φ plane shown in Fig. 2, with the same rotation curve (red dots) and the Sun's orbit (cyan dots) along the L z -level (dashed line).The equidistant levels are used to represent the amplitude of the vertical oscillation of the stellar orbits (see top panel to associate the absolute z max -values), with the initial conditions z = 0 and V z = 20 km s −1 .Top panel: the amplitude of the vertical oscillation, z max , calculated along the rotation curve, as a function of R.

Fig. 5 .
Fig.5.Same as in Fig.3, except the projections were calculated at conditions p R = 0 and ṗR > 0 on the z-V z plane.Red, cyan and green dots show the projections of three periodic orbits on the Poincaré map corresponding to the 2/1, 4/1 and 6/1 resonances, respectively.The large island of the orbits oscillating around the center at z = 0 kpc and V z = 103.0km s −1 (black dot) is related to the strong 1/1 resonance (see later Fig.7).

Fig. 6 .
Fig. 6.Top: Amplitude of the vertical oscillation of stars as a function of the initial Galactocentric distance.The initial values of V φ are chosen along the Sun's L z -level, while z = 0 and V Z = 20 km s −1 .Bottom: Dynamical spectrum of the stellar orbits calculated along the Sun's L zlevel.The nominal positions of some resonances are shown by vertical dashed lines.

Fig.
Fig.8.The R-V z plane of initial conditions chosen along the Sun's L zlevel, with p R = 0 and z = 0.The equidistant levels of the maximal vertical deviation from the equatorial plane, z max , are shown by continuous curves; its value is mainly increasing with the increasing V z (see beside panel to associate the absolute z max -values).The red curve shows the location of the circular orbits with L z = 1995 kpc km s −1 as a function of the initial V z .The structures, which are characteristic of the e-z resonances, appears on both sides of the red curve, expanding their domains with increasing eccentricities.The blue curve shows loci of the very strong 1/1 resonance.The projection of the Sun's orbit is shown by cyan dots; its maximal z-amplitude is around 0.85 kpc.Beside panel: z max -values as a function of the initial V z , for circular orbits (red curve) and orbits with the initial R = 13 kpc (black curve).

Table 1 .
Physical and structural parameters of the disk components of the axisymmetric Galactic model