© ESO 2020

1. Introduction

The manifold theory of a spiral structure in barred galaxies (Romero-Gomez et al. 2006; Voglis et al. 2006) predicts bisymmetric spirals emanating from the end of galactic bars as a result of the outflow of matter connected with the unstable dynamics around the bar’s Lagrangian points L₁ and L₂ (see also Danby 1965). The manifold theory has been expressed in two versions, namely the “flux-tube” version (Romero-Gomez et al. 2006, 2007; Athanassoula et al. 2009a,b; Athanassoula 2012) and the “apocentric manifold” version (Voglis et al. 2006; Tsoutsis et al. 2008, 2009; Harsoula et al. 2016). In both versions, the orbits of stars along the manifolds are chaotic, thus the manifolds provide a skeleton of orbits that support “chaotic spirals” (Patsis 2006). Furthermore, the theory predicts that the orbital flow takes place in a direction preferentially along the spirals. This is in contrast to standard density wave theory (Lin & Shu 1964; see Binney & Tremaine 2008), which predicts a regular orbital flow forming “precessing ellipses” (Kalnajs 1973) that intersect the spirals. This difference has been proposed as an observational criterion to distinguish chaotic (manifold) from regular (density wave) spirals (Patsis 2006).

Among a number of objections to manifold theory (see introduction in Efthymiopoulos et al. 2019 for a review as well as Font et al. 2019; Díaz-García et al. 2019 for more recent references), a common one stems from the long-recognized possibility that the bar and the spirals could rotate at different pattern speeds (Sellwood & Sparke 1988; see Binney 2013; Sellwood 2014 for reviews). Since it is only possible to define the unstable equilibria L₁ and L₂ when the potential is static in a frame corotating with the bar, manifold spirals emanating from L₁ and L₂ are also necessarily static in the same frame, hence, they should corotate with the bar. This prediction seems hard to reconcile with either observations (Vera-Villamizar et al. 2001; Boonyasait et al. 2005; Patsis et al. 2009; Meidt et al. 2009; Speights & Westpfahl 2012; Antoja et al. 2014; Junqueira et al. 2015; Speights & Rooke 2016) or simulations (Sellwood & Sparke 1988; Little & Carlberg 1991; Rautiainen & Salo 1999; Quillen 2003; Minchev & Quillen 2006; Dubinski et al. 2009; Quillen et al. 2011; Minchev et al. 2012; Baba et al. 2013; Roca-Fàbrega et al. 2013; Font et al. 2014; Baba 2015). On the other hand, Efthymiopoulos et al. (2019) recently found empirically that the manifold spirals, which are computed in an N-body simulation by momentarily “freezing” the potential and making all calculations in a frame that rotates with the instantaneous pattern speed of the bar, reproduce the time-varying morphology of the N-body spirals rather well. This occurs despite the fact that multiple patterns are demonstrably present in the latter simulation. Such a result points toward the question of whether manifold theory is possible to generalize under the presence of more than one patterns in the disk that rotate at different speeds.

We hereafter present such a generalization of the manifold theory under the explicit assumption of multiple pattern speeds. In particular, we consider models of barred galaxies in which the disk potential at time t, considered in cylindrical coordinates (ρ, ϕ) in a frame corrotating with the bar, has the form

$$\begin{array}{cc}\hfill V_{\mathrm{disc}}(\rho ,\varphi ,t)=& V_0(\rho )+V_{\mathrm{bar}}(\rho ,\varphi )\hfill \\ \hfill & +V_2(\rho ,\varphi -({\mathrm{\Omega }}_2-{\mathrm{\Omega }}_{\mathrm{bar}})t)\hfill \\ \hfill & +V_3(\rho ,\varphi -({\mathrm{\Omega }}_3-{\mathrm{\Omega }}_{\mathrm{bar}})t)+\cdots \hfill \end{array}$$(1)

In such a model, the bar rotates with pattern speed Ω_bar, while Ω₂, Ω₃, etc., are the pattern speeds (incommensurable between themselves and with Ω_bar) of additional non-axisymmetric perturbations, modeled by the potentials V₂, V₃, etc. The latter can be secondary spiral, ring, or bar-like modes, assumed to be of smaller amplitude than the principal bar. Our main result is as follows. Through Hamiltonian perturbation theory, we demonstrate that spiral-like invariant manifolds exist in the above generalized potential given by Eq. (1). These manifolds emanate from special orbits, which can be regarded as continuations of the unstable Lagrangian equilibria of the potential V₀ + V_bar after “turning on” the terms V₂, V₃, etc. Specifically, adding one more term V₂ with the frequency Ω₂, we can prove the existence of two periodic solutions in the bar’s rotating frame, each with a period equal to π|Ω₂ − Ω_bar|⁻¹. These orbits are hereafter denoted as GL₁ and GL₂ (standing for “generalized L₁ and L₂”); they form epicycles of size O(V₂) with a center near the bar’s end, and they reduce to the usual Lagrangian points L₁ and L₂ when V₂ goes to zero. In the same way, adding two terms V₂ and V₃ with incommensurable frequencies Ω₂ and Ω₃ allows one to prove the existence of two quasi-periodic orbits (also denoted GL₁ and GL₂) that reduce to the points L₁ and L₂ in the limit when both V₂ and V₃ go to zero. Each of the orbits GL₁ and GL₂ then appears as an epicyclic oscillation with the two frequencies |Ω₂ − Ω_bar|, |Ω₃ − Ω_bar|, thus forming a Lissajous figure around L₁ or L₂. One can continue in the same way by adding more frequencies. The key result, shown in Sect. 2 below, is that, independently of the number M of assumed extra frequencies, the orbital phase space in the neighborhood of the generalized orbits GL₁ and GL₂ admits a decomposition into a center ^M + 1× saddle linearized dynamics (see Gómez et al. 2001). Hence, the orbits GL₁ and GL₂ possess stable and unstable manifolds, which generalize the manifolds of the points L₁ and L₂ of the pure bar model. In particular, the unstable manifolds of the orbits GL₁ and GL₂ support trailing spirals and ring-like structures. In fact, these manifolds have a similar morphology as the manifolds of the L₁ and L₂ points, but they are no longer static in the frame that corotates with the bar. In physical terms, the manifolds of the GL₁ and GL₂ orbits adapt their form in time periodically or quasi-periodically in order to follow the additional patterns present in the disk. An explicit numerical example of this behavior is given in Sect. 3, referring to a Milky-Way-like model in which a bar and spirals rotate at different pattern speeds. In this example, we explicitly computed the orbits GL₁ and GL₂ as well as the manifolds emanating from them. Remarkably, despite using only a coarse fitting approach, the manifolds provide a good fit to the model’s imposed spirals. Since the latter have a relative rotation with respect to the bar, one has to test this fitting at different phases of the displacement of the spirals with respect to the bar’s major axis. We find that the fitting is good at nearly all phases except when close to ±π/2. A possible connection of this effect with the phenomenon of recurrent spirals is discussed.

The paper is structured as follows. Section 2 gives the general theory on the existence of the generalized unstable Lagrangian orbits GL₁ and GL₂ and their manifolds under the presence of multiple pattern speeds. Section 3 presents our numerical example in which the manifolds are constructed under a different pattern speed of the bar and the spirals. Section 4 gives the summary of the results and conclusions. Mathematical details on the Lie method used for the series computations are given in the Appendix A.

2. Theory

The Hamiltonian in the disk plane in a galactic model with the potential (1) can be written as

$$\begin{array}{c}\hfill H=H_0+H_1,\end{array}$$(2)

where H₀ is the axisymmetric plus the bar Hamiltonian:

$$\begin{array}{c}\hfill H_0=\frac{1}{2}(p_{\rho }^2+\frac{p_{\varphi }^2}{2{\rho }^2})-{\mathrm{\Omega }}_{\mathrm{bar}}{p}_{\varphi }+V_0(\rho )+V_{\mathrm{bar}}(\rho ,\varphi ),\end{array}$$(3)

and H₁ = V₂ + V₃ + ⋯. The pair (ρ, ϕ) are the test particle’s cylindrical coordinates in a frame that rotates with angular speed Ω_bar, while (p_ρ, p_ϕ) are the values of the radial velocity and angular momentum per unit mass of the particle in the rest frame. The dependence of the Hamiltonian H on time can be formally removed by introducing extra action-angle pairs. Setting the angles ϕ₂ = (Ω₂ − Ω_bar)t, ϕ₃ = (Ω₃ − Ω_bar)t, etc., with conjugate dummy actions I₂, I₃, etc., we arrive at the extended Hamiltonian

$$\begin{array}{cc}\hfill H=& \frac{1}{2}(p_{\rho }^2+\frac{p_{\varphi }^2}{2{\rho }^2})-{\mathrm{\Omega }}_{\mathrm{bar}}{p}_{\varphi }+V_0(\rho )+V_{\mathrm{bar}}(\rho ,\varphi )\hfill \\ \hfill & +({\mathrm{\Omega }}_2-{\mathrm{\Omega }}_{\mathrm{bar}})I_2+({\mathrm{\Omega }}_3-{\mathrm{\Omega }}_{\mathrm{bar}})I_3+\dots \hfill \\ \hfill & +V_2(\rho ,\varphi -{\varphi }_2)+V_3(\rho ,\varphi -{\varphi }_3)+\dots ,\hfill \end{array}$$(4)

