Globally optimal stretching foliations of dynamical systems reveal the organizing skeleton of intensive instabilities
GGlobally optimal stretching foliations of dynamical systemsreveal the organizing skeleton of intensive instabilities
Sanjeeva BalasuriyaSchool of Mathematical Sciences, University of Adelaide, Adelaide SA 5005, AustraliaErik M. BolltDepartment of Electrical and Computer Engineering and C S the Clarkson Center forComplex Systems Science, Clarkson University, Potsdam, New York 13699January 27, 2021 Abstract
Understanding instabilities in dynamical systems drives to the heart of modern chaostheory, whether forecasting or attempting to control future outcomes. Instabilities in thesense of locally maximal stretching in maps is well understood, and is connected to the con-cepts of Lyapunov exponents/vectors, Oseledec spaces and the Cauchy–Green tensor. In thispaper, we extend the concept to global optimization of stretching, as this forms a skeletonorganizing the general instabilities. The ‘map’ is general but incorporates the inevitabilityof finite-time as in any realistic application: it can be defined via a finite sequence of dis-crete maps, or a finite-time flow associated with a continuous dynamical system. Limitingattention to two-dimensions, we formulate the global optimization problem as one over arestricted class of foliations, and establish the foliations which both maximize and minimizeglobal stretching. A classification of nondegenerate singularities of the foliations is obtained.Numerical issues in computing optimal foliations are examined, in particular insights intospecial curves along which foliations appear to veer and/or do not cross, and foliation be-havior near singularities. Illustrations and validations of the results to the H´enon map, thedouble-gyre flow and the standard map are provided.
Keywords : Lyapunov vector , finite-time flow, punctured foliation
Mathematics Subject Classification : a r X i v : . [ n li n . C D ] J a n Graphical Abstract
Sanjeeva Balasuriya, Erik Bollt
Highlights
1. Understanding the organizing skeleton of instability for orbits must be premised on anal-ysis of globally optimal stretching.2. Provides the theory to obtain foliation for globally optimizing stretching for any two-dimensional map (analytically specified, derived from a finite-time flow or a sequence ofmaps, and/or given via data);3. Classifies singularities and provides insight and solutions to spurious artefacts emergingwhen attempting to numerically determine such a foliation;4. Establishes connections with a range of well-established methods: locally optimizingstretching, Cauchy–Green eigenvalues and singularities, Lyapunov exponents, Lyapunovvectors, Oseledec spaces, and variational Lagrangian coherent structures.2
Introduction
A central topic of dynamical systems theory involves analysis of instabilities, since this is thecentral ideas behind the possibility of forecast time horizon, or even of ease of control of futureoutcomes. The preponderance of work has involved analysis of local instability, whether by theHartman-Grobman theorem and center manifold theorem [1] for periodic orbits and similarlyfor invariant sets [2]. For general orbits, local instability is characterized by Oseledec spaces[3] which are identified via Lyapunov exponents [4] and Lyapunov vectors [5, 6]. Via thesetechniques, locally optimizing stretching due to the operation of a map from subsets of R n tosubsets of R n is well-understood. Computing the map’s derivative matrix at each point is allowsfor computation of Oseledecs/Lyapunov information: its singular values and correspondingsingular vectors are respectively associated with stretching rates and relevant directions in thedomain, and its (scaled) operator norm is the classical Lyapunov exponent of the orbit beginningat that point.In this paper, we assert that understanding the global dynamics—how a system organizesorbits—is related to a global view of instabilities. The related organizing skeleton of orbitsmust therefore be premised on analysis of globally optimal stretching . Here, orbits will be inrelation to two-dimensional maps which can be derived from various sources: a finite sequence ofdiscrete maps, or a flow occurring over a finite time period. The latter situation is particularlyrelevant when seeking regions in unsteady flows which remain ‘coherent’ over a given time period[7]. In all these cases, we emphasize that we are not seeking to understand stretching in theinfinite-time limit—which is the focus in many classical approaches [3, 2]—but rather stretchingassociated with a one-step map derived from any of these approaches. From the applicationsperspective, the one-step map would be parametrized by the discrete or continuous time overwhich the map operates, and this number would of necessity be finite in any computationalimplementation.When additionally seeking global optimization , the first issue is defining what this meanswith respect to a bounded open domain on which the map operates. In Section 3, we posethis question as an optimization over foliations , but need to restrict these foliations in a cer-tain way because they would generically have singularities. We are able to characterize therestricted foliations of optimal stretching (minimal or maximal) in a straightforward geomet-ric way, while establishing connections to well-known local stretching optimizing entities. Weprovide a complete classification of the nondegenerate singularities using elementary argumentsin Section 4, thereby easily identifying 1- and 3-pronged singularities as the primary scenarios.We argue in Section 5 the inevitability of a ‘branch cut’ phenomenon if attempting to com-pute these restricted foliations using a vector field; this will generically possess discontinuitiesacross one-dimensional curves which we can characterize. Other computational ramificationsare addressed in Section 6, which includes issues of curves stopping abruptly when coming inhorizontally or vertically, and veering along spurious curves. We are able to give explicit insightsinto the emergence of these issues as a result of standard numerical implementations, and wesuggest an alternative integral-curve formulation which avoids these difficulties. In Section 7,we demonstrate computations of globally optimal restricted foliations for several well-knownexamples: the H´enon map [8], the Chirikov (standard) map [9], and the double-gyre flow [4],each implemented over a finite time. The aforementioned numerical issues are highlighted inthese examples. 3 Globally optimizing stretching
Let Ω be a bounded two-dimensional subset of R consisting of a finite union of connected opensets, each of whose closure has at most a finite number of boundary components. So Ω may, forexample, consist of disconnected open sets and/or entities which are topologically equivalent tothe interior of an annulus. We will use ( x, y ) to denote points in Ω. Let F be a map on Ω to R which is given componentwise by F ( x, y ) = (cid:18) u ( x, y ) v ( x, y ) (cid:19) . (1) Hypothesis 1 (Smoothness of F ) . Let the map F ∈ C (Ω) . Physically, we note that F can be generated in various ways. It can be simply one iterationof a given map, multiple (finitely-many) iterations of a map, or even the application of a finite sequence of maps. It can also be the flow-map generated from a nonautonomous flow in two-dimensions over a finite time. In this sense, F encapsulates the fact that finiteness is inevitablein any numerical, experimental or observational situation, while allowing for both discrete andcontinuous time, as well as nonautonomy. The time over which the system operates can bethought of as a parameter which is encoded within F , and its effect can be investigated ifneeded by varying this parameter.The relative stretching of a tiny line (of length δ >
0) placed at a point ( x, y ) in Ω, with anorientation given by θ ∈ [ − π/ , π/
2) due to the action of F isΛ( x, y, θ ) = lim δ → (cid:107) F ( x + δ cos θ, y + δ sin θ ) − F ( x, y ) (cid:107) δ . This is the magnitude of F ’s directional derivative in the θ direction. It is clear thatΛ( x, y, θ ) := (cid:13)(cid:13)(cid:13)(cid:13) ∇ F ( x, y ) (cid:18) cos θ sin θ (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) = (cid:13)(cid:13)(cid:13)(cid:13)(cid:18) u x ( x, y ) u y ( x, y ) v x ( x, y ) v y ( x, y ) (cid:19) (cid:18) cos θ sin θ (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) . (2)We refer to Λ( x, y, θ ) in (2) as the local stretching associated with a point ( x, y ) ∈ Ω; note thatthis also depends on a choice of angle θ in which an infinitesimal line is to be positioned. Ifwe take the supremum over all θ ∈ [ − π/ , π/
2) of the right-hand side of (2), we would get theoperator (matrix) norm (cid:107) ∇ F (cid:107) , computable for example via Cauchy–Green tensor C ( x, y ) := [ ∇ F ( x, y )] (cid:62) ∇ F ( x, y ) . (3)Thus, our development has close relationships to well-established methods related to the Cauchy–Green tensor, finite-time Lyapunov exponents, and methods for determining Lagrangian coher-ent structures, which we describe in more detail in D. However, at this stage our local stretchingdefinition in (2) is θ -dependent. Definition 1 (Isotropic and remaining sets) . The isotropic set I ⊂ Ω is defined by I := (cid:26) ( x, y ) ∈ Ω : ∂ Λ( x, y, θ ) ∂θ = 0 (cid:27) , (4) and the remaining set is Ω := Ω \ I . (5)4he isotropic set I consists of points at which the local stretching does not depend ondirectionality of a local line segment. Given the smoothness we have assumed in F , I must bea ‘nice’ closed set; it cannot, for example, be fractal. In general, I may be empty, equal to Ω,or consist of a mixture of finitely many isolated points and closed regions of Ω.We are seeking a partition of Ω into a family of nonintersecting curves, such that globalstretching is optimized in a way to be made specific. Since the local stretching at points in I is impervious to the directionality of lines passing through them, these families of curves onlyneed be defined on Ω = Ω \ I , with the understanding that this has nonempty interior. In moreformal language, we need to think of singular codimension-1 foliations on Ω, whose singularitiesare restricted to I . We codify this in terms of the required geometric properties of the familyof curves: Definition 2 (Restricted foliation) . A restricted foliation , f , on Ω consists of a family of curvesdefined in the remaining set Ω such that(a) The curves of f (‘the leaves of the foliation’) are disjoint;(b) The union of all these curves covers Ω ;(c) The tangent vector varies in a C -smooth fashion along each curve. Our definition is consistent with the local properties expected from a formal definition offoliations on manifolds [10], but bears in mind that Ω is not a manifold because of the omissionof the closed set I from Ω. We remark that if I consists of a finite number of points, our restrictedfoliation definition is equivalent to that of a ‘punctured foliation’ [11] on Ω, where the puncturesare at the points in I . This turns out to be a generic expectation for I , and we will examinethis (both theoretically and numerically) in more detail later.The properties of Definition 2 ensure that every restricted foliation f is associated with aunique C -smooth angle field on the remaining set Ω in the following sense. Given a point( x, y ) ∈ Ω , there exists a unique curve from f which passes through it. The tangent line drawnat this point makes an angle θ f with the positive x -axis. This angle can always be chosenuniquely modulo π , from the set [ − π/ , π/ θ f = − π/
2, while horizontallines have θ = 0. Thus, every foliation induces a unique angle field θ f : Ω → [ − π/ , π/ π ). The angle field must be C -smooth to complement the continuous variation in thetangent spaces of f ’s leaves. Conversely, suppose a C -smooth angle field θ f : Ω → [ − π/ , π/ π ) is given. Given an arbitrary point ( x α , y α ) ∈ Ω , the existence of solutions to thedifferential equation (sin θ f ( x, y )) x. − (cos θ f ( x, y )) y. = 0passing through the point ( x α , y α ) ensures that there is an integral curve of the form g α ( x, y ) = 0,in which g α is C -smooth in both arguments. This is possible for each and every ( x α , y α ) ∈ Ω ,and uniqueness ensures that the curves g α ( x, y ) = 0 do not intersect one another. Moreover,Ω is spanned by (cid:83) α { ( x, y ) : g α ( x, y ) = 0 } because Ω = (cid:83) α { ( x α , y α ) } , ensuring that there isa curve passing through every point ( x α , y α ). Hence, this process generates a unique restrictedfoliation f on Ω .We are now in a position to define the global stretching which we seek to optimize. Definition 3 (Global stretching) . Given any restricted foliation f , we define the global stretch-ing on Ω as the local stretching integrated over Ω, i.e.,Σ f := (cid:90)(cid:90) Ω Λ ( x, y, θ f ( x, y )) x. y. + (cid:90)(cid:90) I Λ ( x, y, (cid:5) ) x. y. , (6)in which θ f is the angle field induced by a choice of restricted foliation f .5otice that the integral over the full domain Ω has been split into one over Ω (on which f and thus θ f is well-defined) and over I (over which the directionality has no influence on Λ,and has thus been omitted). Thus, any understanding of foliation leaves on I is irrelevant tothe global stretching, motivating our definition of restricted foliation defined only on Ω .As central premise of this work, we seek restricted foliations f which optimize (maximize,as well as minimize) Σ f . Partitions of Ω which are extremal in this way represent the greatestinstability or most stability associated with the dynamical system, and so orbits associatedwith these are distinguished for their corresponding difficulties in forecasting, or alternatively,relative coherence. Before we state the main theorems, some notation is needed. On Ω, wedefine the C -smooth functions φ ( x, y ) = u x ( x, y ) + v x ( x, y ) − u y ( x, y ) − v y ( x, y ) ψ ( x, y ) = u x ( x, y ) u y ( x, y ) + v x ( x, y ) v y ( x, y ) (8)in terms of the partial derivatives u x , u y , v x and v y of the mapping F . First, we show theconnection between zero level sets of φ and ψ and the isotropic set I . Lemma 1 (Isotropic set) . The isotropic set I defined in (4) can be equivalently characterizedby I := { ( x, y ) ∈ Ω : φ ( x, y ) = 0 and ψ ( x, y ) = 0 } , (9) Proof.
See A.We reiterate from this recharacterization of I that generically, it will consist of finitely manypoints (at which the curves φ ( x, y ) = 0 intersect the curves ψ ( x, y ) = 0), but may contain curvesegments (if the two curves are tangential in a region), or areas (if both φ and ψ are zero intwo-dimensional regions). Even for the generic case (finitely many isolated points), we will seethat I will strongly influence the nature of the optimal foliations in Ω .Next, we define the angle field θ + : Ω → [ − π/ , π/
2) by θ + ( x, y ) := 12 ˜tan − ( ψ ( x, y ) , φ ( x, y )) (mod π ) , (10)in terms of the four-quadrant inverse tangent function ˜tan − (˜ y, ˜ x ) (sometimes called atan2 incomputer science applications, which assigns the angle in [ − π, π ) associated with the quadrant in(˜ x, ˜ y )-space when computing tan − (˜ y/ ˜ x )). We also define the angle field θ − : Ω → [ − π/ , π/ θ − ( x, y ) = π − ( ψ ( x, y ) , φ ( x, y )) (mod π ) , (11)and observe that θ + ( x, y ) − θ − ( x, y ) = − π π ) . (12) Lemma 2 (Equivalent characterizations of angle fields, θ ± ) . On Ω , θ ± ∈ [ − π/ , π/ arerepresentable as θ + ( x, y ) := tan − − φ ( x, y ) + (cid:112) φ ( x, y ) + ψ ( x, y ) ψ ( x, y ) (mod π ) (13) and θ − ( x, y ) := tan − − φ ( x, y ) − (cid:112) φ ( x, y ) + ψ ( x, y ) ψ ( x, y ) (mod π ) . (14)6 roof. See B.
Remark 1 (Removable singularities at ψ = 0 and φ (cid:54) = 0) . While it appears that points where ψ = 0 but φ (cid:54) = 0 are not in the domain of θ + as written in (13) and (14), these turn out tobe removable singularities, and thus can be thought of in the sense of keeping φ constant andtaking the limit ψ →
0. More specifically, this implies that θ + ( x, y ) (cid:12)(cid:12)(cid:12) ψ =0 = (cid:26) − π/ φ <
00 if φ > . (15)With this understanding of dealing with the removable singularities, we will simply view (13)as being defined on Ω . Similarly, θ − ( x, y ) (cid:12)(cid:12)(cid:12) ψ =0 = (cid:26) φ < − π/ φ > . (16) Remark 2 (Smoothness of θ ± in Ω ) . Subject to the removable singularity understandingof Remark 1, θ + and θ − are C -smooth in Ω , and thereby respectively induce well-definedfoliations f + and f − on Ω .While theoretically, the alternative expressions in (13)-(14) for θ ± are equivalent to thedefinitions in (10)-(11), practically in fact, which of these is chosen will cause differences whenperforming numerical optimal foliation computations. We will highlight similarities and differ-ences between their usage in Section 6, and demonstrate these issues numerically in Section 7.We can now state our first main result: Theorem 1 (Stretching Optimizing Restricted Foliation - Maximum (SORF max )) . The re-stricted foliation f + which maximizes the global stretching (6) is that associated with the anglefield θ + . The corresponding maximum of the global stretching (6) is Σ + = (cid:90)(cid:90) Ω (cid:34) | ∇ u | + | ∇ v | (cid:112) φ + ψ (cid:35) / x. y. . (17) Proof.
See C.
Remark 3 (Lyapunov exponent field) . The integand of (17) is the Λ field associated withmaximizing stretching, and is given byΛ + ( x, y ) = (cid:34) | ∇ u | + | ∇ v | (cid:112) φ + ψ (cid:35) / = (cid:107) ∇ F ( x, y ) (cid:107) . (18)This is (a scaled version of) the standard Lyapunov exponent field . We avoid a time-scaling heresince, for example, F may be derived from a sequence of application of various forms of maps(indeed, any sequential combination of discrete maps and continuous flows). Neither will wetake a logarithm, since we do not necessarily want to think of the stretching field as an exponentbecause the finite ‘amount of time’ associated with F depends on its discrete/continuous nature,which is flexible in our implementation. 7 emark 4 (Stretching on the isotropic set I ) . The value of the global stretching restricted to I (i.e., the second integral in (6)) is, from (17), (cid:90)(cid:90) I (cid:34) | ∇ u | + | ∇ v | (cid:35) / x. y. = 12 (cid:90)(cid:90) I (cid:107) ∇ F (cid:107) Frob x. y.= 12 (cid:90)(cid:90) I (cid:110) Tr (cid:104) ∇ F ( ∇ F ) (cid:62) (cid:105)(cid:111) / x. y. , = 12 (cid:90)(cid:90) I { Tr [ C ( x, y )] } / x. y. , (19)expressed in terms of the Frobenius norm (cid:107) (cid:5) (cid:107) Frob or trace Tr [ (cid:5) ] of the Cauchy–Green tensor (3).Similar to the maximizing result, we also have the minimal foliation:
Theorem 2 (Stretching Optimizing Restricted Foliation - Minimum (SORF min )) . The re-stricted foliation f − which minimizes the global stretching (6) is that associated with the anglefield θ − . The corresponding minimum of the global stretching (6) is Σ − = (cid:90)(cid:90) Ω (cid:34) | ∇ u | + | ∇ v | − (cid:112) φ + ψ (cid:35) / x. y. . (20) Proof.
See C.
Corollary 1 (SORF max and SORF min are orthogonal) . If any curve from SORF max intersectsa curve from SORF min in Ω , then it does so orthogonally.Proof. The SORF max and SORF min curves are respectively tangential to the angle fields θ + and θ − , which are known to be orthogonal by (12).There is clearly a strong interaction between local properties and quantities related to global stretching optimization. We summarize some properties below. We do not discuss them indetail, but provide additional explanations in D. Remark 5 (Maximal or minimal local stretching) . (a) Given a point ( x, y ) ∈ Ω , if we pose the question of determining the orientation of aninfinitesimal line positioned here in order to experience the maximum stretching, then thisis at an angle θ + .(b) The local maximal stretching associated with choosing the angle of orientation θ + isexactly the operator norm of the gradient of the map F , which is expressible in terms ofthe Cauchy–Green tensor (3).(c) The above quantity is associated with the Lyapunov exponent field , given in (18), whichis defined on all of Ω despite having the above interpretation only on Ω .(d) In Ω , the SORF max leaves (curves) lie along streamlines of the eigenvector field of theCauchy–Green tensor corresponding to the larger eigenvalue. This eigenvector field canalso be thought of as the Lyapunov or Oseledec vector field associated with F .(e) If the question is instead to find the orientation of an infinitesimal line positioned at ( x, y )in order to experience the minimum stretching, then the angle of this line is θ − . Comparethis statement to observation (a) together with Corollary 1.(f) In Ω , Eq. (5), the SORF min leaves lie along streamlines of the eigenvector field of theCauchy–Green tensor corresponding to the smaller eigenvalue.(g) The set I corresponds to points in Ω at which the two eigenvalues of the Cauchy–Greentensor coincide. 8igure 1: Topological classification of nondegenerate singularities with respect to SORF max or-min (a) a 1-pronged (intruding) point, and (b) a 3-pronged (separating) point. See Property1. Compare to Fig. 13. The previous section’s optimization ignored the isotropic set I , since the local stretching within I was independent of direction. In this section, we analyze the topological structure of ouroptimal foliations near generic points in I , which can be thought of as singularities with respectto optimal foliations. By Lemma 1, these are points where both φ and ψ are zero. Definition 4 (Nondegenerate singularity) . If a point p ∈ I is such that (cid:104) ∇ φ × ∇ ψ (cid:105) p (cid:54) = or equivalently det ∂ ( φ, ψ ) ∂ ( x, y ) (cid:12)(cid:12)(cid:12) p (cid:54) = 0 , (21) then p is a nondegenerate singularity . Since by Hypothesis 1 both φ and ψ are C -smooth in Ω, their gradients are well-defined onΩ. Nondegeneracy precludes either φ or ψ possessing critical points at p ; thus, we cannot getself-intersections of either φ = 0 or ψ = 0 contours at p , have local extrema of φ or ψ at p , orhave a situation where φ or ψ is constant in an open neighborhood around p . Nondegeneracyalso precludes φ = 0 and ψ = 0 contours intersecting tangentially at p (although we will be ableto make some remarks about this situation later). Thus, at nondegenerate points p , the curves φ = 0 and ψ = 0 intersect transversely. We explain in E how we obtain the following complete classification for nondegenerate singularities, as illustrated in Fig. 1: Property 1 (1- and 3-pronged singularities) . Let p ∈ I be a nondegenerate singularity, andlet ˆ k be the unit-vector in the + z -direction (i.e., ‘pointing out of the page’ for a standard right-handed Cartesian system). Then, • If p is right-handed , i.e., if (cid:104) ∇ φ × ∇ ψ (cid:105) p · ˆ k = det ∂ ( φ, ψ ) ∂ ( x, y ) (cid:12)(cid:12)(cid:12) p > , (22) then p is a (an ‘intruding point’), with nearby foliation of both f + and f − topologically equivalent to Fig. 1(a); and • If p is left-handed , i.e., if (cid:104) ∇ φ × ∇ ψ (cid:105) p · ˆ k = det ∂ ( φ, ψ ) ∂ ( x, y ) (cid:12)(cid:12)(cid:12) p < , (23) then p is a (a ‘separating point’), with nearby foliation of both f + and f − topologically equivalent to Fig. 1(b). max near p when transversality is relaxed (see Efor explanations of these structures). The intrusions/separations occur in opposite directions for the two orthogonal foliations f ± . We use the ‘1-pronged’ and ‘3-pronged’ terminology from the theory of singularities of mea-sured foliations [12, 13]. We also note that in the case of all singularities being nondegenerate,the curves on Ω may be thought of as a punctured foliation [11, e.g.] on Ω. These two singu-larities also correspond to the index of the foliation being +1 / − / f − is similar to that of f + as illustrated inFig. 1. To see why this is so, imagine reflecting these curves about the vertical line going through p . This generates an orthogonal set of curves, which are the complementary (orthogonal)foliation. Thus, f + and f − have the same topology near p .At the next-order of degeneracy, we will have φ = 0 and ψ = 0 contours continuing tobe curves, but now intersecting at p nontangentially . In that case, it turns out that Fig. 2gives the possible topologies for SORF max , which are explained in detail in E. If p is not anisolated point in I , then many other possibilities exist. The SORF min in the mildly degeneratesituations in Fig. 2 represent curves which are orthogonal to the pictured ones, by Corollary 1.Their topology will be identical. We have determined slope fields θ + and θ − corresponding to maximizing and minimizing the global stretching. By Remark 5, maximizing the local stretching at a point in Ω also results inan angle corresponding to θ + . Such local stretching is well-studied; it is related to the Lyapunovexponent, and the directions are associated with Lyapunov vectors [5] or
Oseledec spaces [3].Additionally, the direction associated with θ + can be characterized in terms of the eigenvectorassociated with the larger Cauchy–Green eigenvalue. See D for a more extensive discussion ofthese connections.Here, we analyze the vector fields associated with θ ± in some detail, using the behavior inthe ( φ, ψ )-plane introduced in the previous section. The main observation is that, generically,it is not possible to express a C -vector field on the closure of Ω from the θ ± angle fields.This has implications in numerically computing curves in the optimal foliations, where we giveinsight into spurious effects that arise. 10he θ + field in Ω is given by (10). To determine a curve from the SORF max , we need topick an initial point in Ω , and evolve it according to ‘the’ vector field generated from θ + . Asimple possibility would be to take the (unit) vector field w + ( x, y ) := (cid:18) cos [ θ + ( x, y )]sin [ θ + ( x, y )] (cid:19) , (24)in which θ + is computed from (10). In evolving trajectories associated with this vector field—i.e., in determining streamlines of (10)—one can of course multiply w + by a scalar function m ( x, y ), which simply changes the parametrization along the trajectory/streamline. As verifiedin D, (24) is indeed the eigenvector associated with the larger eigenvalue of the Cauchy–Greentensor at ( x, y ), with the understanding that it can be multiplied by a nonzero scalar. The factthat the eigenvector at each point is unique, modulo a constant multiple, is of course directlyrelated to these observations.Exactly the same arguments hold when attempting to compute the SORF min : from theangle field θ − we can construct the vector field w − ( x, y ) := (cid:18) cos [ θ − ( x, y )]sin [ θ − ( x, y )] (cid:19) , (25)where θ − is defined from (11). Property 2 (Generating foliation curves using vector fields) . If generating a SORF max orSORF min curve in Ω , we can in general find solutions to dds xy = w ( x ( s ) , y ( s )) ; x (0) y (0) = x y , (26) where s is the parameter along the curve and ( x , y ) ∈ Ω , and we can choose a Lyapunovvector field in the form w ( x, y ) = m ( x, y ) w ± ( x, y ) (27) for a suitable scalar function m . If we use m ≡ , the parametrization s along the trajectory is exactly the arclength.However, more general scalar functions m can be used in (26), reflecting the fact that the vectorfields which generate the foliations are actually direction fields , and thus can be multipliedat each point by a scalar. The only restrictions are (i) m can never be zero, because if itis, we introduce a spurious fixed point in the system (26) which ‘stops’ the curve, and (ii) m is sufficiently smooth to ensure that the equation (26) has unique C -smooth solutions.From the perspective of a SORF curve, making a choice of the function m simply adjusts theparametrization along the curve. Notice that if we flip the sign of m we would be going alongthe curve in the opposite direction.To understand the generation of curves from (27), it helps to think of the mapping from Ωto ( φ, ψ )-space, illustrated in Fig. 3. We have already characterized an important subset of Ωin relation to this mapping: the isotropic set I is the kernel of this mapping (by Lemma 1). Itsimage is denoted by I (cid:48) , the origin in ( φ, ψ )-space.Another important set that we require is Definition 5 (Branch cut) . The branch cut B is the set of points ( x, y ) ∈ Ω such that B := { ( x, y ) ∈ Ω : φ ( x, y ) < ψ ( x, y ) = 0 } . (28)11igure 3: The map from Ω to ( φ, ψ )-space, illustrating the sets I (cid:48) and B (cid:48) to which the sets I and B map. In red, we have stated the value of the field θ + in (10) in each quadrant.Figure 4: Vector field of (26) using w + , near a nondegenerate singularity p , with the branchcut B shown in green: (a) right-handed p and (b) left-handed p .The image B (cid:48) of the branch cut is also shown in Fig. 3 as the negative φ -axis. In each ofthe four quadrants of Fig. 3, we have carefully stated the value of the θ + field in terms of the standard inverse tangent function. We focus here near a nondegenerate singularity p , where the φ = 0 and ψ = 0 contours must cross p transversely, given that the Jacobian determinant of( φ, ψ ) with respect to ( x, y ) is nonzero. The axis-crossings in Fig. 3 will have the same topologyas these contours if the determinant is positive (the map is orientation-preserving).The relevant set B in Ω , near p , must therefore have the structure as seen in Fig. 4(a).Consider a small circle around p as drawn in Fig. 4(a), and indicated via arrows the directionsof the vector field w + along it. The reasons for these directions stems directly from Fig. 3; weneed to take the cosine (for the x -component) and the sine (for the y -component) of the anglefield defined therein. While w + must vary smoothly along the circle, it exhibits a discontinuityacross the branch cut B , because the angle has rotated around from − π/ π/
2. Clearly,the same behavior occurs for left-handed p : in this case we need to consider Fig. 3 with the ψ -axis flipped (this orientation-reversing case is indeed pictured in Fig. 13(b)). Once again, itis the φ − axis to which the branch cut B ∈ Ω gets mapped. The intuition of Fig. 4 gives us a theoretical issue related to using a vector field to find curves: Theorem 3 (Impossibility of continuous Lyapunov vector field) . If there exists at least onenondegenerate singularity p ∈ Ω , then no nontrivial scalar function m in (26) exists such thatthe right-hand side (i.e., vector field associated with the angle field θ + ) is a C -smooth nonzerovector field in Ω . The same conclusion holds for vector fields generated from θ − .Proof. See F. 12
Computational issues of finding foliations
In the previous section, we have outlined a theoretical concern in defining a vector field for com-puting optimal foliations. We show here related numerical issues which emerge when attemptingto compute foliating curves.First, we remark that using a vector field to generate curves of streamlines of eigenvectorfields of a tensor—which as seen here are equivalent to SORF max and SORF min curves—isstandard practice. Numerical issues in doing so have been observed previously, and ad hoc remedies proposed: • In generating trajectories following ‘smooth’ fields from grid-based data, one suggestedapproach is to keep checking the direction of the vector field within each cell a trajectoryventures into, and then flip the vector field at the bounding gridpoints to all be in thesame direction before interpolating [16]. • In dealing with points at which the eigenvector field is not defined, an approach is tomollify the field by multiplying with a sufficiently smooth field which is zero at suchpoints (e.g., the square of the difference in the two eigenvalues [15]).Our Theorem 3 gives explicit insights into the nature of both these issues. Both ad hoc numericalmethods relate to choosing the function m (respectively as ±
1, or a smooth scalar field which iszero at singularities). In either case, actual behavior near the singularities gets blurred by thisprocess.The branch cut near singularities also leads to more subtle—and apparently hitherto uniden-tified in the literature of following streamlines of tensor fields—issues when performing numericalcomputations. In G, we explain why the following occur.
Property 3 (Numerical computation of optimal foliations using vector fields) . Suppose wenumerically compute a SORF max (resp. SORF min ) curve by using (26) with m = 1 and thevector field w + (resp. w − ), by allowing the parameter s to evolve in both directions. Then(a) SORF max curves will not cross a one-dimensional part of B vertically, and may also veeralong B even though B may not be a genuine SORF max curve;(b) SORF min curves will not cross a one-dimensional part of B horizontally, and may alsoveer along B even though B may not be a genuine SORF min curve.These problems are akin to branch splitting issues arising when applying curve continuationmethods in instances such as bifurcations [17]. Is it possible to choose a function m which isnot identically 1 to remove these difficulties? The proof of Theorem 3 tells us that the answeris no. Either the branch cut gets moved to a different curve connected to p across which thereis a similar discontinuity, or it gets converted to a curve which has spurious fixed points (i.e., acenter manifold curve) because m = 0 on it. In either case, the numerical evaluation will giveproblems.Thus, there are several numerical issues in computing foliations using the vector fields w ± .Lemma 2 suggests a straightfoward alternative method for numerically computing such curvesin generic situations, while systematically avoiding all these issues. LetΦ − := { ( x, y ) : φ ( x, y ) < ψ ( x, y ) = 0 } andΦ + := { ( x, y ) : φ ( x, y ) > ψ ( x, y ) = 0 } ;13hese are points mapping to the ‘negative φ -axis’ and the ‘positive φ -axis’ (see Figs. 3 and 13),and we also note that Φ − = B . In seeking the maximizing foliation, we define on Ω \ Φ − , h + ( x, y ) = (cid:40) − φ ( x,y )+ √ φ ( x,y )+ ψ ( x,y ) ψ ( x,y ) if ψ ( x, y ) (cid:54) = 00 if ψ ( x, y ) = 0 and φ ( x, y ) > . (29)This is essentially the function tan θ + as defined in (13), and is C in Ω \ Φ − because ofRemark 1. The reason for not defining h + on Φ − is because the relevant tangent line becomesvertical. Hence we define its reciprocal, C on Ω \ Φ + , by h + ( x, y ) := (cid:40) φ ( x,y )+ √ φ ( x,y )+ ψ ( x,y ) ψ ( x,y ) if ψ ( x, y ) (cid:54) = 00 if ψ ( x, y ) = 0 and φ ( x, y ) < . (30)The minimizing foliation is associated with the angle field θ − . Thus we define on Ω \ Φ + , h − ( x, y ) := (cid:40) − φ ( x,y ) − √ φ ( x,y )+ ψ ( x,y ) ψ ( x,y ) if ψ ( x, y ) (cid:54) = 00 if ψ ( x, y ) = 0 and φ ( x, y ) < , (31)which gives the slope field associated with θ − , and on Ω \ Φ − its reciprocal h − ( x, y ) := (cid:40) φ ( x,y ) − √ φ ( x,y )+ ψ ( x,y ) ψ ( x,y ) if ψ ( x, y ) (cid:54) = 00 if ψ ( x, y ) = 0 and φ ( x, y ) > . (32) Property 4 (Foliations as integral curves) . Within Ω , a SORF max curve can be determinedby taking an initial point ( x , y ) and then numerically following dydx = h + ( x, y ) if (cid:12)(cid:12) h + ( x, y ) (cid:12)(cid:12) ≤ dxdy = h + ( x, y ) if else , (33)where we keep switching between the equations depending on the size of | h + | . This generatesa sequence ( x i , y i ) to numerically approximate an integral curve. Similarly, a SORF min curvecan be determined in Ω as integral curves of dydx = h − ( x, y ) if (cid:12)(cid:12) h − ( x, y ) (cid:12)(cid:12) ≤ dxdy = h − ( x, y ) if else . (34)Property 4 is an attractive alternative which avoids issues related to the branch cut andvector field discontinuities. Moreover, it is directly expressed in terms of the functions φ and ψ via the straightforward definitions of h ± and h ± . The switching between the dy/dx and dx/dy forms avoids the infinite slopes which may result if only one of these forms is used. Thus, wecan follow a particular curve as it meanders around Ω , having vertical and horizontal tangents,and also crossing branch cuts, with no problem. We will demonstrate applications of the theory to several maps F , generated from severalapplications of discrete maps, and from sampling flows driven by unsteady velocities. Theexamples include situations which are highly disordered (e.g., maps known to be chaotic underrepeated iterations, flows known to possess chaos over infinite times). Moreover, the maps F need not be area-preserving.In order to retain sufficient resolution to view relevant features in the many subfigures thatwe present in this Section, we will dispense with axes labels when these are self-evident: x willbe the horizontal axis and y the vertical as per standard convention.14 .1 H´enon map As our first example, consider the H´enon map, which is defined by [8] H ( x, y ) = (cid:18) − ax + ybx (cid:19) on Ω = R , and where we make the classical parameter choices a = 1 . b = 0 .
3. We choose F to be four iterations of the H´enon map, i.e., F = H . Fig. 5 demonstrates the computedfoliations and related graphs. The stretching field Λ + is first displayed in Fig. 5(a). In Fig. 5(b),we show the zero contours of φ and ψ . In this case, there are no nice transversalities. Indeed,there are several regions of almost tangencies, and the fact that several of the zero contoursalmost coincide in the two outer streaks in the figure, indicate that degenerate foliations areto be expected in their vicinity. The ‘squashing together’ that is occurring here is because weare at an intermediate stage in which initial conditions are gradually collapsing to the H´enonattractor. The vector fields w ± , computed using (24) and (25) and shown in Figs. 5(c,d)display discontinuities, which impact the computation of the SORF curves in (e) and (f). Theseare obtained by seeding 300 initial locations randomly in the domain, and then computingstreamlines generated from (26) with m = 1 in forward, as well as backward, s . Since the φ and ψ fields have large variations at small spatial scales because of the chaotic nature of the map,finding the branch cut B (where where ψ = 0 and φ <
0) as obtained from (28) requires care.We assess each gridpoint, and color it in (in green) if it has a different sign of ψ in comparisonto any of its four nearest neighbors, and the φ value at this point is negative. The lowermostpanel overlays the (green) set B on the SORF curves, indicating why some of the apparentbehavior in (e) and (f) is not representative of the true foliation; the center vertical line in (f),for example, occurs because of Property 3(b), while the SORF max (resp. SORF min ) curves stopabruptly on B if crossing vertically (resp. horizontally).On the other hand, Fig. 5(b) indicates that the zero contours of φ and ψ almost coincide ontwo curves: ‘outer’ and ‘inner’ parabolic shapes. These are also identified as part of the branchcut set B because ψ ≈ φ is slightly negative here. These curves are ‘almost’ a curveof I , and we see accumulation of SORF max curves towards these, indicating—at this level ofresolution—potential degeneracy of the foliation. We zoom in to this in Fig. 6. In conjunctionwith the explanations in Fig. 13, what occurs here is that the inner green line in Fig. 6(a) musthave a slope field which is − π/ − = B with respect to Fig. 6), while on the innerpink line it should be − π/ − in Fig. 13(a)). The extreme closeness of thecontours means that a very sharp change in direction must be achieved in a tiny region, whichthen visually appears as a form of degeneracy.This example highlights an important computational issue which is very general: eventhough relevant foliations will exist, in order to resolve them, one needs a spatial resolutionwhich can resolve the spatial changes in the φ and ψ fields. As an example of when F is generated from a finite-time flow, let us consider the flow mapfrom time t = 0 to 2 generated from the differential equation ddt xy = − πA sin [ πg ( x, t )] cos [ πy ] πA cos [ πg ( x, t )] sin [ πy ] ∂g∂x ( x, t ) , (35)in which g ( x, t ) := ε sin ( ωt ) x +[1 − ε sin ( ωt )] x and Ω = (0 , × (0 , F = H : (a) the logarithm of the maximumstretching field Λ + , (b) zero contours of φ and ψ , (c) vector field w + generated from (24),(d) vector field w − generated from (25), (e) SORF max by implementing vector field in (c),(f) SORF min by implementing vector field in (d), (g) SORF max with branch cut (green), (h)SORF min with branch cut (green). 16igure 6: Zooming in to an area associated with the map F = H (a) the zero contours of φ and ψ , (b) the SORF max , and (c) the SORF min . A = 1, ω = 2 π and ε = 0 .
1, and the optimal reduced foliations are demonstrate in Fig. 7.Fig. 7(a) is a classical figure in this context: the logarithm of the field Λ + ; if divided bythe time-of-flow 2, this is the finite-time Lyapunov exponent field. Fig. 7(b) indicates the φ = 0 and ψ = 0 contours, with their intersections defining I . We use the ‘standard’ w ± unitversions, Eq. (24), to generate the vector fields in (c) and (d), and the corresponding SORFs aredetermined in (e) and (f). Figs. 7(g) and (h) overlay the branch cuts (green), which are partsof the green curves in Fig. 7(b) at which φ <
0. As expected, the SORF max curves fail to crossthe branch cut vertically, as do the SORF min curves horizontally. Moreover, foliation curveswhich do get pushed in towards the branch cuts tend to meander along them, giving an impactof spurious accumulations. We zoom in towards one of these regions in Fig. 8; the SORF max curves requirements of having slopes − π/ π/
2) on Φ − (resp. Φ + ) result in abruptcurving. The accumulation is not exactly to Ψ − , but rather to a curve which is very close,as seen in Fig. 8(b). Thus, it is not true that there is a one-dimensional part of the isotropicset I along here. The geometric insights of the previous sections allows us to understand andinterpret these issues, while appreciating how resolution may give misleading visual cues.In Fig. 9, we zoom in to two difference locations, chosen by zeroeing in to two differentintersection points of the zero φ and ψ -contours. The top panels illustrate the SORF max (left)and the SORF min (right) curves at the same location. The theory related to 1-pronged intrudingpoints is well-demonstrated, with there being two such points adjacent to each other. The twoorthogonal families ‘reverse’ the locations of the singularities for the maximizing and minimizingfoliations, and the branch cut (green) forms vertical/horizontal barriers as appropriate. Incontrast, the bottom figures are of a 3-pronged separating point; again, the numerics validatethe theory. The Chirikov (also called ‘standard’) map is defined on the doubly-periodic domain Ω = [0 , π ) × [0 , π ) by [9] C k ( x, y ) = (cid:18) x + y + k sin x (mod 2 π ) y + k sin x (mod 2 π ) (cid:19) . We choose F = C nk , that is, n iterations of the Chirikov map for a given value of the param-eter k . Increasing k increases the disorder of the map, as does having n large. (The mapis a classical example of chaos, with Ω consisting of quasiperiodic islands in a chaotic sea,where ‘chaos/chaotic’ must be understood in the limit n → ∞ .) In more disorderly situations,increasingly fine resolution is required to reveal the structures that we have defined.Relevant computations for k = 2 and n = 4 are shown in Fig. 10. There are significantregions where the behavior is quite orderly. There is ‘greater disorder’ in the region foliated withlarge values of Λ + in (a)—indeed, this region is associated with the ‘chaotic sea’ when the map isiterated many more times—with the outer parts of low Λ + being associated with quasiperiodic17igure 7: Optimal foliation computations for the double-gyre flow: (a) The logarithm of thefield Λ + , (b) zero contours of φ and ψ , (c) vector field w + generated from (24), (d) vector field w − generated from (25), (e) SORF max by implementing vector field in (c), (f) SORF min byimplementing vector field in (d), (g) SORF max with branch cut (green), (g) SORF min withbranch cut. 18igure 8: Zooming in to near an ‘accumulating’ SORF max from Fig. 7: (a) the relevant zerocontours of φ and ψ , and (b) the SORF max .Figure 9: Zooming in to the SORF max (left) and SORF min (right) in the double-gyre. The topand bottom panels correspond to different locations, respectively near two adjacent intruding(1-pronged) points, and a separating (3-pronged) point. The branch cut is shown in green.Compare to Fig. 1 and Property 1. 19igure 10: Optimal foliation computations for the Chirikov map F = C : (a) the logarithm ofthe field Λ + , (b) zero contours of φ and ψ , (c) mboxSORF max with branch cut (green), (d)SORF min with branch cut (green).Figure 11: A degenerate singularity of the map F = C , shown zoomed-in: (a) the zero contoursof φ and ψ , (b) SORF max , and (c) SORF min .islands and hence order. All features mentioned in previous examples are reiterated in thepictures. Moreover, the SORF min foliation somewhat mirrors the structure expected fromclassical Poincar´e section numerics.If we instead consider k = 1 and n = 2, an interesting degenerate singularity (correspondingto the ψ = 0 contour crossing exactly a saddle point of φ ) is displayed in Fig. 11. The singularityin the S SORF max foliation (b) appears like a degenerate form of a separating point, if thinkingin terms of curves coming from above. However, if viewed in terms of curves coming in frombelow, it appears as an intruding point with a sharp (triangular) end. The SORF min conformsto this, having elements of a separating point, and an intruding point, as well. (The numericalissue of SORF min not crossing B horizontally is displayed in Fig. 11(c); in reality, the SORF min curves should connect smoothly across.)Next, we demonstrate in Fig. 12, using F = C , the efficacy of using the integral-curve forms(33) and (34) of the foliations, rather than using a vector field. The ln Λ + field in Fig. 12(a)has several sharp ridges; these are well captured by locations where the φ and ψ zero-contoursin Fig. 12(b) coincide. The SORF max/min foliations in (b) and (c) are computed respectively20sing the vector fields w ± as in previous situations, and exhibit the usual issues when crossing B . In contrast, the lower row is generated by using the integral-curve forms (33) and (34),where we have once again started from 300 random initial conditions. For each initial condition( x , y ), we define the next point ( x , y ) on a SORF max curve by x = x + h + ( x , y ) δy where δy > y -direction, and dx/dy is based on (33). Similarly, y = y + h + ( x , y ) δx using (33), and where δx > x -direction.This initializes the process. Next, we check the value of h + ( x , y ), thereby deciding which ofthe equations in (33) to implement. If the dy/dx equation, we take x = x + sign ( x − x ) δx ,and thus find y using the ODE solver. Having now obtained ( x , y ), we again use the last twopoints to make decisions on which of the two equations to use, and continue in this fashion for apredetermined number of steps. Next, we go back to ( x , y ) and now set x = x − h + ( x , y ) δy and y = y − h + ( x , y ) δy , thereby going in the opposite direction. Having initiated thisprocess, we can then continue this curve using the same continuation scheme. The SORF min are obtained similarly, using the two equations in (34). There is sensitivity in the process tolocations where φ and ψ change rapidly (they are each of the order 10 in this situation), and inparticular where zeros are near. The resolution scales δx and δy need to be reduced sufficientlyto not capture spurious effects. Notice that there are no branch-cut problems in the resultingfoliations obtained using the integral-curve approach, since we do not have to worry about adiscontinuity in a vector field. Neither are there any abrupt stopping of curves. In this paper, we have examined the issue of determining foliations which globally maximizeand minimize stretching associated with a two-dimensional map, where the map can be definedin terms of a finite sequence of discrete maps, or a finite-time flow of a differential equation.Our formulation establishes a connection to the well-known local optimizing issue, and providesnew insights into the resulting foliations and their singularities. In particular, an easy criterionfor classifying the nature of generic singularities is expressed. Some numerical artefacts arisingwhen computing these foliations in standard ways are characterized in terms of a ‘branch cut’phenomenon, and a methodology of avoiding these is developed. We have expressed connec-tions with a range of related and highly studied concepts (Cauchy–Green tensor, Lyapunovvectors, singularities of vector fields), and demonstrated computations in both discretely- andcontinuously-derived maps.We expect these results to help researchers interpret, and improve, numerical calculationsin related situations. In particular, misinterpretations of numerics can be mitigated via theunderstandings presented here. Regions of high sensitivity towards spatial resolutions are alsoidentifiable in terms of the near-zero sets of the φ and ψ functions.We wish to highlight from our numerical results the role of SORF min restricted foliations asbeing effective demarcators of complication flow regimes. These curves—observable for examplein blue in Figs. 5, 7, 10 and 12—indicate curves along which there is minimal stretching. Conse-quently, there is maximal stretching in the orthogonal direction to these curves. This indicatesthat the SORF min curves are barriers in some senses: disks of initial conditions positioned onsuch a curve experience sharp stretching orthogonal to them. That is, initial conditions onone side of such a curve get separated quickly from those on the other side, with the curvepositioned optimally to maximize the separation. Our methodology enables this intuitive ideato be put into a global optimizing foliation framework. Looking at this another way, the denseregions of the SORF min (blue) foliations in Figs. 5, 7, 10 and 12 are reminiscent of separationcurves which attempt to demarcate chaotic from regular regions. We emphasize, though, that‘chaotic’ has no proper meaning in the finite-time context since it must be understood in terms21igure 12: Comparison between using the integral-curve forms (33) and (34) and the vectorfield forms for F = C : (a) ln Λ + field, (b) zero contours of φ and ψ , (c) SORF max using thevector field (24), (d) SORF min using the vector field (25), (e) SORF max using the integral curveform (33), and (f) SORF min using the form (34).22f infinite-time limits; in this case, the separation one may try to obtain is between more ‘dis-orderly’ and ‘orderly’ regions. The ambiguity of defining these is reflected in the Figures, inwhich the SORF min foliation nonetheless identifies coherence-related topological structures inΩ which are strongly influenced by the nature of the singularities in the foliation.Note that the interaction of φ = 0 and ψ = 0 level sets as seen in Fig. 5(b) bear a strikingresemblance to Figures regarding zero angle between stable and unstable foliations of Lyapunovvectors such as in Fig. 1 for the H´enon map from [18] that was part of a search for primaryheteroclinic tangencies when developing symbolic dynamic generating partitions of the Henonmap, [19, 20, 21, 22]. Indeed this analysis likely bears a relationship, in that in a infinite timelimit, the Lyapunov vectors suggested come to the same point as those much earlier storiesunderlying the topological dynamics of smooth dynamical systems. What is clear in the finitetime discussion here is that when we see a coincidence between the stretching and folding, thatin successively longer time windows, these properties repeat in progressively smaller regions. Assuggested by Fig. 5, e.g. (h), any point of tangency would in turn be infinitely repeated in thelong time limit. The perspective of this current work may further understanding of what hasalways been the intricate topic of why and how hyperbolicity is lost in nonuniformly hyperbolicsystems wherein seemingly paradoxically, errors can grow along the directions related to stablemanifolds, such as highlighted by Fig. 5 in [23]. Acknowledgements:
SB acknowledges with thanks partial support from the AustralianResearch Council via grant DP200101764. EB acknowledges with thanks the Army ResearchOffice (N68164-EG) and also DARPA.
A Proof of Lemma 1
Given a general point ( x, y ) ∈ Ω , let θ ∈ [ − π/ , π/ x, y, θ ) = (cid:113) ( u x cos θ + u y sin θ ) + ( v x cos θ + v y sin θ ) . where the ( x, y )-dependence on u x , u y , v x and v y has been omitted from the right-hand side forbrevity. Hence,Λ = u x + v x − u y − v y θ + ( u x u y + v x v y ) sin 2 θ + u y + v y + u x + v x . Using the definitions for the functions φ and ψ from (7) and (8),Λ = φ cos 2 θ + ψ sin 2 θ + | ∇ u | + | ∇ v | . (36)Given the linear independence of the sine and cosine functions, the value of Λ at ( x, y ) isindependent of θ if and only if φ and ψ are both zero. Thus, the isotropic set is characterizedas the intersection of the zero sets of the functions φ and ψ . B Proof of Lemma 2
We begin with (10), and obtain (13). Assuming for now that both φ and ψ are not zero, we usethe double-angle formula to obtain 2 tan θ + − tan θ + = tan 2 θ + = ψφ . max near a nondegenerate singularity: (a) Value of θ + ∈ [ − π/ , π/
2) in ( φ, ψ )-space using (10), (b) as in (a), but shown in a left-hand system, (c) and (d) qualitative slopefields for (a) and (b); (e) 1-pronged ‘intruding point’ associated with the structure (c); (f) 3-pronged ‘separating’ point associated with the structure (d); (g) intruding point when axes aretilted; (h) separating point when axes are tilted. Compare to Fig. 1 and Property 1.24olving the quadratic for tan θ + , we see thattan θ + = − ± (cid:112) ( ψ/φ ) + 1 ψ/φ = − φ ± (cid:112) φ + ψ ψ (37)We now need to choose the sign in this expression, bearing in mind the usage of the four-quadrant inverse tangent as used in (10). The four quadrants here are in the ( φ, ψ )-space,which is indicated in Fig. 13(a). If φ > ψ >
0, this implies that 2 θ + is in the firstquadrant, and thus so is θ + . This means that tan θ + >
0, and consequently the positive signmust be chosen. If φ > ψ <
0, 2 θ + is in fourth quadrant, or 2 θ + ∈ ( − π/ , θ + <
0, and so the positive sign must be chosen in (37) to ensure that the division by ψ < φ < ψ >
0, 2 θ + ∈ ( π/ , π ), and θ + ∈ ( π/ , π/ θ + > φ < ψ <
0, 2 θ + ∈ ( − π, − π/
2) and θ + ∈ ( − π/ , − π/ θ + < θ + = − φ + (cid:112) φ + ψ ψ , whence (13) when neither φ nor ψ is zero.Next, we rationalize the fact that (13) arises from (10) even if one or the other of φ or ψ is zero. The arguments to follow are equivalent to considering the four emanating axes inFig. 13(a). If φ = 0 and ψ (cid:54) = 0, (10) tells us that 2 θ + = ( π/
2) sign ( ψ ) and thus tan θ + =tan( π/
4) sign ( ψ ) = sign ( ψ ). This is consistent with what (13) gives when φ = 0 is inserted. If ψ = 0 and φ (cid:54) = 0, (10), which tells us that 2 θ + = − π if φ <
0, or 2 θ + = 0 if φ >
0. Thus if ψ = 0, θ + = − π/ φ <
0, and θ + = 0 if φ >
0. This verifies that (13) is equivalent to (10) inΩ . Now, θ − in (11) is defined specifically to be orthogonal to θ + . There is only one angle in[ − π/ , π/
2) which obeys this condition. It is straightforward to verify from (13) and (14) that (cid:0) tan θ + (cid:1) (cid:0) tan θ − (cid:1) = − . Thus, θ − as defined in (14) is at right-angles to θ + as defined in (13), which has beenestablished to be equivalent to (10). C Proofs of Theorems 1 and 2
First, we tackle Theorem 1, related to maximizing the global stretching. Let f be a restrictedfoliation on Ω, and θ f be the unique angle field in Ω associated with it. From (36) from theproof of Lemma 1, we have that the local stretching Λ at a point ( x, y ) ∈ Ω related to theangle θ f obeysΛ = (cid:112) φ + ψ (cid:34) φ (cid:112) φ + ψ cos 2 θ f + ψ (cid:112) φ + ψ sin 2 θ f (cid:35) + | ∇ u | + | ∇ v | (cid:112) φ + ψ (cid:2) cos 2 θ + cos 2 θ + sin 2 θ + sin 2 θ f (cid:3) + | ∇ u | + | ∇ v | (cid:112) φ + ψ cos (cid:2) (cid:0) θ + − θ f (cid:1)(cid:3) + | ∇ u | + | ∇ v | θ + = θ + ( x, y ) satisfiescos 2 θ + = φ (cid:112) φ + ψ and sin 2 θ + = ψ (cid:112) φ + ψ . (39)25hus, tan 2 θ + = ψ/φ . If applying the inverse tangent to determine 2 θ + from this, we need totake the two equations (39) into account in choosing the correct branch. This clearly dependson the signs of φ and ψ , which is automatically dealt with if the four-quadrant inverse tangentis used. Consequently, (39) implies that θ + ( x, y ) = 12 ˜tan − ( ψ ( x, y ) , φ ( x, y )) , which is chosen modulo π because of the premultiplier of 1 / π ). Thus, θ + as defined here is identical to that given in (10), which by Lemma 2is equivalent to (13).Next, given that the cosine function is always between − (cid:34) − (cid:112) φ + ψ + | ∇ u | + | ∇ v | (cid:35) / ≤ Λ ≤ (cid:34)(cid:112) φ + ψ + | ∇ u | + | ∇ v | (cid:35) / , and consequently the global stretching (6) satisfiesΣ f ≥ (cid:90)(cid:90) Ω (cid:34) − (cid:112) φ + ψ + | ∇ u | + | ∇ v | (cid:35) / x. y. and (40)Σ f ≤ (cid:90)(cid:90) Ω (cid:34)(cid:112) φ + ψ + | ∇ u | + | ∇ v | (cid:35) / x. y. . (41)for any choice of foliation.Let f + be the foliation identified with the angle field θ + ( x, y ) at every location in Ω .Inserting this into (38) renders the cosine term 1, and thus the right-hand side of (41) is achievedfor this foliation. There can be no foliation with gives a larger value of Σ f . This foliation isequivalent to pointwise maximizing Λ in Ω .Can there be a different acceptable foliation, ˜ f , which also attains this maximum value forΣ f (i.e., that Σ ˜ f = Σ f + )? If so, there must be a point (˜ x, ˜ y ) ∈ Ω where the induced slopes θ ˜ f and θ + of the two different foliations are different. Given that foliations must be smooth,this implies the presence of an open neighborhood N ε (with positive measure) around this pointsuch that cos 2 (cid:16) θ ˜ f − θ + (cid:17) < − ε , for any given ε >
0. Thus the integrated local stretching in N ε for ˜ f is strictly less than that of f + . Since it is not possible to obtain a greater integratedstretching outside of N ε (because f + , by forcing the cosine term to take its maximum possiblevalue, cannot be bettered), this would imply that the integrated stretching of ˜ f over Ω isstrictly less than that of f + . Given that the contribution to the integral in I is independentof the foliation, this provides a contradiction. Therefore, the foliation f + , corresponding to thechoice of angle field θ + as given in (10), maximizes Σ f , and is uniquely defined in Ω .The proof of Theorem 2 related to minimizing the global stretching is similar. We use (40),which corresponds to choosing θ f such that the term cos 2 ( θ + − θ f ) is always −
1. This tellsus that θ f must be chosen perpendicular to θ + . Thiis is exactly the characterization used todetermine θ − in (11), and the equivalence to (14) has been established in Lemma 2. D Local stretching connections related to Remark 5
Given a location ( x, y ), suppose we wanted to determine the direction (encoded by an angle θ )to place an infinitesimal line segment such that it stretches the most under F . From (2), we26eed to solve sup θ (cid:13)(cid:13)(cid:13)(cid:13) ∇ F ( x, y ) (cid:18) cos θ sin θ (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) := (cid:107) ∇ F (cid:107) , where the right-hand side is the operator norm of ∇ F . This is computable by the square-rootof the larger eigenvalue of [ ∇ F ] (cid:62) ∇ F , i.e., of the Cauchy–Green tensor C as defined in (3).Given the map (1), since ∇ F = (cid:18) u x u y v x v y (cid:19) , it is clear that the Cauchy–Green strain tensor (as defined in (3)) is C := [ ∇ F ] (cid:62) ∇ F = (cid:18) u x + v x u x u y + v x v y u x u y + v x v y u y + v y (cid:19) . Accordingly, the eigenvectors λ of the Cauchy–Green tensor (i.e., the singular values of ∇ F )obey λ − (cid:16) | ∇ u | + | ∇ v | (cid:17) λ + (cid:104) ( u x + v x )( u y + v y ) − ( u x u y + v x v y ) (cid:105) = 0 , and thus λ = | ∇ u | + | ∇ v | ± (cid:112) φ + ψ (42)by using the definitions for φ and ψ in (7) and (8).We assume that φ and ψ are not simultaneously 0 (in our framework, that we are not in I ).Clearly, the larger value of λ is obtained by taking the positive sign, and the square-root of thisis the matrix norm of C . This gives exactly the pointwise maximized local stretching of Λ asdefined in (38), which satisfies (cid:0) Λ + (cid:1) = | ∇ u | + | ∇ v | (cid:112) φ + ψ . The quantity Λ + defined above (and also given in the main text as (18)), is related to the finite-time Lyapunov exponent or simply the Lyapunov exponent . We note that for defining θ + (for optimizing stretching) we required that ( x, y ) (cid:54) = I , but Λ + can be thought of as a field onall of Ω. Obtaining the eigenvector of the Cauchy–Green tensor C corresponding to the λ associ-ated with (18) is somewhat unpleasant. However, our equation for θ + in (13) indicates thateigenvector—modulo a nonzero scaling—can be written as˜ w + = (cid:18) ψ − φ + (cid:112) φ + ψ (cid:19) , as long as this value is not zero (which is when ψ = 0 and φ >
0, in which case w + = (1 0) (cid:62) ).Tedious calculations reveal that C ˜ w + = (cid:18) u x + v x ψψ u y + v y (cid:19) (cid:18) ψ − φ + (cid:112) φ + ψ (cid:19) = . . . = (cid:0) Λ + (cid:1) ˜ w + , verifying that our expression does indeed give the relevant eigenvector. The situation of ψ = 0and φ > w + = (1 0) (cid:62) , we once again get C ˜ w + = (cid:32) φ + | ∇ u | + | ∇ v | (cid:33) ˜ w + = (cid:0) Λ + (cid:1) ˜ w + . w + of C (or a scalar multiple of it) is only defined on Ω . In the literature,this is variously referred to as the Lyapunov [5, 6] or
Oseledec [3] vector field, related to thelocal direction (in the domain of F ) in which the stretching due to the application of F will bethe most. If F were a flow map derived from a flow over a finite-time, then these would dependboth on the initial time t and a time t at the end. In other words, F would be the flow mapfrom time t to t . In this situation, the variation of the vector field with respect to both t and t is to be noted.The smaller eigenvalue of the Cauchy–Green tensor is obtained by taking the negative signin (42), which gives (cid:0) Λ − (cid:1) = | ∇ u | + | ∇ v | − (cid:112) φ + ψ . This is clearly the local stretching minimizing choice, corresponding to choosing θ = θ − (i.e.,making the cosine term equal to − w − can be verified (asabove) to be in the direction specified by θ − . However, given that (cid:112) φ + ψ (cid:54) = 0, we havedistinct eigenvalues for the symmetric matrix C , and thus the two eigenvectors must be orthog-onal by standard spectral theory. Hence we can easily conclude that θ − corresponds to ˜ w − , theeigenvector of C corresponding to the smaller eigenvalue.The situation in which the eigenvalues of C coincide corresponds to ‘singularities,’ in par-ticular because this means that an orthogonal eigenbasis may not exist. This can only occurwhen the eigenvalues are repeated, and from (42) this occurs only when φ + ψ = 0. Thus, both φ and ψ must be zero. Thus corresponds exactly to the isotropic set I , in Definition 1 andLemma 1.We note that Haller [24] uses streamlines of the eigenvector fields from the Cauchy–Greentensor in his theories of variational Lagrangian coherent structures, looking for example forcurves to which there is extremal attraction or repulsion due to a flow over a given time period.Our foliations obtained here, corresponding to globally maximizing and minimizing stretching,are generated from fibers of the same fields. Therefore, our insights into singularities andbranch-cut discontinuities are therefore relevant to these approaches as well. E Singularity classification
This section provides explanations for the nondegenerate singularity classification of Property 1.Given the transverse intersection of the φ = 0 and ψ = 0 contours at a singularity p , we examinenearby contours not in standard ( x, y )-space, but in ( φ, ψ )-space, in which p is at the origin. Theangle fields θ ± are the defining characteristics of the foliation, and thus we show in Fig. 13(a)a schematic of the maximizing angle field θ + . A nonstandard labelling of the φ and ψ axes isused here because the relative orientations of the positive axes φ + and ψ + (the directions inwhich φ > ψ > φ − and ψ − is related to whether p is right-or left-handed. Thus, Fig. 13(a) corresponds to p being right-handed. The slope fields andexpressions indicated are based on the four-quadrant inverse tangent (10), expressed in termsof the regular inverse tangent in each quadrant. We also express the values of θ + on each of theaxes in Figs. 13(a), along which θ + is seen to be constant.In Figs. 13(c), just below, we indicate the angle field θ + by drawing tiny lines which havethe relevant slope. What happens when we ‘connect these lines’ to form a foliation is shownunderneath in Figs. 13(e). The foliation bends around the origin (shown as the blue point p ), effectively rotating around it by π . However, it must be cautioned that while Fig. 13(e)seems to indicate that the fracture ray lies along φ + , this is in general not the case. The anglefields shown in Figs. 13(c) and (e) display directions in physical (Ω) space, in which the φ = 0and ψ = 0 contours intersect in some slanted way. We show one possibility in Fig. 13(g), in28igure 14: SORF max near p when transversality is relaxed: (a), (b) and (c) show differentpossibilities for axes to intersect, and the corresponding SORF max topologies are illustrated inFig. 2.which the fracture ray will be approximately from the northwest. We identify p in this case an intruding point or a 1 -pronged singularity . The nearby SORF max curves rotate by π around it.In the right-hand panels of Fig. 13 we examine the other possibility of p being left-handed.This is achieved in Fig. 13(b) by simply flipping the ψ − and ψ + axes, and retaining the informa-tion that we have already determined in Fig. 13(a). The corresponding slope field is displayedin Fig. 13(d). The fracture ray (also along the φ + -axis in this case) now separates out curvescoming from the right, rather than causing them to turn around the origin. Fig. 13(f) demon-strates this behavior, obtained by connecting the angle fields into curves. There are two otherfracture rays generated by this process of separation, because curves in the φ − region are forcedto rotate away from the origin without approaching it. Fig. 13(h) is an orientation-preservingrotation of the axes in Fig. 13(f), which highlights that the directions of the three fracturerays are based on the orientations of the axes in physical space. Based on the topology of thefoliation, when p is left-handed, we thus have a separating point or 3 -pronged singularity .Suppose next that the nondegeneracy of p is relaxed mildly by allowing the φ = 0 and ψ = 0contours (both still considered to be one-dimensional) to intersect tangentially at p . To achievethis, imagine bending the ψ -axis in Figs. 13(a) and (c) so that it becomes tangential to the φ -axis, but the axes still cross each other. This degenerate situation is shown in Fig.14(a), andwe note that the orientation of the axes remains right-handed despite the tangency. Connectingthe angle field lines gives the relevant topological structure of Fig. 2(a). The topology is veryclose to the nondegenerate intruding point, but there is an accumulation of curves towards thefracture ray from one side. It is easy to verify (not shown) that there is no change in thistopology if the tangentiality shown in Fig. 14(a) goes in the other direction, with ψ + becomingtangential to φ + and ψ − to φ − . Fig. 14(b) examines the impact on the degenerate left-handedsituation; Fig. 2(b) indicates that the fracture ray acquires a similar one-sided accumulationeffect, while the remainder of the portrait remains essentially as it was. So this is a degenerateseparation point. Finally, in Fig. 14(c) we consider the case where the tangentiality is such thatthe φ - and ψ -axes do not cross one another. In this case, drawing connecting curves reveals thatthe topology is a combination of degenerate intruding and separating points, and is illustratedin Fig. 2(c). Testing the other possibilities (interchanging the ψ − and ψ + axes locations, anddoing the same analysis with them below the φ -axis) yields no new topologies. One way torationalize this is that the relative (degenerate) orientation between the negative axes and thatbetween the positive axes is in this case exactly opposite; one is as if there is a right-handedorientation, while the other is left-handed. 29 Proof of Theorem 3
We have established via Fig. 4 that if there exists a nondegenerate singularity p , then w + is notcontinuous across the branch cut B . This vector field is ‘the’ Lyapunov vector field, generatedfrom the eigenvector field corresponding to the larger eigenvalue of the Cauchy–Green tensorfield, where this is well-defined (i.e., in Ω ). However, a vector field associated with the anglefield θ + is not unique, as is reflected in the presence of the arbitrary function m in (26). Thenonuniqueness is equivalent to the potential of scaling Lyapunov vectors in a nonuniform wayin Ω \ I , by multiplying by a nonzero scalar. The question is: is it possible to remove thediscontinuity that w + has across B by choosing a scaling function m ?From Fig. 4, we argue that the answer is no. Imagine going around the black dashed curve, C , and attempting to have w + be continuous while doing so. Since w + has a jump discontinuityacross B , it will therefore be necessary to choose m to have the opposite jump discontinuity for m w + to be smooth. So m must jump from +1 to − w + is continuous on C \ B , to retain this continuity m must also remain continuous along C \ B . This implies that m must cross zero at some point in C \ B . Doing so would renderthe Lyapunov vector w + invalid. We have therefore established Theorem 3 using elementarygeometric means. We remark that this theorem is analogous to the classical “hairy ball” theoremdue to Poincar´e [25]. G Branch cut effects on computations If p is a nondegenerate singularity, then the vector field of (26) with m = 1 and the choice ofthe positive sign (SORF max ) will locally have the behavior as shown in Fig. 4. Now, in general,in finding a SORF max which passes through ( x , y ), we can implement (26) for the choice of m = 1, in both directions (increasing and decreasing s ), thereby obtaining the curve whichcrosses the point. An equivalent viewpoint is that we implement (26) with m = 1, and s > m = − s > m = +1 (globally) and w + to generate a SORF max curve, the vector fieldin Fig. 4(a) must be followed. However, it is clear that anything approaching the branch cut B gets pushed away in the vertical direction. Thus, SORF max curves near B will in general bedifficult to find.The solution appears to be to set m = −