Issue 
A&A
Volume 597, January 2017



Article Number  A39  
Number of page(s)  17  
Section  Galactic structure, stellar clusters and populations  
DOI  https://doi.org/10.1051/00046361/201628895  
Published online  20 December 2016 
Modelling resonances and orbital chaos in disk galaxies
Application to a Milky Way spiral model
Instituto de Astronomia, Geofísica e Ciências Atmosféricas, USP, Rua do Matão 1226, 05508090 São Paulo, Brazil
email: tatiana@astro.iag.usp.br
Received: 10 May 2016
Accepted: 1 September 2016
Context. Resonances in the stellar orbital motion under perturbations from the spiral arm structure can play an important role in the evolution of the disks of spiral galaxies. The epicyclic approximation allows the determination of the corresponding resonant radii on the equatorial plane (in the context of nearly circular orbits), but is not suitable in general.
Aims. We expand the study of resonant orbits by analysing stellar motions perturbed by spiral arms with Gaussianshaped groove profiles without any restriction on the stellar orbital configurations, and we expand the concept of Lindblad (epicyclic) resonances for orbits with large radial excursions.
Methods. We define a representative plane of initial conditions, which covers the whole phase space of the system. Dynamical maps on representative planes of initial conditions are constructed numerically in order to characterize the phasespace structure and identify the precise location of the corotation and Lindblad resonances. The study is complemented by the construction of dynamical power spectra, which provide the identification of fundamental oscillatory patterns in the stellar motion.
Results. Our approach allows a precise description of the resonance chains in the whole phase space, giving a broader view of the dynamics of the system when compared to the classical epicyclic approach. We generalize the concept of Lindblad resonances and extend it to cases of resonant orbits with large radial excursions, even for objects in retrograde motion. The analysis of the solar neighbourhood shows that, depending on the current azimuthal phase of the Sun with respect to the spiral arms, a star with solar kinematic parameters (SSP) may evolve in dynamically distinct regions, either inside the stable corotation resonance or in a chaotic zone.
Conclusions. Our approach contributes to quantifying the domains of resonant orbits and the degree of chaos in the whole Galactic phasespace structure. It may serve as a starting point to apply these techniques to the investigation of clumps in the distribution of stars in the Galaxy, such as kinematic moving groups.
Key words: galaxies: spiral / Galaxy: kinematics and dynamics
© ESO, 2016
1. Introduction
A general consensus on the spiral arm structure of the Galaxy has not yet been reached. Questions like the physical nature of the arms, their exact location, the existence and position of resonances, their role in the evolution of the galactic disk, and the connection of the spiral arms with the bar have originated papers with quite diverse views. One of the main questions which divides the astronomical community is the lifetime of the spiral structure. Many authors consider that the spiral arms of a granddesign spiral are longlived, with lifetimes of a few billion years, and rotate like a rigid body (e.g. Bertin & Lin 1996, and references therein). The first model of spiral arms proposed by Lin & Shu (1964), in which the stars were treated as a fluid, adopted this view. However, the more recent “longlived arms” models are based on a very different understanding of the nature of the arms. They focus on stellar orbits, and have adopted a concept introduced by Kalnajs (1973), who considers that the arms are places where neighbouring galactic orbits of stars become close one to the other, resulting in regions of high stellar density or elongated spiralshaped potential wells.
Considerable effort has been devoted to the selfconsistency of models (e.g. Contopoulos & Grosbol 1986; Amaral & Lépine 1997; Pichardo et al. 2003; Martos et al. 2004; Junqueira et al. 2013, among others) that we briefly explain below. Series of closed stellar orbits are calculated in a potential which is the sum of the axisymmetric potential of the disk plus an imposed spiralshaped perturbation. These orbits naturally tend to present concentrations (high proximity) in some regions and, therefore, to produce new potential perturbations. If these newly created perturbations have a shape similar to the one that was originally imposed then the organization of the orbits will persist for long periods, so that the spiral structure is longlived or selfconsistent. The fact that the above studies show that such selfconsistent solutions do exist justifies that it is an acceptable approximation to impose a spiralshaped perturbation in order to perform studies of stellar orbits without worrying about the origin of this constant perturbation any more.
Most of the models aiming to describe stellar orbits adopt a perturbation potential given by a cosine law in the azimuthal direction, as explained in more detail in this paper. This comes from a long tradition and for simplicity. However, this potential is too smooth, with broad maxima and minima, while the density of orbits that are obtained vary sharply. Junqueira et al. (2013) proposed a new expression in which the arms are Gaussianshaped wells, or grooves, in the azimuthal direction. This model results in a much better selfconsistency since the width of the grooves derived from the stellar density is equal to that of the imposed perturbation. In the present work we adopt the potential of Junqueira et al. (2013), rescaled to new parameters R_{0} and V_{0}.
Another line of research which is current in the literature, namely Nbody numerical simulations (e.g. Sellwood & Carlberg 1984; Sellwood & Kahn 1991, among others), is quite distinct and competes with the one that we adopt in this paper. In general, these studies find that spiral arms are transient, appearing and disappearing in a recurrent way. In these studies, the positions of resonances are not well defined as they can move with time, and different parts of the structure can have, for instance, different corotation radii.
We believe, however, that there is evidence that the corotation resonance usually stays at a same radius for a few billion years. This is suggested, for instance, by the step in metallicity at corotation in our Galaxy (see e.g. Lépine et al. 2011, their Fig. 4), and the breaks in the gradients of metallicity at corotation, in external galaxies (Scarano & Lépine 2013).
Regarding the dynamical modelling of stellar orbits in spiral galaxies, most of the studies analyse the dynamics of stars on the equatorial plane of galaxy models when subjected to small azimuthal perturbations (such as central bars and spiral arms). Orbits in barred galaxies are reviewed in Contopoulos & Grosbøl (1989) and Contopoulos et al. (1996). Recently, numerical studies of chaos in longlived spiral and barred galaxies have been performed in Pichardo et al. (2003, 2004), Chakrabarty & Sideris (2008), Contopoulos (2009), Patsis (2012), and Moreno et al. (2015), among others. Prospects of observing chaos in disk galaxies are summarized in Grosbøl (2002, 2003, 2009). Also, many recent studies have focused on the phasespace structure of the solar neighbourhood (Dehnen 2000; Quillen 2003; Chakrabarty 2004, 2007; Pichardo et al. 2004; Chakrabarty & Sideris 2008; Antoja et al. 2008, 2009, 2011; Moreno et al. 2015).
The present work, which follows the approach described above, is a theoretical and numerical study of dynamical consequences of a longlived spiral structure, such as the presence of resonances and regions of chaotic orbits, in a galaxy described by the spiral potential perturbation proposed by Junqueira et al. (2013). We obtain dynamical maps and dynamical power spectra for the orbits restricted to the equatorial plane, which are tools widely used in celestial mechanics (e.g. Michtchenko et al. 2002; FerrazMello et al. 2005), and compare the results with analytical results for the axisymmetric potential. In particular, we show that our approach describes the positions of the resonance chains in the whole phase space with great precision, expanding the concept of Lindblad (epicyclic) resonances for orbits with large radial excursions, and for the case of retrograde orbits.
The model that we generate is useful in order to understand the observed features of our Galaxy, at least in the range of radii where the adopted perturbation seems to be similar to the observed spiral arms. This range should exclude the inner 3 kpc, since the adopted perturbation does not consider the presence of a bar, and possibly also the outer regions of the Galaxy, since we do not know up to what distance the potential perturbation that we adopt is a good description of the real one. We emphasize that the present study deals only with resonance effects induced by the spiral perturbation. Therefore, any possible resonant structures influenced by the central bar will not be considered in the present analysis.
The outline of the paper is as follows. In Sect. 2, using observational data, we develop a rotationcurve model to obtain the axisymmetric gravitational potential. In Sect. 3 we include the perturbing potential of spiral galaxies in the Hamiltonian model, and in Sect. 4 we briefly introduce the analytical and numerical techniques used in this study. The stationary solutions of the Hamiltonian are obtained and analysed in Sect. 5, where we also introduce the concept of spiral branches. Section 6 is devoted to the detailed analysis of the topology of the Hamiltonian on the representative plane. Section 7 presents dynamical maps of the phase space of the system under study, together with dynamical power spectra, which allow us to identify the main Lindblad resonances. The comparison of our results with those obtained via the epicyclic approximation is also done in this section. In Sect. 8, we briefly discuss the influence of the corotation dip in the rotation curve on the stellar dynamics. Finally, the discussions and conclusions are presented in Sect. 9, while Appendices A and B describe two additional points, namely the comparison between the Gaussian and cosine profiles of the spiral arms and the detailed description of the tools and methods employed in this study.
2. Rotation curve and axisymmetric potential
In this paper, we consider a realistic rotationcurve model of the Milky Way based on published observational data. For the Galactocentric distance of the Sun, we adopt R_{0} = 8.0 kpc, which is based on the statistical analysis performed by Malkin (2013) on several R_{0} measurements published in the literature. The circular velocity at R_{0}, which is the velocity V_{0} of the local standard of rest (LSR), is chosen to satisfy the relation V_{0} = R_{0}Ω_{⊙}−v_{⊙}, where Ω_{⊙} = 30.24 km s^{1} kpc^{1} is the angular rotation velocity of the Sun (Reid & Brunthaler 2004), and v_{⊙} = 12.24 km s^{1} is the peculiar velocity of the Sun in the direction of Galactic rotation (Schönrich et al. 2010). These values result in V_{0} = 230 km s^{1}.
For the rotation curve, we use the tangentpoint data in H I from Burton & Gordon (1978) and Fich et al. (1989), and the COline tangentpoint data from Clemens (1985). We also use data of maser sources associated with highmass starforming regions, obtained from Table 1 in Reid et al. (2014). From these compiled data, the rotation velocities and Galactic radii were calculated using the Galactic constants (R_{0}, V_{0}) adopted in this work.
Figure 1 shows the rotation curve of the Galaxy built up with the observational data described above: red points represent the H I and CO tangentpoint data, while blue points correspond to the maser sources data. In order to obtain a realistic model for the rotation curve, we fit the observational data by a convenient expression given by the sum of three exponentials in the form (1)with the factors that multiply the exponentials given in units of kilometers per second, and the factors in the arguments of the exponentials given in kiloparsecs. The numerical values of the model rotation curve in Eq. (1) were obtained after fitting the observational data: we minimized the sum of the squares of the residuals between V_{rot} and the measured rotation velocities, weighted by their respective uncertainties.
The smooth black curve in Fig. 1 represents the rotation velocity given by Eq. (1). It is worth noting a dip in the rotation curve centred at 8.9 kpc with a superKeplerian falloff. This velocity dip (represented by the third term in Eq. (1)) is related to a local minimum followed by a local maximum in the disk’s surface density distribution, as explained in Sect. 8 (see also Barros et al. 2013). However, the surface density in this range always remains positive, which would not be the case if the spherical version of Eq. (1) is inadvertently used.
The radial gradient of the axisymmetric potential Φ_{0}(R) is related to the rotation velocity V_{rot}(R) through (2)where V_{rot} is given by the rotation curve from Eq. (1). To obtain the axisymmetric potential Φ_{0}(R), we use the trapezium rule with adaptive step to solve numerically the integral of Eq. (2). The constant of integration is formally obtained from the limit condition limR → + ∞Φ_{0} = 0. In practice, ∞ can be replaced by a large value of R; for instance, 1000 kpc is found to be a good choice. We note that the maximum range of applicability of our model is constrained by a radius of ~30 kpc, as we show in Sect. 7.
Fig. 1 Rotation curve of the Galaxy. Red points indicate H I and CO tangentpoint data from Burton & Gordon (1978), Clemens (1985), and Fich et al. (1989); blue points indicate masers from highmass starforming regions from Reid et al. (2014). The black curve indicates the fitted rotation curve expressed by Eq. (1). 
3. Model
Junqueira et al. (2013) proposed a new description of the perturbed gravitational potential of spiral galaxies where the spiral arms have Gaussianshaped groove profiles. In that approach, the surface density of a zerothickness disk is represented analytically as the sum of an axisymmetric (unperturbed) surface density Σ_{0}(R) plus a small perturbation Σ_{1}(R,ϕ), which describes the spiral pattern in a rotating frame with angular speed Ω_{p}. The azimuthal coordinate in the rotating frame is ϕ = θ−Ω_{p}t, where θ is the angular coordinate with respect to the inertial frame.
The Hamiltonian which describes the stellar dynamics on the equatorial plane under perturbations of the spiral galaxy potential is written as (3)with ℋ_{0} and Φ_{1} being the unperturbed and perturbation components, respectively. The momenta p_{r} and J_{ϕ} = p_{θ} are the linear and angular momenta per unit mass, respectively. It is worth noting that J_{ϕ} is measured with respect to the inertial frame, but it is also a conjugate momentum to the canonical coordinate ϕ of the rotating frame: .
The onedegreeoffreedom unperturbed Hamiltonian is given by Jacobi’s integral (4)where Φ_{0}(R) is the galactic axisymmetric potential obtained in Sect. 2. In the expression above, the first term defines the kinetic energy T of a star and the second is a gyroscopic term.
Adopted spiral arms parameters.
The perturbation component of the Hamiltonian (3) is defined in Junqueira et al. (2013) as (5)where ζ_{0} is the perturbation amplitude, is the scale length of the spiral, σ is the width of the Gaussian profile in the galactocentric azimuthal direction, and f_{m}(R) is the shape function given by (6)where m is the number of arms fixed at 2 in our model, i is the pitch angle, R_{i} is an arbitrary radius chosen to adjust the phase of the spirals (here we adopt R_{i} = R_{0}), and γ is a phase angle, which does not influence the dynamics and can be initially fixed at 0. It is worth noting that the rotation curve of the Galaxy was fitted without analysing the mass components that explain it, and a kind of axisymmetric average of the spiral structure is one of the mass components that are naturally included in the fit. When we add the perturbation potential to the axisymmetric potential (Eq. (4)), the effect of the arms is counted a second time. The function of Eq. (5) corresponds to a positive density at all points, so that we are increasing the total mass of the disk. We estimate that the mass of the spiral arms corresponding to the potential of Eq. (5) is of the order of 5% of the mass of the disk. Following Junqueira et al. (2013), we consider that not removing the contribution of the spiral arms from the unperturbed Hamiltonian before going to Eq. (4) is a valid approximation.
The parameters in Eq. (5) are given in Table 1. Their choice was based on the adopted rotation curve with Galactic constants R_{0} = 8 kpc and V_{0} = 230 km s^{1}. Since these constants are different from those used by Junqueira et al. (2013) (R_{0} = 7.5 kpc, V_{0} = 210 km s^{1} in that case), the parameters of the axisymmetric potential Φ_{0}(R) are also different from the previous ones. To keep the selfconsistency of the spiral perturbation model of Junqueira et al. (2013) unaltered, the parameters of the perturbed Hamiltonian in Eq. (5) also needed to be updated. To do so, we considered that the ratio between the radial forces due to the spiral perturbation and due to the axisymmetric potential is a quantity that has to be preserved after the rescaling of Φ_{0}(R) with the new pair of constants (R_{0}, V_{0}). In other words, we searched for the correspondence between the ratio η = (∂Φ_{1}/∂R)/(∂Φ_{0}/∂R) calculated with the parameters of the Junqueira et al. model and that calculated with the parameters adopted in the present work. Here, Φ_{1} is the potential due to the perturbation, which is given by the expression in Eq. (5). Figure 2 shows the comparison between the ratio η calculated with the potential model of Junqueira et al. (2013; red dashed curve) and the ratio obtained with our updated models (blue dashed curve), in the direction of the azimuth ϕ = 90° for illustration. To make clear the correspondence between the curves, Fig. 2 shows the ratio η as a function of Galactocentric radius normalized with respect to R_{0}. This correspondence was obtained by setting new values for the amplitude ζ_{0}, the scale length , and the width σ of the spiral potential model, which are presented in Table 1.
The angular speed of the spiral pattern was directly measured by Dias & Lépine (2005), using the birthplaces of samples of observed open clusters. The authors found Ω_{p} = 25 ± 1 km s^{1} kpc^{1} (see also the review by Gerhard 2011, on the values of Ω_{p} estimated in the literature). In this work, we adopt Ω_{p} = 26 km s^{1} kpc^{1}. For the rotation curve of Eq. (1), it results in the corotation radius R_{CR} = 8.54 kpc estimated from the condition V_{rot}(R_{CR}) = R_{CR} Ω_{p}. This value for the corotation radius is in agreement with the ratio R_{CR}/R_{0} = 1.06 as determined in Dias & Lépine (2005).
Fig. 2 Variation of the ratio η = (∂Φ_{1}/∂R)/(∂Φ_{0}/∂R) as a function of normalized Galactic radius R/R_{0} and along the azimuth ϕ = 90°, for the potential model of Junqueira et al. (2013, red dashed curve) and for the potential model with updated parameters adopted in the present work (blue dashed curve). Φ_{0} and Φ_{1} are the unperturbed and perturbed parts of the Galactic potential, respectively. 
Cosine spiral pattern.
We also work in this paper with the widely used cosine perturbation given by (7)Our goal is to compare the dynamics under two different perturbations, Gaussian (5) and cosine, which will be done in Appendix A. We advance that the two potentials do not present significant qualitative differences in the dynamics, except near the main Lindblad resonances and for large radii.
4. Tools and methods
We briefly present in this section the tools employed in the orbital analysis throughout the paper. First we present the analytical tools, and then we describe the numerical techniques for detecting resonances and chaos.
The analytical study of motion is done in the context of the unperturbed Hamiltonian, described by the term ℋ_{0}(R,p_{r},J_{ϕ}) in Eq. (3). In particular, the two independent frequencies are estimated; they are given as (8)where T_{ϕ} and T_{R} are the azimuthal and radial periods in the rotating frame, respectively, and Ω_{ϕ} and Ω_{R} are the corresponding angular frequencies. We then have Ω_{ϕ} = Ω_{θ}−Ω_{p}, where Ω_{θ} = 2π/T_{θ} is the orbital frequency in the inertial frame. Resonances are then given by (9)with j and n coprime integers.
The orbital (θ) and radial (R) periods of a regular bounded orbit are estimated using the axially symmetric potential Φ_{0} (Binney & Tremaine 2008, see also Appendix B.1). In a frame rotating with angular velocity Ω_{p}, the azimuthal coordinate ϕ has a variation Δϕ along one radial period of the orbit. The positions of resonances calculated from the unperturbed problem then occur for (10)(see Appendix B.1 for a more detailed discussion). We note that this condition is valid regardless of the approximation of a nearly circular orbit; our calculations are done for general orbits with large radial span.
Regarding numerical techniques, we utilize the Spectral Analysis Method (SAM; Michtchenko et al. 2002; FerrazMello et al. 2005) in order to analyse the full Hamiltonian (3). First, we integrate numerically the equations of motion of stars in order to obtain the corresponding orbits. Each orbit is then Fourier transformed and the number of frequency peaks above a given threshold is calculated (here we adopt 5% of the amplitude of the highest peak). This number is referred to as spectral number N. The spectral number quantifies the chaoticity of each orbit. Regular orbits are conditionally periodic and therefore have a small number of significant frequency peaks (small N). On the other hand, chaotic orbits are not confined to an invariant torus. Their frequency spectrum is not discrete and is described by large values of N (Powell & Percival 1979).
A dynamical map associates a spectral number N to each initial condition on a representative plane of initial conditions. The result is a colour map, with a greyscale associated with N: lighter regions represent regular orbits, while darker regions correspond to chaotic motion. The construction of the dynamical maps is complemented by calculating a dynamical power spectrum along a family of initial conditions parameterized by one of the coordinates, for instance, the initial radial distance of the star. The spectrum shows the evolution of the main frequencies of the problem and their linear combinations as functions of the chosen parameter. The smooth evolution of frequencies is characteristic of regular motion, while the erratic scattering of the frequency values is characteristic of chaotic motion. The domains where one of the frequencies tends to zero accurately indicate the location of the separatrices between distinct regimes of motion, and resonance islands appear as regions between separatrices. A detailed explanation of the numerical methods used in this paper is given in Appendix B.2 (see also Powell & Percival 1979; Michtchenko et al. 2002; FerrazMello et al. 2005, and references therein).
5. Equations of motion, stationary solutions, and spiral branches
The phase space of the Hamiltonian system under study is fourdimensional. The equations of motion of a star in the gravitational potential given by the Hamiltonian (3) are written as (11)where , as defined in Eq. (2).
We look for stationary solutions of the Hamiltonian (3), which correspond to circular orbits in the inertial frame rotating with angular velocity Ω_{p}. The conditions for a star being at equilibrium in the rotating frame are dp_{R}/ dt = dR/ dt = dJ_{ϕ}/ dt = dϕ/ dt = 0. The last three conditions provide immediately that (12)where ϕ_{0} = ± nπ and n = 0,1,... The symmetry of the problem is 2 π/m, with the number of spiral arms m equal to 2 in our case.
The conditions (12) are visualized on the (X = Rcosϕ, Y = Rsinϕ)plane in Fig. 3, where they are plotted by the red and black lines, for ϕ_{0} = ± π and ϕ_{0} = 0,2π, respectively. In the background of this figure, we plot energy levels of the Hamiltonian (3) with J_{ϕ} given by the last condition of Eq. (12), calculated with the parameters from Table 1, except ζ_{0} = 6300 km^{2} s^{2} kpc^{1} and the pitch angle i = + 14°. We choose a larger value of ζ_{0} in order to enhance the visual effect of the perturbation. Also, the pitch angle i is chosen positive in order to follow the conventional maps of the spiral structure of the Milky Way (e.g. Georgelin & Georgelin 1976; Drimmel & Spergel 2001; Russeil 2003; Vallée 2013; Hou & Han 2014, among others), with the Galactic rotation in the clockwise direction from the viewpoint of an observer located towards the direction of the north Galactic pole. All further calculations in this paper are done with the values of the parameters from Table 1.
Figure 3 shows that the Hamiltonian topology follows the black and red curves: the black branches correspond to the locations associated with the maxima of spiral arm density, while the red branches cross the level of maximum energy. These configurations satisfy the WKB approximation (for a tightly wound spiral pattern; see Binney & Tremaine 2008), which requires that R df_{m}(R)/dR ≫ 1, with f_{m}(R) given in Eq. (6); for our adopted parameters m = 2 and i = −14°, we have R df_{m}(R)/dR = m/ tan(i) ≃ 8. The stationary solutions of the Hamiltonian (3) must belong to the branches mentioned above, which we refer to hereafter as spiral branches.
Fig. 3 Energy levels of the Hamiltonian function (3) on the (X = Rcosϕ, Y = Rsinϕ)plane, for p_{R} and J_{ϕ} fixed at their stationary values given in Eq. (12) and with the parameters from Table 1, except ζ_{0} = 6300 km^{2} s^{2} kpc^{1} and i = + 14°. Black lines correspond to the azimuthal minima of the Hamiltonian and are identified with the physical spiral arms. They correspond to the choice ϕ_{0} = 0, 2π in the second condition of Eq. (12). Red lines correspond to azimuthal maxima of ℋ; they have ϕ_{0} = ± π in the second condition of Eq. (12). Dots represent the stationary solutions of the Hamiltonian system, elliptic (red) and hyperbolic (black). The grey region is of escaping orbits. The escape radius is obtained as R_{esc} = 21.16 kpc. 
The calculation of the solution for stationary R is complicated and requires the implementation of some numerical procedures to resolve the condition dp_{R}/ dt = 0 in Eq. (11). The solutions obtained for equilibrium are shown by four large dots in Fig. 3. The equilibria belong to each one of the spiral branches. The two elliptic fixed points (red) correspond to the global maxima of the Jacobi constant ℋ (3) and lie on the red spiral branches, while the two hyperbolic saddlelike points (black) lie on the black spiral branches. It is worth noting that the stable solutions define the corotation radius. For the parameters from Table 1, the corotation radius is approximately equal to 8.54 kpc and the phase angle is . According to expression (6), the free phase angle γ initially fixed at zero, can now be assumed to be equal to −mϕ_{eq}, placing the stable solutions always on the Yaxis, at ϕ = ± 90°.
Figure 4 shows stellar orbits calculated for some initial conditions along a spiral branch defined by ϕ_{0} = π. The trajectories were obtained by integrating numerically the equations of motion (11) over only a few radial periods T_{R}; they are plotted in the phase subspace R–p_{R}. The location of the corotation radius is shown by the dashed line in Fig. 4. The evolution of the orbits starting on the spiral branch, inside the corotation radius, is bound to this domain, and the amplitude of the Roscillation grows when the initial conditions decrease from the corotation radius. The evolution of the spiralbranch orbits starting outside the corotation radius is analogous: they remain in this region and their amplitudes increase with increasing distance. At a given distance, the orbits become unbounded; this distance defines the escape radius, which is R_{esc} ≅ 21.16 kpc, for the parameters from Table 1. The shaded regions in Fig. 3 are domains of initial conditions leading to escape orbits.
Fig. 4 Examples of stellar orbits calculated with initial conditions along the spiral branch with ϕ_{0} = π (see Fig. 3). 
Outside the spiral branches, the escape velocity can be approximately obtained from the unperturbed problem in Eq. (B.1), at the condition E = 0, as (13)This is a good approximation whenever the amplitude of the azimuthal perturbation is small. In the region surrounding the corotation point, the escape boundary (given by Eq. (13)) is around  V  = 600 km s^{1}.
Figure 5 shows the energy components of the Hamiltonian (3), calculated along the spiral branches given by the conditions (12), as functions of the radius R. We show the kinetic energy T, the axisymmetric potential Φ_{0}, and the Jacobi integral ℋ. The radius value of 21.16 kpc, where the kinetic energy is equal to the modulus of the total gravitational potential energy Φ_{0} + Φ_{1}, defines the escape radius beyond which the spiralbranch stars are escaping.
Fig. 5 Dependence of the energy components of Hamiltonian (3) on the radius R, along the spiral branches (in logarithmic scale): T is kinetic energy, Φ_{0} is the axially symmetric potential, and ℋ is the Jacobi integral. Two spiral branches of the perturbation Φ_{1} are defined by ϕ_{0} = 0 and ϕ_{0} = π (see Eq. (12)). The amplitude of the minima of the potential (ϕ_{0} = 0) falls by two orders of magnitude from its maximum to the escape radius (outer vertical dashed line). The falling is almost exponential after R ≅ 6 kpc. We note that our model does not account for additional structures in the region below R = 3 kpc (inner dashed line), for instance, the Galaxy’s bar. 
It is interesting to observe in Fig. 5 the evolution of the perturbation potential Φ_{1} along the two spiral branches, one defined by ϕ_{0} = 0 and the other by ϕ_{0} = π. Both branches have maxima (of modulus) which lie close to 2 kpc and are three orders of magnitude smaller than the Jacobi integral ℋ. For increasing distances, the strength of the perturbation decreases exponentially.
6. Topology of ℋ on the representative plane and stellar orbits
Fig. 6 Left: topology of ℋ (3) on the R–V_{θ} plane of initial conditions, where V_{θ} = J_{ϕ}/R, p_{R} = 0, and ϕ = 90°. Right: dynamical maps on the same plane calculated for the time series J_{ϕ}(t). The shaded blue regions correspond to domains of initial conditions leading to orbits which go beyond 30 kpc (the region where our model is no longer applicable). The bar relates grey tones and the values of the spectral number N between 1 and 100 in logarithmic scale. All orbits with N> 100 are labelled in black. A typical orbit intersects this plane in two points corresponding to the minimum and maximum values of the coordinate R. 
To visualize the dynamical features of the Hamiltonian system (3), we introduce a representative plane of initial conditions. The space of initial conditions of the twodegreesoffreedom Hamiltonian system given by ℋ is fourdimensional, but the problem can be reduced to the systematic study of initial conditions on the plane, which is a projection of a twodimensional surface embedded in phase space. This plane may be chosen in such a way that all possible configurations of the system are included, and thus all possible regimes of motion of the system under study can be represented on it. Hereafter, we refer to this plane as a representative plane of initial conditions.
In this work, we define the representative plane in the following way. First, we fix the initial values of the momentum p_{R} at zero. Indeed, by definition, all bounded orbits should have two turning points defined by the condition p_{R} = 0. In the following, we fix the initial values of the azimuthal angle ϕ. We know that this angle is generally circulating; it oscillates around 90° (or −90°) when the system is close to the stable stationary solution shown in Fig. 3. In both cases, it goes through 90° (or −90°) for all initial conditions. Hence, without loss of generality, the angular variable ϕ can be initially fixed at 90°.
Now the topology of ℋ can be fully represented on the plane of initial values R–V_{θ}, where the stellar azimuthal velocity is defined as V_{θ} = J_{ϕ}/R. We choose a different, noncanonical variable V_{θ} instead of J_{ϕ} in order to have a better dynamical representation of orbits in terms of observables. The astrometric data from propermotion measurements, along with the lineofsight velocities and distances, can be transformed into V_{θ} data, allowing us to develop a framework in which future observations can be fit in.
It is worth emphasizing that the X–Y plane widely used in the stellar dynamics studies (see Fig. 3) cannot be chosen as a representative plane since, in its construction, the initial values of the angular momentum are fixed at the condition J_{ϕ} = Ω_{p}R^{2}. The Hamiltonian topology, in this case restricted to the stationary values of p_{r} and J_{ϕ}, is equivalent to the analysis of the level curves of the effective potential Φ_{eff} = Φ_{0} + Φ_{1}−Ω_{p}R^{2}/ 2. The energy function is written as , where (see e.g. Binney & Tremaine 2008; Barros et al. 2013). The equations of motion are , with Ω_{p} = Ω_{p}ẑ.
Figure 6 shows, in the left panel, the energy levels of the Hamiltonian on the plane R–V_{θ} (solid grey lines). The level which contains the stationary solution with coordinates R = 8.54 kpc and V_{θ} = 221.8 km s^{1} is shown by a dashed blue line. It separates the whole domain in three dynamically distinct regions, A, B, and C, which is explained below. The straight blue line shows the (R,V_{θ}) coordinates of the spiral branches.
The black curves present locations of the rotation curve V_{rot} (1) for both positive and negative values of V_{θ}. In the unperturbed problem, when J_{ϕ} remains constant in time, the rotation curves define the locations of circular orbits according to Eqs. (11), while the initial conditions outside the rotation curves lead to oscillations around the circular orbits with the corresponding values of J_{ϕ}. Hereafter we refer to this pattern of oscillation as a Rmode of motion. The orbits starting with V_{θ}> 0 (V_{θ}< 0) are prograde (retrograde) orbits and oscillate around the positive (negative) branch of the rotation curve.
When the perturbation due to spiral arms is introduced in the problem, the Rmode remains as a dominating mode of motion, while J_{ϕ} begins oscillating according to Eqs. (11). We refer to this new pattern of oscillation as a J_{ϕ}mode of motion. Owing to this additional mode of motion, the circular orbits of the unperturbed problem become noncircular. The initial conditions, for which the amplitude of the Rmode tends to zero, form a family of periodic orbits; if the perturbation is small, the periodic orbits are nearly circular orbits and their family is located in the vicinity of the positive and negative branches of the rotation curve, depending on the direction of the initial V_{θ}value. Therefore, the perturbed system oscillates in the Rmode around noncircular periodic orbits, with the exception of the equilibrium points.
Thus, a typical stellar motion in the full phase space can be represented as a combination of two independent oscillations, determined by the R and J_{ϕ}modes of motion with characteristic frequencies f_{R} and f_{ϕ}, respectively.
To illustrate the behaviour of a star in the Rmode, we show in Fig. 7 the surfaces of section in the phase subspace R–p_{R} calculated along the fixed value of J_{ϕ} = 1937.92 km s^{1} kpc. This value corresponds to the current coordinates of the Sun, R_{⊙} = R_{0} = 8 kpc and V_{θ, ⊙} = V_{0} + v_{⊙} = 242.24 km s^{1} (V_{0} and v_{⊙} were defined in Sect. 2). The family of orbits is centred around a periodic orbit, with zeroamplitude Rmode oscillation, corresponding to the chosen value of J_{ϕ}. The energy of the periodic orbit is minimal for a given J_{ϕ}, and increases when the oscillation amplitude increases. Knowing the current value of the Sun’s linear momentum, p_{R} = −11 km s^{1}, we plot the orbit of the star with solar kinematic parameters (SSP) and its position by the red curve and a cross symbol, respectively, in Fig. 7. The period of the Rmode oscillation of the SSP (inverse of f_{R}) is approximately 190 Myr. It is worth noting that, for small perturbations when J_{ϕ} ≈ const., the described behaviour in the phase subspace R–p_{R} is determined mainly by the unperturbed term ℋ_{0} in the Hamiltonian (3).
The typical behaviour of a star under the perturbation due to spiral arms is illustrated in the phase subspace ϕ–V_{θ} (the azimuthal velocity V_{θ} = J_{ϕ}/R) in Fig. 8. The top panel shows the surface of section for prograde orbits calculated along the energy level ℋ = −0.2298 kpc^{2} Myr^{2} close to the SSP energy. The bottom panel shows the surface of section for retrograde orbits along the energy level ℋ = −0.20 kpc^{2} Myr^{2}. The conditions used in the construction of the surfaces of section were p_{R} = 0 and ṗ_{R}< 0, which reduce the fourdimensional phase space of the full problem to a twodimensional surface defined by the variables ϕ and V_{θ}. We associate the behaviour of the system shown in Fig. 8 with the J_{ϕ}mode of oscillation of the system, with the corresponding frequency f_{ϕ}.
Since the perturbation determined by the term Φ_{1} in the Hamiltonian (3) is small (see Fig. 5), the variation of the azimuthal velocity is weak. It is enhanced only by the resonances which appear as chains of islands; some of the resonances are identified in Fig. 8, top. The most significant of these is the corotation resonance, which occurs at V_{θ}> 0, when the azimuthal frequency f_{ϕ} defined in Eqs. (8) tends to zero. The estimate of the current ϕvalue of the SSP is discussed in Sect. 7.3. The 1/1 resonance observed in Fig. 8, bottom, for retrograde orbits with V_{θ}< 0 is discussed in the next section.
Finally, the analysis of the topology of the Hamiltonian, shown in Fig. 6, left, allows us to infer several constraints on the stellar orbital evolution under the potential of the model galaxy. For instance, the orbits starting in the domains A and B are characterized by energies below the equilibrium value. Since the Jacobi integral ℋ (3) is conserved during the orbital evolution of the star, all initial conditions from region A lead to either prograde or retrograde orbits (apart from the resonant regions around V_{θ} = 0), which are all confined inside this region, with maximum distance smaller than the corotation radius. In contrast, all stars starting in the region B have prograde orbits and evolve beyond the corotation radius. We note that spiral branches, as well as prograde orbits starting along the rotation curve, belong to regions A and B.
In region C of Fig. 6, left, the stellar orbits have energies higher than the equilibrium energy: the amplitudes of the Rmode of oscillation are small when the star starts close to the families of periodic solutions and rapidly increase with increasing deviation from these families. The only constraint is that prograde orbits remain prograde forever, and also that retrograde orbits remain retrograde (apart from the resonant regions around V_{θ} = 0), since the variation of J_{ϕ} is very small under small perturbations. In particular, the inner turning point R_{1} will have in general  V_{θ}(R_{1})  >  V_{rot}(R_{1})  and the outer turning point R_{2} will have  V_{θ}(R_{2})  <  V_{rot}(R_{2})  for both prograde and retrograde motion. This region is unbounded for V_{θ}  → ± ∞.
Fig. 7 Surfaces of section on the R–p_{R} plane calculated with J_{ϕ} = 1937.92 km s^{1} kpc, which corresponds to the current parameters of the SSP (see text). The orbit of the SSP and its current position are shown by a red curve and a cross symbol, respectively. 
Fig. 8 Top: surfaces of section on the ϕ–V_{θ} plane calculated along the energy level ℋ = −0.2298 kpc^{2} Myr^{2}, close to the SSP energy, and V_{θ}> 0. Several resonances appear as chains of islands (the number of islands is equal to the order of the resonance); some are indicated. The estimated values of ϕ for the SSP (52° or 142°, see Sect. 7.3) place it inside the corotation resonance (see also Fig. 13); the period of the J_{ϕ}mode of oscillation of the SSP is about 2 Gyr. Bottom: same as the top panel, except that ℋ = −0.20 kpc^{2} Myr^{2} and V_{θ}< 0. 
Fig. 9 Top: dynamical power spectrum: the evolution of the proper frequencies f_{R} (red) and f_{ϕ} (black), their harmonics, and linear combinations along the spiral branch with ϕ_{0} = π, in logarithmic scale. The frequencies f_{R} and f_{ϕ} were calculated by analysing the time evolution of R and J_{ϕ}. The smooth evolution of the frequencies is characteristic of regular motion, while the erratic scattering of the points is characteristic of chaotic motion (see Sect. 4 and Appendix B). The main Lindblad resonances (see Sect. 7.1) are indicated by vertical lines and the corresponding ratio. The sign “−” corresponds to the condition Δϕ< 0 in Eq. (10) (see also Appendix B.1). Bottom: frequencies f_{R} (red) and f_{ϕ} (black) calculated via the equations in Appendix B.1, along the curve J_{ϕ}(R) = Ω_{p}R^{2}, with the unperturbed potential Φ_{0} (4). They match the fundamental frequencies of the perturbed problem (top panel), and the resonance radii (grey vertical lines) are very close to the corresponding radii of the perturbed problem. The bottom table shows the radii of resonances in the unperturbed problem. 
7. Dynamical maps, dynamical power spectra, and resonances
The dynamical map of the representative plane R–V_{θ} is shown in Fig. 6, right. During its construction, for each initial condition of the 1000 × 1000grid, the equations of motion (11) were integrated numerically over ~2600 radial epicyclic periods, with the value of R given by the initial conditions of the orbit. The values of p_{R} and ϕ are fixed at 0 and 90°, respectively. Next, the time series of the variable J_{ϕ} were SAManalysed and their spectral numbers N determined. The spectral numbers were then plotted on the representative plane using a greyscale code in logarithmic scale (see colour bar in Fig. 6, right). This scale varies from 0 to 2, corresponding to a spectral number N between 1 and 100. All values N> 100 are also labelled in black.
All initial conditions leading to orbits whose radius exceeds 30 kpc at some moment during the integrations were plotted in blue in Fig. 6, right. This value was chosen because several massive structures, such as the Sagittarius dwarf spheroidal galaxy (Helmi 2004; Law et al. 2009; Deg & Widrow 2013; Ibata et al. 2013) and the Magellanic clouds among other possible candidates, orbit our Galaxy at approximately this distance. Our model does not take into account interactions with these structures and therefore is no longer valid at regions beyond 30 kpc.
Comparing the two graphs in Fig. 6, we note the correlation between the white horizontal strips on the dynamical map (right panel) and the location of the rotation curves in the left panel, for both prograde and retrograde orbits. Indeed, the colour white in the greyscale corresponds to harmonic motion, characterized by a very small value of the spectral number N. This is also an important property of the periodic solutions, for which the amplitude of the Rmode of oscillation tends to zero. We note that there is no family of periodic solutions composed of orbits with zeroamplitude J_{ϕ}mode of oscillations.
On the contrary, the narrow black strips observed on the dynamical map in Fig. 6, right, are characterized by highly nonharmonic and chaotic stellar motions. They are associated with the dynamical phenomenon known as a resonance. A resonance occurs when one of the fundamental frequencies of the system or one of the linear combinations of these frequencies tends to zero. In this case, the topology of the phase space is transformed, giving rise to islands of stable resonant motion which are surrounded by the layers of chaotic motion associated with the separatrix of the resonance (FerrazMello 2007). Some examples of the formation of the chains of islands and the libration of the angle ϕ can be observed in Fig. 8.
The identification of the resonances observed on the dynamical map in Fig. 6, right, is done by applying the dynamical power spectrum approach described in Appendix B.2. The dynamical spectrum obtained shows the evolution of the proper frequencies, f_{R} (red) and f_{ϕ} (black), and their linear combinations along the spiral branch with ϕ_{0} = π, in Fig. 9, top. The dominant feature of the spectrum is the trend of the frequency f_{ϕ} associated with the J_{ϕ}mode of motion towards zerovalues at distances close to the corotation radius, around 8.54 kpc (vertical dashed line). We associate this behaviour with the corotation resonance which occurs at the condition , where is the averaged value of the angular velocity of the star along its trajectory. It is worth emphasizing here that, at the corotation radius, this resonance is crossed by the family of periodic solutions with V_{θ}> 0; thus, the intersection of the spiral branch with the family of periodic solutions gives rise to a stationary configuration of the system. The dynamical map constructed with ϕ = 90° (Fig. 6, right) shows a wide stable region around the corotation point, surrounded by layers of highly unstable motion.
The corotation resonance appears in Fig. 6, right, as the black strip emanating from the stationary solution and crossing the positive halfplane R–V_{θ}. The location of the corotation resonance can be estimated analytically through the unperturbed approximation assuming Δϕ = 0 in Eq. (10) (see also Appendix B.1). Figure 9, bottom, shows the fundamental frequencies f_{R} (red) and f_{ϕ} (black) calculated via Eq. (10) and the equations in Appendix B.1, with the unperturbed potential Φ_{0} (4). The comparison between the fundamental frequencies obtained analytically and those of the perturbed problem obtained numerically shows a good qualitative agreement, at least when the perturbations are small. Obtained this way, the approximate location of the corotation resonance on the representative plane R–V_{θ} is shown by a thick red curve on the positive halfplane in Fig. 6, left.
In addition to the corotation resonance, there is another kind of resonance whose effects can be observed in Fig. 9, top, as discontinuities in the smooth evolution of the frequencies in the dynamical power spectrum. These discontinuities are associated with the initial configurations for which the two fundamental frequencies of the problem become commensurable; in other words, they obey the relation jf_{R}−nf_{ϕ} ≅ 0, where j and n are simple integers. We say that, when the above condition is satisfied, a resonance occurs. We have identified these resonances and marked their positions by vertical lines and the corresponding ratio n/j in Fig. 9. The positions of the resonances, obtained for the unperturbed problem via Eq. (10) and Appendix B.1, are also shown in the bottom panel in Fig. 9, together with the approximated values of the resonance radii. It is worth noting a very good agreement between the results obtained analytically for the axisymmetric problem and numerically for the full problem.
7.1. Comparison with the epicyclic approximation
It is known that the epicyclic approximation allows us to predict the positions of resonances along nearly circular orbits which have small radial deviations (Lindblad 1974; Contopoulos & Grosbol 1986; Binney & Tremaine 2008; Efthymiopoulos 2010, and references therein). Those studies gave rise to the concept of Lindblad resonances (e.g. Binney & Tremaine 2008). Some recent studies base the interpretation of numerical results on the positions of resonances predicted by the epicyclic approximation (e.g. Dehnen 2000; Quillen 2003; Quillen & Minchev 2005; Minchev & Quillen 2006; Antoja et al. 2011; Gómez et al. 2013; Barros et al. 2013; Faure et al. 2014, among others).
In our approach, we expand the definition of these resonances to the whole representative domain of stellar orbits; moreover, since these newly determined resonances are of the same nature, we continue to refer to them as Lindblad resonances. We show in Fig. 10 the dynamical power spectrum of the orbits along V_{θ} = V_{rot}(R), with ϕ = 90°. The top panel shows the full dynamical power spectrum, calculated for the perturbed Hamiltonian. The bottom panel shows the unperturbed prediction for the radial and azimuthal frequencies f_{R} and f_{ϕ}, and the corresponding predictions for the resonance radii. We see that the fundamental frequencies associated with the R and J_{ϕ}modes agree with those obtained in the unperturbed case (compare the top and bottom panels of Fig. 10). The radii of the Lindblad (epicyclic) resonances match the unperturbed prediction, Eq. (B.8), the difference being due to the spiral perturbation.
The difference between the epicyclic (Lindblad) approximation (Fig. 10) and the resonances of the dynamical power spectrum calculated along the spiral branch (Fig. 9) are due to the adopted form of the curve V_{θ}(R) for the initial conditions on the R–V_{θ} dynamical map of Fig. 6. The epicyclic approximation considers V_{θ} given by circular motion in the unperturbed potential, V_{θ}(R) = V_{rot}(R).
Finally, the prediction for the Lindblad resonance locations (in our context) on the representative plane R–V_{θ}, are shown in the left panel in Fig. 6 by red curves and the corresponding ratio n/ j. These predictions were obtained from the unperturbed problem via Eq. (10) and the equations in Appendix B.1. The corresponding epicyclic resonances are precisely the intersections between quasicircular orbits (represented by a white strip around the rotation curve in Fig. 6, right) and the resonance chains, and their predictions from the unperturbed problem are given by the intersection between the rotation curve and the red curves in Fig. 6, left.
Fig. 10 Same as in Fig. 9, except that the frequencies were calculated along the rotation curve (1) with fixed ϕ = 90°. The epicyclic approximation was employed in order to estimate the frequencies analytically. 
7.2. Retrograde orbits
The dynamical portrait of the domain of retrograde orbits can be seen on the negative halfplane in Fig. 6. According to the solutions for the extrema of the Hamiltonian function (3), there is no corotation point in this region (and, consequently, no corotation resonance), which is dominated by the strong 1/1 Lindblad resonance. The position of this resonance can be predicted from the unperturbed potential Φ_{0} analytically, assuming Δϕ = 0 in Eq. (10) and noting that, for V_{θ}< 0, we have Δθ< 0 in Eq. (B.4). The 1/1 Lindblad resonance calculated this way is shown by a thick red curve on the negative halfplane in Fig. 6, left. The 1/1 Lindblad resonance on the dynamical map can be easily identified as a dominating black strip crossing the negative halfplane and intersecting the family of periodic orbits at R ≅ 4.6 kpc.
It is interesting to note in Fig. 6, right, that the 1/1 Lindblad resonance at V_{θ}< 0 seems to be a continuation of the corotation resonance at V_{θ}> 0. To understand the connection between the two resonances, we integrated numerically the equations of motion of the full system (11), setting the perturbation strength ζ_{0} at zero. The dynamical power spectrum of this unperturbed problem was calculated along V_{θ} at the fixed value of the radius R = 7 kpc and is shown in Fig. 11.
Fig. 11 Evolution of the unperturbed frequencies f_{R} (red) and f_{ϕ} (black) along V_{θ}, calculated numerically with the initial values R = 7.0 kpc, p_{R} = 0 km s^{1}, ϕ = 90°, and ζ_{0} = 0 (see Eq. (5)). We note that f_{ϕ} is presented by two branches. The location of the corotation and 1/1 resonances is shown by the vertical lines. 
Two fundamental frequencies obtained are present in the spectrum: the radial frequency f_{R} (red) and the azimuthal frequency f_{ϕ} (black). The evolution of the radial frequency is continuous, similar for prograde and retrograde orbits, with local minima coinciding with the position of periodic orbits corresponding to the given value R = 7 kpc.
The azimuthal frequency f_{ϕ} is presented by two branches whose evolution seems to be independent. The existence of the two independent components in the variation of the angle ϕ can be understood by analysing the corresponding equation of motion in Eqs. (11), which shows the dependence on the two variables, R and J_{ϕ}. The question which arises here is, Which of the two components of f_{ϕ} is fundamental? For V_{θ}> 0, the correct answer would be the lower branch, which tends to zero when the corotation resonance is crossed, at V_{θ} = 273 km s^{1} (the vertical line labelled “corotation” indicates this position). We note that, at the corotation resonance, the other component of f_{ϕ} becomes equal to f_{R} which, by definition, configures the 1/1 Lindblad resonance, as defined in this paper.
The dynamics is similar for the retrograde orbits, although the corotation motion is not expected for V_{θ}< 0 owing to the definition of the azimuthal angle as ϕ = θ−Ω_{p}t. In this case, the lower branch of f_{ϕ} tends to zero when the higher branch crosses the f_{R}family at V_{θ} = −109 km s^{1}. By definition, this event is associated with the 1/1 Lindblad resonance (the vertical line labelled “1/1 resonance” indicates its position). Thus, the connection between the corotation and 1/1 resonances indicates the same dynamical nature of both, which is also supported by the same topology of the surfaces of section in Fig. 8.
Finally, there are other features of retrograde orbits which can be observed on the dynamical map in Fig. 6, right. They are associated with Lindblad resonances of higher order, whose identification requires an additional study employing the models of resonant motion in actionangle variables.
7.3. Dynamical portrait of the solar neighbourhood
The observations place the SSP close to the corotation radius. To visualize the dynamical features of the domain surrounding this point, we construct two dynamical maps. In their construction, we use the 1000 × 1000 grid of initial conditions with p_{R} = 0. We test two critical values of the corotation resonance: ϕ = 90° and ϕ = 0. The map calculated with ϕ = 0 is shown in the left panel in Fig. 12, while the map calculated with ϕ = 90° is in the right panel. The time series of R was used in the construction of both maps.
We also plot the predicted position of the SSP on both planes in Fig. 12. For this, the current initial conditions of the SSP were integrated backwards in order to coincide with each plane at the condition p_{R} = 0.
Comparing the two maps, we note that the current azimuthal position of the SSP has a significant influence on the qualitative behaviour of its trajectory for a long time. In particular, if ϕ is close to zero (left panel), the SSP will be in a region dominated by dynamical instabilities which are associated with the separatrices of the corotation resonance. This means that its behaviour for future times is unpredictable.
On the other hand, if ϕ is close to 90° (right panel), the SSP will evolve inside the stable region of the corotation resonance, as shown in the top panel in Fig. 8. Its motion on the X–Y plane will be confined around the stable fixed point of the Hamiltonian (red point in Fig. 3). Assuming that the Sun is currently located in a void between two physical spiral arms, at radius R_{⊙} = 8 kpc, according to our model it is more likely that the SSP would evolve on a stable orbit inside the corotation resonance. On the contrary, if we assumed its position close to ϕ = 0, the SSP would be located very close to the arm, possibly lying inside of it, and evolving chaotically.
The angle ϕ for the SSP can be estimated as follows. In the observed spiral pattern proposed by Vallée (2013), the two arms nearest to the Sun are the SagittariusCarina arm, which passes about 1 kpc from the SSP at an inner radius, and the Perseus arm, at about 2 kpc from the SSP at an outer radius. Vallée presents a fourarm spiral structure, while in our model we adopt two arms. We consider that one of the arms in our model should coincide with either the observed SagittariusCarina arm or with the Perseus arm. Since both models use approximately the same pitch angle, we can match one of these two arms with one from our model by simple rotation. We take as reference the point at which the corotation circle crosses the arm in the clockwise rotating Galaxy (with positive pitch angle, see Fig. 3), and then calculate its angular position according to the above discussion. With this definition of ϕ, we find ϕ = 128° for the position of the Sun if we adopt the SagittariusCarina as one of the two arms in our model, and ϕ = 38° if the Perseus arm is adopted. However, since we use the negative value for the pitch angle in our calculations (see Table 1), we must reflect these obtained ϕvalues with respect to the Yaxis in Fig. 3, in order to obtain the correct values of ϕ in our model. In this way, the azimuthal position of the SSP will be at ϕ = 52° or ϕ = 142° for the SagittariusCarina and Perseus arms, respectively. We believe that SagittariusCarina is the best choice since it is the most prominent and welldefined arm for most of the tracers.
Fig. 12 Left: dynamical map of the neighbourhood of the corotation point calculated with the variables p_{R} and ϕ fixed at zero. The red cross shows the position of the SSP integrated backwards until p_{R} = 0. Right: same as in the left graph, except ϕ = 90°. Comparing the graphs, we note that the star would be inside the corotation resonance island or in a very chaotic area, depending on the angle ϕ. 
Figure 13 shows three dynamical maps obtained for the initial conditions R_{0},p_{R, ⊙},V_{θ, ⊙} corresponding to the SSP. The top panels show the plane Ω_{p}–ζ_{0}, calculated with the two different values of the azimuthal angle ϕ, 52° and 142°. The red crosses represent the location of the SSP according to the parameters used in our study. The bottom panel shows the plane Ω_{p}–ϕ, calculated with ζ_{0} = 630 km^{2} s^{2} kpc^{1}. We draw two red crosses, corresponding to the two possible values of the SSP coordinate ϕ described above. In all panels, the central lightcoloured region corresponds to the corotation resonance. The greyscale is the same as in Fig. 6. We see that the SSP remains inside the corotation resonance for values of Ω_{p} in the range [24.5–27] km s^{1} kpc^{1}. For ϕ = 52°, this conclusion is valid for any value of ζ_{0} in the plotted range. For ϕ = 142°, however, larger values of ζ_{0} show secondary resonant islands associated with the corotation resonance. These results show that the SSP is still located inside the corotation resonance for small variations of the spiral amplitude ζ_{0}, pattern speed Ω_{p}, and azimuthal angle ϕ. For larger variations of these parameters, though, the SSP can either stay inside the corotation resonance or be in a chaotic zone. The possibility of being in a secondary resonance island associated with corotation is less likely.
8. Corotation dip
In this section we slightly modify the potential studied in the previous sections, varying the parameters of the last term in Eq. (1). Hereafter, we refer to this term as the corotation dip and rewrite it in terms of new parameters, in the form (14)where A_{dip} (in units of km s^{1}), δ_{dip}, and R_{dip} (both in units of kpc) are the amplitude, the halfwidth, and the mean radius of the dip, respectively. The corotation dip corresponds to the smallamplitude valley which can be observed in the adopted form of the rotation curve in Fig. 1, near R = 9 kpc. In the literature, some authors have also included a similar term related to the velocity dip to fit the observed rotation curve of the Galaxy traced by H II regions (e.g. Clemens 1985; Sofue et al. 2009). However, Binney & Dehnen (1997) showed that the form of the rotation curve for R>R_{0} as traced by H II regions is strongly affected by the uncertainties on the distances of these sources. The abovementioned velocity dip could then be an artefact of this method of deriving the rotation curve for R>R_{0}.
Fig. 13 Top: dynamical maps on the parametric plane Ω_{p}–ζ_{0}, defined by the pattern speed and the amplitude of the spiral perturbation (Eq. (5)), respectively. The maps were calculated for the initial conditions of the SSP and ϕ = 52° (left), ϕ = 142° (right), corresponding to the possible SSP azimuthal positions (see Sect. 7.3). The red crosses represent the location of the SSP according to the parameters used in this paper. The central lightcoloured region corresponds to the stable corotation resonance. Bottom: same as in the top panel, except on the parametric plane Ω_{p}–ϕ, defined by the pattern speed and the azimuthal angle, respectively. 
In the present paper, we employ maser sources, with accurately determined distances and velocities, to trace the rotation curve beyond R_{0} (see Sect. 2). The maser sources provide strong evidence for a velocity dip after R_{0}, which makes us believe that this feature is real. Moreover, there is both theoretical and observational evidence for a deficit in the densities of stars and gas in the Galactic disk (in the corotation region). Such a deficit is intrinsically linked to the presence of the rotation velocity dip. For instance, Amôres et al. (2009) observed voids of H I gas density distributed in a ringlike structure with mean radius slightly outside the solar circle. Zhang (1996), based on Nbody simulations, showed that secular processes of energy and angular momentum transfer between stars and the spiral density wave induce the formation of a minimum in the stellar density, centred at the corotation radius. Based on numerical integrations of testparticle orbits in a spiral potential, Lépine et al. (2001), and also Barros et al. (2013) observed the formation of minima of density around the corotation radius of the simulated disks. Sofue et al. (2009) employed a ring density structure in their Galactic disk model to fit the observed dip in the rotation curve at R ~ 9 kpc.
This dip near corotation influences the stability of the orbits inside the corotation region. To show this, we test different values for dip amplitudes and widths in order to construct the parametric plane shown in Fig. 14. For each pair (A_{dip}, δ_{dip}), the corotation radius was calculated looking for equilibrium solutions of the full Hamiltonian (3) described in Sect. 5. The corresponding orbit was integrated numerically and SAManalysed. The spectral number obtained was plotted on the parametric plane using a greyscale code. Light grey regions correspond to stable corotation solutions, whereas black regions correspond to chaotic corotation resonances. The red domain in Fig. 14 corresponds to large instabilities in the stellar motion. This region coincides with the region where κ^{2}(R_{CR}) < 0 (see Eq. (B.7)), which means hyperbolic equilibrium solutions of ℋ_{0}, when the amplitude of the spiral perturbation tends to zero.
Our choice of the parameters A_{dip} = 11.86 km s^{1} and δ_{dip} = 0.48 kpc, marked by a red cross on the parametric plane in Fig. 14, places the corotation solution inside the stable region. In the case when the amplitude of the velocity dip beyond the solar radius R_{0} = 8 kpc is larger or its width is smaller, the corotation region becomes chaotic. It is worth noting, however, that the values of the parameters were determined through the fitting of the observed rotation curve of the Galaxy and using the specific exponential function (14) to describe the corotation dip. Thus, the region around corotation can be considered stable in the framework of our model and for the adopted set of the parameters.
Fig. 14 Top: parametric plane of the amplitude A_{dip} (in km s^{1}) and the thickness δ_{dip} (in kpc) of the corotation dip. The red cross indicates the values used in this paper. Bottom: dependence of R_{CR} on the dip’s amplitude for a fixed value of δ_{dip}. The red region corresponds to strong instabilities in the stellar motion. 
9. Conclusions
In this paper, we presented a new approach to the study of the dynamics of stars in spiral galaxies, based on widely used methods of celestial mechanics (Michtchenko et al. 2002; FerrazMello et al. 2005).
First, we modeled the rotation curve of our Galaxy based on observations of H I, CO, and maser sources. The rotationcurve model provides the axisymmetric component of the Galactic potential in the midplane of the disk. In this way, the potential model relies directly on rotation curve data, without any assumption of a massmodel description of the Galaxy’s components.
Our model is limited by the extent to which the fit to the rotation curve is reliable and by the existence of massive bodies orbiting the Galaxy, and thus will be valid up to a radius of R ~ 30 kpc. The inner region (R< 3 kpc) may also be excluded from the physical point of view, since there is evidence of the presence of a bar in this central region.
We employed a spiral potential described by Gaussianshaped potential wells (Junqueira et al. 2013), and rescaled its parameters in order to be consistent with the parameters of the LSR, adopted in this paper, R_{0} = 8 kpc and V_{0} = 230 km s^{1}. We believe that the Gaussiangrooved model is a more accurate description of the azimuthal profile of the spiral arms than the classical cosine perturbation. Indeed, as shown in Junqueira et al. (2013), a number of observed structures in the Milky Way can be associated with resonant orbits predicted by such a Gaussian spiral potential model.
We then introduced the representative plane R–V_{θ} of initial configurations, with fixed p_{R} = 0 and ϕ = 90°, to explore the stellar dynamics described by the adopted spiralgalaxy model. This plane is representative of the whole phase space of the system and its analysis provides the precise positions of the corotation and Lindblad resonances, as well as the degree of chaoticity of stellar orbits. It is also suitable for the analysis of existent and upcoming data regarding the peaks in the distribution function of stars in the Galactic phase space, such as kinematic moving groups (Antoja et al. 2008, 2010; Smith 2016).
The analysis of the dynamics was done by means of dynamical maps based on the SAM approach (Michtchenko et al. 2002; FerrazMello et al. 2005). A complementary study was done by means of dynamical power spectra of the orbits corresponding to curves V_{θ}(R) on the representative plane, particularly along spiral branches and along the rotation curve. This method allows us to identify the fundamental frequencies of the problem. The comparison with theoretical results in the axisymmetric potential permits us to correctly classify resonances in phase space. The results obtained agree with the predictions of classical Lindblad resonances given by the epicyclic approximation. The concept of Lindblad resonances was then extended over the whole phase space of the system.
We also analysed the solar neighbourhood on the representative plane of initial conditions, showing that a star with solar kinematic parameters is likely to evolve inside the corotation resonance. Depending on its current value of ϕ, it may be either in a stable or in a chaotic zone. The analysis of the corotation dip, present in the adopted model of the rotation curve, shows that our choice of parameters for the rotation curve corresponds to a stable corotation resonance, whereas slightly different parameters could drastically change the dynamics around corotation.
Finally, the approach presented in this work describes with high precision the resonance chains in the whole phase space and provides a more complete dynamical portrait than the classical epicyclic approximation. Although many stars in the solar neighbourhood seem to be on quasicircular orbits, there is a great deal of evidence of noncircular motion; there is also evidence of retrograde motion (Carney 1993), for which our model is more suitable. Thus, our approach may have an impact on the identification of dynamical signatures of resonant orbits and on the degree of chaos in the solar neighbourhood, and on its local phasespace structure (Chakrabarty 2007; Chakrabarty & Sideris 2008). It may also shed light on the origin of kinematic moving groups, believed to be related to resonances (Antoja et al. 2008, 2009, 2011; Moreno et al. 2015; MartinezMedina et al. 2016). These topics, as well as a dynamical analysis of threedimensional potentials, will be the subject of a forthcoming study.
Acknowledgments
This work was supported by the São Paulo State Science Foundation, FAPESP, and the Brazilian National Research Council, CNPq. R.S.S.V. acknowledges the financial support from FAPESP grant 2015/105779. D.A.B. acknowledges support from the Brazilian research agency CAPES, under the program PNPD. This work has made use of the facilities of the Laboratory of Astroinformatics (IAG/USP, NAT/Unicsul), whose purchase was made possible by FAPESP (grant 2009/540064) and the INCTA. We acknowledge the anonymous referee for the detailed review and for the many helpful suggestions which allowed us to improve the manuscript.
References
 Amaral, L. H., & Lépine, J. R. D. 1997, MNRAS, 286, 885 [NASA ADS] [CrossRef] [Google Scholar]
 Amôres, E. B., Lépine, J. R. D., & Mishurov, Y. N. 2009, MNRAS, 400, 1768 [NASA ADS] [CrossRef] [Google Scholar]
 Antoja, T., Figueras, F., Fernández, D., & Torra, J. 2008, A&A, 490, 135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Antoja, T., Valenzuela, O., Pichardo, B., et al. 2009, ApJ, 700, L78 [NASA ADS] [CrossRef] [Google Scholar]
 Antoja, T., Figueras, F., Torra, J., Valenzuela, O., & Pichardo, B. 2010, Lecture Notes and Essays in Astrophysics, 4, 13 [NASA ADS] [Google Scholar]
 Antoja, T., Figueras, F., RomeroGómez, M., et al. 2011, MNRAS, 418, 1423 [NASA ADS] [CrossRef] [Google Scholar]
 Barros, D. A., Lépine, J. R. D., & Junqueira, T. C. 2013, MNRAS, 435, 2299 [NASA ADS] [CrossRef] [Google Scholar]
 Bertin, G., & Lin, C. C. 1996, Spiral Structure in Galaxies: A density wave theory (MIT Press) [Google Scholar]
 Binney, J., & Dehnen, W. 1997, MNRAS, 287, L5 [NASA ADS] [Google Scholar]
 Binney, J., & Tremaine, S. 2008, Galactic Dynamics, 2nd edn. (Princeton, NJ: Princeton Univ. Press) [Google Scholar]
 Burton, W. B., & Gordon, M. A. 1978, A&A, 63, 7 [NASA ADS] [Google Scholar]
 Carney, B. W. 1993, in Galaxy Evolution. The Milky Way Perspective, ed. S. R. Majewski, ASP Conf. Ser., 49, 83 [Google Scholar]
 Chakrabarty, D. 2004, MNRAS, 352, 427 [NASA ADS] [CrossRef] [Google Scholar]
 Chakrabarty, D. 2007, A&A, 467, 145 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
 Chakrabarty, D., & Sideris, I. V. 2008, A&A, 488, 161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Clemens, D. P. 1985, ApJ, 295, 422 [NASA ADS] [CrossRef] [Google Scholar]
 Contopoulos, G. 1973, ApJ, 181, 657 [NASA ADS] [CrossRef] [Google Scholar]
 Contopoulos, G. 2009, Cel. Mech. Dyn. Astron., 104, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Contopoulos, G., & Grosbol, P. 1986, A&A, 155, 11 [NASA ADS] [Google Scholar]
 Contopoulos, G., & Grosbøl, P. 1989, A&ARv, 1, 261 [NASA ADS] [CrossRef] [Google Scholar]
 Contopoulos, G., Voglis, N., & Efthymiopoulos, C. 1996, in Barred Galaxies and Circumnuclear Activity, eds. A. Sandqvist, & P. O. Lindblad (Berlin: Springer Verlag), Lect. Notes Phys., 474, 19 [Google Scholar]
 Daniel, K. J., & Wyse, R. F. G. 2015, MNRAS, 447, 3576 [NASA ADS] [CrossRef] [Google Scholar]
 Deg, N., & Widrow, L. 2013, MNRAS, 428, 912 [NASA ADS] [CrossRef] [Google Scholar]
 Dehnen, W. 2000, AJ, 119, 800 [NASA ADS] [CrossRef] [Google Scholar]
 Dias, W. S., & Lépine, J. R. D. 2005, ApJ, 629, 825 [NASA ADS] [CrossRef] [Google Scholar]
 Drimmel, R., & Spergel, D. N. 2001, ApJ, 556, 181 [NASA ADS] [CrossRef] [Google Scholar]
 Efthymiopoulos, C. 2010, Eur. Phys. J. Special Topics, 186, 91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Faure, C., Siebert, A., & Famaey, B. 2014, MNRAS, 440, 2564 [NASA ADS] [CrossRef] [Google Scholar]
 FerrazMello, S. 2007, Canonical Perturbation Theories – Degenerate Systems and Resonance, Astrophys. Space Sci. Lib., 345 [Google Scholar]
 FerrazMello, S., Michtchenko, T. A., Beaugé, C., & Callegari, N. 2005, in Chaos and Stability in Planetary Systems, eds. R. Dvorak, F. Freistetter, & J. Kurths (Berlin: Springer Verlag), Lect. Notes Phys., 683, 219 [Google Scholar]
 Fich, M., Blitz, L., & Stark, A. A. 1989, ApJ, 342, 272 [NASA ADS] [CrossRef] [Google Scholar]
 Georgelin, Y. M., & Georgelin, Y. P. 1976, A&A, 49, 57 [NASA ADS] [Google Scholar]
 Gerhard, O. 2011, Mem. Soc. Astron. It. Suppl., 18, 185 [Google Scholar]
 Gómez, G. C., Pichardo, B., & Martos, M. A. 2013, MNRAS, 430, 3010 [NASA ADS] [CrossRef] [Google Scholar]
 Grosbøl, P. 2002, Space Sci. Rev., 102, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Grosbøl, P. 2003, in Galaxies and Chaos, eds. G. Contopoulos, & N. Voglis (Berlin: Springer Verlag), Lect. Notes Phys., 626, 201 [Google Scholar]
 Grosbøl, P. 2009, Order and Chaos in Spiral Galaxies Seen through their Morphology, eds. G. Contopoulos, & A. P. Patsis (Berlin, Heidelberg: Springer), 23 [Google Scholar]
 Helmi, A. 2004, MNRAS, 351, 643 [NASA ADS] [CrossRef] [Google Scholar]
 Hou, L. G., & Han, J. L. 2014, A&A, 569, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ibata, R., Lewis, G. F., Martin, N. F., Bellazzini, M., & Correnti, M. 2013, ApJ, 765, L15 [NASA ADS] [CrossRef] [Google Scholar]
 Junqueira, T. C., Lépine, J. R. D., Braga, C. A. S., & Barros, D. A. 2013, A&A, 550, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kalnajs, A. J. 1973, PASA, 2, 174 [Google Scholar]
 Law, D. R., Majewski, S. R., & Johnston, K. V. 2009, ApJ, 703, L67 [NASA ADS] [CrossRef] [Google Scholar]
 Lépine, J. R. D., Mishurov, Y. N., & Dedikov, S. Y. 2001, ApJ, 546, 234 [NASA ADS] [CrossRef] [Google Scholar]
 Lépine, J. R. D., Cruz, P., Scarano, Jr., S., et al. 2011, MNRAS, 417, 698 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, C. C., & Shu, F. H. 1964, ApJ, 140, 646 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Lindblad, P. O. 1974, in The Formation and Dynamics of Galaxies, ed. J. R. Shakeshaft, IAU Symp., 58, 399 [Google Scholar]
 Malkin, Z. 2013, in IAU Symp., ed. R. de Grijs, 289, 406 [Google Scholar]
 MartinezMedina, L. A., Pichardo, B., Moreno, E., Peimbert, A., & Velazquez, H. 2016, ApJ, 817, L3 [NASA ADS] [CrossRef] [Google Scholar]
 Martos, M., Hernandez, X., Yáñez, M., Moreno, E., & Pichardo, B. 2004, MNRAS, 350, L47 [NASA ADS] [CrossRef] [Google Scholar]
 Michtchenko, T. A., Lazzaro, D., FerrazMello, S., & Roig, F. 2002, Icarus, 158, 343 [NASA ADS] [CrossRef] [Google Scholar]
 Minchev, I., & Quillen, A. C. 2006, MNRAS, 368, 623 [NASA ADS] [CrossRef] [Google Scholar]
 Moreno, E., Pichardo, B., & Schuster, W. J. 2015, MNRAS, 451, 705 [NASA ADS] [CrossRef] [Google Scholar]
 Papayannopoulos, T. 1979a, A&A, 77, 75 [NASA ADS] [Google Scholar]
 Papayannopoulos, T. 1979b, A&A, 79, 197 [NASA ADS] [Google Scholar]
 Patsis, P. A. 2012, Int. J. Bifurcation Chaos, 22, 1230029 [NASA ADS] [CrossRef] [Google Scholar]
 Pichardo, B., Martos, M., Moreno, E., & Espresate, J. 2003, ApJ, 582, 230 [NASA ADS] [CrossRef] [Google Scholar]
 Pichardo, B., Martos, M., & Moreno, E. 2004, ApJ, 609, 144 [NASA ADS] [CrossRef] [Google Scholar]
 Powell, G. E., & Percival, I. C. 1979, J. Phys. A Math. Gen., 12, 2053 [Google Scholar]
 Quillen, A. C. 2003, AJ, 125, 785 [Google Scholar]
 Quillen, A. C., & Minchev, I. 2005, AJ, 130, 576 [NASA ADS] [CrossRef] [Google Scholar]
 Reid, M. J., & Brunthaler, A. 2004, ApJ, 616, 872 [NASA ADS] [CrossRef] [Google Scholar]
 Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2014, ApJ, 783, 130 [Google Scholar]
 Russeil, D. 2003, A&A, 397, 133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Scarano, S., & Lépine, J. R. D. 2013, MNRAS, 428, 625 [NASA ADS] [CrossRef] [Google Scholar]
 Schönrich, R., Binney, J., & Dehnen, W. 2010, MNRAS, 403, 1829 (SBD) [NASA ADS] [CrossRef] [Google Scholar]
 Sellwood, J. A., & Carlberg, R. G. 1984, ApJ, 282, 61 [NASA ADS] [CrossRef] [Google Scholar]
 Sellwood, J. A., & Kahn, F. D. 1991, MNRAS, 250, 278 [NASA ADS] [Google Scholar]
 Smith, M. C. 2016, in Astrophys. Space Sci. Lib., 420, eds. H. J. Newberg, & J. L. Carlin, 113 [Google Scholar]
 Sofue, Y., Honma, M., & Omodaka, T. 2009, PASJ, 61, 227 [NASA ADS] [Google Scholar]
 Vallée, J. P. 2013, Int. J. Astron. Astrophys., 3, 20 [Google Scholar]
 Zhang, X. 1996, ApJ, 457, 125 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Comparison with the cosine profile of the spiral arms
