A globally stable attractor that is locally unstable everywhere

We construct two examples of invariant manifolds that despite being locally unstable at every point in the transverse direction are globally stable. Using numerical simulations we show that these invariant manifolds temporarily repel nearby trajectories but act as global attractors. We formulate an explanation for such global stability in terms of the `rate of rotation' of the stable and unstable eigenvectors spanning the normal subspace associated with each point of the invariant manifold. We discuss the role of this rate of rotation on the transitions between the stable and unstable regimes.


I. INTRODUCTION
Invariant manifolds organize the trajectories of a dynamical system and determine in a qualitative manner the behavior of a system. An invariant manifolds can be thought of as unstable if nearby trajectories diverge away from it and stable otherwise. The global stability of an invariant manifold in R n , n ≥ 3, can be assessed using Lyapunov type numbers that were first described by Fenichel We construct two examples of dynamical systems with invariant manifolds that despite being unstable at every point in the normal direction are nevertheless globally stable. The examples are inspired by some recent findings, Ref. 4, on the dynamics of inertial particles modeled using a simplified Maxey-Riley equation, Refs. 5-7. The planar motion of inertial particles in a fluid generates a four dimensional dynamical system with the fluid streamlines forming the invariant manifold for a time independent fluid flow. Though this invariant manifold is globally stable, it was observed that this globally attracting invariant manifold contained sub-domains of local instability. 6,8 In Ref. 4 it was hypothesized that locally unstable subsets of the invariant manifold are globally stable due to the 'rotation' of the stable and unstable eigenvectors of an associated reduced two dimensional system along a particles trajectory. The physical origin of this rotation of the stable and unstable eigenvectors lies in the vorticity field of the fluid.
The two dynamical systems we discuss in this paper are constructed by explicitly using this phenomenon of the rotation of the stable and unstable eigenvectors that span the normal subspace associated with each point of the invariant manifold. We also demonstrate via numerical simulations the relationship between the global stability of the invariant manifold and the rate of rotation of these eigenvectors. The first system has an invariant manifold that is not compact, while the second system has a limit cycle. In both cases the invariant manifolds repel a subset (of nonzero measure) a ptallap@clemson.edu of trajectories in any neighborhood for a brief period of time. However every trajectory eventually converges to the invariant manifold. The concrete examples in this paper demonstrate a novel type of a global attractor that is locally unstable everywhere. It is important to draw attention to past work that demonstrate the existence of limit cycles that are locally unstable but globally (or orbitally) stable. For instance Refs. 9 and 10 discuss dynamical systems with a limit cycle a subset of which is locally unstable and a subset of which is locally stable. Some perturbations normal to a subset of the limit cycle would grow, while in the stable subset of the limit cycle, all perturbations would decay. The limit cycle is however globally or orbitally stable. Hybrid dynamical systems can also demonstrate such limit cycles. For instance Ref. 11 shows that the hybrid dynamical system of the so called passive walker has a locally unstable limit cycle, but one that is globally stable due to a reset map that makes the dynamical system piecewise smooth. In contrast in this paper we construct smooth dynamical systems with invariant sets that are locally unstable everywhere and yet are globally stable. More importantly, the constructive examples in this paper distill the mechanism by which an invariant set that is locally unstable could become globally stable.
The paper is organized as follows. In section II we review the definition of the normal stability of an invariant manifold. This stability is in a local sense and is measured by the normal infinitesimal Lyapunov exponent (NILE), Ref. 2. This indicator of stability has also been referred to as the local Lyapunov exponent, Refs. 3 and 10. In Sec. III we construct a dynamical system in R 3 that has a non compact invariant manifold that is locally unstable everywhere but is globally stable. In Sec. IV we construct a dynamical system in R 3 that has a limit cycle that is locally unstable everywhere but is globally stable.

II. NORMAL LOCAL STABILITY OF AN INVARIANT MANIFOLD
Consider a time invariant dynamical system of the forṁ where n ≥ 2 and let a trajectory of the system, ξ(t; t 0 , ξ 0 ) be defined as the solution to the initial value problem, ξ(t 0 ) = ξ 0 . The flow map, φ t t 0 : R n → R n , associated with this system is defined as the transformation A m dimensional submanifold M ⊂ R n is defined as invariant under the flow if The manifold M is a minimal invariant manifold if no proper subset of M is invariant. Let T ξ M and N ξ M denote the tangent and normal spaces to M at ξ ∈ M, i.e., T ξ R n = T ξ M ⊕N ξ M. We will denote the projection from the tangent space of R n at ξ to N ξ M by Denoting the Jacobian of the flow map φ t+s t by Dφ t+s t , the evolution of normal vectors along a trajectory is given by where u n (t) ∈ N ξ (t ) M and u n (t + s) ∈ N ξ (t + s) M. The instantaneous growth in the norm of normal vectors is captured by the so called normal infinitesimal Lyapunov exponent (NILE), defined as, The normally stable and unstable sets M s ∈ M and M u ⊂ M respectively are defined as The stable sets, M s , are regions of M where every normal perturbation contracts in norm as illustrated in fig

III. A GOBALLY STABLE NON COMPACT INVARIANT MANIFOLD THAT IS LOCALLY UNSTABLE EVERYWHERE
Consider a vector field in R 3 of the form where A(z) is a two dimensional square matrix whose entries are dependent on z. The manifold which is a straight line (the 'z-axis') is invariant. In the simple case where such that λ 1 < 0, λ 2 > 0 and λ 1 + λ 2 < 0, all trajectories eventually diverge away from the z-axis, making the invariant manifold M unstable. We now modify the matrix A(z) in the following way, where R = cos ωz − sin ωz sin ωz cos ωz with ω being a real number. We will denote the two columns of the matrix R by p 1 (z) and p 2 (z). From (10) it is clear that the eigenvalues of A are λ 1 and λ 2 with the corresponding eigenvectors being p 1 (z) and p 2 (z).

A. Local stability of the invariant manifold
The Jacobian corresponding to the linearization of (8) at any point (x, y, z) is where A ij denotes the element of the matrix A on the ith row and jth column. The Jacobian J has a block upper triangular form The eigenvectors and eigenvalues of J can be computed easily due to the block upper triangular form of J. A formal substitution shows that are two of the eigenvectors of J with the corresponding eigenvalues being λ 1 and λ 2 .
An analogous calculation can be performed for the second eigenvector w 2 . The sum of the eigenvalues of J is Where, tr(A) = λ 1 + λ 2 . Hence, the third eigenvalue of J, λ 3 , is indeed ∂g ∂z . We will choose g(z) such that ∂g ∂z 0 for all z.
At any (0, 0, z) ∈ M, the Jacobian is The third eigenvector of the Jacobian, J(x = 0, y = 0, z) is The eigenvectors w 1 , w 2 and w 3 are mutually orthogonal. The vectors w 1 and w 2 span the normal space of M, with w 1 spanning the stable subspace of N ξ M and w 2 spanning the unstable space of N ξ M. The effect of the transformation (10) is such that these vectors undergo a continuous rotation along M as shown in fig. 2. With (15), the Jacobian matrix resulting from the linearization of the dynamical system given by (8) has only real eigenvalues. Furthermore we will choosė where k = 0.01. Therefore | ∂g ∂z | |λ 1 | and | ∂g ∂z | |λ 2 |. This means the dynamics on M itself are slower than the dynamics normal to M, making M a normally hyperbolic invariant manifold (NHIM). We note however that the global attractivity of M can occur even if the dynamics on M are not slow. From the definition of A(z), λ 2 > 0 at all points ξ = (0, 0, z) ∈ M, hence the NILE, σ(ξ) = λ 2 > 0 everywhere on the invariant manifold. This implies that M is locally unstable at all points ξ ∈ M. Any neighborhood of ξ ∈ M has a subset that will be repelled away from M.
The relevance of the NILE as an indicator of local instability can be seen through a direct calculation to show that any neighborhood of z ∈ M contains a subset that is temporarily repelled from the invariant set. A perturbation u = x y is decreasing in norm if the condition d |u | 2 dt < 0 is met. Noting that |u| 2 = u T u andu = RAR T u, the rate of growth of a perturbation can be calculated as follows where B = RAR T . If every perturbation in a neighborhood of z ∈ M decays, then S(z) < 0, which requires that B(z) be negative definite. Since B is obtained from A through rotations, the trace and determinant of B are the same as that of A. This can of course be verified through a direct calculation.
Since we are considering the case where the eigenvalues of A are real and of opposite sign, with the negative eigenvalue having the larger magnitude, we obtain tr(B) = λ 1 + λ 2 < 0 and det(B) = λ 1 λ 2 < 0. Therefore B(z) is not negative definite for any value of z. Therefore any neighborhood of z ∈ M has a subset that can grow in norm, at least temporarily and the z axis is locally unstable everywhere. Suppose one chooses perturbations on a circle centered at some z, u = r(cos θ, sin θ), and denoting S(z, r, θ) = 2u T Bu, where B ij represent the entries in the matrix B. The sign of S can be found numerically for different values of θ for a given z. The locally stable subsets of a circular neighborhood of z are S (z) = {θ: S(z, r, θ) < 0}, while the locally unstable subsets of the circular neighborhood are S + (z) = {θ: S(z, r, θ) > 0}. Figure 3 shows the stable (blue) and unstable (red) subsets of a unit circle for different ranges of z along the invariant manifold. At the origin S(0, r, θ) = 2r 2 (λ 1 cos 2 θ + λ 2 sin 2 θ) and S (0) is the set that satisfies, The stable and unstable subsets on the unit circle at any value of z are obtained through a rigid rotation of these sets when z = 0. The relative measures of the unstable and stable subsets on a unit circle remain constant at any point on the z-axis.
Since the rate of growth of |u| 2 instead of u was computed in (19), the rate of growth of the perturbation itself is given by S 2r . Furthermore normalizing this by the norm of the initial perturbation, gives S 2r 2 . The relationship of S to the NILE can be understood by computing the maximum value of S 2r 2 . By inspecting (21), max( S 2r 2 ) = λ 2 , since λ 1 < 0 < λ 2 . Therefore max( S 2r 2 ) = σ, the NILE.
Here we comment on an alternative way to characterize the local stability of a limit cycle, such as in Ref. 10. A transformation of coordinates that rotates with the eigenvectors of A along M can be performed, to obtain v = R(ωz)u. Differentiating this expression with time one obtains Here we use the fact where Ω is the skew-symmetric matrix In the event thatż is constant, (23) is independent of the flow on z-axis. The stability of the fixed point of (23), u = 0 can be decided by the eigenvalues of C = (AΩż). If both the eigenvalues of C then all perturbations will decay. The eigenvalues of C are negative for a high enough ω. While this transformation to a rotating frame of reference is convenient, the local stability of M cannot be inferred from the eigenvalues of C ifż is not constant. Even whenż is constant, it can be seen that some perturbations can temporarily grow, The right hand side of (24) is the same as that of (21) when ωz = 0. From the previous analysis there is a finite subset of perturbations that grow temporarily. This local instability is in fact independent of ω, but depends only on λ 1 and integration is performed in MATLAB using the solver ode113 which is a multistep variable order Adams Bashworth Moulton solver with an error tolerance of 10 14 . The numerical integration was performed for large time periods ranging from t = 10 5 to t = 10 7 to confirm the eventual non divergence of trajectories.
These numerics are discussed here first for two special cases; a perturbation along locally unstable directions at the origin, with (x(0) = 0.0, y(0) = 0.001, z = 0) and a perturbation along the locally stable direction at the origin with (x(0) = 0.001, y(0) = 0, z = 0). Figure 4(a) shows the trajectory and Fig. 4(b) shows the distance d(t) of the trajectory from the z axis. Since the perturbation is along the locally unstable direction at the origin, the distance of the trajectory from the zaxis initially increases more than six fold. However as the trajectory simultaneously spirals around the zaxis, its distance from the z-axis decreases and the trajectory converges to the z axis.   Fig. 5(b) shows the distance d(t) of the trajectory from the z axis. Since the perturbation is along the locally stable direction at the origin, the distance of the trajectory from the z axis initially decreases. However as the trajectory simultaneously spirals around the z-axis, it begins to be repelled. This repulsion increases the distance of the trajectory from the z axis before it eventually decreases again and converges to zero.
The surprising decay of normal perturbations can be understood through a simpler example; a planar dynamical system with a saddle type fixed point at the origin, where µ 1 < 0 and µ 2 > 0. If the initial conditions are (x 1 (0), x 2 (0)), the distance of the trajectory from the origin is D(t) = x 1 (0) 2 e 2µ 1 t + x 2 (0) 2 e 2µ 2 t . Rescaling the initial conditions as x 2 (0) = x 1 (0), the distance is When = 0, D(t) decays to zero. For a small enough > 0 the distance of the trajectory to the origin first decays, in an interval (0, T ) before increasing monotonically for t > T. The critical value of for which the D(t) is an increasing function at t = 0 can be obtained by setting dD 2 (t) dt | t=0 > 0. This gives the critical value of , (27) is the same as tan θ in (22) and the two equations are equivalent in identifying the locally unstable subsets of a neighborhood of the invariant sets (8) and (25).

If y(t) x(t) < cr then D(t) is a decreasing function. When y(t) x(t) = cr D(t) is a minimum and increases thereafter. It should be noted that in
We will denote the subset of R 2 where D(t) is decreasing, by B s and call it the decay set, The sets B s ⊂ R 2 and S are both trivially related; S is merely the arc of a circle contained in the set B s . The rate of growth of perturbations, S, (20) is graphed in fig. 6 for two sets of parameters. The red dotted graph is for the case λ 1 = 1.1, λ 2 = 0.1 and ω = 25 and the blue dotted graph is for the case λ 1 = 1.1, λ 2 = 0.9 and ω = 98. In both cases λ 1 < 0 and λ 1 > λ 2 . This leads to the arc length of S to be greater than S + and the maximum magnitude of decay to be larger than the maximum magnitude of repulsion. The role of the saddle is played by the z-axis in the dynamical system (8). The difference is that the stable and unstable subspaces undergo a continuous rotation along z axis. A trajectory (x(t), y(t), z(t)) = ξ(t; t 0 , ξ 0 ) moves closer to the invariant manifold when (x(t), y(t)) ∈ B s . The decay set too undergoes a rotation along the z axis. When perturbed trajectories enter this decay set their distance to the zaxis decreases.
The decay set for the system (25) is shown in gray in fig. 7 where we chose µ = 1.1 and µ = 0.1. In fig. 7(a)-(b) the trajectory at (x(t), y(t)) lies in the region where d(t) increases, while in Figs. 7(c)-(d) the trajectory at (x(t), y(t)) lies in the decay set. The initial conditions for this trajectory are the same as those in fig. 4. The distance from the z-axis increases initially and from about t = 56 decreases as shown in fig. 4(b). This is when the trajectory enters the decay set (shown in gray in fig. 7(c)). The increase or decrease of the distance of the trajectory from the zaxis is determined by whether or not the trajectory lies in the decay set as shown in fig. 7.
The temporary growth of perturbations and their eventual decay is observed all along the z axis with the repulsion of experienced by a trajectory varying. Figure. 8 shows the projection on to the x y plane of the evolution of the perturbations at two different values of z(0) away from the origin.
A circle whose radius is equal to the norm of the initial perturbation is also shown to illustrate the evolution of the distance of the trajectory from the zaxis. The stable and unstable subsets of the circle are shown in blue and red respectively. A sample perturbation from each subset is shown in both the figures. When perturbations begin in the stable subset they first decay while perturbations that begin on the unstable subset first increase in norm. The local repulsion of a subset of trajectories and the eventual decay of all perturbations is observed in both cases.
Perturbations are guaranteed to decay only if every trajectory either eventually lies in the stable set or oscillates between the stable and unstable sets but spends more time in the stable set. Transforming the equation of the dynamical system (8) to a polar form, x = r cos θ and y = r sinθ, where r = x 2 + y 2 , one can show thatθ The initial perturbation is shown by a black filled circle. Two perturbations are shown in each case. One perturbation is chosen from unstable subset (dashed red circular arc) and one from the stable subset (dotted blue arc).
We define a new variable β = θ ωz, which represents the relative polar angle of the trajectory with respect to the stable and unstable eigenvectors at some z. Theṅ Equation (30) is independent of r, but depends on z. When ωk is small instantaneous stagnation points (ISP) of (30) plays a dominant role in the evolution of the β. At each z(t), the ISPs are obtained by setting the right hand side of (30) to zero. There are four ISPs of (30). For example at z = 0, fig. 9(a) showsβ by the solid green curve which is zero at four distinct values of θ. Two of these, shown by black solid dots, are stable due to the negative slope. The time evolution of these ISPs are the so called distinguished hyperbolic trajectories and all solutions of (30) converge to a neighborhood of these distinguished trajectories, Ref. 12. Figure 9(a)-(c) shows such a convergence of eight solutions of (30) with distinct values of θ(0) shown by the red filled circles. The growth or decay of perturbations of the dynamical system (8) to (x, y) = (0, 0) is determined by the evolution of the stable ISPs of (30). As the position along M varies, the two stable ISPs could lie either in the stable set S (z)or the unstable set S + (z). The blue (dashed) curve in fig. 9(a)-(c) shows the instantaneous rate of repulsion, S. When the stable ISP lies in the stable set (where the graph of S is negative) a perturbation away from the z axis decays. Numerical simulations show that the stable ISPs lie in the stable set S (z) for a longer duration of time than in the set S + (z), Fig. 9(d). Sinceż is periodic, with period π, the graph of z isp (t) in Fig. 9(d) is also periodic. The repulsion or decay rate of a trajectory is repeated along the z axis. Furthermore the highest instantaneous rate of repulsion of perturbations is smaller than the highest instantaneous rate of attraction if λ 1 + λ 2 < 0, see for example fig. 6. This leads to the eventual decay of all perturbations. Figure 9 also explains the weaker repulsion of perturbations seen in fig. 8 compared to the case shown in fig. 4. The observed repulsion is the least in the neighborhood of z = (2n + 1) π 2 for any integer n. This is because the stable ISPs lie in the stable set S in a large neighborhood of z = π 2 , fig. 9(d). Since convergence of any solution of (30) to one of the stable ISPs occurs very rapidly, any perturbation whether it begins in the stable or the unstable set, decays rapidly, experiencing only a small transient repulsion. Repulsion of perturbations around the z axis is the highest in the neighborhood of z = nπ, since the stable ISPs lie in the unstable set S + (z) in a large neighborhood of z = nπ, see fig. 9(d). When ωk λ 1 −λ 2 2 , instantaneous stagnation points do not exist for some or all values of z and the relative angle β is driven to oscillate with a large amplitude (≈ωk). In this case a trajectory, (x(t), y(t), z(t)) moves into and out of the decay set S (z) almost periodically. However since the rate of decay in the decay set is larger than the rate of repulsion, the perturbation eventually decays with oscillations. An example of such a case is shown in fig. 12(b), where the norm of a perturbation decays to zero with oscillations.

C. Parametric dependence of global stability of the invariant manifold
It should be obvious that if ω = 0 then the zaxis is both locally and globally unstable. Any normal perturbation such that y(0) 0 would lead to the trajectory diverging from the z axis. Therefore the global stability of the z axis certainly depends on the rate of the rotation of the stable and unstable eigenvectors along the z axis. Starting from ω = 0, we performed numerical simulations for increasing values of ω to determine the critical value of the rate of rotation at which the z axis become globally stable. These simulations show that the transition of the z-axis to a globally stable manifold occurs in four stages. In each of these stages the behavior of the function d(t) the distance of a trajectory from the z-axis is qualitatively distinct.
For small values of ω any trajectory diverges away from z axis and d(t) is a monotonically increasing function, see fig. 10(a). In the second stage the distance of a trajectory undergoes small oscillations with the mean distance increasing. Figure 10(b) shows these oscillations in the log 10 (d(t)). We identified the local maxima of the function log 10 (d(t)) and find that these maxima are linearly increasing. As ω increases further, the growth in these local maxima of log 10 (d(t)) decreases. At a certain critical value of the rate of rotation, ω cr ≈ 22.4686859, we find that the distance function, d(z), is periodic, see fig. 10(c). Numerical integration of (8) for very large periods of time yields a trajectory that is bounded for perturbations across a range of magnitudes. Figures 10(c) and 10(d) show the periodic variation in the distance of a trajectory from the zaxis when the norm of the initial perturbation is 10 8 and 10 3 respectively. The time period between consecutive maxima in the graphs is T = 222.1442 accurate up to three decimal places. Identifying the local maxima of the distance of this trajectory, shows that these maxima are nearly constant with variation on the order of 10 13 and 10 9 . The numerics therefore indicate bounded trajectories whose distance from the zaxis varies periodically. For ω > ω cr normal perturbations from the z axis decay with the distance decreasing with oscillations. Increasing the value of ω even up to 20000 showed that all normal perturbations do not decay monotonically. This is to be expected since any perturbation in the locally repelling direction has to first increase before decaying. The value of ω cr depends on λ 1 and λ 2 for a fixed vector fieldż. The critical angular velocity, ω cr , increases with λ 2 and the magnitude of λ 1 . This relationship is shown in fig. 11. When λ 1 + λ 2 is closer to zero, the temporary repulsion of the trajectories from the zaxis is more dramatic, as shown in fig. 12. The z axis loses its global stability if λ 1 + λ 2 > 0.

IV. A GOBALLY STABLE LIMIT CYCLE THAT IS LOCALLY UNSTABLE EVERYWHERE
The invariant manifold of the dynamical system (1) is unbounded. We construct a dynamical system in R 3 that has an invariant limit cycle. The equations for this dynamical system are arrived by forcing the normal stable and unstable eigenvectors at each point on the limit cycle to undergo a rotation. We express these equations in cylindrical coordinates, where the rotation of the stable and AIP Advances 7, 125012 (2017) unstable eigenvectors is easily seen. The equations are where As in the earlier example, we choose λ 1 < 0 and λ 2 > 0. The set is an invariant manifold (a limit cycle) for the system (31). The limit cycle repels almost all normal perturbation for an intermediate period of time. If the rate of rotation, ω is sufficiently high, these repelled trajectories eventually converge to the limit cycle. Figure 13(a) shows a trajectory with initial conditions (r(0) = 1.1, y(0) = 0.1, θ = 0). The trajectory spirals around the limit cycle and converges to it eventually. The distance of the trajectory to the limit cycle does not decrease monotonically. It initially decreases, then increases and eventually decreases as shown in fig. 13(b). The initial conditions chosen for this example are generic, containing perturbations from the limit cycle in both the stable and unstable direction. . In (a) the distance grows more than 10 times its initial value and in (b) it grows about 1.4 times its initial value before decaying to zero. As in the previous example the global stability of the limit cycle depends on the value of ω for a fixed set of λ 1 and λ 2 . In this case too, there are four different regimes of ω where the behavior of trajectories is qualitatively different. Exploring a range of values of ω, we find that for small values of ω, trajectories that begin close to the limit cycle diverge with the distance to the limit cycle, d(t), increasing monotonically. As the value of ω increases, this divergence is such that d(t) decreases at regular intervals of time, but on the average increases, as in fig. 10(b). At a critical value of the rate of rotation, ω cr ≈ 22.4686859 trajectories that begin close to the limit cycle are bounded but do not converge to the limit cycle, fig. 13(a). The distance of such a generic trajectory, (r(0) = 1.1, y(0) = 0.1, θ(0) = 0), is shown in fig. 13(b) which oscillates periodically. The trajectory itself moves through two neighborhoods around the limit cycle, one where the distance from the limit cycle increases and one where it decreases. Numerical simulations suggest that these neighborhoods of the limit cycle where distance decays are nearly independent of the initial conditions.

V. CONCLUSION
We have shown that the invariant manifold of a dynamical system in R n , n = 3 can be such that it is locally unstable everywhere, but still act as a global attractor. These examples can be trivially extended to higher dimensions, n > 3, by increasing the dimension of the invariant manifold. It turns out that a simplified form of the Maxey-Riley equation for the motion of inertial particles in a fluid is in fact a non trivial extension of the examples considered here to R 4 , Ref. 4. The findings in this paper clearly attribute the cause of the global stability of the invariant manifold, M, to the rate of rotation, ω, of the stable and unstable eigenvectors spanning the normal subspace at each point on M. Our findings in this paper should be of interest from both a theoretical perspective and from the point of view of applications such as those in micro scale fluid flows where globally stable invariant regions of a fluid are sought to be created that act as particle traps.