Edge manifold as a Lagrangian coherent structure in a high-dimensional state space

Dissipative dynamical systems characterized by two basins of attraction are found in many physical systems, notably in hydrodynamics where laminar and turbulent regimes can coexist. The state space of such systems is structured around a dividing manifold called the edge, which separates trajectories attracted by the laminar state from those reaching the turbulent state. We apply here concepts and tools from Lagrangian data analysis to investigate this edge manifold. This approach is carried out in the state space of autonomous arbitrarily high-dimensional dissipative systems, in which the edge manifold is reinterpreted as a Lagrangian coherent structure (LCS). Two different diagnostics, ﬁnite-time Lyapunov exponents and Lagrangian descriptors, are used and compared with respect to their ability to identify the edge and their scalability. Their properties are illustrated on several low-order models of subcritical transition of increasing dimension and complexity, as well on well-resolved simulations of the Navier-Stokes equations in the case of plane Couette ﬂow. They allow for a mapping of the global structure of both the state space and the edge manifold based on quantitative information. Both diagnostics can also be used to generate efﬁcient bisection algorithms to approach asymptotic edge states, which outperform classical edge tracking.


I. INTRODUCTION
Many deterministic physical systems can operate in two different regimes depending on initial conditions.This property calls naturally for a geometrical partition of the associated state space in terms of basins of attractions [1].The identification of the boundaries between basins is crucial for the cartography of the state space, for prediction, and for control.Besides such basin boundaries, being invariant sets of the system also supports their own specific dynamics.Some physical situations where the precise dynamics on the basin boundaries matter include climate dynamics [2], endothermic chemical reactions barriers [3], synchronization of phase oscillators [4], the instability of accretion disks [5], laser dynamics in modulated optical cavities [6], magnetic reconnection [7], free fall of objects in a gravity field [8], chaotic plasma devices [9], drift-wave turbulence [10], and many others.However the main illustration for the study of such a mixed-state space comes from hydrodynamics, more particularly the century-old problem of transition from laminar to turbulence [11].
In the field of hydrodynamics, most incompressible viscous fluid flows near solid boundaries can undergo transition to turbulence in a subcritical manner [12].The flow configurations of interest include the flow inside circular or rectangular pipes [13], between two plates in motion [14], inside a rotor-stator [15], or the flow developing on a semiinfinite flat plate [16].The laminar state in the examples above is time independent, characterized by spatial symmetries, and free of fluctuations, while the turbulent state is, depending on the level of modeling, arbitrarily complicated.For the parameters of interest here, the laminar state is linearly stable and coexists with the turbulent set.Some specific trajectories reach neither the laminar state nor the turbulent regime and have been labeled edge trajectories [17].Based on computational attempts to identify edge trajectories, the set of such trajectories has been conjectured to be a differentiable manifold of codimension 1 in the state space, called either edge of chaos, laminar-turbulent boundary, boundary of turbulence, or simply the edge [18][19][20].At sufficiently high Reynolds numbers, trajectories in the turbulent state are sustained: Both the turbulent and the laminar state are attractors of the system.As far as the Navier-Stokes equations are concerned, this bistable situation was first approached computationally in Refs.[21,22], in which the authors used a simple bisection method to identify the invariant dynamics on the edge, corresponding to a traveling wave solution (a result revised later).The analysis of this orbit and of its instability yields useful information for the understanding of the turbulent dynamics itself, notably regarding bursting dynamics, interpreted as a homoclinic connection between the traveling wave on the edge and a more complicated turbulent-like regime [23].
The bistable case is, however, only an ideal case, because the edge simply coincides with the intersection of the closures of the two basins of attraction.Following Ref. [24], we will refer to it as a strong edge.Ample experimental and numerical evidence has shown that for sufficiently low Reynolds numbers, the turbulent regime is only transient.The notion of basin boundary becomes more fragile and incompatible with the temporally asymptotic notions of stability: The laminar state becomes formally the only attractor of the system.The notion of basin boundary is, however, routinely generalized to such leaky cases by focusing on the values of the lifetimes associated with each initial condition, at least if the time of residence in the immediate vicinity of turbulent set is long enough compared to the typical timescales to be specified.We will refer in a such case to a weak edge.As a result, if the turbulent dynamics is transiently chaotic [25], the basin boundary can display a transversely fractal structure, as a consequence of the lemma [26].Edge trajectories remain smooth in all cases and stand out as trajectories with infinite lifetimes.Several algorithms have been constructed in order to approximate such accessible orbits [27], all based on a maximization of the lifetime [18,19,[28][29][30].In practice, the algorithms do not differ much from the standard shooting method.Most importantly, they rely on quantitative criteria globally measured in the state space to distinguish whether a given trajectory has, at a given time, safely entered the basin of attraction of one of the two states.These criteria have so far taken the form of numerical bounds on a given scalar observable.For instance, the kinetic energy of the departure from the laminar state is low when the laminar is approached and higher in the turbulent regime itself.Even bounds-based methods, for some flow cases like the Blasius boundary layer flow, have recently shown limitations for the tracking of edge trajectories beyond certain horizon times [31].Also, even in the general cases the identification of relevant bounds for a well-chosen observable implies a very good knowledge of the system that is not always available.Moreover, the outcome of a classical bisection process is a list of Booleans: Little information about neighboring trajectories is saved apart from the "laminar" or "turbulent" label at large enough times.Critically, the very notion of edge is defined so far both by the algorithmic way to identify the edge and by the infinite-time notions of stability found in the usual textbooks.A different paradigm, based on quantitative criteria, optimal quantities, and a finite-time framework, is hence necessary for a more transparent and universal characterization of edge manifolds.
Beyond the Boolean bounds-based definition of the edge manifold itself, other important properties can be quantified.The long-time outcome of classical bisection itself reveals a third asymptotic regime, unstable by construction but specific to the dynamics within that invariant manifold.This asymptotic regime, the edge state, is the relative attractor on the codimension one manifold.Depending on the flow case and the symmetries of the system (discrete and continuous), an equilibrium, a relative equilibrium (a traveling wave), a periodic orbit, a relative periodic orbit, a torus or even a chaotic attractor can be edge states.In each of these cases, the invariant manifold appears by definition as the set of initial conditions that converge in forward time to the edge state: defining the edge manifold as the stable manifold of the edge state.Complications do arise in some cases, seen in practice, where more than one distinct edge state exists for the same parameter values [20,32].The edge can then be globally interpreted as the closure of the union of such stable manifolds.
A few specific parts of the edge manifold are also important: By assuming a given distance defined in the entire state space, the minimal seed is defined loosely as a state on the edge, minimizing the distance to the laminar state [33,34].The notion of minimal seed (and the associated minimal distance) is useful in assessing the stability of the laminar state with respect to finite-amplitude disturbances.Recently, it was demonstrated that the edge state and the turbulent state often arise from the same saddle-node bifurcation followed by further bifurcations [35].Even the leaky property of the turbulent state, when present, emerges in boundary crises involving the collision between the turbulent state and the edge manifold as parameters are varied [36,37].This makes the identification of edge states instrumental for the determination of the full, exact bifurcation diagram of the system.
Interpreting the edge manifold as a stable manifold of some object to be found suggests bridges to other areas of dynamical systems.Many different tools have been developed and analyzed recently to identify the locally most repelling (resp., attracting) material surfaces.This forms the concept of Lagrangian coherent structure (LCS), where Lagrangian refers to the tracking of individual trajectories [38].LCSs correspond often in applications to stable (and unstable) manifolds of fixed points, as in Fig. 1(a).The original frame for which they have been developed is also inspired by hydrodynamics: It deals with the finite-time transport of fluid particles by known time-dependent fluid flows [40,41].Other applications, for instance, chemical reactions, have also been considered [3].Mathematically, the class of systems under study is radically different: The dynamics are always conservative, there are no attractors, the trajectories can be computed in either forward or backward time, and the dimension of the state space (which coincides with the physical space) is very low, typically two.Besides, the system is in general nonautonomous due to the time dependence of the vector field.In most applications of Lagrangian concepts, the state space associated with Lagrangian tracers (governed by an equation ẋ = v) and the physical space coincide.In the Eulerian point of view relevant to the transition problem, both spaces differ radically in their dimension.Despite these differences, the goal of this article is to demonstrate that the toolbox developed over the years for the study of LCSs can be used for the investigation and the numerical identification of the edge manifold in transitional flows according to the analogy expressed in Fig. 1.It is in line with the point of view used first by Ref. [42] for multidimensional dissipative systems.Dissipative PDEs formally have a state space of infinite dimension [in practice, their numerical discretization still yields a finite dimension of O (10 5 ) or more].The scalability property of the tools considered is hence crucial for the feasibility of the whole method.However, we restrain the analysis to autonomous velocity fields, this being the generic situation in all the hydrodynamic models considered here.We focus in particular Sketched attracting and repelling LCS (blue and red, respectively) identified with the unstable (resp., stable) manifold of a saddle point (green).(b) Sketch of the state space for a bistable hydrodynamical system, where the stable manifold of the edge state (green) is the basin boundary (red), whereas the unstable manifold (blue) leads to either the laminar or the turbulent attractor.The 3D figures represent flow fields from actual numerical simulations of plane Couette flow.Contours of streamwise velocity in the flow (red to blue) and an isovalue λ 2 criterion (green) used for vortex identification.mainly on two distinct tools popular in the recent literature on LCS: finite-time Lyapunov exponents (FTLEs) [38] and Lagrangian Descriptors (LDs) [43,44].We give key examples of applications of these tools in subcritical shear flow models of increasing complexity and dimensionality, following a hierarchy which ranges from a two-dimensional state space to infinite dimensions in the incompressible Navier-Stokes equations.Using the presented framework, we suggest and test the use of these diagnostics to improve edge tracking algorithms.

II. DEFINITIONS AND LCS DIAGNOSTICS
We consider here a general dynamical system defined on a forward-invariant subset A ⊂ R n , governed by where f : A × [t 0 , t 0 + τ ] → R n is a sufficiently smooth vector field and τ > 0. Let F t t 0 be the flow map which maps an initial position x 0 at time t 0 to its position at time t.When we linearize (1) around the trajectory, x(t ) ≡ x(t; t 0 , x 0 ) reads where z is the linearized solution and ∇ x f (x(t ), t ) is the Jacobian of the vector field f .The deformation gradient ∇F t t 0 corresponds to the fundamental solution matrix of the equation of variations so that z(t; t 0 , z 0 ) = ∇F t t 0 (x 0 )z 0 (4) in the interval t ∈ [t 0 , t 0 + τ ].

A. Finite-time Lyapunov exponents
Finite-time Lyapunov exponents are popular tools for the identification of LCS (cf.Ref. [38] and the references in Sec.4.2 therein).Let C t t 0 be the (positive definite) Cauchy-Green tensor where (•) * represents the transpose.Let λ i i = 1, ..., n denote the eigenvalues of C t t 0 and ξ i denote its associated eigenvectors such that and where ||•|| denotes the L 2 norm.The maximum of all factors by which small perturbations vectors are stretched over the time interval (t 0 , t 0 + τ ) around an initial condition x 0 is given by √ λ 1 (t 0 , t 0 + τ, x 0 ).The ith finite-time Lyapunov exponent of the system at position x 0 is given by for i = 1, . . ., n.The ridges in the field of the largest FTLE at time t 0 + τ (labeled simply λ t 0 +τ t 0 ) can be used as a diagnostic of hyperbolic LCS [38,45,46].Reference [45] provides rigorous theorems that establish a precise connection between hyperbolic LCS and FTLE ridges when further conditions on the rate-of-strain tensor are satisfied along the ridges.The Cauchy-Green tensor has size n × n, and its entries are usually computed using second-order centered finite differences [38].Hyperbolic LCSs refer to attracting and repelling distinguished invariant manifolds [40], with forward-time FTLEs relevant for the identification of repelling LCS and backwardtime FTLEs relevant for the identification of attracting LCS [47].
The computational cost of FTLEs following the above method increases rapidly with the state space dimension and is unavoidably costly for large dimensions.FTLEs can, however, be computed in arbitrarily high dimension using the recent algorithm based on optimally time-dependent (OTD) modes, which are determined from a minimization principle [48].Under generic conditions, these OTD modes converge exponentially fast to the eigendirections of the Cauchy-Green tensor associated with the largest eigenvalues, i.e., the dominant finite-time Lyapunov exponents [49].In practice for the present LCS diagnostic, only the largest exponent needs to be computed, for which other methods are available.

B. Lagrangian descriptors
Lagrangian descriptors (LDs) are a more recent diagnostic for Lagrangian coherence which does not require a differentiation of the flow map with respect to the initial condition.This diagnostic was introduced in Refs.[43,50] and further developed in Ref. [51].LDs are based on the integration of a given observable along trajectories.The original quantity of interest is In Eq. ( 8), the observable g is taken as where the f i 's are the components of the velocity field f , and p ∈ (0, 1] and τ ∈ R + are two parameters.We focus on the definition (8) used in the literature [52,53], although other alternative definitions have been suggested.It is convenient to split Eq. ( 8) into its forward and backwards contributions [3,52]: Due to the dissipative nature of Eq. ( 1), numerical backwards integration is ruled out for stability issues and only M + can be considered here.Since our focus is on stable manifolds rather than unstable ones, it is sufficient to focus on the computation of M + .LDs have been used to identify boundaries between qualitatively different dynamics, based on abrupt variations of M + .Alternatively, it is useful to quantify the abrupt changes of M + and thus to consider the Euclidian norm of its gradient Many other diagnostics for LCS have been suggested; see Ref. [46] for a recent comparative review.Alongside the diagnostics based on a scalar field, such as FTLEs and LDs, other approaches based on transfer operators or dynamic Laplace operators seek coherent structures by formulating rigorous mathematical coherence principles (see, e.g., Refs.[46,54]).These diagnostics generally display limited scalability properties and are not considered here.

III. LCS IDENTIFICATION OF THE EDGE
In this section, we demonstrate that the edge is highlighted as an LCS for several nonlinear models of increasing complexity, from two-dimensional models to thousands of degrees of freedom in the Navier-Stokes equations.We show that the edge can be effectively identified in state space even when its identification using global observables can be challenging.

A. Hierarchy of low-order shear flow models
The hierarchy of low-order models on which the LCS indicators are tested is based on a Galerkin truncation of the Navier-Stokes equations, in the spirit of the derivation of the Lorenz model [55].We assume that the bold vector X ∈ R 3 represents the position in the physical space coordinates (which differ from the state space coordinates) and that the velocity field v b is the stable base-flow solution, i.e., a steady solution of the governing PDEs.Let u = v − v b be the perturbation velocity to the base flow, not to be confused with the tangent field f .At any time, each model assumes that u can be written as with n the state space dimension and ûi , i = 1, .., n a basis of predetermined vector fields.The vector x(t ) ∈ R n contains all the amplitudes a i , i = 1, ..., n.We consider that for all models, Eq. ( 1) can be written under the generic form with x ∈ R n , L being an n × n linear operator, and N being a quadratic form containing the nonlinear terms.As explained in Ref. [56], all models of subcritical transition consistent with the original PDEs are subject to two constraints: (i) L is a nonnormal operator (LL * = L * L) with a stable eigenspectrum and (ii) the nonlinear terms do not contribute to the change in energy, i.e., N (x), x = 0 for all x, where •, • defines an inner product.

B. Dauchot-Manneville model (DM2D)
As an illustrative case, we consider the simple twodimensional model introduced by Dauchot and Manneville [57], henceforth referred to as DM2D.Although introduced originally in the context of hydrodynamic stability, its typical phase portrait (see Fig. 2) appears in many different applications, e.g., in chemistry for potential barriers [3], in ecology for the competition of two species [58], or in mechanics for the free fall of objects [8].There is no chaos in this model, only fixed points as attractors and one saddle point as edge state.The matrix L and the vector N terms read, respectively, . s 1 and s 2 are two negative parameters, so that the laminar state ) is linearly stable.The main control parameter is the discriminant = 1 − 4s 1 s 2 , for which x L is the only fixed point if < 0. Two additional fixed points x E and x T appear in a saddle-node bifurcation for 0, given by For 0 < 1, the model is bistable: It features two welldefined basins of attraction.They are separated by a smooth edge manifold = W s ({x E }) of the strong type.The saddle point x E is the edge state, whereas x L and x T are attractors interpreted as the laminar and turbulent state, respectively.Figure 2 shows a phase portrait of the model for the parameters s 1 = −0.1875and s 2 = −1.Two LCS diagnostics, FTLEs and LDs, are applied to this model with a time horizon fixed to τ = 60.For the FTLE field, a clear ridge is visible in Fig. 3  of Fig. 2.There is a nondifferentiable minimum of M + across , as in former applications of LDs to nondissipative systems.Isocontours of its gradient norm B happens to be convenient for plotting purposes [Fig.3(c)].B is computed using centered finite differences: where ε i is a small variation along each state space direction.B is shown in Fig. 3(c).The isolevels of B highlight the dip of M + map in Fig. 3(b) as a singular feature that coincides with .The variations of M + and B along the edge are quantified in Fig. 3(d).There is a nonsmoooth minimum (actually a zero) of M + = 0 coinciding with the saddle fixed point x E .The nonsmooth minimum shows a discontinuity at the same location for B. The minimal seed does not show any distinctive feature along .The influence of p 1 has been well studied in Ref. [51].It is mainly concerned with type of singularity.p = 1 favors a linear behavior for M + ; however, the ridge in B remains.A 1D cut along the LD maps is shown in Figs.10(a) and 11(b) for the discussion in Sec.IV B. The ridge in the FTLE (or LD) maps becomes thinner with increasing horizon time τ .From each of these two diagnostics alone, the edge manifold is hence identified as a repelling LCS once the horizon time is long enough.