which yields the same equations of motion as the Hamiltonian (2).

The Hamiltonian H₀ gives rise to the two Lagrangian equilibrium points: $L_1=({\rho }_{\text{L1}},{\varphi }_{\text{L1}},p_{\rho ,\text{L1}}=0,p_{\varphi ,\text{L1}}={\Omega }_{\text{bar}}{\rho }_{\text{L1}}^2)$ and $L_2=({\rho }_{\text{L2}},{\varphi }_{\text{L2}},p_{\rho ,\text{L2}}=0,p_{\varphi ,\text{L2}}={\Omega }_{\text{bar}}{\rho }_{\text{L2}}^2)$ such that ∂H₀/∂ρ = ∂H₀/∂ϕ = ∂H₀/∂p_ρ = ∂H₀/∂p_ϕ = 0 at the points L₁ and L₂. Focusing on, say, L₁, and by defining δρ = ρ − ρ_L1, δϕ = ϕ − ϕ_L1, and J_ϕ = p_ϕ − p_ϕ, L1, the Hamiltonian H₀ can be expanded around the phase-space coordinates of the point L₁. This yields H₀ = const + H_0, 2 + H_0, 3 + …, where H_0, 2, H_0, 3, … are quadratic, cubic, etc. in the variables (δρ, δϕ, p_ρ, J_ϕ). Then, with a standard procedure (see the Appendix A), we can define a linear transformation

$$\begin{array}{c}\hfill \left(\begin{array}{c}\delta \rho \\ \delta \varphi \\ p_{\rho }\\ J_{\varphi }\end{array}\right)=\mathcal{A}·\left(\begin{array}{c}U\\ Q\\ V\\ P\end{array}\right),\end{array}$$(5)

where 𝒜 is a 4 × 4 matrix with constant entries such that the quadratic part of the Hamiltonian H₀ takes a diagonal form in the new variables (U, Q, V, P)

$$\begin{array}{c}\hfill H_{0,2}=\lambda UV+\frac{\kappa }{2}(Q^2+P^2),\end{array}$$(6)

with λ and κ real constants. The matrix 𝒜 satisfies the symplectic condition 𝒜 ⋅ 𝒥 ⋅ 𝒜^T = 𝒜^T ⋅ 𝒥 ⋅ 𝒜 = 𝒥, where 𝒥 is the 4 × 4 fundamental symplectic matrix. The constants λ and κ are related to the eigenvalues of the variational matrix

$$\begin{array}{c}\hfill \mathcal{M}=\mathcal{J}·\left(\begin{array}{c}\frac{{\partial }^2{H}_0}{\partial {\rho }^2}\frac{{\partial }^2{H}_0}{\partial \rho \partial \varphi }\frac{{\partial }^2{H}_0}{\partial \rho \partial p_{\rho }}\frac{{\partial }^2{H}_0}{\partial \rho \partial p_{\varphi }}\\ \frac{{\partial }^2{H}_0}{\partial \varphi \partial \rho }\frac{{\partial }^2{H}_0}{\partial {\varphi }^2}\frac{{\partial }^2{H}_0}{\partial \varphi \partial p_{\rho }}\frac{{\partial }^2{H}_0}{\partial \varphi \partial p_{\varphi }}\\ \frac{{\partial }^2{H}_0}{\partial p_{\rho }\partial \rho }\frac{{\partial }^2{H}_0}{\partial p_{\rho }\partial \varphi }\frac{{\partial }^2{H}_0}{\partial p_{\rho }^2}\frac{{\partial }^2{H}_0}{\partial p_{\rho }\partial p_{\varphi }}\\ \frac{{\partial }^2{H}_0}{\partial p_{\varphi }\partial \rho }\frac{{\partial }^2{H}_0}{\partial p_{\varphi }\partial \varphi }\frac{{\partial }^2{H}_0}{\partial p_{\varphi }\partial p_{\rho }}\frac{{\partial }^2{H}_0}{\partial p_{\varphi }^2}\end{array}\right),\end{array}$$(7)

evaluated at the point L₁, via the relations λ_1, 3 = ±λ, λ_2, 4 = ±iκ. Furthermore. the columns of the matrix 𝒜 are derived from the unitary eigenvectors of ℳ (see Appendix A). Finally, the constant κ is equal to the epicyclic frequency at the distance ρ_L1, namely ${\kappa }^2={\partial }^2{V}_0/\partial {\rho }_{\text{L1}}^2+3p_{\varphi .\text{L1}}^2/{\rho }_{\text{L1}}^4$, assuming that V₀(ρ) in Eq. (1) represents the entire disk’s axisymmetric potential term, that is, ⟨V_bar⟩ = 0, where the average is taken with respect to all angles at fixed ρ.

The Hamiltonian H_0, 2 in Eq. (6) describes the linearized dynamics around L₁: the harmonic oscillator part (κ/2)(Q² + P²) describes epicyclic oscillations with the frequency κ, while the hyperbolic part λUV implies an exponential dependence of the variables U and V on time, namely U(t) = U(0)e^λt, V(t) = V(0)e^−λt. The linearized phase space can be decomposed in three subspaces, namely the invariant plane $E_{\mathrm{L}1}^{\mathrm{C}}=(Q,P)$, called the “linear center manifold”; as well as the axes U, called the linear unstable manifold $E_{\mathrm{L}1}^{\mathrm{U}}$; and V, called the linear stable manifold $E_{\mathrm{L}1}^{\mathrm{S}}$ of the point L₁. The linearized equations of motion yield independent motions in each of the spaces $E_{\mathrm{L}1}^{\mathrm{C}}$, $E_{\mathrm{L}1}^{\mathrm{U}}$, and $E_{\mathrm{L}1}^{\mathrm{S}}$. Those on $E_{\mathrm{L}1}^{\mathrm{U}}$ describe orbits that recede exponentially fast from L₁. A simple analysis shows that the outflow defined by such orbits has the form of trailing spiral arms. Basic theorems on invariant manifolds (Grobman 1959; Hartman 1960) predict that the invariant subspaces $E_{\mathrm{L}1}^{\mathrm{C}}$, $E_{\mathrm{L}1}^{\mathrm{U}}$, and $E_{\mathrm{L}1}^{\mathrm{S}}$ of the linearized model are continued as invariant sets $W_{\mathrm{L}1}^{\mathrm{C}}$, $W_{\mathrm{L}1}^{\mathrm{U}}$, and $W_{\mathrm{L}1}^{\mathrm{S}}$, respectively, in the full nonlinear model given by the Hamiltonian H₀. In particular, the linear unstable manifold $E_{\mathrm{L}1}^{\mathrm{U}}$ is tangent, at the origin, to the unstable manifold $W_{\mathrm{L}1}^{\mathrm{U}}$ of the full model. The latter is defined as the set of all initial conditions tending asymptotically to L₁ when integrated backward in time. In the forward sense of time, these orbits form an outflow, which deviates exponentially from L₁. This outflow forms trailing spiral arms or ring-like structures, which can be visualized either as “flux tubes” (Romero-Gomez et al. 2006) or as “apocentric manifolds” Voglis et al. (2006). For more details and precise definitions (see Efthymiopoulos 2010; Efthymiopoulos et al. 2019, and the references therein).

We now extend the previous notions from the Hamiltonian H₀ to the full model of Eq. (4). To this end, we first consider the canonical transformation (5), which is defined through the coefficients of the second order expansion around the coordinates of L₁ of the H₀ part of the Hamiltonian (i.e., the first line in Eq. (4)). In substituting this transformation to the full Hamiltonian, and by omitting constants, the Hamiltonian takes the form:

$$\begin{array}{cc}\hfill H=& \lambda UV+\frac{\kappa }{2}(Q^2+P^2)\hfill \\ \hfill & +{\displaystyle \sum _{s=3}^{\infty }\sum _{\begin{array}{c}{k}_1,l_1,k_2,l_2\ge 0\\ k_1+k_2+l_1+l_2=s\end{array}}{h}_{k_1,k_2,l_1,l_2}{U}^{k_1}{Q}^{k_2}{V}^{l_1}{P}^{l_2}}\hfill \\ \hfill & +({\mathrm{\Omega }}_2-{\mathrm{\Omega }}_{\mathrm{bar}})I_2+({\mathrm{\Omega }}_3-{\mathrm{\Omega }}_{\mathrm{bar}})I_3+\dots \hfill \\ \hfill & +{\displaystyle \sum _{s=1}^{\infty }\sum _{\begin{array}{c}{k}_1,l_1,k_2,l_2\ge 0\\ k_1+k_2+l_1+l_2=s\end{array}}{\mathrm{\Psi }}_{(2),k_1,k_2,l_1,l_2}({\varphi }_2)U^{k_1}{Q}^{k_2}{V}^{l_1}{P}^{l_2}}\hfill \\ \hfill & +{\displaystyle \sum _{s=1}^{\infty }\sum _{\begin{array}{c}{k}_1,l_1,k_2,l_2\ge 0\\ k_1+k_2+l_1+l_2=s\end{array}}{\mathrm{\Psi }}_{(3),k_1,k_2,l_1,l_2}({\varphi }_3)U^{k_1}{Q}^{k_2}{V}^{l_1}{P}^{l_2}+\dots }\hfill \end{array}$$(8)

