Attraction-Based Computation of Hyperbolic Lagrangian Coherent Structures

Recent advances enable the simultaneous computation of both attracting and repelling families of Lagrangian Coherent Structures (LCS) at the same initial or final time of interest. Obtaining LCS positions at intermediate times, however, has been problematic, because either the repelling or the attracting family is unstable with respect to numerical advection in a given time direction. Here we develop a new approach to compute arbitrary positions of hyperbolic LCS in a numerically robust fashion. Our approach only involves the advection of attracting material surfaces, thereby providing accurate LCS tracking at low computational cost. We illustrate the advantages of this approach on a simple model and on a turbulent velocity data set.


Introduction
Hyperbolic Lagrangian Coherent Structures (LCS) in a two-dimensional unsteady flow are locally most repelling or most attracting material lines over a given finite time interval I = [t 1 , t 2 ] of interest [8]. Both mathematical methods and intuitive diagnostic tools have been developed to locate LCS in finite-time unsteady velocity fields with general time dependence (see [5] for a recent review.) Most computational approaches to LCS seek their initial or final positions as curves of initial conditions that lead to locally maximal trajectory separation in forward or backward time. This repulsion-based approach requires two numerical runs: one forward-time run that renders the time-t 1 position of forward-repelling LCS, and one backward-time run that reveals the time-t 2 position of forward-attracting LCS. Determining the positions of these material surfaces at an intermediate time t accurately, however, comes at high computational cost: it requires the accurate numerical advection of curves that are unstable in the time direction of advection (see Fig. 1.1, as well as the discussion in [3]).
A recent computational advance is offered by [3], showing how both repelling and attracting LCS can be simultaneously obtained either at t 1 or t 2 from a single numerical run. This approach renders an attracting LCS at a time t ∈ [t 1 , t 2 ] as the advected image of the initial LCS position at time t 1 . Similarly, the time-t position of a repelling LCS can be obtained by backward-advecting its position from time t 2 to t. Both of these computations track attracting material surfaces, and hence are numerically robust. However, they involve the advection of LCS from two different initial times, and hence are necessarily preceded by two separate numerical advections of a dense enough grid of initial conditions. Altogether, therefore, the Exact Numerical Figure 1.1. Forward advection of a classic stable manifold (a repelling LCS over finite times). Inaccuracies in determining the initial location of the LCS lead to exponentially growing errors and accumulation along the unstable manifold even if numerical errors were fully absent in the advection.
computational cost of constructing both repelling and attracting LCS at arbitrary times t ∈ [t 1 , t 2 ] has remained relatively high.
Here we propose a new computational strategy for two-dimensional incompressible flows. Our strategy builds on results from [3], [9] and [13,11], enabling the reconstruction of all hyperbolic LCS for arbitrary times t ∈ [t 1 , t 2 ] in a numerically robust fashion. This approach involves a single integration of trajectories from a full numerical grid, followed by the advection of select attracting material segments from local extrema of the singular value field of the flow gradient. This procedure yields substantial savings in computational time, as well as increased numerical accuracy in LCS detection and tracking. We demonstrate these advantages on a simple analytical flow example and on a direct numerical simulation of two-dimensional turbulence.
This paper is organized as follows. In Section 2, we fix our notation and recall relevant findings from [9] on the singular value decomposition of the linearized flow map. In Section 3 we present our attraction-based approach to hyperbolic LCS in the context of the recent geodesic theory of LCS [6,2,1]. In Section 4, we provide a proof of concept in the autonomous Duffing oscillator and compare our approach to previous ones in a simulation of two-dimensional turbulence, before concluding in Section 5.

