© 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 grand-design spiral are long-lived, 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 “long-lived 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 spiral-shaped potential wells.

Considerable effort has been devoted to the self-consistency 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 spiral-shaped 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 long-lived or self-consistent. The fact that the above studies show that such self-consistent solutions do exist justifies that it is an acceptable approximation to impose a spiral-shaped 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 Gaussian-shaped wells, or grooves, in the azimuthal direction. This model results in a much better self-consistency 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), re-scaled to new parameters R₀ and V₀.

Another line of research which is current in the literature, namely N-body 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 co-rotation radii.

We believe, however, that there is evidence that the co-rotation resonance usually stays at a same radius for a few billion years. This is suggested, for instance, by the step in metallicity at co-rotation in our Galaxy (see e.g. Lépine et al. 2011, their Fig. 4), and the breaks in the gradients of metallicity at co-rotation, 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 long-lived 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 phase-space 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 long-lived 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; Ferraz-Mello 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 rotation-curve 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 co-rotation 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 rotation-curve model of the Milky Way based on published observational data. For the Galactocentric distance of the Sun, we adopt R₀ = 8.0 kpc, which is based on the statistical analysis performed by Malkin (2013) on several R₀ measurements published in the literature. The circular velocity at R₀, which is the velocity V₀ of the local standard of rest (LSR), is chosen to satisfy the relation V₀ = R₀Ω_⊙−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₀ = 230 km s^-1.