The functions Ψ_{(2),k1, k2, l1, l2}(ϕ₂), Ψ_{(3),k1, k2, l1, l2}(ϕ₃), etc. are trigonometric in the angles ϕ₂, ϕ₃, etc., while h_{k1, k2, l1, l2} are constants. The Lagrangian point L₁ has coordinates U = V = Q = P = 0. This is no longer an equilibrium solution of the full Hamiltonian; with Hamilton’s equations, one obtains, in general, $\dot{U}\ne 0$, $\dot{V}\ne 0$, $\dot{Q}\ne 0$, Ṗ ≠ 0 for U = V = Q = P = 0, provided that at least one of the functions Ψ_{(2),k1, k2, l1, l2}(ϕ₂), Ψ_{(3),k1, k2, l1, l2}(ϕ₃), etc., are different from zero for k₁ + k₂ + l₁ + l₂ = 1. However, the existence of an equilibrium solution of the Hamiltonian (8) can be proven using perturbation theory. In particular, as a consequence of a theorem proven in Giorgilli (2001), there is a near-to-identity canonical transformation

$$\begin{array}{cc}& (U,Q,V,P,{\varphi }_2,{\varphi }_3,\dots ,I_2,I_3,\dots )\hfill \\ \hfill & \to (\xi ,q,\eta ,p,{\theta }_2={\varphi }_2,{\theta }_3={\varphi }_3,\dots ,J_2,J_3,\dots )\hfill \end{array}$$(9)

with

$$\begin{array}{cc}\hfill U& =\xi +F_U(\xi ,q,\eta ,p;{\varphi }_2,{\varphi }_3,\dots )\hfill \\ \hfill Q& =q+F_Q(\xi ,q,\eta ,p;{\varphi }_2,{\varphi }_3,\dots )\hfill \\ \hfill V& =\eta +F_V(\xi ,q,\eta ,p;{\varphi }_2,{\varphi }_3,\dots )\hfill \\ \hfill P& =p+F_P(\xi ,q,\eta ,p;{\varphi }_2,{\varphi }_3,\dots )\hfill \end{array}$$(10)

where the functions F_U, F_Q, F_V, and F_P are polynomial series of second or higher degree in the variables (ξ, q, η, p) and trigonometric in the angles ϕ₂, ϕ₃, etc., such that, in the variables (ξ, q, η, p, θ₂, θ₃, …, J₂, J₃, …) the Hamiltonian (8) takes the form

$$\begin{array}{cc}\hfill H=& \lambda \xi \eta +\frac{\kappa }{2}(q^2+p^2)\hfill \\ & +{\displaystyle \sum _{s=3}^{\infty }\sum _{\begin{array}{c}{k}_1,l_1,k_2,l_2\ge 0\\ k_1+k_2+l_1+l_2=s\end{array}}{g}_{k_1,k_2,l_1,l_2}{\xi }^{k_1}{q}^{k_2}{\eta }^{l_1}{p}^{l_2}}\hfill \\ & +({\mathrm{\Omega }}_2-{\mathrm{\Omega }}_{\mathrm{bar}})J_2+({\mathrm{\Omega }}_2-{\mathrm{\Omega }}_{\mathrm{bar}})J_3+\cdots \hfill \\ & +{\displaystyle \sum _{s=0,2,3,\dots }\sum _{\begin{array}{c}{k}_1,l_1,k_2,l_2\ge 0\\ k_1+k_2+l_1+l_2=s\end{array}}{\mathrm{\Phi }}_{(2),k_1,k_2,l_1,l_2}({\varphi }_2){\xi }^{k_1}{q}^{k_2}{\eta }^{l_1}{p}^{l_2}}\hfill \\ & +{\displaystyle \sum _{s=0,2,3,\dots }\sum _{\begin{array}{c}{k}_1,l_1,k_2,l_2\ge 0\\ k_1+k_2+l_1+l_2=s\end{array}}{\mathrm{\Phi }}_{(3),k_1,k_2,l_1,l_2}({\varphi }_3){\xi }^{k_1}{q}^{k_2}{\eta }^{l_1}{p}^{l_2}+\cdots }\hfill \end{array}$$(11)

The functions Φ_{(2),k1, k2, l1, l2}(ϕ₂), Φ_{(3),k1, k2, l1, l2}(ϕ₃), etc., are again trigonometric in the angles ϕ₂, ϕ₃, etc., while g_{k1, k2, l1, l2} are constants.

The formal difference between the Hamiltonians (8) and (11) is the lack, in the latter case, of polynomial terms that are linear in the variables (ξ, q, η, p). As a consequence, the point ξ = q = η = p = 0 is an equilbrium point of the system as transformed to the new variables, since from Hamilton’s equations for the Hamiltonian (11) one has $\dot{\xi }=\dot{q}=\dot{\eta }=\dot{p}=0$ if ξ = q = η = p = 0. Using the transformation (10) the equilibrium solution can be represented in the original variables as a “generalized L₁ solution”:

$$\begin{array}{cc}\hfill U_{\mathrm{GL}1}& =F_u(0,0,0,0;{\varphi }_2,{\varphi }_3,\dots )\hfill \\ \hfill Q_{\mathrm{GL}1}& =F_Q(0,0,0,0;{\varphi }_2,{\varphi }_3,\dots )\hfill \\ \hfill V_{\mathrm{GL}1}& =F_v(0,0,0,0;{\varphi }_2,{\varphi }_3,\dots )\hfill \\ \hfill P_{\mathrm{GL}1}& =F_P(0,0,0.0;{\varphi }_2,{\varphi }_3,\dots ).\hfill \end{array}$$(12)

Through the linear transformation (5), we also obtain the functions δρ_GL1(ϕ₂, ϕ₃, …), δϕ_GL1(ϕ₂, ϕ₃…), p_ρ, GL1(ϕ₂, ϕ₃…), and J_ϕ, GL1(ϕ₂, ϕ₃…). Since ϕ₂ = (Ω₂ − Ω_bar)t, ϕ₃ = (Ω₃ − Ω_bar)t, etc., the above functions determine the time-dependence of the original phase-space coordinates of the generalized trajectory GL₁. The trajectory depends trigonometrically on the phases ϕ₂, ϕ₃, etc., hence it depends on time through the frequencies |Ω₂ − Ω_bar|, |Ω₃ − Ω_bar|, etc. In particular, the trajectory GL₁ is a periodic orbit (“a 1-torus”) when there is one extra pattern speed. This generalizes to a Lissajous-like figure (M-torus) when there are M >  1 extra pattern speeds, etc. Through the linear part of Hamilton’s equations for the Hamiltonian (11), we find that the equilibrium point (ξ, q, η, p) = (0, 0, 0, 0) is simply unstable, that is, the variational matrix has one pair of real eigenvalues equal to ±λ and one pair of imaginary eigenvalues equal to ±iκ). By also taking into account the frequencies |Ω₂ − Ω_bar|, |Ω₃ − Ω_bar|, etc., the complete phase space in the neighborhood of the solution GL₁ can be decomposed into a center^M + 1 × saddle topology (Gómez et al. 2001).

First, this means, in particular, that the phase-space invariant subset ${\mathcal{W}}_{\mathrm{GL}1}^{\mathrm{C}}$ defined by the condition ξ = η = 0 is invariant under the flow of the Hamiltonian (11). It is hereafter called the center manifold of the orbit GL₁. Its dimension is 2 + 2M, where M is the number of additional frequencies. By the structure of Hamilton’s equations, ${\mathcal{W}}_{\mathrm{GL}1}^{\mathrm{C}}$ is a normally hyperbolic invariant manifold (NHIM; see Wiggins 1994).

Second, the set ${\mathcal{W}}_{\mathrm{GL}1}^{\mathrm{U}}$ of all initial conditions tending asymptotically to the generalized orbit GL₁ in the backward sense of time is the unstable manifold of the orbit GL₁. Basic theorems of dynamics (Grobman 1959; Hartman 1960) guarantee that such an invariant manifold exists, and at the origin it is tangent to the linear unstable manifold ${\mathcal{E}}_{\mathrm{GL}1}^{\mathrm{U}}$, which coincides with the axis ξ with q = η = p = 0. Both ${\mathcal{E}}_{\mathrm{GL}1}^{\mathrm{U}}$ and ${\mathcal{W}}_{\mathrm{GL}1}^{\mathrm{U}}$ are one-dimensional. The product of ${\mathcal{W}}_{\mathrm{GL}1}^{\mathrm{U}}$ with the angles ϕ₂ = (Ω₂ − Ω_bar)t, ϕ₃ = (Ω₃ − Ω_bar)t, etc., defines the “generalized unstable tube manifold” of the orbit GL₁, denoted hereafter as ${\mathcal{W}}_{\mathrm{GL}1}^{\mathrm{TU}}$.