Set-up
Consider a smooth, two-dimensional vector field v(x, t), defined over a finite interval I := [t 1 , t 2 ] and over spatial locations x ∈ D ⊂ R 2 . The trajectories generated by v(x, t) satisfy the ordinary differential equation The t 1 -based flow map is denoted by F t2 t1 : x 1 → x 2 , mapping initial values x 1 from time t 1 to their position at time t 2 along the corresponding solution of (2.1). We recall that the flow map is as smooth in x 1 as is v in x.
At any x 1 ∈ D, the deformation gradient DF t2 t1 (x 1 ) is a matrix that admits a singular value decomposition (SVD) of the form The numbers σ f 2 , σ f 1 are the singular values of DF t2 t1 ; the columns of Ξ (i.e., ξ 2 and ξ 1 ) are the right singular vectors of DF t2 t1 ; the columns of Θ (i.e., θ 2 and θ 1 ) are the left singular vectors of DF t2 t1 . From (2.2), we see that We recall that the singular values σ f 2 (x 1 ) and σ f 1 (x 1 ) measure infinitesimal stretching and compression along the trajectory starting from x 1 . Furthermore, the unit vectors ξ 2 (x 1 ) and ξ 1 (x 1 ) are the tangent vectors pointing to the directions of strongest stretching and compression under the linearized flow DF t2 t1 (x 1 ). If the velocity field is incompressible, i.e., ∇ x · v(x, t) ≡ 0, then det (DF ) = σ f 1 σ f 2 = 1, and consequently As a result, local maxima of σ f 2 (locally strongest-stretching points) coincide with local minima of σ f 1 (locally strongest-compressing points). At any point x 1 ∈ D, the average exponential rate of largest stretching over the time interval [t 1 , t 2 ] of length T = t 2 − t 1 is defined as , which is referred to as the (forward) finite-time Lyapunov exponent (FTLE). In the incompressible case, Eq. (2.3) shows that the FTLE can equally well be considered as a measure of the strongest local compression at x 1 .
For the backward flow from t 2 to t 1 , the backward deformation gradient is given by The singular values of DF t1 t2 are therefore given by , and the backward right singular vectors are given by θ 1 and θ 2 , the strongest-and weakest-stretching directions at x 2 in backward time. 1 Eq. (2.3) shows that the maximal (minimal) singular value of the linearized flow map is equal to the maximal (minimal) singular value of the linearized inverse flow map. Thus, local maxima of σ f 2 are mapped bijectively to local maxima of σ b 2 by the flow map. We summarize the equivalences of local extrema of the forward and backward singular value fields as follows: (2.6) In terms of the backward FTLE field, we recover [7, Prop. 2] for the incompressible case: In summary, as argued in [9], the SVD of DF t2 t1 yields complete forward and backward stretch information from a uni-directional flow computation.

