Detecting invariant manifolds using hyperbolic Lagrangian coherent structures

IAA-AAS-DyCoSS1 -08-06 DETECTING INVARIANT MANIFOLDS USING HYPERBOLIC LAGRANGIAN COHERENT STRUCTURES Daniel P´ erez∗ , Gerard G´ omez† , Josep J. Masdemont‡ Using as reference test model the Planar Circular Restricted Three Body Problem, this paper explores its Lagrangian Coherent Structures, as well as its Hyperbolic Lagrangian Coherent Structures. The purpose is to identify stable and unstable manifolds acting as separatrices between orbits with different qualitative behaviour and, therefore, relevant to the dynamics of the problem. Particular attention is given to the manifolds associated to the collinear libration points and to the practical stability regions around the triangular equilibrium points. INTRODUCTION Lagrangian Coherent Structures (LCS), introduced by G. Haller et al.5for the study of dynamical systems, give a methodology to identify the boundaries between regions in the configuration space with orbits that have different dynamical behaviour. LCS are usually computed by means of Finite Time Lyapunov Exponents (FTLE) at time T, just looking at the FTLE scalar-field. The value of the FTLE at ~x gives an idea of the behaviour of the orbits around ~x: if the FTLE is small then the orbits in a neighbourhood of ~x will be close at time T; however, if FTLE is high then the image under the flow, up to t=T, of points close to ~x will have different behaviours and usually some of them will tend to depart from the others. In particular a high value of a FTLE can be an indicator of the existence of an invariant manifold. According to Haller,5LCSs are defined as follows. Let D ⊂ Rnbe an open set and fthe vector field: f:D×R−→ Rn (~x, t)7−→ f(~x, t). Given a dynamical system defined by the ordinary differential equation: ˙ ~x(t) = f(~x(t), t), whose associated flow will be denoted by φ φ:R×R×D −→ D (t, t0, ~x0)7−→ φ(t;t0, ~x0), ∗Departament de Matem` atica Aplicada i An` alisi, Universitat de Barcelona, Gran Via 585, 08007 Barcelona, Spain †IEEC & Departament de Matem` atica Aplicada i An` alisi, Universitat de Barcelona, Gran Via 585, 08007 Barcelona, Spain ‡IEEC & Departament de Matem` atica Aplicada I, ETSEIB, Universitat Polit` ecnica de Catalunya, Diagonal 647, 08028 Barcelona, Spain 1 where ~x0denotes the initial condition at t=t0. The value of the FTLE at ~x from t=t0up to t=Tis defined as: σT t0(~x) := 1 |T|ln      dφt0+T t0(~x) dx     2 ,(1) this is, the FTLE is the normalised square of the largest eigenvalue of: dφt0+T t0(~x) dx !∗ · dφt0+T t0(~x) dx !≡Dφt0+T t0(~x)∗·Dφt0+T t0(~x). A linear analysis shows that the eigenvector associated to the largest eigenvalue of (Dφt0+T t0(~x))∗· (Dφt0+T t0(~x)), which is usually called the Cauchy-Green tensor, gives the maximum expansion direction in a neighbourhood of ~x. For a two-dimensional dynamical system (n= 2) the LCS with initial and final times t0and T, respectively, is an injective curve ~c =~c(s)⊂ D ⊂ R2parametrised by s∈(a, b)⊂Rsuch that verifies the following two conditions (in the n-dimensional case, the LCS will be a codimension 1 manifold): 1. ~c 0(s)is parallel to the gradient of the FTLE field at ~c(s), i.e: ~c 0(s)k ∇σT t0(~c(s)) ∀s∈(a, b). 2. The unit normal vector ~n(s)to ~c(s)verifies: ~nTΣ~n = min kuk=1 ~uTΣ~u < 0, for all s∈(a, b), where Σis the Hessian of the FTLE field defined by: Σ = Σ(~x, t0, T) = d2σT t0(~x) d~x2. Some examples show that not all the LCS are related with invariant manifolds. To deal with this situation, G. Haller6introduced the so called Hyperbolic LCS (HLCS). For this purpose, consider a smooth curve M(t0)⊂ D at t=t0, which is advected by the flow into a time-evolving material line defining the so called material surface M(t)⊂ D. To measure how strongly repelling the material line M(t0)is, at each point x0∈ M(t0)we select a unit normal vector ~n0to M(t0)(~n0∈N~x0M(t0)) and follow its evolution under the linearised flow given by Dφt0+t t0(~x). Denote by ρt t0(~x0,~n0)the length of the normal component to M(t)of Dφt0+t t0(~x0)~n0 ρt t0(~x0,~n0) = ~nt, Dφt0+t t0(~x0)~n0, where ~ntis the unit normal vector to M(t)at ~xt=φt0+t t0(~x0).ρt t0(~x0,~n0)is called the normal repulsion rate of M(t)along the orbit φ(t;t0, ~x0). If is larger (smaller) than one, then M(t)has been overall repelling (attracting) between t0and t0+talong the trajectory starting at ~x0. In order to deal with any possible tangential growth within M(t)larger than the growth normal to M(t), the so called repulsion ratio is introduced as νt t0(~x0,~n0) = min e0∈Tx0M(t0)~nt, Dφt0+t t0(~x0)~n0 kDφt0+t t0(~x0)~e0k, 2 where Tx0M(t0)denotes the set of unit tangent vectors to M(t0)at ~x0. With the above definitions, M(t)is said to be normally hyperbolic over [t0, t0+T]if there exists constants a, b > 0such that for all ~x0∈ M(t0)and ~n0∈N~x0M(t0)we have: •ρt0+T t0(~x0,~n0)≥eaT , •νt0+T t0(~x0,~n0)≥ebT . Ahyperbolic LCS is a normally hyperbolic material surface M(t)such that its normal repulsion rate admits a point-wise non-degenerate maximum value along M(t)for all locally C1-close material surfaces. Haller6gives also a characterisation of the hyperbolic LCS by means of the eigenvalues λ1≤ λ2≤... ≤λn, and their associated eigenvectors ξ1, ξ2, ..., ξn, of the Cauchy-Green positive definite matrix Dφt0+T t0(~x0)∗Dφt0+T t0(~x0). The result is summarised in the following: Theorem: A compact material surface M(t)⊂ D over the interval [t0, t0+T]is a hyperbolic LCS if and only if the following conditions are verified for all ~x0∈ M(t0): 1. λn−1(~x0, t0, T)6=λn(~x0, t0, T)>1, 2. ξn(~x0, t0, T)⊥Tx0M(t0), 3. h∇λn(~x0, t0, T), ξn(~x0, t0, T)i= 0, 4. The matrix L(~x0, t0, T)defined by L=       ∇2C−1[ξn, ξn, ξn, ξn] 2λn−λ1 λ1λnhξ1,∇ξnξni. . . 2λn−λn−1 λn−1λnhξn−1,∇ξnξni 2λn−λ1 λ1λnhξ1,∇ξnξni2λn−λ1 λ1λn. . . 0 . . .. . ..... . . 2λn−λn−1 λn−1λnhξn−1,∇ξnξni0. . . 2λn−λn−1 λn−1λn        where Cdenotes the Cauchy-Green tensor at ~x0, is positive definite for all ~x0∈ M(t0). This result gives an algorithm to detect hyperbolic LCS: Step 1 Compute the two highest eigenvalue fields λnand λn−1as well as ξnfor all ~x0∈ D. Step 2 Determine the solution set Z Z={~x0∈ D | h∇λn(~x0, t0, T), ξn(~x0, t0, T)i= 0} Step 3 Identify repelling weak LCS (WLCS) at t=t0as the subsets of ZW LCS ⊂Zsuch that: 1. λn−16=λn>1. 2. ξn⊥Tx0ZW LCS Step 4 Identify repelling LCS as the (n−1)-dimensional surfaces ZLCS ⊂ZW LCS such that the matrix Lis positive definite. Step 5 Repeat steps 1-4 backwards in time (from t0+Tto t0) to obtain attracting LCS. Step 6 Verify the robustness of the LCS found (see Haller6for details). 3 COMPUTATION OF LCS The first step to be done for the computation of LCS in a domain, is to set a grid of points in it. Then, for every point ~x in the grid, the Cauchy-Green tensor C=Dφt0+T t0(~x)∗·Dφt0+T t0(~x)= J∗·J, together with its eigenvalues and eigenvectors, must be computed. Concerning these two points, some technical details should be decided. The first one has to do with the computation of the matrix C, which can be done, at least, in two different ways: •Using numerical differentiation for the computation of the derivatives that appear in its components. In this way, the matrix Jbecomes: J=         φt0+T t0(~x+∆x1)−φt0+T t0(~x−∆x1) 2δx1[1] . . . φt0+T t0(~x+∆xn)−φt0+T t0(~x−∆xn) 2δxn[1] . . ..... . . φt0+T t0(~x+∆x1)−φt0+T t0(~x−∆x1) 2δx1[n] . . . φt0+T t0(~x+∆xn)−φt0+T t0(~x−∆xn) 2δxn[n]         , where the subindex [i], i = 1, . . . , n denotes the i-th component of the associated vector, and ∆xi= (0, . . . , δxi,...,0). If we compute the advected grid (i.e. compute the values φt0+T t0(~x)for all the points ~x of the grid) then the matrix Jcan be determined without extra integrations, except for some points in the border of the region explored. Note that in the advected grid there are all the values φt0+T t0(~x ±∆xi)that are needed to compute J. •Integrating the differential equations of the dynamical system together with the associated variational equations. In this way, a linear system of n×ndifferential equations, for the components of the matrix J, is added to the original one. In exchange to this increment of equations usually more precision is obtained. This has been the option selected for the computations shown in this paper. An important issue that must be considered is that the components of Jcan be very large and, therefore, the accurate computation of its eigenvalues and eigenvectors be difficult. To overcome this, is convenient to split the total time interval [t0, t0+T]in Nintervals with endpoints at t0, t1, . . . , tN−1, tN=t0+Tas well as power’s method for the computation of the eigenvalues and eigenvectors. In this way the matrix Jcan be written as the composition of Nmatrices J=Jt0→tN=JtN−1→tNJtN−1→tN−2·. . . ·Jt1→t2·Jt0→t1, so we can apply power’s method step by step. The gradient of the greatest eigenvalues field ∇λn(~x0, t0.T)should be computed at all the points of the grid (Step 2). Since the eigenvalues are computed in a grid of points the values of ∇λn(~x0, t0.T) can be easily obtained using numerical derivatives with the values of λn(~x0, t0.T)at the points of the grid, again except at some border points. Another difficulty that arises is that, if the radius of the grid size is not very small, only a reduced set of points of Zis determined, even using a relatively large threshold for the orthogonality condition required in Step 2. To overcome this problem we have implemented a procedure to find, in 4 an iterative way, a set covering of Z. The procedure recalls the two-step algorithm implemented in GAIO (see Dellnitz et al.2) to find cover sets of attractor in dynamical systems. Each iteration of the method consists in two steps: division and selection. In the division the the cells of the recovering are divided in smaller ones. In the selection step it is verified which ones of the new cells contain (or can contain) the desired object. We start with rectangle R0(defined by its centre c0= (c0 1, . . . , c0 n)and its radius r0= (r0 1, . . . , r0 n)) where we look for the points of Z. Note that both, the centre and the radius, are vectors, therefore, R0={x= (x1, . . . , xn)|c0 i−r0 i≤xi≤c0 i+r0 i, i = 1, . . . , n}. Using R0as initial rectangle, the two steps of the iterative procedure are: 1. Division: Divide the rectangle in two smaller rectangles with radius (ri+1 j) and centre (ci+1 j,±) given by: ri+1 j=(ri j 2if i=j(mod n) ri jotherwise , ci+1 j,±=(ci j±ri+1 jif i=j(mod n) ci jotherwise . 2. Selection: In this step we must determine if the divided rectangles contain or not points of Z. The selection criteria varies along the iterative procedure. In the first selection steps we ask for weak conditions while in the last steps we ask for a stronger one according to the following. A rectangle Rjis selected if: •exist at least a point x∈Rjsuch that it verifies the first condition of the hyperbolic LCS characterisation (for values of jless than a fixed given one). •exist at least a point x∈Rjsuch that it verifies the first and third conditions of the hyperbolic LCS characterisation (for values of jgreater than a fixed given one). Since it is not possible to verify all the points in Rj, a grid is defined inside the rectangle. Then, inside the grid we do the same kind of computations as for the non GAIO-like procedure. This two selection styles are included in order to capture as many LCS as possible. The second style is the one which allows the determination of the set Zgiven in step 2. But in the initial iterates the grid is not enough fine to detect the orthogonality condition. Therefore the algorithm eliminates some cells that contains points in Z. To avoid this situation the first style is included. Even with the usage of this two conditions some rectangles are discarded although they contain a portion of Z. To avoid that the semi-selected rectangles are introduced. A rectangle is said semi-selected if it does not verifies the selection condition but it is besides one verifying it.Then a rectangle is said selectable (i.e. the selection step should be applied on it) if: •it is contained in a selected or semi-selected rectangle in the previous iterate, or •it is beside a selected rectangle in the actual iterate. So at the end of an iterate all the selected rectangles will be surrounded by semi-selected ones. Doing this it is possible to add some of the rectangles which are discarded in some previous steps. 5 As final technical remark, due to the great amount of numerical integrations to be performed during the procedure, is important to use a fast numerical integrator. In the examples below, a Taylor method has been chosen (see Jorba et al.7), which is both accurate and fast for the differential equations considered. LCS IN THE RESTRICTED THREE BODY PROBLEM The Circular Restricted Three Body Problem (CRTBP) describes the motion of a body with negligible mass under the gravitational attraction of two primaries, with normalised masses µand 1−µ, that rotate around their common centre of masses following circular orbits. In the planar case and using a synodical reference frame, the motion of the negligible mass follows is given by (see Szebehely9): (¨x−2 ˙y= Ωx ¨y+ 2 ˙x= Ωy, where Ω(x, y) = x2+y2 2+1−µ √(x−µ)2+y2+µ √(x+1−µ)2+y2+µ(1−µ) 2. This system has a first integral given by: C(x, y, ˙x, ˙y) = −( ˙x2+ ˙y2) + 2Ω(x, y). In the computations that follow, we have used the mass ratio of the Earth-Moon system µ= 0.0121506683. It is well known that the planar CRTBP has three collinear equilibrium points of the center × saddle type. Due to this fact, and according to Lyapunov centre theorem, it is possible to find periodic orbits (corresponding to the centre part) around them. Additionally, due to the saddle part, these periodic orbits have invariant manifolds associated, which can be computed with high precision either using numeric or semi-analytic procedures. We will use this last approximation as a test for the invariant manifold identification using LCS. Invariant manifolds and LCS associated to periodic orbits around L1 First we remark that the phase space of the CRTBP has dimension four so, in order to visualise the results of the preceding section, we must do some reduction. This is done using the Jacobi first integral (we fix its value to a constant) and considering the intersections of the orbits with a fixed hyperplane. This has already been used by Gawlik et al.3for the elliptic CRTBP. There, they introduce a mixed position-velocity space using xand ˙xas the two dimension space, fixing y= 0 and determining ˙yaccording to the value of the Jacobi constant chosen. Other strategies can be used to to reduce the dimension of the system, for instance Tricoche et al.10 use as reduction method the Periapsis Poincar´ e Mapping, and compute the LCS for that particular map. In the computations that follow, we have fixed the hyperplane x=x0(for some x0fixed and usually small) as surface of section, and (y, ˙y) are the two free variables of the system. For a fixed value of the Jacobi constant, C0, one can determine ˙xfrom: C(x0, y, ˙x, ˙y) = C0.(2) An advantage of using the x= 0 hyperplane is that the hyperbolic invariant manifolds associated to the periodic orbits considered are transversal to it. Hence, one can determine their intersections with this hyperplane without problems (without tangencies) as is shown in Fig. 1. 6 Figure 1. For C(x, y, ˙x, ˙y)=3.1638 and x= 0 (top) and x= 0.2(bottom), (y, ˙y) FTLE field together with the first two intersections of the unstable invariant manifold of the Lyapunov periodic orbit about L1(left hand side figures). In the right hand side plots, only the intersections of the manifolds with the two hyperplanes have been represented. Removing some spurious regions Figure 1 shows, in the (y, ˙y) space, the regions with high values of the FTLE field as well as the first two intersections with x= 0 and x= 0.2of the invariant unstable manifold of the periodic orbit associated to L1for C0= 3.1638 = C(L1)−0.035. There is a good agreement between the intersections of the manifolds with some of the regions with high FTLE values. The computation of the FTLE has been done using a final time T=−10 and a grid of 200 ×200 points. In the same figure, some other regions with high FTLE values appear that can not be identified with intersections of the unstable manifolds under consideration. We call them false positives. More false positives become evident if we move the intersecting hyperplane x0. For instance, Figure 1 7 (bottom right) shows the FTLE for the hyperplane x0= 0.2with the same values of the Jacobi constant, C0, the final time Tand the grid dimensions of the top figure. This false positives are due to the procedure used for the computation of the FTLE. Recall that the FTLE give an indication of the variations (sensitivity) of the final conditions of orbits at t=T with respect to their initial conditions. Large values for this sensitivity may be due to the fact that an hyperbolic manifold is acting as a separatrix of two close initial conditions or to some other effect. In the case of the false positives detected in Figure 1 they are due to close approaches of the massless particle to the large primary. In Figure 2 we have represented both the distance of the massless particle to the large primary at t=T(lower curve) together with the FTLE values (for the same value of T) of orbits with initial conditions along the line (x, y, ˙x, ˙y) = (0.2, y, ˙xc,−7y−3), where the value of ˙xchas been chosen such that the Jacobi constant of all them is equal to C0= 3.1638. From this figure it clearly follows that the false positives are associated to close approaches to the large primary. To remove these false positives, one option is to compute the FTLE at the apoapsis t=Ta. To do this, we first integrate the equations of motion, together with their variational equations, from time t0up to time T. Once this time is reached we continue the integration until the next apoapsis of the orbit. Observe that now the final time of the computation is not constant for all the points of the grid. 0 0.5 1 1.5 2 2.5 -0.32 -0.3 -0.28 -0.26 -0.24 -0.22 y Figure 2. FTLE (upper curve) and final distance to the big primary (lower curve) for initial conditions of the form (x, y, ˙x, ˙y) = (0.2, y, ˙xc,−7y−3) Therefore two similar orbits can look different at time Tdue to the fact that one of them is close to one of the primaries while the other is still far. The first orbit will be, after some time, far from the primary but not at time T, so its seems convenient to allow variations of the final time. Implementing the computation of the FTLE with the apoapsis strategy we get Figure 3, there we can see that the false positives appearing in Figure 1 have almost disappeared. FTLE vs invariant manifolds Once the false positives have been eliminated, we can try to apply the detection algorithm of the LCS. Besides the computation of the highest eigenvalues, the first thing to do is to compute the set 8 Figure 3. y-˙yapoapsis-FTLE field with x= 0.0(left) and x= 0.2(right) for C(x, y, ˙x, ˙y)=3.1638, together with the first two intersections of the invariant manifold with the hyperplane x=x0. Zof Step 2 of the proposed algorithm. Observe that this set is defined by means of an equality so, in the numerical computations, it is substituted by: Ztol ={~x0∈ D | h∇λn(~x0, t0, T), ξn(~x0, t0, T)i<tol}, where tol is a certain tolerance. The set Ztol determined using the computations done to get the FTLE field represented in Figure 1 taking tol=0.1 does not detect the invariant manifolds. If we add the first condition of Step 3 (λn>1) then it is clearly seen (in Figure 4(top) that the points that fulfil both conditions remain on one of the intersections of the manifold, nevertheless there are still very few points verifying both conditions. To enlarge the set of points in Ztol, in principle, we could reduce the radius of the grid in order to have more points available for the required computations. Figure 4 shows the set Z0.1in a smaller region with the same number of points in the grid but with a smaller grid radius (consequently the grid is finer). Observe that now there are more points verifying both conditions and, in fact, we are able to detect more intersections of the manifolds, which are associated to other branches of the invariant manifolds or with subsequent intersections of the ones that already appeared in the preceding figures. Therefore, in principle, the computation of the manifolds with a certain accuracy will require a very fine grid. Nevertheless, if we use this strategy there will be a huge set of points that do not verify the required conditions but still will be computed. Refined numerical determination of hyperbolic LCS To overcome the difficulties associated with the use of a very large set of test points, most of them not in Ztol, the already explained algorithm, that recalls the one used in GAIO2for the computation of attractors in dynamical systems, has been used. Therefore, the regions close to the set Ztol are computed with high precision (i.e. with a fine grid) while the points far from Ztol are computed with low precision. 9