Third, similar definitions remain valid for the stable manifold ${\mathcal{W}}_{\mathrm{GL}1}^{\mathrm{S}}$ and stable tube manifold ${\mathcal{W}}_{\mathrm{GL}1}^{\mathrm{TS}}$ of the orbit GL₁. These represent the sets of orbits tending asymptotically to the orbit GL₁ in the forward sense of time.

As in the standard manifold theory of spirals, the basic objects giving rise to spirals are the generalized unstable tube manifolds ${\mathcal{W}}_{\mathrm{GL}1}^{\mathrm{TU}}$ and ${\mathcal{W}}_{\mathrm{GL}2}^{\mathrm{TU}}$ of the orbits GL₁ and GL₂, respectively. A basic argument allows one to show that the projections of ${\mathcal{W}}_{\mathrm{GL}1}^{\mathrm{TU}}$ and ${\mathcal{W}}_{\mathrm{GL}2}^{\mathrm{TU}}$ on the configuration space are trailing spirals that emanate from the neighborhood of the bar’s Lagrangian points L₁ and L₂, but with a position and shape that vary in time quasi-periodically. The variation is small and characterized by as many frequencies as the additional pattern speeds. The argument is as follows: instead of the transformation (10), one can formally compute a standard Birkhoff transformation (see Efthymiopoulos 2012) of the form

$$\begin{array}{cc}\hfill U& ={\xi }_B+B_u({\xi }_B,q_B,{\eta }_B,p_B;{\varphi }_2,{\varphi }_3,\dots )\hfill \\ \hfill Q& =q_B+B_Q({\xi }_B,q_B,{\eta }_B,p_B;{\varphi }_2,{\varphi }_3,\dots )\hfill \\ \hfill V& ={\eta }_B+B_v({\xi }_B,q_B,{\eta }_B,p_B;{\varphi }_2,{\varphi }_3,\dots )\hfill \\ \hfill P& =p_B+B_P({\xi }_B,q_B,{\eta }_B,p_B;{\varphi }_2,{\varphi }_3,\dots )\hfill \end{array}$$(13)

such that the Hamiltonian (8) expressed in the new variables (ξ_B, q_B, η_B, p_B) becomes independent of the angles ϕ₂, ϕ₃, etc., and it takes the form (apart from a constant)

$$\begin{array}{cc}\hfill H=& \lambda \xi {\eta }_B+\frac{\kappa }{2}(q_B^2+p_B^2)\hfill \\ & +{\displaystyle \sum _{s=1}^{\infty }\sum _{\begin{array}{c}{k}_1,l_1,k_2,l_2\ge 0\\ k_1+k_2+l_1+l_2=s\end{array}}{f}_{k_1,k_2,l_1,l_2}{\xi }_B^{k_1}{q}_B^{k_2}{\eta }_B^{l_1}{p}_B^{l_2}.}\hfill \end{array}$$(14)

Contrary to the normalization leading to the Hamiltonian (11), the Birkhoff normalization that leads to the Hamiltonian (14) is not guaranteed to converge (see Efthymiopoulos 2012), thus it cannot be used to theoretically demonstrate the existence of the manifolds ${\mathcal{W}}_{\mathrm{GL}1}^{\mathrm{TU}}$ and ${\mathcal{W}}_{\mathrm{GL}2}^{\mathrm{TU}}$. For practical purposes, however, the Birkhoff normalization can proceed up to an exponentially small remainder, hence the Hamiltonian (14) approximates the dynamics with an exponentially small error. In this approximation, the coefficients f_{k1, k2, l1, l2} are of the order of the amplitude of the extra patterns if k₁ + k₂ + l₁ + l₂ ≤ 2, while one has f_{k1, k2, l1, l2} = h_{k1, k2, l1, l2} + h.o.t if k₁ + k₂ + l₁ + l₂ >  2. Hence, the resulting Hamiltonian is dominated by the bar terms. The equilibrium solutions representing the generalized Lagrangian equilibria GL₁ and GL₂ can be computed as the (nonzero) roots of Hamilton’s equations ${\dot{\xi }}_B={\dot{q}}_B={\dot{\eta }}_B={\dot{p}}_B=0$. The key remark is that since Hamiltonian (14) no longer depends on time, the unstable tube manifolds ${\mathcal{W}}_{B,\mathrm{GL}1}^{\mathrm{TU}}$ and ${\mathcal{W}}_{\mathrm{B},\mathrm{GL}2}^{\mathrm{TU}}$ remain unaltered in time when regarded in the variables (ξ_B, q_B, η_B, p_B). Then, due to the transformation (13), the manifolds as expressed in the original variables have a dependence on the angles ϕ₂, ϕ₃, etc., implying a dependence on time through M independent frequencies. Physically, the manifolds are subject to small oscillations (of order maxV_i, i = 2, 3, …) with respect to a basic static shape, which is given by their time-invariant form in the variables (ξ_B, q_B, η_B, p_B). Hence, the manifolds yield spirals with a pattern exhibiting quasi-periodic oscillations around the basic spiral patterns induced by the manifolds of the pure bar model.

3. Application in a Milky-Way-type model

We now apply the above mentioned theory in the case of a Milky-Way-type galactic model, assuming a different pattern speed for the bar and for the spiral arms. We emphasize that this is not intended as a modeling of the real spiral structure in the Milky Way, but only as a “proof of concept” of the possibility of manifold spirals to support structures with more than one pattern speed.

3.1. Potential

In our model we use a variant of the Galactic potential proposed in Pettitt et al. (2014). This consists of the components listed below.

The axisymmetric component is a superposition of a disk and halo components, V_ax(ρ, z) = V_d(ρ, z)+V_h(r), where r = (x² + y² + z²)^1/2. The disk potential has the Miyamoto–Nagai form (Miyamoto & Nagai 1975)

$$\begin{array}{c}\hfill V_{\mathrm{d}}=\frac{-G{M}_{\mathrm{d}}}{\sqrt{{\rho }^2+{(a_{\mathrm{d}}+\sqrt{z^2+b_{\mathrm{d}}^2})}^2}},\end{array}$$(15)

where M_d = 8.56 × 10¹⁰ M_⊙, a_d = 5.3 kpc and b_d = 0.25 kpc. The halo potential is a γ-model (Dehnen 1993) with parameters as in Pettitt et al. (2014)

$$\begin{array}{c}\hfill V_{\mathrm{h}}=\frac{-G{M}_{\mathrm{h}}(r)}{r}-\frac{-G{M}_{\mathrm{h},0}}{\gamma r_{\mathrm{h}}}{[-\frac{\gamma }{1+{(r/r_{\mathrm{h}})}^{\gamma }}+ln{(1+\frac{r}{r_{\mathrm{h}}})}^{\gamma }]}_{\mathrm{r}}^{r_{\mathrm{h},max}},\end{array}$$(16)

where r_h, max = 100 kpc, γ = 1.02, M_h, 0 = 10.7 × 10¹⁰ M_⊙, and M_h(r) is the function

$$\begin{array}{c}\hfill M_{\mathrm{h}}(r)=\frac{M_{\mathrm{h},0}{(r/r_{\mathrm{h}})}^{\gamma +1}}{1+{(r/r_{\mathrm{h}})}^{\gamma }}·\end{array}$$(17)

The bar potential is as in Long & Murali (1992), that is

$$\begin{array}{c}\hfill V_{\mathrm{b}}=\frac{G{M}_{\mathrm{b}}}{2a}ln\left(\frac{x-a+T_{-}}{x+a+T_{+}}\right)\end{array}$$(18)

with $T_{\pm }=\sqrt{[{(a\pm x)}^2+{\mathit{y}}^2+{(b+\sqrt{c^2+z^2})}^2]}$, M_b = 6.25 × 10¹⁰ M_⊙, a = 5.25 kpc, b = 2.1 kpc, and c = 1.6 kpc. The values of a and b set the bar’s scale along the major and minor axes in the disk plane (x and y, respectively), while c sets the bar’s thickness in the z-axis (see Gerhard 2002; Rattenbury et al. 2007; Cao et al. 2013). These values where chosen so as to bring the bar’s corotation, for Ω_bar = 45 km s⁻¹ kpc⁻¹, to the value as specified by the L_1, 2 points’ distance from the center, R_L1, 2 = 5.4 kpc. Assuming corotation to be at 1.2–1.3 times the bar’s length, the latter turns out to be about 4 kpc with the adopted parameters.

Regarding the spiral arms, we use a variant of the logarithmic spiral arms model adopted in Pettitt et al. (2014). The spiral potential reads (Cox & Gómez 2002)

