The edge as a Lagrangian Coherent Structure in a high-dimensional state space

Dissipative dynamical systems characterised 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 automous arbitrarily high-dimensional dissipative systems, in which the edge manifold is re-interpreted as a Lagrangian Coherent Structure (LCS). Two different diagnostics, finite-time Lyapunov exponents and Lagrangian Descriptors, are used and compared with respect to their ability to identify the edge and to 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 flow. 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 efficient 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, prediction, as well as for control. Besides such basin boundaries, being invariant sets of the system, also support 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], synchronisation 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] 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 [9].
In the field of hydrodynamics, most incompressible viscous fluid flows near solid boundaries, can undergo transition to turbulence in a subcritical manner [10]. The flow configurations of interest include the flow inside circular or rectangular pipes [11], between two plates in motion [12], inside rotor-stators flow [13] or the flow developing on a semi-infinite flat plate [14]. The laminar state in the examples above is time-independent, characterised by spatial symmetries, and free of fluctuations, while the turbulent state is, depending on the level of modelling, arbitrarily complicated. For the parameters of interest here, the laminar state is linearly stable and coexists for the parameters with the turbulent set. Some specific trajectories neither reach the laminar state nor the turbulent regime, but instead stay arbitrary close to their respective basins of attraction,. They have been labelled edge trajectories [15]. Based on computational attempts to identify edge trajectories, the set of such trajectories has been conjectured to be a differentiable manifold of codimension one in the state space, called either edge of chaos, laminar/turbulent boundary, boundary of turbulence or simply the edge [16][17][18]. At high enough 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 [19,20] who used a simple bisection method to identify the invariant dynamics on the edge, corresponding to a travelling wave solution (a result revised later). The analysis of this orbit and of its instability yielded useful information for the understanding of the turbulent dynamics itself, notably regarding bursting dynamics, interpreted as a homoclinic connection between the travelling wave on the edge and a more complicated turbulent-like regime [21].
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 [22] we will refer to it as a strong edge. Ample experimental and numerical evidence has shown that for low enough 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 generalised to such leaky cases by focusing on the values of the lifetimes associ- ated 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 [23] the basin boundary can display a transversely fractal structure, as a consequence of the Lambda-lemma [24]. 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 [25], all based on a maximisation of the lifetime [16,17,[26][27][28]. 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 invariably take 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 [29]. 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 neighbouring 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 usual textbooks. A different paradigm, based on quantitative criteria, optimal quantities and on a finite-time framework, is hence necessary for a more transparent and universal characterisation 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 travelling 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 states exist for the same parameter values [18,30]. 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 minimising the distance to the laminar state [31,32]. 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 [33]. 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 [34,35]. 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 analysed 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 [36]. LCS 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 [37,38]. Other applications, for instance chemical reactions, have also been considered [3]. Mathematically yet, 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 non-autonomous 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 equatioṅ x = 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 LCS can be exported to the realm of high-dimensional dissipative systems. More specifically it can be used for the investigation of the edge manifold of transitional flows according to the conceptual analogy expressed in Fig. 1. 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, being this the generic situation in all the hydrodynamic models considered herein. We focus in particular mainly on two distinct tools popular in the recent literature on LCS: finite-time Lyapunov Exponents (FTLEs) [36,39] and Lagrangian Descriptors (LDs) [40,41]. We give key examples of application 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 developed framework we suggest and test the use of these new diagnostic to improve edge tracking algorithms.

II. DEFINITIONS AND LCS DIAGNOSTICS
We consider here a general dynamical system defined on a subset A ⊂ R n , governed bẏ where f : A × [t 0 , t 0 + τ ] → R n is a sufficiently smooth vector field and τ > 0. Let F t t0 be the flow map which maps an initial position x 0 at time t 0 to its position at time t. The system linearised around the trajectory where ∇F t t0 is the deformation gradient.

A. Finite-Time Lyapunov exponents
Finite-time Lyapunov exponents have been popularised as a tool for the identification of LCS in [37,42]. Let C t t0 be the (positive definite) Cauchy-Green tensor where (·) * represents the Hermitian transpose. Let λ i i = 1, ..., n denote the eigenvalues of C t t0 and ξ i its associated eigenvectors such that and where ||·|| denotes the L 2 -norm. By defining τ = t − t 0 > 0, the average growth rate 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 i th finite-time Lyapunov exponents 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 + τ (labelled simply λ t 0 t 0 + τ ) can be used as a diagnostic of hyperbolic LCS [36,43,44]. The Cauchy-Green tensor has size n × n, and its entries are usually computed using second-order centered finite-differences [36]. Hyperbolic LCS refer to attracting and repelling distinguished invariant manifolds [37], with forward time FTLEs relevant for the identification of repelling LCS, and backward time FTLEs relevant for the identification of attracting LCS [45]. The computational cost of FTLEs following the above method increases rapidly with the state space dimension and is hopeless 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 [46]. 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 [47]. In practice for the present LCS diagnostic, only the largest exponent needs to be computed, for which other methods are possible too.

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 [40,48] and its theoretical framework further developed in [49]. LDs are based on the integration of a given observable along trajectories. The original quantity of interest is In Eq. (7), 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 (7) used in the literature [50,51] although other alternative definitions have been suggested. It is convenient to split Eq. (7) into its forward and backwards contribution [3,50]: Due to the dissipative nature of Eq. (1), numerical backwards integration is ruled out for stability issues, 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 (x 0 , t 0 , τ ) + . LDs have been used to identify boundaries between qualitatively different dynamics, based on abrupt variations of M + . The abrupt change also means that the derivative of M + transverse to the boundaries is discontinous. These singular features in LD plots are often connected to the stable and unstable manifolds of saddle points. In such a case this leads to a finite-time redefinition of the manifold according to W s (x 0 , t) = argmin M + (x 0 , t 0 , τ ) [50,52]. However, based on our investigation, we observe that this definition does not generalise to more complex dynamics on Σ. Alternatively, it is useful to quantify the abrupt changes of M + and thus to define the L 2 -norm of its gradient Many other diagnostics for LCS have been suggested, see [44] 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. [44,53]). 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 [54]. We assume that the bold vector X ∈ R 3 represents the position in the physical space coordinates (which differ from the state space coordinates), 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 confuse 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 an n×n linear operator and N a quadratic form containing the nonlinear terms. As explained in [55], all models of subcritical transition consistent with the original PDEs are subject to two constraints: i) L is a non-normal 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.
As an illustrative case, we consider the simple twodimensional model introduced by [56], 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 [57] 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 lam- 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 x 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. Fig. 2 shows a phase portrait of the model for the parameters s 1 = −0.1875 and 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(a). Fig. 3(d) shows the largest FTLE along part of Σ -normalised to 1 at the start of the curve-, indicating a rather smooth and uniform value of the FTLE along the ridge. The application of LDs is assessed by the contours of M + in Fig. 3(b), which highlight a state space structure perfectly conform to the phase portrait of Fig. 2. There is a non-differentiable minimum of M + across Σ, as in former applications of LDs to non-dissipative systems. Isocontours of its gradient norm B can be convenient for plotting purposes (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 exclusively the dip of M + map of Fig. 3(b) that coincides with Σ. The variations of M + and B along the edge are quantified in Fig. 3(d). There is a non-smoooth minimum (actually a zero) of M + = 0 coinciding with the saddle fixed point x E . The non-smooth 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 [49]. It is mainly concerned with type of singularity. p = 1 favours a linear behaviour for M + , however the ridge in B remains. A 1D cut along the LD maps is shown in Fig.  9(a) and Fig. 10(b) for the discussion in Section 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 non-trivial attractors, we consider the model introduced by [58], 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 [11,59]. The derivation of the model is similar to that of [55] and the details of the derivation are introduced in [60]. The matrix L is such that and the vector N reads where k i , i = 1, . . . , 6 and σ i i = 0, . . . , 6 are parameters of the model [61]. The choice of parameters for the present study shows a bistable system, and thus with a hard edge. The bifurcation diagram in [58] shows that for the working parameters there is a stable laminar fixed point and a torus to be interpreted as 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 Σ, asymptoting to the UPO (blue) (we assume by default that the manifold Σ is orientable). The LCS diagnostics for LM6D are illustrated in an  The next example in the model hierarchy is a ninedimensional model suggested by [23], hereafter referred to as MFE9D. It possesses again a linearly stable laminar state x L , but unlike the models in Section III B and Section III C no turbulent-like attractor is present. Instead the turbulent state appears as a chaotic non-attracting set [11,62], the systems posses therefore a weak edge. This model thus contains properties common to real subcritical fluid systems [11,63]. The parameter defining the state space in the MFE9D is the Reynolds number, set to R = 400 in the current study. For the chosen set of parameters [64] the dynamics of the model is described in [23], and the edge state corresponds to an UPO [16].
The LCS diagnostics are applied in MFE9D for a horizon time of τ = 800. The results are shown in For both LCS diagnostics two different regions, respectively smooth and speckled, emerge in Fig. 5. They correspond respectively to state space regions with two different behaviours, either uneventful relaminarisation or visiting the chaotic saddle. The FTLE field is plotted in Fig. 5(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 1 x 2 , a ridge emerges in Fig. 5(b). The performances of the LDs are illustrated for p = 0.5 in Fig. 5(c) using isovalues of the gradient norm B. In this figure the ridge also stands out as a continuous line enclosing the uniform, smooth region and separating it from the speckled region. The existence of this sudden transition in the LCS diagnostics is highlighted in Fig. 5(d) by plotting the diagnostics along a 1D cut of state space. 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 [16] for the same model. 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 made finite-dimensional by the presence of physical viscosity and by the numerical discretisation. The resulting system is still of such huge dimension -here of order O(10 5 )-that scalable methods are required 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 parametrised by a non-dimensional Reynolds number R proportional to the plate velocities. All the quantites are made non-dimensional with the half gap between the plates h and the plate velocity U w , the Reynolds number is R = U w h/ν = 400, with ν 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 (17) and where p = p(X, Y, Z, t) is the pressure field [65]. pCf is the archetype of a subcritical system, where for all values of R 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 Ref. [66,67] 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 R, L X , L Y and L Z . However there are robust common features to all of them from a structural point of view: 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 (R=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.  6(a). The structure of its stable manifold, however, and of the edge manifold in general, is however poorly understood despite some recent progresses [33,68]. For this parameters, as in the MFE9D model, the turbulent trajectories are supertransients [62] and there is no turbulent attractor. An observable characterising the laminar-turbulent transition is is the streamwise vorticity squared |ω X | 2 averaged over the computational domain a(t) 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 threedimensional 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. By tuning a real parameter λ, one expects to identify a change in behaviour of both FTLEs and LDs at λ=1. In order to extend the definition (7) to the Navier-Stokes case, we first need to specify which projection system is used to define the corresponding high-dimensional state space. There is no unique choice. 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 M + along a trajectory is where τ >0. Other physically meaningful formulations are possible, e.g. using 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 ∂ λ M + (V 0 (λ), t 0 , τ ). 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 parametrised by a coordinate λ i . The expression for B becomes The expression for B is computed numerically using second order centered finite differences as in Eq. (16). 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 byproduct of the calculation of Optimally Time-Dependent (OTD) modes [46]. 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 directy computed under mild conditions [47,69]. The present simulations involve r = 4 OTD modes. 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 [47]. The four corresponding initial conditions are the edge state point u E , and the three linear optimal perturbations (LOPs) computed for the wavenumbers (k X , k Z )=(0, 2π/L z ), (2π/L x , 0) and (2π/L x , 2π/L z ). For each wavenumber (k X , k Z ), the corresponding LOP is defined by [10] as the initial condition u 0 ) on the unit sphere maximising the energy gain G τ , where where ||·|| is the usual L 2 -norm and F t0+τ t0 is the propagator associated with the inearised Navier-Stokes operator [10] linearised around U L , which is independent of t 0 . The value of τ used here is the one that gives the largest G τ over all initial conditions. Interestingly the quantity (log G τ ) /2τ can be interpreted as the leading FTLE at  time τ around the fixed point U L at time t 0 , in the subspace spanned by the corresponding eigenvector [70]. It is known that eigenvalue crossing can affect the resulting OTD subspace [47] and therefore the reduced-order FTLEs. In order 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 Fig. 7 and 8 we show landscapes of the two indicators in a one-dimensional cut and in a two-dimensional slice of the state space, respectively. Fig. 7 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 behaviour for λ < 1, where trajectories approach the laminar state, while for λ > 1 the landscape looks jagged in a way similar to the MFE9D in Fig. 5(d). The LD landscape also displays the same properties, however the values of B in 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 [16,28]. 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 both normalised so that ||u i || = 1, i = 1, 2. 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 7(a). As for u 2 we chose the LOP corresponding to the wavenumbers (k X , k Z ) = (1, 1). The velocity fields for the initial conditions u i i = 1, 2 are shown in Fig. 6. The resulting LCS diagnostics for τ = 300 are shown in Fig. 8. In each subfigure, two very different regions appear, separated by a smooth boundary. It has been monitored that the zone containing the laminar state corresponds to rapidly relaminarising trajectories, whereas the speckled zone on the other side of the boundary contains trajectories visiting the tur- bulent state. In other words, for both diagnostics the basin boundary is successfully identified in finite time.
For both indicators, some 'holes' appear inside the turbulent basin, characterised by low values typical of the laminar basin. No strong difference between FTLEs and LDs appears in the comparison between the two subfigures. These holes correspond to initial conditions which relaminarise after a turbulent transient shorter than τ . A qualitative comparison with lifetime maps as used in [16] and [33] suggests again comparable features.

IV. EDGE TRACKING REVISITED
The previous sections have demonstrated that LCS diagnostics are efficient to capture stable manifolds in their globality (the edge), we now show that these diagnostics can be adapted in order to identify the relative attractors sitting on it, 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 lowdimensional chaotic sets. In the hydrodynamic literature, edge states have proven crucial to unfold bifurcation diagrams [33,34], for dynamic control [67] and because of their role in turbulence nucleation [71,72]. Until now edge states have been routinely identified numerically using standard bisection methods coupled with prior knowledge of lower and upper bounds for the location of the edge states based on an appropriate observable [16,19]. The present section introduces a new class of bisection method based on LCS diagnostics, and we will contrast the related concepts of global vs. local bisection methods in terms of their ability to identify the edge state, their scalability and their associated computational cost.

A. Global vs. 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 F t0+τ t0 (x 0 ) remains in (α A , α B ) for τ as large as possible. Such a method warrants that F t0+τ t0 (x 0 ) converges towards 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 a hard and a soft edge manifold.
Global methods can however fall short whenever no available global observable can elucidate on which side of the edge the different trajectories evolve [29,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 : • No preliminary choice of observable is needed based on physical intuition • No prior knowledge of bounds α A and α B is needed • The approach to edge tracking is not binary, edge trajectories are labelled as optimisers of an LCSbased observable.
• The edge tracking has a defined 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.
• 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 however to be nuanced in practice: the shorter the time horizon τ the more local the approach.
For larger values of τ , non-convergent trajectories are not bound to local neighbourhoods and can explore remote parts of the state space. The analysis of the LCS diagnostics in the previous section suggests long time horizons τ for sharper idenfication 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 by 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 optimiser 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 maximising λ t0+τ t0 (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 Section III point towards 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. 9 where one-dimensional landscapes of M + are displayed. FTLE landscapes are shown for comparison in Fig. 10(b), 4(c) and 5(c).
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 from Eq. (15). For the sake of generality, the LD-based edge tracking is yet 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 maximisation is always   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 (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. (16) using p = 0.5. Fig. 10(a) shows the minimum distance D min to the edge state for the classical, the FTLE-based and the LDbased 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. Fig.  10(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 . Fig. 10(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 [43,49]. As a consequence, the trajectory in the neighbourhood 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 ∇F t t0 . The function M + is hence trivially generalised to arbitrary dimension according to (7), as demonstrated in (19). The full computation of all FTLEs, on the other hand, becomes unfeasible in higher dimension since ∇F t t0 is of prohibitive size n × n, as pointed out in [47]. 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 [47]. For r=1 this amounts to O(n), the computation of the largest FTLE is hence based on a scalable algorithm as well.

V. CONCLUSION, DISCUSSION AND OUTLOOKS
Motivated by the physical problem of subcritical laminar-turbulent transition to turbulence in hydrodynamics, we have investigated the established notion of edge manifold which charts the state space into different basins of attraction [11]. The main mathematical contribution of this study is the generalisation of the concept of Lagrangian Coherent Structure (LCS) to arbitrarily high-dimensional dissipative dynamical systems and PDEs. The edge manifold divides the state space in the same way as a hyperbolic repelling Lagrangian coherent structure, and hence a direct relation between Lagrangian coherent structures 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 generalisation is that only forward-time 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) [36] and forward-time Lagrangian Descriptors (LDs) [40], 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 [56], ii) a bistable six-dimensional model where the edge state is an unstable periodic orbit and turbulence corresponds to a stable periodic orbit [60], iii) a nine-dimensional model where the edge state is a periodic orbit and turbulence is represented by a chaotic saddle [23], and eventually iv) well-resolved Navier-Stokes simulations of a periodic cell of plane Couette flow, understood as a O(10 5 )-dimensional dynamical system [67]. 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 and a finite-time formalism [75]. 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 either FTLEs or 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. A hitherto unreported drawback of LDs, worth a more systematic characterisation, is that their ability to highlight the edge manifold depends on the dynamical nature of the invariant sets of the system. The FTLEs represent a computationally more expensive tool. They have a solid mathematical ground, are objective diagnostics [44] and, although their interpretation as LCS is mainly heuristic [43], they have a more direct interpretation in terms of material barriers. 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 optimisers 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 [22,29,73] using this new framework is strongly encouraged.