Forward and backward Geodesic Theory of Hyperbolic LCS
The following definition recalls the hyperbolic LCS candidates obtained from two-dimensional geodesic LCS theory.
Definition (Shrink and stretch lines, [2,3]). We call a smooth curve γ a forward ( or backward) shrink line, if it is pointwise tangent to the ξ 1 (or θ 2 ) field. Similarly, we call γ a forward ( or backward) stretch line, if it is pointwise tangent to the ξ 2 (or θ 1 ) field.
Shrink and stretch lines are solutions of a variational principle put forward in [1] for LCS. This principle stipulates as a necessary condition that the time t 1 positions of hyperbolic LCS must be stationary curves of the averaged Lagrangian shear [1]. This variational principle leads to the result that time t 1 positions of hyperbolic LCS are necessarily null-geodesics of an appropriate Lorentzian metric associated with the deformation field [1]. This prompts us to refer to the underlying approach as geodesic LCS theory.
Away from points where σ f 2 = σ f 1 at t = t 1 and σ b 2 = σ b 1 at t = t 2 , both the initial and the final flow configuration is foliated continuously by mutually orthogonal forward and backward shrink and stretch lines. As discussed in [3,9], the following equivalence relations hold: The forward shrink and stretch lines provide candidate curves for the positions of repelling and attracting LCS at time t 1 . To find the positions of actual hyperbolic LCS as centerpieces of observed tracer deformation, we seek members of these two line families that evolve into locally most attracting or repelling material lines over the time interval [t 1 , t 2 ].
To this end, we follow [13,11] to require a sufficient condition that hyperbolic LCS must satisfy. Specifically, the time t 1 positions of forward repelling LCS are shrink lines that intersect local maxima of σ f 2 ; the time t 1 positions of forward attracting LCS are stretch lines that intersect local minima of σ f 1 . The time t 2 positions of backward-repelling and backward-attracting LCS are defined analogously using the backward singular-value fields σ b 2 and σ b 1 . By the equivalences detailed above, forward-attracting LCS, as evolving material lines, coincide with backward repelling LCS. Similarly, forward-repelling LCS, as evolving material lines, coincide with backward-attracting LCS.
If the vector field v(x, t) is incompressible, then the relation (2.3) forces local maxima of σ f 2 to coincide with local minima of σ f 1 . As a consequence, both forward-repelling and forward-attracting LCS intersect the maxima of σ f 2 at time t 1 . This fact will simplify our upcoming computational algorithm considerably for incompressible flows.
As noted earlier, reconstructing a full forward-attracting LCS as a material line involves advecting its time t 1 position under the flow map. This is a self-stabilizing numerical process, as it tracks an attracting surface. In contrast, reconstructing a forward-repelling LCS from its time t 1 position by flow advection is an unstable numerical process. Indeed, the smallest initial errors in identifying the LCS position are quickly amplified, as shown in Fig. 1.1.
Relations (2.6) and (3.1), however, allow us to compute the forward-repelling LCS equivalently as backward-attracting LCS. Specifically, forward-repelling LCS positions at a time t ∈ [t 1 , t 2 ] can be equivalently obtained from advection under the backward flow map F t t2 . The curves to be advected under F t t2 are just the backward stretch lines running through local minima of σ b 1 . By (2.5), however, local minima of σ b 1 are just the images of local maxima of σ f 2 under the flow map F t2 t1 . The computation of stretch lines still involves the integration of direction fields, for which orientation issues have to be resolved (see [14,2]). A new feature we introduce here is to advect short line segments (tangents) as opposed to whole stretch lines running through the appropriate extrema of the singular value fields. This idea exploits the tangentially stretching and normally attracting nature of stretch lines, saves on computational cost, and produces highly accurate results, as we demonstrate later. We summarize our attraction-based LCS construction in Fig. 3.1 for the case of incompressible flows. For compressible flows, forward-attracting LCS at time t 1 are still constructed from local minima of σ f 1 , but backward-attracting LCS at t 2 are constructed from advected local maxima of σ f 2 , which generally differ from advected local minima of σ f 1 .
Numerical implementation. Here we summarize the computational steps resulting from our previous considerations, assuming a forward-time advection of the chosen numerical grid.
(1) Compute flow map and its linearization: We solve the ODE (2.1) from a sufficiently dense grid of initial conditions to obtain a discrete approximation to the flow map F t2 t1 . We also obtain a numerical approximation to the linearized flow map DF t2 t1 at the grid points by one of four methods: (i) solving the equation of variations associated with (2.1), (ii) finite-differencing F t2 t1 along the grid, (iii) finite-differencing on a smaller auxiliary grid [2], (iv) via convolution with Gaussian kernels [12]. . Here A denotes a short vector parallel to ξ 2 at a local maximum x * 1 of the σ f 2 (x) field . Similarly, R denotes a short vector parallel to θ 1 at the point Recall that both forward-repelling and forward-attracting LCS intersect the maxima of σ f 2 at time t 1 in case of incompressibility.
values σ f i as well as the right-and left-singular vector fields ξ i and θ i , respectively. The singular values σ b i are obtained directly from the relation (2.5).
(3) Select seeding points for LCS: We need to identify points of strongest attraction, i.e. local minima of σ f 1 at the initial time and local minima of σ b 1 at the final time. While the first are identified directly, the latter are advected images of local maxima of σ f 2 under the flow map F t2 t1 . In the incompressible case, the points of strongest attraction coincide with local maxima of σ f 2 and their advected images under F t2 t1 , respectively. As in [11], we start by sorting all local maxima in ascending order by the values of σ f 1 or descending order by the values of σ f 2 . We then pick the first point p 1 in the ordered list and discard all local extrema in a small neighborhood of p 1 . From the remaining points on the list, we pick the first point p 2 and discard extrema in a small neighborhood of p 2 , and so on. This procedure filters out local extrema in noisy singular value fields. (4) Compute hyperbolic LCS: For any time t ∈ [t 1 , t 2 ] of interest, we use the flow map F t t1 to advect short line segments tangent to ξ 2 (p i ) at the points p i identified in the previous step. The resulting set of curves form the time t positions of attracting LCS. In the incompressible case, we use the flow map F t t2 to advect short line segments tangent to θ 1 (F t t2 (p i )) at the points F t2 t1 (p i ). Recall that the characteristic stretching directions for the backward flow are obtained from the forward time computation in step (2) due to Eq. (2.4). The resulting set of curves form the time t positions of repelling LCS. For the advection of line segments, the use of an adaptive integration scheme may be necessary. This is to fill emerging large gaps between adjacent points due to stretching, and to mitigate the possibly high curvature in the tracked material curve (see, e.g., [10] This example has already been used to illustrate shrink and stretch line context by [3] locally around the origin, showing the convergence of forward-and backward maximal stretch directions to the unstable and stable subspaces, respectively. In our present computations, we use the times t 1 = 0 and t 2 = T = 2.5. In Fig. 4.1, we compare the results from the earlier numerical LCS detection scheme used in [6] to our approach described in Section 3. While the left plot shows all structures to highlight the homoclinic loop, the middle plot shows that the shrink line deviates from the loop visibly at the first turn. In contrast, the backward-advected line segment stays close to the loop. The right plot shows that at the origin, both the shrink line and the advected stretch line indicate consistently the direction of strongest attraction. Fig. 4.2 gives further quantitative evidence that the backward-advected backward stretch line gives a better approximation to the actual repelling LCS position at time t 1 than the direct computation of this LCS position from forward shrink lines. Even in this simple example, therefore, the actual evolution of a shrink line and a backward-advected backward stretch line are noticeable different, although they should theoretically be identical. The root cause is numerical errors in the singular vector computation, as well as the limited ability of the discrete numerical grid to approximate a repelling LCS (local stable manifold) as a continuous curve. The error is initially invisible, but starts to accumulate rapidly during integration of the ξ 2 (θ 1 ) field and advection. 4.2. Two-dimensional turbulence. As a second example, we consider the twodimensional Navier-Stokes equations where the unsteady velocity field v(x, t) is defined on the two-dimensional domain U = [0, 2π] × [0, 2π] with doubly periodic boundary conditions. As in [3,4], we use a standard pseudo-spectral method with 512 modes in each direction, and 2/3 de-aliasing to solve the above Navier-Stokes equations with viscosity ν = 10 −5 on the time interval [0, 100]. The flow integration is then carried out over the interval t ∈ [50, 100], in which the turbulent flow has fully developed, by a fourthorder Runge-Kutta method with variable step-size. The initial condition v 0 is the instantaneous velocity field of a decaying turbulent flow. The external force f is random in phase and band-limited, acting on the wave-numbers 3.5 < k < 4.5.  Fig. 4.4). Analytically, these curves should coincide. In some locations, they indeed agree well, but in other locations, the discrepancy is dramatic (see the close-up view in the right panel of Fig. 4.4). This is the consequence of the effect illustrated in Fig. 1.1, showing the clear advantage of our method over the forward-time tracking of a repelling LCS.

Conclusions
We have proposed a paradigm shift in the detection of hyperbolic Lagrangian Coherent Structures (LCS). Instead of detecting initial positions of LCS as curves of maximal forward repulsion, we seek them as backward-advected locations of maximal backward attraction. While these two approaches are theoretically equivalent, the latter approach (developed here) eliminates an inherent numerical instability of the former approach (used in prior work). We have demonstrated that our attraction-based approach leads to substantial improvements in accuracy and computational cost.
We have discussed our approach in the framework of the geodesic theory [6,3,1], because this theory allows for the explicit computation of hyperbolic LCS as  parametrized curves. The proposed focus on attraction, however, automatically extends to potential future refinements in LCS computations.
The advection of identified hyperbolic LCS in the stable time direction is a simple idea, but relies heavily on the notion of a forward-time attracting LCS, which has been proposed only recently [3]. We have combined this notion with the SVD of the deformation gradient and with the seeding of straight line segments at points of locally strongest attraction to obtain a dynamically consistent and numerically robust approach to compute LCS. Extensions of these ideas to higher dimensions are possible and will be communicated elsewhere.