$$\begin{array}{cc}\hfill V_{\mathrm{sp}}=& -4\pi \frac{G{h}_z{d}_{0}C}{\mathit{KD}}{\left[\mathrm{sech}\left(\frac{\mathit{Kz}}{\beta }\right)\right]}^{\beta }F(\rho )exp(-\frac{\rho -{\rho }_0}{R_{\mathrm{s}}})\hfill \\ \hfill & \times cos\left(N[\varphi -({\mathrm{\Omega }}_{\mathrm{sp}}-{\mathrm{\Omega }}_{\mathrm{bar}})t-\frac{ln(\rho /{\rho }_0)}{tan(\alpha )}]\right),\hfill \end{array}$$(19)

where N is the number of spiral arms and

$$\begin{array}{cc}& K=\frac{N}{\rho |sin(\alpha )|},\hfill \end{array}$$(20)

$$\begin{array}{cc}& D=\frac{1+K{h}_z+0.3{(K{h}_z)}^2}{1+0.3K{h}_z},\mathrm{and}\hfill \end{array}$$(21)

$$\begin{array}{cc}& \beta =K{h}_z(1+0.4K{h}_z).\hfill \end{array}$$(22)

The function F(ρ) plays the role of a smooth envelope that determines the radius beyond which the spiral arms are important. We adopt the form F(ρ) = b − carctan((R_s0 − ρ)/kpc), with R_s0 = 6 kpc, b = 0.474, and c = 0.335. The values of the remaining constants are: N = 2, α = −13°, h_z = 0.18 kpc, R_s = 3 kpc, ρ₀ = 8 kpc, and C = 8/3π. The spiral amplitude is determined by setting the value of the density d₀. We consider three values, namely d₀ = A_s0 × 10⁸ M_⊙ kpc⁻³, with A_s0 = 1.5, 3, and 4, called the weak, intermediate, and strong spirals, respectively. These values were chosen so as to yield spiral Q-strength values that are consistent with those reported in the literature for a mild bar (see Buta et al. 2009). Our basic model is the intermediate one, but as shown below, there are only small variations to the basic manifold morphology in any of these three choices since the manifolds’ shape is mostly determined by the bar. Finally, for the spiral pattern speed we adopt the value Ω_spiral  =  20 km s⁻¹ kpc⁻¹, which is different from the bar pattern speed Ω_bar  =  45 km s⁻¹ kpc⁻¹ (Gerhard 2011; Bland-Hawthorn & Gerhard 2016).

Figure 1 shows the rotation curve arising from the axisymmetric components as well as the azimuthally averaged part of the bar’s potential; the corresponding component is equal to zero for the spirals. The model is close to the “maximum disk”, that is, the rotation curve up to ∼10 kpc is essentially produced by the components of the disk and bar alone. On the other hand, Fig. 2 shows an isodensity color map of the projected surface density $\sigma (x,\mathit{y})={\int }_{-\infty }^{\infty }\rho (x,\mathit{y},z)\mathrm{d}z$ in the disk plane, where the density ρ is computed from Poisson’s equation ∇²V = 4πGρ for the potential V = V_d + V_bar + V_sp. The fact that the spiral potential has a nonzero relative pattern speed in the bar’s frame results in a time-dependent spiral pattern in the disk plane. However, it is well known (Sellwood & Sparke 1988) that, under reasonable assumptions for the bar and spiral parameters, such a time dependence results in a morphological continuity, at most time snapshots, between the end of the bar and the spiral arms. In order to numerically test the manifold theory, we choose below four snapshots to be characteristic, corresponding to the times t = 0, T/4, T/2, and 3T/4 in Eq. (19), where T = π/|Ω_sp − Ω_bar|. It is important to note that since the imposed spiral potential has only cos2ϕ and sin2ϕ terms, the spiral patterns shown in Fig. 2 are repeated periodically with period T. Defining the “phase” of the spirals at a radial distance ρ as the angle ϕ_s(ρ) where the spiral potential is minimum, given by

$$\begin{array}{c}\hfill {\varphi }_{\mathrm{s}}(\rho ,t)=({\mathrm{\Omega }}_{\mathrm{sp}}-{\mathrm{\Omega }}_{\mathrm{bar}})t+(1/tan(a))ln(\rho /{\rho }_0)(\mathrm{mod}2\pi ),\end{array}$$(23)

Fig. 1. Rotation curve (black) corresponding to the potential V = Vh + Vd + ⟨Vb⟩, where ⟨Vb⟩ is the m = 0 (average with respect to all azimuths) part of the bar’s potential Vb. The contribution of each component is shown with a different color.

Fig. 2. Color map of the surface density Σ(ρ, ϕ) corresponding to the potential V = Vd + Vb + Vsp (see text), as viewed in the bar’s rotating frame, at four different snapshots, namely t = 0 (top left), t = T/4 (top right), t = T/2 (bottom left), and t = 3T/4 (bottom right), where T = π/|Ωsp − Ωbar|. Since Ωsp <  Ωbar, the spirals have a relative clockwise angular displacement in time with respect to the bar. However, the morphological continuity between the bar and spirals is retained in all of these snapshots.

we characterize below the relative position of the spirals with respect to the bar by the angle ϕ_s(ρ₀, t), which is a periodic function of time. In physical terms, the angle ϕ_s measures the angular distance between the point L₁ (L₂), which lies in the semi-plane x >  0 (x <  0), and the point of local minimum with respect to ϕ of the spiral potential at ρ = ρ₀, which lies in the semi-plane y ≤ 0 (y ≥ 0).

Regarding the relative bar and spiral contributions to the nonaxisymmetric forces, Fig. 3 allows one to estimate the relative importance of the bar’s and spirals’ nonaxisymmetric force perturbation by showing the corresponding Q-strengths as functions of the radial distance ρ in the disk. The Q-strength at fixed ρ (e.g., Buta et al. 2009) is defined for the bar as

$$\begin{array}{c}\hfill Q_{\mathrm{b}}(r)=\frac{F_{\mathrm{b},\mathrm{t}}^{max}(\rho )}{⟨F_{\mathrm{r}}(\rho )⟩},\end{array}$$(24)

Fig. 3. Bar, spiral, and total Q-strengths (Qb, Qs and Qtotal respectively) as functions of the radius ρ in the model including the potential terms Vd, Vb, and Vsp for the intermediate spiral model (see text).

where $F_{\mathrm{b},\mathrm{t}}^{max}$ is the maximum, with respect to all azimuths ϕ, tangential force generated by the potential term V_b at the distance ρ, while ⟨F_r(ρ)⟩ is the average, with respect to ϕ, radial force at the same distance generated by the potential V_d + V_b + V_sp. The bar yields a Q-value Q_b ≈ 0.25 in its inner part, which falls to Q_b ≈ 0.15 to 0.10 in the domain outside the bar where the manifolds (and spirals) develop, that is, 5 kpc < ρ <  10 kpc. The spirals, in turn, yield a maximum Q_s around ρ ≈ 7 kpc, which is equal to Q_s ≈ 0.08 in the intermediate model, turning to 0.04 or 0.11 in the weak and strong models, respectively. Thus, the total Q-strength is about 0.15–0.2 in the domain of interest.

3.2. Manifold spirals

A useful preliminary computation pertains to the form of the apocentric manifolds in the above models in two particular cases: (i) a pure bar case, and (ii) a bar and spiral case, assuming, however, that the spirals rotate with the same pattern speed as the bar. The corresponding results are shown in Fig. 4. It is noteworthy that even the pure bar model yields manifolds that support a spiral response (Figs. 4a and c). In addition, the manifolds induce a R₁-type ring-like structure, which is reminiscent of pseudorings (see Buta 2013 for a review), that is, rings with a diameter that is comparable to the bar’s length and a spiral-like deformation with respect to a symmetric shape on each side of the bar’s minor axis. Now, by adding the spiral term, with the same pattern speed as the bar, these structures are considerably enhanced (Figs. 4b and d). The most important effect is on the pseudo-ring structure, which is now deformed in order to support the imposed spirals over a large extent. It is of interest to follow in detail how the intricate oscillations of the manifolds result in supporting the imposed spiral structure. Figure 4b gives the corresponding details. We note that the manifolds emanating from the point L₁ (blue) initially expand outward, yielding spirals with a nearly constant pitch angle. However, after half a turn, the manifolds turn inward and move toward the neighborhood of the point L₂. While approaching this point, the manifolds develop oscillations, known in dynamics as the “homoclinic oscillations” (see Contopoulos 2002 for a review). As a result, the manifolds form thin lobes. In Figs. 4a and b, we mark the tips of the first four lobes with the numbers 1 to 4, and label these lobes accordingly. Focusing on Fig. 4b, we note that lobe 1 is in the transient domain between the spirals and the Lagrangian points. However, lobe 2 of the manifold emanating from L₁ supports the spiral arm originating from the end of the bar at L₂, and, conversely, lobe 2′ of the manifold emanating from L₂ supports the spiral arm originating from L₁. We call this phenomenon a bridge (see also Efthymiopoulos et al. 2019) and mark the corresponding parts of the manifolds with B and B′. One can check that this phenomenon is repeated for higher order lobes of the manifolds. Thus, in Fig. 4b, lobe 3 supports the outer part of the pseuroding, which is assosiated with the spiral originating from L₂, while, in the same way, lobe 3′ supports the spiral originating from L₁. Furthermore, between lobes 2 and 3, a gap is formed (marked G), which separates the pseudoring from the outer spiral (and similarly for the gap G′ formed between lobes 2′ and 3′). On the other hand, lobe 4 returns to support the spiral originating from L₁. Higher order lobes repeat the same phenomenon, but their succession becomes more and more difficult to follow, as shown in Fig. 4d. One can remark that the manifolds mostly support the spiral geometry in the outer parts of the pseudorings. In fact, in the pure bar model, we again have the appearance of manifold oscillations, which lead to lobes, a bridge, and gaps (Figs. 4a and c), but now the ring part is only mildly deformed and clearly separated from the outer lobes that support the spirals.