Many works treat spiral arms as perturbations with the azimuthal dependence in the cosine form, first introduced by Lin & Shu (1964) and then largely adopted in the literature (e.g. Contopoulos 1973; Papayannopoulos 1979a,b; Chakrabarty & Sideris 2008; Daniel & Wyse 2015, and references therein).
Figure A.1 shows the R–V_{θ} dynamical map calculated with the perturbation Φ_{1,cos}(R,ϕ) given in Eq. (7), for the same parameters listed in Table 1. We see that the resonance chains remain at approximately the same radii (as expected from the nonperturbed calculations) as those shown in Fig. 6; however, the shape and degree of chaos of the resonances is different: while the Gaussiangrooved spiral perturbation (5) gives strips that are more highly defined for the resonances, the cosine perturbation shows high chaoticity for smallradii resonances. Outside corotation, however, the grooveshaped perturbations give rise to more chaotic orbits.
Particular attention must be paid to the inner 2/1 Lindblad resonance (ILR). While the ILR for the cosine perturbation shows strong chaotic behaviour, as we see in Fig. A.1, the grooveshaped spirals seem to give rise to robust resonance islands, as we see for the ILR in Fig. 6. We can conclude that cosine perturbations tend to generate a wider chaotic region near the ILR, especially close to the rotation curve. The width of the ILR is also greater than in the Gaussianshaped spiral potential, as we can see by comparing Figs. 6 and A.1.
The other resonances of the system are more pronounced in the Gaussiangrooved case. These results may be related to the differences in the azimuthal profiles of the two spiral potentials. But they may also be related to the differences in the radial dependence of the width of the spirals in the two models: the Gaussianshaped spiral arms are modelled with a constant width along the Galactic disk, whereas in the cosine perturbation it increases with radius.
Fig. A.1 Dynamical map of the representative plane of initial conditions, for the cosine perturbation given in Eq. (7) calculated for the time series J_{ϕ}(t) and the fixed angle ϕ = 90°. The blue shaded regions correspond to the initial conditions of the orbits which go beyond 30 kpc, regions where we consider that our model is no longer applicable. 
Appendix B: Analytical and numerical techniques
Appendix B.1: Motion in the unperturbed system: azimuthal and radial frequencies
For completeness, in this appendix we present the equations used in the analytical predictions of resonances. We also compare the analytical predictions with the predictions from the epicyclic approximation, showing how quantities calculated from orbits which are far from circular are approximated by the corresponding quantities predicted for nearly circular orbits. The epicyclic approximation is not valid if deviations from circular motion become nonnegligible, as we see for instance in Figs. 6 and 9.
The orbital (θ) and radial (R) periods of a regular bounded orbit can be estimated using the axially symmetric potential Φ_{0} (Binney & Tremaine 2008). The angular momentum J_{ϕ} and the energy (B.1)are integrals of motion which can be calculated from the orbit’s initial conditions. It can be shown that (Binney & Tremaine 2008, Eq. (3.17)) (B.2)where the turning points R_{1} and R_{2} are the solutions of Eq. (B.1) with Ṙ = 0. We also have that the azimuthal period T_{θ}, in the inertial frame, is (Binney & Tremaine 2008, Eq. (3.19)) (B.3)where the azimuthal angle variation Δθ during this radial oscillation is (Binney & Tremaine 2008, Eq. (3.18b)) (B.4)In the rotating frame, Δϕ = Δθ−Ω_{p}T_{R}. We then have (B.5)Resonances in this approximation then occur for (B.6)with n,j coprime integers and Δθ and T_{R} given by the equations above.
It is usual to approximate the above quantities by their linearized values along the rotation curve, that is: Ω_{θ} ≈ Ω (i.e., Ω_{ϕ} ≈ Ω−Ω_{p}) and Ω_{R} ≈ κ, where the angular velocity of circular orbits Ω(R) is given by Ω(R) = V_{rot}/R and κ(R) is the (radial) epicyclic frequency, defined by the equation ṗ_{R} = −κ^{2}(R_{c})·(R−R_{c}) for each R_{c} (Binney & Tremaine 2008). It is then given by (Binney & Tremaine 2008, Eq. (3.79a)) (B.7)The “epicyclic approximation” is therefore defined by the linearized problem: Ω_{ϕ} = Ω−Ω_{p} and Ω_{R} = κ. The resonances which appear in this approximation are known as Lindblad resonances. They are given by (B.8)with j and n coprime integers. Assuming Δϕ> 0 (see Eq. (B.6)), we have (B.9)where in the last term we approximate the frequencies by the corresponding nearly circular (epicyclic) motion. It is worth remarking that this approximation is only valid for orbits with small deviations from the rotation curve. Therefore, the resonance radii given by the epicyclic approximation (B.8) need not be equal to the resonance radii obtained from the full numerical, perturbed calculations, or to the axisymmetric estimations (B.6) if the corresponding orbit is not nearly circular. On the other hand, the unperturbed approximation (see Eqs. (B.2)–(B.6)) describes the positions of resonances with good accuracy, as shown in the text.
Appendix B.2: Numerical techniques: spectral analysis method
The stellar motion in the full galactic potential (given by Eq. (3)) is analysed by purely numerical techniques applying the Spectral Analysis Method (Michtchenko et al. 2002; FerrazMello et al. 2005), frequently used in planetary and asteroidal sciences.
The SAM is used to distinguish between regular and chaotic domains on the phase space of Hamiltonian systems and is based on the wellknown features of power spectra (Powell & Percival 1979). It involves two main steps. The first step is the numerical integration of the equations of motion defined by the Hamiltonian (3), while the second step consists of the spectral analysis of the output of the numerical integrations. The series giving the time variation of stellar orbital elements (e.g. the canonical phasespace coordinates) are Fouriertransformed using a standard fast Fourier transform (FFT) algorithm.
The Fourier transforms of the output make it possible to distinguish between regular and irregular motions. Indeed, regular motions are conditionally periodic and any orbital element ele(t) depends on time as a function (B.10)where f is a frequency vector whose components are the fundamental frequencies of motion and k ∈ Z^{D}, where D is the dimension of the configuration space (or the number of degrees of freedom) of the problem. When the independent frequencies are constant in time, the power spectrum of ele, obtained from its Fourier transform, consists of the lines associated with the independent frequencies, whose number is equal to the number of degrees of freedom of the dynamical system, and with those corresponding to higher harmonics and linear combinations of the independent frequencies. The number of peaks that are above an arbitrarily defined “noise” level is defined as the spectral number N and is associated with this ele. In this paper, we define the noise level as 5% of the largest peak in the spectrum.
Chaotic motions are no longer conditionally periodic and the fundamental frequencies of the system vary in time. The power spectrum of chaotic motion is not discrete, showing broadband components; therefore, the spectral number N will be characterized by very high values in this case.
Thus, the spectral number N can be used to quantify the chaoticity of the system in the following way: small values of N correspond to regular motion, while large values of N indicate the onset of chaos. Using this criterion, we construct dynamical maps of the system under study.
Dynamical maps associate a spectral number N of a given orbital element to each point on the plane of initial conditions by means of the spectral analysis of the corresponding orbit described above. They are useful tools in the characterization of the phase space and identification of the principal dynamical mechanisms acting on stellar motion. Once N is determined for all initial conditions on a grid covering the plane of initial conditions, we plot it using a greyscale that varies logarithmically from white (N = 1) to black (N maximum). Since high values of N indicate the onset of chaos, the grey tones indicate different degrees of stochasticity of solutions with initial conditions starting on the plane: lighter regions correspond to regular motion and darker tones indicate chaotic motion.
The spectral analysis allows an efficient identification of the main oscillations contained in the trajectory of the star. Power spectra are plots of the amplitude of the Fourier transform against frequency. In order to see how the spectra change when initial conditions vary, we adopt a dynamic power spectrum, which is constructed plotting the frequencies of the peaks on the power spectra as functions of one parameter describing a family of solutions. Over the domains of regular motion, the proper frequencies vary continuously when the initial conditions are gradually varied. When resonances are approached, the topology of the phase space suffers qualitative transformations and the frequency evolution shows a discontinuity characterized by the erratic scatter of values when chaotic layers associated with separatrices are crossed. In systems with two degrees of freedom, where the chaoticity may be studied with the help of surfaces of section (Poincaré maps), dynamic power spectra allow us to understand the dynamics in areas where the maps show intricate features or where the features are too thin to be visible.
It is worth emphasizing that the SAM is simple to implement. Its performance, tested on the standard map problem with known solutions, is robust (for details, see FerrazMello et al. 2005, Appendix A). Based on a solid criterion, the method is infallible in distinguishing regular from chaotic motion. The resolution of the SAM is defined by the chosen grid: for smaller spacings between the grid nodes, the structure of the dynamical map exhibits finer details. Also, since the independent modes of motion of the problem under study may show distinct behaviour, it is advisable to calculate and visualize the spectral number N for all variables of the problem.
The spectral number N also depends on the integration time span; in other words, an orbit classified as regular can appear as chaotic if a larger time span is used in the integrations. Indeed, if the diffusion rate of the main frequencies is below the Fourier transform resolution (defined by the time span), the spectral analysis methods are unable to detect chaos. Thus, the chosen total integration time should be large enough to allow the chaos generated by Lindblad resonances to be distinguished. On the other hand, the dependence of N on the noiselevel is not critical, and the choice of 5% of the largest peak in the spectrum yields good results.
All Tables
All Figures
Fig. 1 Rotation curve of the Galaxy. Red points indicate H I and CO tangentpoint data from Burton & Gordon (1978), Clemens (1985), and Fich et al. (1989); blue points indicate masers from highmass starforming regions from Reid et al. (2014). The black curve indicates the fitted rotation curve expressed by Eq. (1). 