C. Lebovitz-Mariotti model (LM6D)
Increasing the complexity of the system to six dimensions and nontrivial attractors, we consider the model introduced in the work of Lebovitz and Mariotti [59], from this point called LM6D.The model was originally suggested to illustrate concepts of complex boundaries and boundary collapse inspired by transient chaos in hydrodynamics [13,60].The derivation of the model is similar to that of Ref. [56] and the details of the derivation are introduced in Ref. [61].The matrix L is such that where k i , i = 1, . . ., 6 and σ i i = 0, . . ., 6 are parameters of the model [62].The choice of parameters for the present study shows a bistable system, and thus with a hard edge.The bifurcation diagram in Ref. [59] shows that for the working parameters there is a stable laminar fixed point and a torus to be interpreted as the turbulent state.The edge manifold separates both basins and contains the lower branch solution originating from a saddle-node bifurcation.This solution is an unstable periodic orbit (UPO) and corresponds to the edge state.Typical trajectories of LM6D are shown in Fig. 4(d).
The figure shows a phase portrait using three variables x 1 , x 3 , and x 5 .The three plotted trajectories are respectively below (purple), above (orange), and on , tending asymptotically to the UPO (blue) (we assume by default that the manifold is orientable).The LCS diagnostics for LM6D are illustrated in an x 4 x 5 plane where the rest of the variables have initial values x 1 = 0.3, x 2 = x 3 = 0, x 6 = 0.1.The results are shown in Figs.4(a), 4(b) and 4(c) for both FTLEs and LDs for a common time horizon τ = 500.A clear ridge is present, consistently with the trajectories approaching different attractors in Fig. 4(d).As for FTLEs, the area on the left of the ridge in Figs.4(a) and 4(c) is smoother due to the trajectories all approaching the same fixed point, while the right side displays oscillations depending on which part of the torus the trajectory has reached at the final time.In order to enhance the ridge in the FTLE map, its gradient is shown in Fig. 5.The LD map is computed for p = 0.5 and shows isovalues of B in Fig. 4(b), where a ridge also emerges, as further confirmed in Fig. 4(c).These results show that the edge manifold is again highlighted as a repelling LCS in a nonchaotic case when neither the edge state nor the turbulent attractor are fixed points.As we shall see in the next examples, this generalization meets its limits in the presence of chaos.