Fig. 4. Panel a: apocentric invariant manifolds in the pure bar model with one pattern speed. The manifolds emanating from the points L1 and L2 are plotted in blue and red, respectively. Panel b: apocentric manifolds if we add the spiral potential (intermediate case), assuming that the spirals rotate with the same pattern speed as the bar. The addition of the spiral term enhances the structures described as “lobes”, “bridges”, and “gaps” (see text). Panels c and d: same as in (a) and (b), but with the manifolds computed over a larger length. The black spiral curves correspond to the minima of the imposed spiral potential, given by Eq. (19).

We now examine how these morphologies are altered, if, instead, we assume that the spirals rotate with a different pattern speed than the bar. The computation of the manifolds in this case can be carried out with the same steps as described in Sect. 2. For the computation of the initial diagonalizing transformation matrix 𝒜 (Eq. (5)) as well as the canonical transformation (10), we proceed as described in the Appendix A. In particular, we used the Lie series method in order to perform all series computations. These series allowed us to compute initial conditions for the periodic orbits GL₁ and GL₂ (Eq. (12)). Finally, we numerically refined the latter computation using Newton–Raphson to obtain the periodic orbits with many significant figures. More specifically, since the potential depends periodically on time (with period T = π/(Ω_bar − Ω_sp)), we consider a stroboscopic map

$$\begin{array}{c}\hfill (x(0),y(0),p_{\rho }(0),p_{\varphi }(0))\to (x(T),y(T),p_{\rho }(T),p_{\varphi }(T)),\end{array}$$(25)

which maps any initial condition at the time t = 0 to its image at the time t = T under the full numerical equations of motion without any approximation. Then, the periodic orbits GL₁ and GL₂ are fixed points of the above map. As shown in Fig. 5, the periodic orbits GL₁ and GL₂ found by the above mentioned method form epicycles around the Lagrangian points L₁ and L₂ of the pure bar model. However, the orbits GL₁ and GL₂ should not be confused with the epicyclic Lyapunov orbits PL₁ and PL₂ used in past manifold calculations in models with one pattern speed (Voglis et al. 2006). In particular, the orbits PL₁ and PL₂ exist as a family of orbits in a fixed bar model, whose size depends continuously on the value of the Jacobi energy E_J >  E_L1. Under specific conditions, the orbits PL_1, 2 can be generalized to 2D-tori in the case of one extra pattern speed. However, this generalization requires the use of Kolmogorov-Arnold-Moser theory (Kolmogorov 1954; Arnold 1963; Moser 1962), which is beyond the scope of this paper. On the contrary, in the two-pattern speed case, for a fixed choice of the potential V₂ = V_sp (Eq. (19)) and Ω₂ = Ω_sp, there exist unique GL₁ and GL₂ orbits, which generalize the unique Lagrangian points of the corresponding pure bar model. In fact, the orbits of Fig. 5 have a relative size of the order of the ratio of the m = 2 Fourier amplitudes of the bar and of the spiral potential at the radius ρ = ρ_L1, 2. This is about 0.5 kpc, 1.2 kpc, and 1.5 kpc in the weak, intermediate, and strong spiral case, respectively.

Fig. 5. Periodic orbit GL1 in the weak (left), intermediate (center), and strong (right) spiral cases. The size of the orbit increases with the spiral amplitude.