In the text 
Fig. 2 Variation of the ratio η = (∂Φ_{1}/∂R)/(∂Φ_{0}/∂R) as a function of normalized Galactic radius R/R_{0} and along the azimuth ϕ = 90°, for the potential model of Junqueira et al. (2013, red dashed curve) and for the potential model with updated parameters adopted in the present work (blue dashed curve). Φ_{0} and Φ_{1} are the unperturbed and perturbed parts of the Galactic potential, respectively. 

In the text 
Fig. 3 Energy levels of the Hamiltonian function (3) on the (X = Rcosϕ, Y = Rsinϕ)plane, for p_{R} and J_{ϕ} fixed at their stationary values given in Eq. (12) and with the parameters from Table 1, except ζ_{0} = 6300 km^{2} s^{2} kpc^{1} and i = + 14°. Black lines correspond to the azimuthal minima of the Hamiltonian and are identified with the physical spiral arms. They correspond to the choice ϕ_{0} = 0, 2π in the second condition of Eq. (12). Red lines correspond to azimuthal maxima of ℋ; they have ϕ_{0} = ± π in the second condition of Eq. (12). Dots represent the stationary solutions of the Hamiltonian system, elliptic (red) and hyperbolic (black). The grey region is of escaping orbits. The escape radius is obtained as R_{esc} = 21.16 kpc. 