For the rotation curve, we use the tangent-point data in H I from Burton & Gordon (1978) and Fich et al. (1989), and the CO-line tangent-point data from Clemens (1985). We also use data of maser sources associated with high-mass star-forming 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₀, V₀) 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 tangent-point 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 super-Keplerian fall-off. 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 Φ₀(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 Φ₀(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. 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 tangent-point data from Burton & Gordon (1978), Clemens (1985), and Fich et al. (1989); blue points indicate masers from high-mass star-forming 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 Gaussian-shaped groove profiles. In that approach, the surface density of a zero-thickness disk is represented analytically as the sum of an axisymmetric (unperturbed) surface density Σ₀(R) plus a small perturbation Σ₁(R,ϕ), which describes the spiral pattern in a rotating frame with angular speed Ω_p. The azimuthal coordinate in the rotating frame is ϕ = θ−Ω_pt, 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 ℋ₀ and Φ₁ 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 one-degree-of-freedom unperturbed Hamiltonian is given by Jacobi’s integral (4)where Φ₀(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.

The perturbation component of the Hamiltonian (3) is defined in Junqueira et al. (2013) as (5)where ζ₀ 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₀), 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₀ = 8 kpc and V₀ = 230 km s^-1. Since these constants are different from those used by Junqueira et al. (2013) (R₀ = 7.5 kpc, V₀ = 210 km s^-1 in that case), the parameters of the axisymmetric potential Φ₀(R) are also different from the previous ones. To keep the self-consistency 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 re-scaling of Φ₀(R) with the new pair of constants (R₀, V₀). In other words, we searched for the correspondence between the ratio η = (∂Φ₁/∂R)/(∂Φ₀/∂R) calculated with the parameters of the Junqueira et al. model and that calculated with the parameters adopted in the present work. Here, Φ₁ 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₀. This correspondence was obtained by setting new values for the amplitude ζ₀, 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 co-rotation radius R_CR = 8.54 kpc estimated from the condition V_rot(R_CR) = R_CR Ω_p. This value for the co-rotation radius is in agreement with the ratio R_CR/R₀ = 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/R0 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 ℋ₀(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 Φ₀ (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; Ferraz-Mello 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; Ferraz-Mello 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 four-dimensional. 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 ϕ₀ = ± 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 ϕ₀ = ± π and ϕ₀ = 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 ζ₀ = 6300 km² s^-2 kpc^-1 and the pitch angle i = + 14°. We choose a larger value of ζ₀ 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 pR and Jϕ fixed at their stationary values given in Eq. (12) and with the parameters from Table 1, except ζ0 = 6300 km2 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 Resc = 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 saddle-like points (black) lie on the black spiral branches. It is worth noting that the stable solutions define the co-rotation radius. For the parameters from Table 1, the co-rotation 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 Y-axis, at ϕ = ± 90°.

Figure 4 shows stellar orbits calculated for some initial conditions along a spiral branch defined by ϕ₀ = π. 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 co-rotation radius, is bound to this domain, and the amplitude of the R-oscillation grows when the initial conditions decrease from the co-rotation radius. The evolution of the spiral-branch orbits starting outside the co-rotation 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 co-rotation 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 Φ₀, 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 Φ₀ + Φ₁, defines the escape radius beyond which the spiral-branch 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 Φ₁ along the two spiral branches, one defined by ϕ₀ = 0 and the other by ϕ₀ = π. 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, pR = 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 two-degrees-of-freedom Hamiltonian system given by ℋ is four-dimensional, but the problem can be reduced to the systematic study of initial conditions on the plane, which is a projection of a two-dimensional 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, non-canonical variable V_θ instead of J_ϕ in order to have a better dynamical representation of orbits in terms of observables. The astrometric data from proper-motion measurements, along with the line-of-sight 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_ϕ = Ω_pR². 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 = Φ₀ + Φ₁−Ω_pR²/ 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 R-mode 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 R-mode 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 non-circular. The initial conditions, for which the amplitude of the R-mode 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 R-mode around non-circular 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 R-mode, 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₀ = 8 kpc and V_{θ, ⊙} = V₀ + v_⊙ = 242.24 km s^-1 (V₀ and v_⊙ were defined in Sect. 2). The family of orbits is centred around a periodic orbit, with zero-amplitude R-mode 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 R-mode 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 ℋ₀ 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² 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² Myr^-2. The conditions used in the construction of the surfaces of section were p_R = 0 and ṗ_R< 0, which reduce the four-dimensional phase space of the full problem to a two-dimensional 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 Φ₁ 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 co-rotation 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 co-rotation radius. In contrast, all stars starting in the region B have prograde orbits and evolve beyond the co-rotation 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 R-mode 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₁ will have in general | V_θ(R₁) | > | V_rot(R₁) | and the outer turning point R₂ will have | V_θ(R₂) | < | V_rot(R₂) | for both prograde and retrograde motion. This region is unbounded for |V_θ | → ± ∞.

Fig. 7 Surfaces of section on the R–pR 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 kpc2 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 co-rotation 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 kpc2 Myr-2 and Vθ< 0.

Fig. 9 Top: dynamical power spectrum: the evolution of the proper frequencies fR (red) and fϕ (black), their harmonics, and linear combinations along the spiral branch with ϕ0 = π, in logarithmic scale. The frequencies fR 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 fR (red) and fϕ (black) calculated via the equations in Appendix B.1, along the curve Jϕ(R) = ΩpR2, 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 × 1000-grid, 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 SAM-analysed 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 R-mode of oscillation tends to zero. We note that there is no family of periodic solutions composed of orbits with zero-amplitude J_ϕ-mode of oscillations.

On the contrary, the narrow black strips observed on the dynamical map in Fig. 6, right, are characterized by highly non-harmonic 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 (Ferraz-Mello 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 ϕ₀ = π, 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 zero-values at distances close to the co-rotation radius, around 8.54 kpc (vertical dashed line). We associate this behaviour with the co-rotation 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 co-rotation 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 co-rotation point, surrounded by layers of highly unstable motion.

The co-rotation resonance appears in Fig. 6, right, as the black strip emanating from the stationary solution and crossing the positive half-plane R–V_θ. The location of the co-rotation 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 Φ₀ (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 co-rotation resonance on the representative plane R–V_θ is shown by a thick red curve on the positive half-plane in Fig. 6, left.

In addition to the co-rotation 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 quasi-circular 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 half-plane in Fig. 6. According to the solutions for the extrema of the Hamiltonian function (3), there is no co-rotation point in this region (and, consequently, no co-rotation resonance), which is dominated by the strong 1/1 Lindblad resonance. The position of this resonance can be predicted from the unperturbed potential Φ₀ 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 half-plane 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 half-plane 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 co-rotation 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 ζ₀ 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 fR (red) and fϕ (black) along Vθ, calculated numerically with the initial values R = 7.0 kpc, pR = 0 km s-1, ϕ = 90°, and ζ0 = 0 (see Eq. (5)). We note that fϕ is presented by two branches. The location of the co-rotation 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 co-rotation resonance is crossed, at V_θ = 273 km s^-1 (the vertical line labelled “co-rotation” indicates this position). We note that, at the co-rotation 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 co-rotation motion is not expected for V_θ< 0 owing to the definition of the azimuthal angle as ϕ = θ−Ω_pt. 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 co-rotation 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 action-angle variables.

7.3. Dynamical portrait of the solar neighbourhood

The observations place the SSP close to the co-rotation 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 co-rotation 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 co-rotation 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 co-rotation 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 co-rotation 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 Sagittarius-Carina 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 four-arm 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 Sagittarius-Carina 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 co-rotation 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 Sagittarius-Carina 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 Y-axis 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 Sagittarius-Carina and Perseus arms, respectively. We believe that Sagittarius-Carina is the best choice since it is the most prominent and well-defined arm for most of the tracers.

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

Figure 13 shows three dynamical maps obtained for the initial conditions R₀,p_{R, ⊙},V_{θ, ⊙} corresponding to the SSP. The top panels show the plane Ω_p–ζ₀, 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 ζ₀ = 630 km² 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 light-coloured region corresponds to the co-rotation resonance. The greyscale is the same as in Fig. 6. We see that the SSP remains inside the co-rotation 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 ζ₀ in the plotted range. For ϕ = 142°, however, larger values of ζ₀ show secondary resonant islands associated with the co-rotation resonance. These results show that the SSP is still located inside the co-rotation resonance for small variations of the spiral amplitude ζ₀, pattern speed Ω_p, and azimuthal angle ϕ. For larger variations of these parameters, though, the SSP can either stay inside the co-rotation resonance or be in a chaotic zone. The possibility of being in a secondary resonance island associated with co-rotation is less likely.

8. Co-rotation 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 co-rotation 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 half-width, and the mean radius of the dip, respectively. The co-rotation dip corresponds to the small-amplitude 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₀ as traced by H II regions is strongly affected by the uncertainties on the distances of these sources. The above-mentioned velocity dip could then be an artefact of this method of deriving the rotation curve for R>R₀.

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 light-coloured region corresponds to the stable co-rotation 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₀ (see Sect. 2). The maser sources provide strong evidence for a velocity dip after R₀, 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 co-rotation 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 ring-like structure with mean radius slightly outside the solar circle. Zhang (1996), based on N-body 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 co-rotation radius. Based on numerical integrations of test-particle 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 co-rotation 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 co-rotation influences the stability of the orbits inside the co-rotation 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 co-rotation radius was calculated looking for equilibrium solutions of the full Hamiltonian (3) described in Sect. 5. The corresponding orbit was integrated numerically and SAM-analysed. The spectral number obtained was plotted on the parametric plane using a greyscale code. Light grey regions correspond to stable co-rotation solutions, whereas black regions correspond to chaotic co-rotation resonances. The red domain in Fig. 14 corresponds to large instabilities in the stellar motion. This region coincides with the region where κ²(R_CR) < 0 (see Eq. (B.7)), which means hyperbolic equilibrium solutions of ℋ₀, 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 co-rotation solution inside the stable region. In the case when the amplitude of the velocity dip beyond the solar radius R₀ = 8 kpc is larger or its width is smaller, the co-rotation 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 co-rotation dip. Thus, the region around co-rotation 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 Adip (in km s-1) and the thickness δdip (in kpc) of the co-rotation dip. The red cross indicates the values used in this paper. Bottom: dependence of RCR 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; Ferraz-Mello et al. 2005).

First, we modeled the rotation curve of our Galaxy based on observations of H I, CO, and maser sources. The rotation-curve model provides the axisymmetric component of the Galactic potential in the mid-plane of the disk. In this way, the potential model relies directly on rotation curve data, without any assumption of a mass-model 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 Gaussian-shaped potential wells (Junqueira et al. 2013), and re-scaled its parameters in order to be consistent with the parameters of the LSR, adopted in this paper, R₀ = 8 kpc and V₀ = 230 km s^-1. We believe that the Gaussian-grooved 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 spiral-galaxy model. This plane is representative of the whole phase space of the system and its analysis provides the precise positions of the co-rotation 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; Ferraz-Mello 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 co-rotation resonance. Depending on its current value of ϕ, it may be either in a stable or in a chaotic zone. The analysis of the co-rotation dip, present in the adopted model of the rotation curve, shows that our choice of parameters for the rotation curve corresponds to a stable co-rotation resonance, whereas slightly different parameters could drastically change the dynamics around co-rotation.

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 quasi-circular orbits, there is a great deal of evidence of non-circular 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 phase-space 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; Martinez-Medina et al. 2016). These topics, as well as a dynamical analysis of three-dimensional 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/10577-9. 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/54006-4) and the INCT-A. We acknowledge the anonymous referee for the detailed review and for the many helpful suggestions which allowed us to improve the manuscript.

References

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 non-perturbed calculations) as those shown in Fig. 6; however, the shape and degree of chaos of the resonances is different: while the Gaussian-grooved spiral perturbation (5) gives strips that are more highly defined for the resonances, the cosine perturbation shows high chaoticity for small-radii resonances. Outside co-rotation, however, the groove-shaped 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 groove-shaped 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 Gaussian-shaped spiral potential, as we can see by comparing Figs. 6 and A.1.

The other resonances of the system are more pronounced in the Gaussian-grooved 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 Gaussian-shaped 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 non-negligible, 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 Φ₀ (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₁ and R₂ 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, Δϕ = Δθ−Ω_pT_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 = −κ²(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; Ferraz-Mello 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 well-known 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 phase-space coordinates) are Fourier-transformed 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 Ferraz-Mello 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 noise-level 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 tangent-point data from Burton & Gordon (1978), Clemens (1985), and Fich et al. (1989); blue points indicate masers from high-mass star-forming 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/R0 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 pR and Jϕ fixed at their stationary values given in Eq. (12) and with the parameters from Table 1, except ζ0 = 6300 km2 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 Resc = 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, pR = 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–pR 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 kpc2 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 co-rotation 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 kpc2 Myr-2 and Vθ< 0.

In the text

Fig. 9 Top: dynamical power spectrum: the evolution of the proper frequencies fR (red) and fϕ (black), their harmonics, and linear combinations along the spiral branch with ϕ0 = π, in logarithmic scale. The frequencies fR 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 fR (red) and fϕ (black) calculated via the equations in Appendix B.1, along the curve Jϕ(R) = ΩpR2, 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 fR (red) and fϕ (black) along Vθ, calculated numerically with the initial values R = 7.0 kpc, pR = 0 km s-1, ϕ = 90°, and ζ0 = 0 (see Eq. (5)). We note that fϕ is presented by two branches. The location of the co-rotation and 1/1 resonances is shown by the vertical lines.
In the text

Fig. 12 Left: dynamical map of the neighbourhood of the co-rotation point calculated with the variables pR and ϕ fixed at zero. The red cross shows the position of the SSP integrated backwards until pR = 0. Right: same as in the left graph, except ϕ = 90°. Comparing the graphs, we note that the star would be inside the co-rotation 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 light-coloured region corresponds to the stable co-rotation 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 Adip (in km s-1) and the thickness δdip (in kpc) of the co-rotation dip. The red cross indicates the values used in this paper. Bottom: dependence of RCR 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