The computation of the unstable manifolds of the orbits GL₁ and GL₂ is now straightforward. By focusing on GL₁, for example, we first computed the 4 × 4 variational matrix Λ of the mapping (25) evaluated at the fixed point of the periodic orbit GL₁. The matrix Λ satisfies the symplecticity condition Λ ⋅ 𝒥 ⋅ Λ^T = Λ^T ⋅ 𝒥 ⋅ Λ = 𝒥, and it has two real reciprocal eigenvalues λ₁, λ₂ = 1/λ₁, with |λ₁|> 1, and two complex congugate ones with unitary measure λ_3, 4 = e^±iωT for some positive ω. By denoting $e_{\mathrm{GL}1}^{\mathrm{U}}$ as the unitary eigenvector of Λ, associated with the eigenvalue λ₁, we then consider a small segment divided in 10⁵ initial conditions of the form (ρ_i, 0, ϕ_i, 0, p_ρ, i, 0, p_ϕ, i, 0, i = 1, …, 100 000 defined by (ρ_i, 0, ϕ_i, 0, p_ρ, i, 0, p_ϕ, i, 0) = (ρ_GL1 + δρ_i, 0, ϕ_GL1 + δϕ_i, 0, p_ρ, GL1 + δp_ρ, i, 0, p_ϕ, GL1 + δp_ϕ, i, 0) where $(\delta {\rho }_{i,0},\delta {\varphi }_{i,0},\delta p_{\rho ,i,0},\delta p_{\varphi ,i,0})=(i/100000)\times \mathrm{\Delta }S\times e_{\mathrm{GL}1}^{\mathrm{U}}$, with ΔS = 0.001. Propagating all of these orbits forward in time yields an approximation of the unstable flux-tube manifold ${\mathcal{W}}_{\mathrm{GL}1}^{\mathrm{TU}}$ (see Sect. 2).

In contrast to what happens in the one-pattern speed model, under the presence of the second pattern speed, the projection of the flux-tube manifolds ${\mathcal{W}}_{\mathrm{GL}1}^{\mathrm{TU}}$ in the disk plane varies in time. In order to efficiently visualize how the manifolds develop in space and time, in the following plots, we use an “apocentric double section” of the manifolds, denoted ${\mathcal{W}}_{\mathrm{GL}1}^{\mathrm{AU}}(t_{\mathrm{s}})$, which depends on a chosen value of the “section time” t_s. The computation of the apocentric double section for a given time t_s includes: keeping track of all the points of the tube manifolds generated by the above initial conditions; retaining those points that correspond to integration times t = nT + t_s ± ΔT, with n = 0, 1, 2, etc., and ΔT small (ΔT = 0.1T in all our calculations); and, at the same time, satisfying the apocentric condition ṗ_ρ ≃ 0 with an accuracy defined by $|\dot{\rho }|<|\ddot{\rho }|\mathrm{d}t$ where $|\ddot{\rho }|$ is the measure of the radial acceleration at the evaluation point, and dt = 0.001T is the integration timestep. This representation allows one to obtain the intersections of the manifolds with an apocentric surface of a section (see Efthymiopoulos 2010 for a discusion of how the apocentric manifolds compare with the full flux-tube manifolds). However, it also allows one to capture the dependence of the form of the manifolds on time, through the chosen value of t_s.

Figure 6 shows the main result: the manifolds ${\mathcal{W}}_{\mathrm{GL}1}^{\mathrm{AU}}(t_{\mathrm{s}})$ (blue points) and ${\mathcal{W}}_{\mathrm{GL}2}^{\mathrm{AU}}(t_{\mathrm{s}})$ (red points) computed as above, are shown at four different times t_s, namely t_s = 0, T/4, T/2 and 3T/4, corresponding to the same snapshots as in Fig. 2. The spiral phase ϕ_s(ρ₀, t_s) has the values 0, −π/4, −π/2, and −3π/4, respectively. The black-dotted curves that are superposed to the manifolds correspond to the maxima of the surface density in the annulus 6 kpc <  ρ <  15 kpc, as found from the data of Fig. 2. These figures periodically repeat after the time t_s = T.

Fig. 6. “Double section” apocentric manifolds ${\mathcal{W}}_{\mathrm{GL}1}^{\mathrm{AU}}(t_{\mathrm{s}})$ (blue) and ${\mathcal{W}}_{\mathrm{GL}2}^{\mathrm{AU}}(t_{\mathrm{s}})$ (red) are plotted at four different times ts as indicated in each panel. The black-dotted curves mark the local maxima of the surface density σ(ρ, ϕ, t), corresponding to the potential V = Vd(ρ)+Vb(ρ, ϕ)+Vsp(ρ, ϕ, t) at the times t = ts. These maxima are plotted in the domain 6 kpc < ρ <  15 Kpc, where the imposed spirals have a significant amplitude.

The key result from Fig. 6 is now evident: The spiral maxima rotate clockwise with respect to the bar (with angular velocity equal to 2π/(Ω_bar − Ω_sp). The manifolds ${\mathcal{W}}_{\mathrm{GL}1,2}^{\mathrm{AU}}$ adapt their form to the rotation of the spiral maxima, thus acquiring a time-varying morphology. In particular, the manifolds always form bridges and gaps, thus supporting a pseudo-ring as well as an outer spiral pattern. The spiral-like deformation of the pseudo-ring is most conspicuous at t_s = 0, corresponding to a spiral phase ϕ_s(ρ₀) = 0, and it remains large at the times t_s = T/4 and 3T/4, that is, at the spiral phases ϕ_s(ρ₀) = − π/4 and −3π/4. At all of these phases, the spiral maxima at ρ = ρ₀ remain close to the bar’s major axis, thus the manifolds tend to take a form that is similar to the one of Fig. 4d (in which ϕ_s(ρ₀) = 0 always since we set Ω_bar = Ω_sp). On the other hand, at t_s = T/2, (ϕ_s(ρ₀) = − π/2), the spiral maxima at ρ = ρ₀ are displaced by an angle π/2 with respect to the bar’s horizontal axis. Then, the manifolds yield more closed pseudo-rings, and they temporarily stop supporting the imposed spirals. In comparing the three phases ϕ_s(ρ₀) = 0, −π/4, and −3π/4, we find that the agreement between the manifolds and imposed spirals is best at the phases ϕ_s(ρ₀) = 0 and −3π/4, while the manifolds mostly support the imposed spiral in their pseudo-ring part at ϕ_s(ρ₀) = − π/4.

Altering the spirals’ amplitude (Figs. 7 and 8) makes no appreciable difference to the above mentioned picture. The main noticed difference regards the thickness of the manifolds’ lobes, which increases with the imposed spiral amplitude since, in general, the manifolds make larger oscillations near the bridges when the nonaxisymmetric perturbation increases. This also means that the trajectories supporting these spirals are more chaotic.

Fig. 7. Same as in Fig. 6, but for the weak spiral model.

Fig. 8. Same as in Fig. 6, but for the strong spiral model.

3.3. Discussion

As a comment on the loss of support of the manifolds to the imposed spiral maxima near ϕ_s(ρ₀) = − π/2, we remark that under the scenario in which the manifolds provide the backbone that supports chaotic spirals, a temporary loss of support implies that the spiral response to the manifolds should have its minimum strength when the spirals have a relative phase ±π/2 with respect to the bar’s major axis. Since the bar-spiral relative configuration (and the manifolds’ shape) is repeated periodically, with period T = π/(Ω_bar − Ω_sp), we conclude that under the manifold scenario, the amplitude in the response spiral should exhibit periodic time variations, with a period equal to T, that is, the manifolds support “recurrent spirals” with the above periodicity. The appearance of recurrent spirals in multipattern speed N-body models is well known (see Sellwood & Wilkinson 1993, Sellwood 2003). The manifold theory provides a specific prediction about the period of the recurrence, which is testable in such experiments by the time-Fourier analysis of the nonaxisymmetric patterns. On the other hand, the picture presented above is still “static”, in the sense that it does not take into account phenomena that alter the imposed nonaxisymmetric modes in time. Such phenomena are nonlinear interactions between distinct modes, and the enhancement or decay of the spirals, which is associated with disk instabilities (e.g., swing amplification) or with dissipation mechanisms (e.g., disk heating at resonances or gas phenomena). In all such circumstances, the manifolds provide a way to understand the behavior of chaotic trajectories beyond the bar. Thus, a full exploration of the connection between manifolds and collective disk phenomena is proposed for further study.

4. Conclusions

In the present study, we examine the possibility that manifold spirals in barred galaxies are consistent with the presence of multiple pattern speeds in the galactic disk. In Sect. 2, we detail the main theory and in Sect. 3 we provide numerical examples of such manifold spirals. Our main conclusions are as follows.

1. In the case of one pattern speed, the basic manifolds are those generated by the unstable manifolds of the Lagrangian points L₁ and L₂. In the case of multiple pattern speeds, it can be established theoretically (see Sect. 2) that, while Lagrangian equilibrium points no longer exist in the bar’s rotating frame, such points are replaced by “generalized Lagrangian orbits” (the orbits GL₁ and GL₂), which play a similar role in dynamics. These orbits are periodic, with a period equal to π/|Ω_sp − Ω_bar|, if there is one spiral pattern rotating with speed Ω_sp different from Ω_bar. If there is more than one extra pattern, with speeds Ω₂, Ω₃, etc., the generalized orbits GL₁ and GL₂ perform epicycles around the Lagrangian points L₁ and L₂ of the pure bar model with, in general, incommensurable frequencies |Ω_i − Ω_bar|, i = 1, 2, etc. Furthermore, in all cases the orbits GL₁ and GL₂ are simply unstable; a fact implying that they possess unstable manifolds 𝒲_GL1 and 𝒲_GL2. When the extra patterns have a small amplitude with respect to the bar’s amplitude, perturbation theory establishes that the manifolds 𝒲_GL1 and 𝒲_GL2 undergo small time variations (with the same frequencies |Ω_i − Ω_bar|, i = 1, 2, etc.), but, generally, their form only exhibits a small deformation with respect to the manifolds 𝒲_L1 and 𝒲_L2 of the pure bar model. Thus, the manifolds 𝒲_GL1 and 𝒲_GL2 support trailing spiral patterns.

2. In Sect. 3 we explore a simple bar-spiral model for a galactic disk with parameters relevant to Milky-Way dynamics. In this model we construct manifold spirals in both cases of a unique pattern speed (Ω_sp = Ω_bar) or two distinct pattern speeds (Ω_sp <  Ω_bar). The pure bar model already generates manifolds that support spiral patterns as well as an inner ring around the bar. Imposing further spiral perturbations on the potential mostly generates a deformation of the manifolds, with the ring evolving to a spiral-like “pseudo-ring”. The spiral and ring structures generated by the manifolds connect to each other through bridges (see Figs. 4–8). This implies that, after a bridge, the manifold emanating from the neighborhood of the bar’s L₁ (L₂) point supports the spiral arm associated with the bar’s end near the L₂ (L₁) point. From the point of view of dynamics, these connections are a manifestation of homoclinic chaos.

3. We find that the manifold theory gives good fit to at least some part of to the imposed spirals in both the single and multiple pattern speed models. Focusing on numerical examples in which the spiral and bar pattern speeds satisfy Ω_sp <  Ω_bar, the main behavior of the manifold spirals can be characterized in terms of the (time-varying) phase ϕ_s(ρ₀) (Eq. (23)). The manifolds support the imposed spirals over all the latter’s length at phases ϕ_s(ρ₀) = 0 or −3π/4, and they mostly support pseudo-ring like spirals near the phase −π/4 (the phase ϕ_s is negative since the spirals have a retrograde relative rotation with respect to the bar). On the other hand, the manifolds deviate from the imposed spirals near the phase −π/2. Both the manifolds’ shape and the imposed bar-spiral relative configuration are repeated periodically, with period T = π/(Ω_bar − Ω_sp). Thus, we argue that the temporary loss of support of the manifolds to the imposed spirals suggests a natural period for recurrent spirals, equal to T.

In summary, our analysis shows that manifold spirals in galactic disks are, in general, consistent, with the presence of multiple pattern speeds. Nevertheless, the manifold spirals in this case oscillate in time, thus, they support the imposed spirals along a varying length, which fluctuates from small to almost complete, depending on the relative phase of the spirals with respect to the bar. The manifolds also produce ring and pseudo-ring structures, which are morphologically connected to the spirals via the phenomenon of bridges (Sect. 3). These features are present in real galaxies (Buta 2013), but testing their connection to manifolds in specific cases of galaxies requires a particular study.

Acknowledgments

We acknowledge support by the Research Committee of the Academy of Athens through the Grant 200/895. C. E. acknowledges useful discussions with Dr. E. Athanassoula.

References

Appendix A: Series construction

Starting from the Hamiltonian (4), we implemented the method of composition of the Lie series in order to arrive at Hamiltonian (11), by using the following steps listed below.

(1) Expansion: We computed the Lagrangian points L₁ and L₂ of the Hamiltonian (3). Selecting point L₁, for instance, with coordinates (ρ_L1, ϕ_L1, 0, p_ϕ, L1), we expanded the full Hamiltonian (4) in a polynomial series in the variables δρ = ρ − ρ_L1, δϕ = ϕ − ϕ_L1, and J_ϕ = p_ϕ − p_ϕ, L1 (see Sect. 2). In our computer-algebraic implementation, all expansions were carried up to a maximum truncation order N_t, set as N_t = 10.

(2) Diagonalization: From the quadratic part H_0, 2 of the Hamiltonian, we computed the variational matrix ℳ at L₁, as in Eq. (7), as well as its eigenvalues λ_1, 3 = ±λ, λ_2, 4 = ±iκ, with λ, κ >  0, and associated eigenvectors e_i, i = 1, …, 4. Each eigenvector has four components, thus it can be written as a 4 × 1 column vector. We then formed the 4 × 4 matrix ℬ = (c₁e₁, c₂e₂, c₁e₃, c₂e₂) with the unspecified coefficients c₁ and c₂. Thus, the matrix ℬ contains the four vectors as its columns (multiplied by the c_i’s). Applying the symplectic condition ℬ^T ⋅ 𝒥 ⋅ ℬ = 𝒥, where 𝒥₄ is the 4 × 4 fundamental symplectic matrix

$$\begin{array}{c}\hfill {\mathcal{J}}_4=\left[\begin{array}{cc}0& I_2\\ -I_2& 0\\ \end{array}\right],\end{array}$$(A.1)

with I₂ equal the 2 × 2 identity matrix, yields two independent equations allowing one to specify the coefficients c₁ and c₂, and hence all the entries of the constant matrix ℬ. This matrix 𝒜 in the transformation (5) is then given by 𝒜 = ℬ ⋅ 𝒞 where

$$\begin{array}{c}\hfill \mathcal{C}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& \frac{1}{\sqrt{2}}& 0& \frac{-i}{\sqrt{2}}\\ 0& 0& 1& 0\\ & -\frac{i}{\sqrt{2}}& 0& \frac{1}{\sqrt{2}}& 0\\ \end{array}\right].\end{array}$$(A.2)

(3) Normalization using Lie series: We used the Lie method of normal form construction (see Efthymiopoulos 2012 Sect. 2.10 for a tutorial) in order to pass from Hamiltonian (8) to Hamiltonian (11). Briefly, we considered a sequence of canonical transformations (U^{(r − 1)}, Q^{(r − 1)}, V^{(r − 1)}, P^{(r − 1)}) → (U^(r), Q^(r), V^(r), P^(r)), with r = 1, 2, …, N_t, where (U⁽⁰⁾, Q⁽⁰⁾, V⁽⁰⁾, P⁽⁰⁾) ≡ (u, Q, v, P) and (U^(N_t), Q^(N_t), V^(N_t), P^(N_t)), ≡ (ξ, q, η, p) defined through suitably defined generating functions χ₁, χ₂, …, χ_Nt, through the recursive relations

$$\begin{array}{c}\hfill (U^{(r-1)},Q^{(r-1)},V^{(r-1)},P^{(r-1)})=exp\left({\mathcal{L}}_{{\chi }_{\mathrm{r}}}\right)(U^{(r)},Q^{(r)},V^{(r)},P^{(r)}),\end{array}$$(A.3)

where ℒ_χr denotes the Poisson bracket operator ℒ_χr ⋅ ={⋅, χ_r}, and $exp({\mathcal{L}}_{{\chi }_{\mathrm{r}}})=\sum _{j=0}^{\infty }{\mathcal{L}}_{{\chi }_{\mathrm{r}}}^j/j!$ truncated at order N_t. Once the involved generating functions χ_rr = 1, 2, …, N_t are specified, Eq. (A.3) allows one to define the transformation of Eq. (10), and hence the periodic orbit GL₁ through Eq. (12).

It still needs to be determined how to compute the functions χ_r. This is accomplished via a recursive algorithm, allowing one to transform the original Hamiltonian H⁽⁰⁾ ≡ H, with H given by Eq. (8) to its final form H^(N_t) given by Eq. (11). We consider the rth normalization step, and give explicit formulas in the case of one extra pattern speed in which we have one extra angle ϕ₂ = ϕ_s (generalization to M extra pattern speeds is straightforward). The Hamiltonian has the form

$$\begin{array}{c}\hfill H^{(r-1)}=Z_2+Z_3+\dots +Z_{r+1}+R_{r+2}^{(r-1)}+\dots +R_{N_{\mathrm{t}}+2}^{(r-1)},\end{array}$$(A.4)

where (i) subscripts refer to polynomial order in the variables (U^{(r − 1)}, Q^{(r − 1)}, V^{(r − 1)}, P^{(r − 1)}), and (ii) the terms Z_i, i = 2, …, r + 1 are “in normal form”, that is, they do not contain any monomials linear in (u^{(r − 1)}, v^{(r − 1)}). The remainder term $R_{r+2}^{(r-1)}$ has the form

$$\begin{array}{cc}\hfill R_{r+2}^{(r-1)}=& \sum _m\sum _{k_1+k_2+l_1+l_2=r+2}{a}_{k_1,k_2,l_1,l_2}^{(r-1)}\hfill \\ \hfill & \times e^{im{\varphi }_{\mathrm{s}}}{(U^{(r-1)})}^{k_1},{(Q^{(r-1)})}^{k_2},{(V^{(r-1)})}^{l_1},{(P^{(r-1)})}^{l_2}\hfill \end{array}$$

with k₁, k₂, l₁, l₂ ≥ 0, and m integer. Then, the generating function χ_r is given by

$$\begin{array}{cc}\hfill {\chi }_{\mathrm{r}}=& {\displaystyle \sum _{k_1+l_1=1}{a}_{k_1,k_2,l_1,l_2}^{(r-1)}}\hfill \\ \hfill & \times \frac{e^{im{\varphi }_{\mathrm{s}}}{(U^{(r)})}^{k_1},{(Q^{(r)})}^{k_2},{(V^{(r)})}^{l_1},{(P^{(r)})}^{l_2}}{(l_1-k_1)\nu +i[(k_2-l_2)\kappa +m({\mathrm{\Omega }}_{\mathrm{s}}-{\mathrm{\Omega }}_{\mathrm{bar}})]}·\hfill \end{array}$$(A.5)

With the above rule, the Hamiltonian takes a normal form up to the terms of polynomial order r + 2, namely

$$\begin{array}{cc}\hfill H^{(r)}& =exp({\mathcal{L}}_{{\chi }_{\mathrm{r}}})H^{(r-1)}\hfill \\ \hfill & =Z_2+Z_3+\dots +Z_{r+2}+R_{r+3}^{(r)}+\dots +R_{N_{\mathrm{t}}+2}^{(r)}.\hfill \end{array}$$(A.6)

Hence, repeating the procedure N_t times leads to the Hamiltonian (11).

All Figures

	Fig. 1. Rotation curve (black) corresponding to the potential V = Vh + Vd + ⟨Vb⟩, where ⟨Vb⟩ is the m = 0 (average with respect to all azimuths) part of the bar’s potential Vb. The contribution of each component is shown with a different color.
In the text

Fig. 2. Color map of the surface density Σ(ρ, ϕ) corresponding to the potential V = Vd + Vb + Vsp (see text), as viewed in the bar’s rotating frame, at four different snapshots, namely t = 0 (top left), t = T/4 (top right), t = T/2 (bottom left), and t = 3T/4 (bottom right), where T = π/|Ωsp − Ωbar|. Since Ωsp <  Ωbar, the spirals have a relative clockwise angular displacement in time with respect to the bar. However, the morphological continuity between the bar and spirals is retained in all of these snapshots.

In the text

	Fig. 3. Bar, spiral, and total Q-strengths (Qb, Qs and Qtotal respectively) as functions of the radius ρ in the model including the potential terms Vd, Vb, and Vsp for the intermediate spiral model (see text).
In the text

Fig. 4. Panel a: apocentric invariant manifolds in the pure bar model with one pattern speed. The manifolds emanating from the points L1 and L2 are plotted in blue and red, respectively. Panel b: apocentric manifolds if we add the spiral potential (intermediate case), assuming that the spirals rotate with the same pattern speed as the bar. The addition of the spiral term enhances the structures described as “lobes”, “bridges”, and “gaps” (see text). Panels c and d: same as in (a) and (b), but with the manifolds computed over a larger length. The black spiral curves correspond to the minima of the imposed spiral potential, given by Eq. (19).

In the text

	Fig. 5. Periodic orbit GL1 in the weak (left), intermediate (center), and strong (right) spiral cases. The size of the orbit increases with the spiral amplitude.
In the text

Fig. 6. “Double section” apocentric manifolds ${\mathcal{W}}_{\mathrm{GL}1}^{\mathrm{AU}}(t_{\mathrm{s}})$ (blue) and ${\mathcal{W}}_{\mathrm{GL}2}^{\mathrm{AU}}(t_{\mathrm{s}})$ (red) are plotted at four different times ts as indicated in each panel. The black-dotted curves mark the local maxima of the surface density σ(ρ, ϕ, t), corresponding to the potential V = Vd(ρ)+Vb(ρ, ϕ)+Vsp(ρ, ϕ, t) at the times t = ts. These maxima are plotted in the domain 6 kpc < ρ <  15 Kpc, where the imposed spirals have a significant amplitude.

In the text

	Fig. 7. Same as in Fig. 6, but for the weak spiral model.
In the text

	Fig. 8. Same as in Fig. 6, but for the strong spiral model.
In the text