In the text 
Fig. 4 Examples of stellar orbits calculated with initial conditions along the spiral branch with ϕ_{0} = π (see Fig. 3). 

In the text 
Fig. 5 Dependence of the energy components of Hamiltonian (3) on the radius R, along the spiral branches (in logarithmic scale): T is kinetic energy, Φ_{0} is the axially symmetric potential, and ℋ is the Jacobi integral. Two spiral branches of the perturbation Φ_{1} are defined by ϕ_{0} = 0 and ϕ_{0} = π (see Eq. (12)). The amplitude of the minima of the potential (ϕ_{0} = 0) falls by two orders of magnitude from its maximum to the escape radius (outer vertical dashed line). The falling is almost exponential after R ≅ 6 kpc. We note that our model does not account for additional structures in the region below R = 3 kpc (inner dashed line), for instance, the Galaxy’s bar. 

In the text 
Fig. 6 Left: topology of ℋ (3) on the R–V_{θ} plane of initial conditions, where V_{θ} = J_{ϕ}/R, p_{R} = 0, and ϕ = 90°. Right: dynamical maps on the same plane calculated for the time series J_{ϕ}(t). The shaded blue regions correspond to domains of initial conditions leading to orbits which go beyond 30 kpc (the region where our model is no longer applicable). The bar relates grey tones and the values of the spectral number N between 1 and 100 in logarithmic scale. All orbits with N> 100 are labelled in black. A typical orbit intersects this plane in two points corresponding to the minimum and maximum values of the coordinate R. 

