  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.