D. Moehlis-Faisst-Eckhardt model (MFE9D)
The next example in the model hierarchy is a ninedimensional model suggested by Moehlis, Faisst, and Eckhardt [25], hereafter referred to as MFE9D.It possesses again a linearly stable laminar state x L , but unlike the models in Secs.III B and III C no turbulent-like attractor is present.Instead, the turbulent state appears as a chaotic nonattracting set [13,63], and the systems possess therefore a weak edge.This model thus contains properties common to real subcritical fluid systems [13,64].The parameter defining the state space in the MFE9D is the Reynolds number, set to Re = 400 in the current study.For the chosen set of parameters [65], the dynamics of the model is described in Ref. [25], and the edge state corresponds to an UPO [18].
The LCS diagnostics are applied in MFE9D for a horizon time of τ = 800.The results are shown in Fig. 6 using a projection on the x 2 x 3 plane defined by x 1 = 0.7066, x 4 = 0.01, x 5 = x 6 = x 7 = x 8 = x 9 = 0.The laminar state lies in the bottom left corner of the figure at x 2 = x 3 = 0.
For both LCS diagnostics, two different regions, respectively smooth and speckled, emerge in Fig. 6.They correspond respectively to state space regions with two different behaviors, either uneventful relaminarization or visiting the chaotic saddle.
The FTLE field is plotted in Fig. 6(a).It shows a sudden transition between the smooth and speckled regions both in terms of values reached and, unlike the results for the previous models, in terms of fluctuation level.No proper ridge of the FTLE emerges, however by plotting the norm of the gradient of the FTLE with respect to the variables x 2 x 3 , a ridge stands out in Fig. 6(b).The performances of the LDs are illustrated for p = 0.5 in Fig. 6(c) using isovalues of the gradient norm B. In this figure, the ridge also stands out around smooth region and separating it from the speckled region.The existence of this sudden transition in the LCS diagnostics is highlighted in Fig. 6(d) by plotting them along a 1D cut of state space.Figure 6(d) shows several abrupt changes of B within the speckled region; however, the edge can still be identified as the first abrupt change encountered when coming from the smooth side.Consequently, before using LDs to highlight the edge as an LCS some preliminary knowledge about the qualitative behavior of the model is required.We attribute again the speckledness property to the presence of chaos in the turbulent basin.Strong visual analogies exist between the results from the LCS diagnostics and the lifetime plots in Ref. [18] for the same model.