In the text 
Fig. 7 Surfaces of section on the R–p_{R} plane calculated with J_{ϕ} = 1937.92 km s^{1} kpc, which corresponds to the current parameters of the SSP (see text). The orbit of the SSP and its current position are shown by a red curve and a cross symbol, respectively. 

In the text 
Fig. 8 Top: surfaces of section on the ϕ–V_{θ} plane calculated along the energy level ℋ = −0.2298 kpc^{2} Myr^{2}, close to the SSP energy, and V_{θ}> 0. Several resonances appear as chains of islands (the number of islands is equal to the order of the resonance); some are indicated. The estimated values of ϕ for the SSP (52° or 142°, see Sect. 7.3) place it inside the corotation resonance (see also Fig. 13); the period of the J_{ϕ}mode of oscillation of the SSP is about 2 Gyr. Bottom: same as the top panel, except that ℋ = −0.20 kpc^{2} Myr^{2} and V_{θ}< 0. 

In the text 
Fig. 9 Top: dynamical power spectrum: the evolution of the proper frequencies f_{R} (red) and f_{ϕ} (black), their harmonics, and linear combinations along the spiral branch with ϕ_{0} = π, in logarithmic scale. The frequencies f_{R} and f_{ϕ} were calculated by analysing the time evolution of R and J_{ϕ}. The smooth evolution of the frequencies is characteristic of regular motion, while the erratic scattering of the points is characteristic of chaotic motion (see Sect. 4 and Appendix B). The main Lindblad resonances (see Sect. 7.1) are indicated by vertical lines and the corresponding ratio. The sign “−” corresponds to the condition Δϕ< 0 in Eq. (10) (see also Appendix B.1). Bottom: frequencies f_{R} (red) and f_{ϕ} (black) calculated via the equations in Appendix B.1, along the curve J_{ϕ}(R) = Ω_{p}R^{2}, with the unperturbed potential Φ_{0} (4). They match the fundamental frequencies of the perturbed problem (top panel), and the resonance radii (grey vertical lines) are very close to the corresponding radii of the perturbed problem. The bottom table shows the radii of resonances in the unperturbed problem. 

In the text 
Fig. 10 Same as in Fig. 9, except that the frequencies were calculated along the rotation curve (1) with fixed ϕ = 90°. The epicyclic approximation was employed in order to estimate the frequencies analytically. 

In the text 
Fig. 11 Evolution of the unperturbed frequencies f_{R} (red) and f_{ϕ} (black) along V_{θ}, calculated numerically with the initial values R = 7.0 kpc, p_{R} = 0 km s^{1}, ϕ = 90°, and ζ_{0} = 0 (see Eq. (5)). We note that f_{ϕ} is presented by two branches. The location of the corotation and 1/1 resonances is shown by the vertical lines. 

In the text 
Fig. 12 Left: dynamical map of the neighbourhood of the corotation point calculated with the variables p_{R} and ϕ fixed at zero. The red cross shows the position of the SSP integrated backwards until p_{R} = 0. Right: same as in the left graph, except ϕ = 90°. Comparing the graphs, we note that the star would be inside the corotation resonance island or in a very chaotic area, depending on the angle ϕ. 

In the text 
Fig. 13 Top: dynamical maps on the parametric plane Ω_{p}–ζ_{0}, defined by the pattern speed and the amplitude of the spiral perturbation (Eq. (5)), respectively. The maps were calculated for the initial conditions of the SSP and ϕ = 52° (left), ϕ = 142° (right), corresponding to the possible SSP azimuthal positions (see Sect. 7.3). The red crosses represent the location of the SSP according to the parameters used in this paper. The central lightcoloured region corresponds to the stable corotation resonance. Bottom: same as in the top panel, except on the parametric plane Ω_{p}–ϕ, defined by the pattern speed and the azimuthal angle, respectively. 

In the text 
Fig. 14 Top: parametric plane of the amplitude A_{dip} (in km s^{1}) and the thickness δ_{dip} (in kpc) of the corotation dip. The red cross indicates the values used in this paper. Bottom: dependence of R_{CR} on the dip’s amplitude for a fixed value of δ_{dip}. The red region corresponds to strong instabilities in the stellar motion. 

In the text 
Fig. A.1 Dynamical map of the representative plane of initial conditions, for the cosine perturbation given in Eq. (7) calculated for the time series J_{ϕ}(t) and the fixed angle ϕ = 90°. The blue shaded regions correspond to the initial conditions of the orbits which go beyond 30 kpc, regions where we consider that our model is no longer applicable. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.