E. Navier-Stokes equations: Plane Couette flow (pCf)
Eventually, we demonstrate the relevance of the previous Lagrangian diagnostics to highlight the edge manifold in a very high-dimensional system governed by the incompressible Navier-Stokes equations.Formally this system is of infinite dimension; however, it is made finite dimensional by the presence of physical viscosity and numerical discretization.The resulting system is still of such huge dimension-here of order O(10 5 )-that scalable methods are necessary for the diagnostics.We focus here on plane Couette flow (pCf), where a viscous fluid is sheared between two plates moving with opposite velocities in a direction called X .pCf is parametrized by a nondimensional Reynolds number Re proportional to the plate velocities.All the quantities are made nondimensional with the half-gap between the plates h and the plate velocity U w , and the Reynolds number is Re = U w h/ν = 400, with ν being the kinematic viscosity.The Navier-Stokes equations for the perturbation velocity field to the laminar state U L , u = (u X , u Y , u Z ), read where p = p(X, Y, Z, t ) is the pressure field [67].pCf is the archetype of a subcritical system, where for all values of Re large enough, a turbulent state T coexists with a linearly stable laminar state L. The laminar state has a velocity field U L = ye X , where e X is the longitudinal unit vector and a homogeneous pressure field.More precisely, we consider as in Refs.[66,68] that the fluid particles move inside a numerical domain of dimensions L X × L Y × L Z (in units of h), made artificially periodic in both in-plane directions X and Z and bounded in the wall-normal direction Y .The exact dynamics of the edge state depends on the parameters Re, L X , L Y , and L Z .However, there are robust common features to all of them: The presence of streamwise streaks (transverse spatial modulations of the streamwise velocity field, weakly modulated in the longitudinal direction, associated with longitudinal vortices).For the parameters under study [Re = 400, (L X , L Y , L Z ) = (5.513,2, 3.770)], the edge state is known and consists of a time-periodic flow field with period T p = 85.5 [66].Its spatial structure is illustrated in Fig. 7(a).The structure of its stable manifold, however, and of the edge manifold in general, is poorly understood despite some recent progress [35,69].For this parameters, as in the MFE9D model, the turbulent trajectories are supertransients [63] and there is no turbulent attractor.
An observable characterizing the laminar-turbulent transition is the streamwise vorticity squared |ω X | 2 averaged over the computational domain a(t and V is the volume of the computational domain.In order to identify the edge manifold as an LCS, we use an arbitrary point on the edge state as a reference.That state space point corresponds to a flow field in the three-dimensional physical space (X, Y, Z ).We will use the perturbation velocity field notation u E for such a point.Note that u E is the perturbation to U L and the full velocity field is V = U L + λu E with λ = 1 corresponding to the edge state.By tuning a real parameter λ, one expects a change in behavior of both FTLEs and LDs at λ = 1.In order to extend the definition (8) to the Navier-Stokes case, we first need to specify which projection system is used to define the corresponding high-dimensional state space.The choice is not unique.For simplicity, we choose to rely on the values u k (X i ) of the velocity at each discrete grid point, where the index (k = 1, 2, 3) refers to the velocity component (k = X, Y, Z) and i = 1, ..., N, where N = N X × N Y × N Z is the total number of grid points.
The expression for the Lagrangian Descriptor scalar valued M + along a trajectory is where τ > 0. Other physically meaningful formulations are possible, e.g., a velocity-vorticity formulation and/or a spectral Galerkin decomposition.If we denote by γ the state space coordinate along which the bisection is performed, the gradient of M + in that direction is given by If B is evaluated in a d-dimensional region of the state space (with d N), its partial derivatives need to be evaluated only along these d directions, each direction being parametrized by a coordinate γ i .The expression for B becomes B is computed numerically using second order centered finite differences as in Eq. ( 17).In all simulations, ε = 10 −9 and p = 0.5.
The direct computation of all FTLEs using classical methods is indubitably out of reach in the present system because of the huge value of N. The so-called reduced-order FTLEs can, however, be accessed as a by-product of the calculation of optimally time-dependent (OTD) modes [48].In this approach, a finite number r of these time-dependent vectors can be evolved in time together with the main trajectory.They have the property of spanning a reduced r-dimensional space which approximates the tangent space, from which the finite-time exponents can be directly computed under mild conditions [49,70].Although only the leading FTLE is of interest here, it is recommended to compute them with r > 1 to avoid spurious results linked with eigenvalue crossings [49].The present simulations involve r = 4 OTD modes.The four corresponding initial conditions are the edge state point u E and the three linear optimal perturbations (LOPs) computed for the wave numbers (k X , k Z ) = (0, 2π/L z ), (2π/L x , 0), and (2π/L x , 2π/L z ).For each wave number (k X , k Z ), the corresponding LOP is defined in Ref. [12] as the initial condition (u 0 ) on the unit sphere maximizing the energy gain G τ , where FIG. 8. One-dimensional landscape of Lagrangian diagnostics for several values of λ in the neighborhood of the edge at λ = 1.Plane Couette flow, with the same parameters as Ref. [66].
where || • || is the usual L 2 norm and F t 0 +τ t 0 is the propagator associated with the Navier-Stokes operator [12] linearized around U L , which is independent of t 0 .The value of τ used here is the one that gives the largest value of G τ over all initial conditions.Interestingly, the quantity (log G τ )/2τ can be interpreted as the leading FTLE at time t 0 over a time horizon τ around the fixed point U L at time t 0 , in the subspace spanned by the corresponding eigenvector [71].It is known that eigenvalue crossing can affect the resulting OTD subspace [49] and therefore the reduced-order FTLEs.In order to converge to the relevant subspace, the OTD modes are first evolved for 200 time units along the edge trajectory and then used as an initial condition for the computation of the reduced-order FTLEs.
In Figs. 8 and 9, we show landscapes of the two indicators in a one-dimensional cut and in a two-dimensional slice of the state space, respectively.Figure 8 displays the 1D landscape of LCS diagnostics as λ is varied, for a time horizon τ = 300, for both the leading FTLE and B. The FTLEs display a smooth behavior for λ < 1, where trajectories approach the laminar state, while for λ > 1 the landscape looks jagged in a way similar to the MFE9D in Fig. 6(d).The LD landscape also displays the same properties; however, the values of B in the laminar basin are considerably lower than for FTLEs, resulting in a flatter landscape.In both cases, the crossing of the edge manifold at λ = 1 corresponds to a steep increase of the quantity plotted, this property being better visible in the case of B. This steepness property has its analog in lifetimes studies in former works dedicated to chaotic saddles [9,18,30].The geometry of the state space can be further explored by considering a two-dimensional slice spanned by two perturbation fields u 1 and u 2 .These two fields are normalized so that u 1 = u 2 = 1.An initial perturbation u 0 to the laminar state can be defined as where A is interpreted as an amplitude.We chose u 1 = u E , such that for α = 0 the same subset of state space is explored as in Fig. 8 between FTLEs and LDs appears in the comparison between the two subfigures.These holes correspond to initial conditions which relaminarize after a turbulent transient shorter than τ .A qualitative comparison with lifetime maps as used in Refs.[18] and [35] suggests again similar features.

IV. EDGE TRACKING REVISITED
The previous sections have demonstrated that LCS diagnostics are efficient to capture stable manifolds in their globality (the edge), and we now show that these diagnostics can be adapted in order to identify the relative attractors sitting on them, i.e., edge states.In dissipative systems of interest, edge states are typically invariant sets of low dimension, e.g., fixed points, limit cycles, or low-dimensional chaotic sets.In the hydrodynamic literature, edge states have proven crucial to unfold bifurcation diagrams [35,36], for dynamic control [68] and because of their role in turbulence nucleation [16,72].Until now, edge states have been routinely identified numerically using standard bisection methods coupled with prior knowledge of the location of the edge states [18,21].The present section introduces another class of bisection methods based on LCS diagnostics.We will contrast the related concepts of global versus local bisection methods in terms of their ability to identify the edge state, their scalability, and their associated computational cost.

A. Global versus local methods
The standard bisection method used for edge tracking relies on some preliminary knowledge of the part of state space where the edge state resides.A scalar observable a(t ) (typically the perturbation kinetic energy) is chosen such that the edge state lies entirely in an interval a ∈ (α A , α B ) whereas the other states (L and T for simplicity) lie outside it.Provided this constraint is met, the initial condition x 0 is adjusted recursively so that F t 0 +τ t 0 (x 0 ) remains in (α A , α B ) for τ as large as possible.Such a method warrants that F t 0 +τ t 0 (x 0 ) converges toward the edge state as τ → ∞.We label this first bisection method as global, since it uses information from the whole state space to define .It has proven useful in the presence of both hard-and soft-edge manifolds.
Global methods can, however, fall short whenever no available global observable can elucidate on which side of the edge the different trajectories evolve [31,73], when one of the attractors undergoes a local bifurcation affecting the values of the bounds [74] or simply when information on the bounds α A,B is not available.
As shown in the examples above, the LCS diagnostics are based on the knowledge of the tangent vector and the Cauchy-Green tensor, which are local in state space and rely only on a finite-time description.Since they are sufficient to highlight the edge manifold, we expect these diagnostics to yield revised local bisection methods by opposition to the global ones.The use of a local method for bisection offers the following advantages: (1) No preliminary choice of observable is needed based on physical intuition.
(2) No prior knowledge of bounds α A and α B is needed.(3) The approach to edge tracking is not binary; edge trajectories are labeled as optimizers of an LCS-based observable.(4) The edge tracking has a predefined time horizon τ chosen as an intrinsic parameter.The values of τ are limited by numerical accuracy; however, local edge tracking can be restarted from any location, making the algorithm iterative and eventually convergent.
(5) LCS-based methods allow for a quick identification of the regions of interest.
Although the distinction between local and global methods is pedagogically convenient, the designation "local method" needs to be nuanced in practice: the shorter the time horizon τ the more local the approach.For larger values of τ , nonconvergent trajectories are not bound to local neighborhoods and can explore remote parts of the state space.The analysis of the LCS diagnostics in the previous section suggests long time horizons τ for a sharper identification of the edge.

B. Comparison of different methods for edge tracking
LCS-based bisection relies as usual on the iterative process of straddling the edge trajectory locating the right initial condition along an arbitrary state space direction.Unlike with the classical bounds-based method, the revised edge trajectory emerges now as the optimizer of a given functional rather than from the difficultly quantifiable "neither laminar, nor turbulent" definition.For a given time horizon τ > 0, FTLE-based bisection seeks the initial condition x 0 maximizing λ t 0 +τ t 0 (x 0 ).LD-based bisection seeks either the minimum of M + (x 0 ) or the maximum of the gradient norm B(x 0 ), depending on the dynamical nature of the edge state.The various examples in Sec.III point toward the following phenomenology: When the edge state is a fixed point, it is sufficient to use minima of M + ; however, this criterion needs to be replaced by maxima of B, which is computationally more costly, for any edge state dynamics of higher dimensionality.These different cases are contrasted in Fig. 10  The simplest case for the comparison of the different edge tracking methods is the two-dimensional Dauchot-Manneville model (DM2D), for which the edge state x E there has an analytical expression in Eq. ( 16).For the sake of generality, the LD-based edge tracking is not defined using M + but rather using its gradient B, whereas the FTLE-based method is based on the estimation of the largest FTLE.The present methods rely, for the proof of concept, on simple algorithms to locate the maximum of the associated fields along an arbitrary one-dimensional line L in state space (in practice the line L : {x 1 = 0} was selected).The maximization is always initiated on the "smooth" side, i.e., within the basin of attraction of x L .
The maximization over the line L of the LCS diagnostic [λ τ (x 0 ) or B(x 0 )] generates a sequence of new initial conditions x (k)  0 , k = 0, 1, 2, . . . on the line L. Convergence is satisfied if the sequence x (k) 0 approaches asymptotically to some x * 0 ∈ .Since the goal is to identify x E , a convergence distance D min , based on the minimal Euclidian distance along the resulting trajectory to x E , is preferred.It is defined as where the minimum is taken over the time interval (t 0 , t 0 + τ ].The largest FTLE is computed using finite differences with ε = 10 −9 .As for the Lagrangian descriptors, the gradient B = ||∂M + /∂x 0,2 || is computed following Eq.( 17) using p = 0.5.Figure 11(a) shows the minimum distance D min to the edge state for the classical, the FTLE-based and the LD-based edge tracking algorithms.In the classical edge tracking, τ > 0 is the time it takes for a single trajectory closest to the edge to reach the bounds, and it varies from iteration to iteration.In the local edge tracking algorithms, however, τ > 0 is a prescribed time horizon.Figure 11(a) shows that, for equivalent D min , the so-called local methods both require a lower value of τ compared to the classical method.
A direct mutual comparison of the local methods suggests that FTLE-based edge tracking requires time horizons τ one order of magnitude larger than the LD-based edge tracking to achieve equivalent D min .Figure 11(a) also illustrates the bottom value reached by D min ≈ 10 −9 .This can be understood as follows: FTLEs, unlike LDs, do not display any ridge near hyperbolic saddle points if the dynamics is governed by a linear system [45,51].As a consequence, the trajectory in the neighborhood of the saddle point x E must "feel" nonlinear effects for FTLEs to be able to identify the edge as a LCS.This implies that longer times are needed for the bisection based on FTLEs compared to that based on LDs.Similarly convergence is expected to be faster in the LD case, at least when the edge state corresponds to a fixed point.Longer times can be requested in cases with more complex attractors, although this demands a more systematic investigation.
The convergence properties of these two LCS-based edge tracking algorithms in higher dimension are expected to proceed along the same lines because of the strong scalability properties of the diagnostics themselves.In the case of LDs, the function M + is computed directly along the trajectories, without using the Jacobian operator ∇ x f (x(t ), t ).The function M + is hence trivially generalized to arbitrary dimension according to (8), as demonstrated in (20).The full computation of all FTLEs, on the other hand, becomes unfeasible in higher dimension since ∇F t t 0 is of prohibitive size n × n, as pointed out in Ref. [49].The computation of r OTD modes, with 1 r n, can, however, grant access to the r leading FTLEs at a cost n × (r + 1) + r 2 [49].For r = 1, this amounts to O(2n), and the computation of the largest FTLE is hence based on a scalable algorithm as well.

V. CONCLUSION, DISCUSSION, AND OUTLOOK
Motivated by the physical problem of subcritical laminarturbulent transition to turbulence in hydrodynamics, we have investigated the established notion of edge manifold which charts the state space into different basins of attraction [13].
The main contribution of this study is the recognition of stable manifolds of edge states in autonomous systems of arbitrarily high dimensions as Lagrangian coherent structures (LCSs).The edge manifold divides the state space in the same way as a hyperbolic repelling Lagrangian coherent structure, and hence a direct relation between LCS and exact coherent states can be established.The same mathematical toolbox used to highlight an LCS can then be used to highlight the edge manifold, no matter how high the dimension of the underlying state space.The price to pay for this generalization is that only forwardtime operators can be considered, which prevents one from considering attracting LCS sets such as unstable manifolds.An illustration based on two LCS diagnostics, finite-time Lyapunov exponents (FTLEs) [38] and forward-time Lagrangian descriptors (LDs) [43], has revealed the edge manifold in several autonomous dynamical systems of increasing dimension and complexity.By order of dimension, these models are (i) a didactic bistable two-dimensional model with only fixed points [57], (ii) a bistable six-dimensional model where the edge state is an unstable periodic orbit and turbulence corresponds to a stable periodic orbit [61], (iii) a nine-dimensional model where the edge state is a periodic orbit and turbulence is represented by a chaotic saddle [25], and eventually (iv) wellresolved Navier-Stokes simulations of a periodic cell of plane Couette flow, understood as an O(10 5 )-dimensional dynamical system [68].For all these models, both LCS diagnostics display steep variations at the location of the edge manifold, which can be monitored using a single scalar quantity.Our findings indicate that having some prior knowledge of the system is helpful to properly identify the edge as an LCS and that LDs are particularly sensitive to this.Approaching a model with unknown dynamical characteristics should be done in the following way: (1) the LCS diagnostic should be applied for the parameter region of interest (2) regions on both sides of a ridge in the FTLE or the boundaries between smooth or speckled regions of B should be systematically explored.After testing different values of τ , our exploration suggests that the edge manifold emerges already for τ = O(1) but stands out more dramatically for longer time horizons, a feature worth a more quantitative study.
Furthermore, the same toolbox allows, when it is unknown, to identify the edge state (the relative attractor on the edge manifold) using iterative bisection algorithms based on the LCS diagnostics, i.e., local measures in state space rather than global considerations as in previous works.The associated algorithms have been tested on the Dauchot-Manneville model in two dimensions using both FTLEs and LDs: Both outperformed the classical approaches, with a strong advantage for the LD-based methods in cases where the edge state is a fixed point both in terms of scalability and computational cost.However, edge states are generally not fixed points.
As pointed out in Ref. [46], LDs are not objective and their connection to LCS is unclear.This has lead to some controversy regarding the interpretation of the LD maps [51,75,76].A hitherto unreported drawback of LDs, worth a more systematic characterization, is that their ability to highlight the edge manifold depends on the dynamical nature of the invariant sets of the system.Consequently, their low computational cost comes with uncertainty about how to interpret the plots in an system where the nature of the attractors is unknown.The FTLEs represent a computationally more expensive tool.They have a solid mathematical background, are objective diagnostics [46], and with additional constraints can be precisely connected to hyperbolic LCSs [38].Furthermore, they have a more direct interpretation in terms of material barriers.
The objectivity of the FTLEs makes them independent of the frame of reference, displaying the same values on the ridges and highlighting the edge as an LCS independently of the chosen frame.LDs are not objective and changing the frame of reference would imply a change in the values of the LD maps and even their nature.A relevant example is that of an edge state consisting of a traveling wave (TW) in physical space.Choosing a fixed frame of reference would yield a different result for the LD landscape than considering a frame moving with the velocity of the TW.
The present work paves the road for the design of effective local manifold tracking methods which do not rely on measurable global properties.The local edge tracking algorithms do not rely on a Boolean approach, but can be defined as optimizers of a scalar observable.We point out that the identification of local edge properties is not limited to bistable systems, and therefore as future work a revisit of manifold-tracking methods in different situations such as [24,31,73] using this framework is strongly encouraged.

FIG. 1 .
FIG. 1. Sketch of the analogy between Lagrangian coherent structures in conservative systems and the state space geometry of subcritical laminar-turbulent transition in hydrodynamics.(a) Lagrangian patterns in nature: wind map from March 2, 2020, at a height of 850 hPa [39].Sketched attracting and repelling LCS (blue and red, respectively) identified with the unstable (resp., stable) manifold of a saddle point (green).(b) Sketch of the state space for a bistable hydrodynamical system, where the stable manifold of the edge state (green) is the basin boundary (red), whereas the unstable manifold (blue) leads to either the laminar or the turbulent attractor.The 3D figures represent flow fields from actual numerical simulations of plane Couette flow.Contours of streamwise velocity in the flow (red to blue) and an isovalue λ 2 criterion (green) used for vortex identification.

FIG. 4 .
FIG. 4. Lagrangian analysis of the LM6D model.(a) FTLE for τ = 500.(b) B for p = 0.5 and τ = 500.(c) Landscapes of LCS diagnostics from the maps above along a line with constant x 4 = 0.4.Both FTLE and B are normalized with respect to their maximum in the plotted interval.(d) 3D phase portrait of three different trajectories, where thicker lines indicate the approached attractor; see text for details.

FIG. 6 .
FIG. 6. Lagrangian analysis of the MFE9D.(a) FTLE for τ = 800.(b) Norm of the gradient of the FTLE.(c) B for p = 0.5 and τ = 800.(d) Landscape of the Lagrangian diagnostics for τ = 2000 along a line of x 3 in the neighborhood of the edge for fixed x 2 = 0.15.

FIG. 7 .
FIG. 7. Flow fields corresponding to (nonrescaled) initial conditions in Eq. (23) for plane Couette flow.(a) u 1 edge state in Ref. [66].Contours of u Z = ±0.045 in red and blue respectively.Black lines: contours of u X = −0.3 on cross-flow planes at X = const.(b) u 2 LOP for wave numbers (2π/L X , 2π/L Z ).Contours u Z = ±0.45 in red and blue respectively.Flow along the x direction.
FIG. 9. 2D maps of LCS diagnostics for plane Couette flow, same parameters as Ref.[66], τ = 300.Left: largest reduced-order FTLE.Right: isovalues of B for p = 0.5.The boundary between the speckled and the smooth area corresponds to the finite time approximation of .The orange dot indicates a projection of the edge state in the α-A state portrait.

FIG. 11 .
FIG.11.Results for edge tracking along x 2 in the DM2D model, for a fixed x 1 = 0. Left: Euclidian distance to the edge state x E depending on the objective time τ of the bisection algorithm employed.Tolerance for algorithm convergence 10 −9 .Right: LCS diagnostic landscape for x 2 in the neighborhood of .Both FTLE and B are normalized with respect to their maximum in the plotted interval.
where one-dimensional landscapes of M + are displayed.FTLE landscapes are shown for comparison in Figs.11(b), 4(c), and 6(c).