Hyperbolic Covariant Coherent Structures in two dimensional flows

A new method to describe hyperbolic patterns in two dimensional flows is proposed. The method is based on the Covariant Lyapunov Vectors (CLVs), which have the properties to be covariant with the dynamics, and thus being mapped by the tangent linear operator into another CLVs basis, they are norm independent, invariant under time reversal and can be not orthonormal. CLVs can thus give a more detailed information on the expansion and contraction directions of the flow than the Lyapunov Vector bases, that are instead always orthogonal. We suggest a definition of Hyperbolic Covariant Coherent Structures (HCCSs), that can be defined on the scalar field representing the angle between the CLVs. HCCSs can be defined for every time instant and could be useful to understand the long term behaviour of particle tracers. We consider three examples: a simple autonomous Hamiltonian system, as well as the non-autonomous"double gyre"and Bickley jet, to see how well the angle is able to describe particular patterns and barriers. We compare the results from the HCCSs with other coherent patterns defined on finite time by the Finite Time Lyapunov Exponents (FTLEs), to see how the behaviour of these structures change asymptotically.


I. INTRODUCTION
In the paradigm of chaotic advection, the trajectories of passive tracers can be complex even when the velocity field of the flow is simple. This is the case, for example, for time dependent two-dimensional flows or even steady three-dimensional flows, like the celebrated ABC flow 1 . However, even flows with complicated time dependent structure allow for the formation of coherent patterns that influence the evolution of tracers. These structures are common in nature, appearing both at short and long time scales as well as small and large spatial scales. Remarkable examples of these structures are eddies and jets in the ocean and atmosphere, the Gulf Stream current, and ring clouds 2 . These patterns can influence, for example, the evolution of nutrients as well as oil spills and other pollutants. Furthermore, these coherent structures could act as local inhibitor for the energy transfer between scales 3 .
They also appear in other planets and other astrophysical systems, such as the jets on the surface of Jupiter, Saturn and other gaseous planets, and in the solar photospheric flows 4 .
These structures, often referred as Lagrangian Coherent Structures (LCSs), can shed light on the mixing and transport properties of a particular system on a finite time interval. The term "Lagrangian" in their name is motivated by the fact that they evolve as material lines with the flow 2 . A particular kind of LCS is called Hyperbolic LCS (HLCS), and can be seen as the locally most attracting or repelling material lines that characterize the dynamical system over a finite time interval. The term hyperbolic is just an analogy with the stability of fixed points in dynamical systems, since usually one wants to study non-autonomous systems for which entities such as fixed points or stable and unstable manifolds are not defined. Moreover these systems are studied for finite time.
Early attempts to detect HLCSs were based on Finite Time Lyapunov Exponents (FTLEs), that measure the rate at which initial conditions (or, equivalently, tracer particles) separate locally after a given interval of time [5][6][7][8][9][10][11][12] . One of the first rigorous definitions of the LCSs was based on ridges of the FTLEs 10 . FTLEs and the so defined LCSs were used to study for example Lagrangian dynamics of atmospheric jets 13 and oceanic stirring by mesoscale eddies [14][15][16][17][18][19] , describing for example the chaotic advection emerging by mixed layer instabilities and its sensitivity to the vertical shear 20 . A similar method to detect coherent structures uses the Finite Size Lyapunov Exponents (FSLEs), that represents the of evolution [45][46][47] . The tracer particles could be maximally repelled for a short time, but on long time they could have a different behaviour. Is it natural to wonder what happen asymptotically to the tracer particles? Is it possible to find some coherent structures that suggest the asymptotic behaviour of the tracer particles?
In this work we thus propose an alternative method to detect coherent patterns emerging in chaotic advection, which is based on the Covariant Lyapunov Vectors (CLVs). CLVs were first introduced by Oseledec 48 and Ruelle 49 , but for a long time they have received very little attention due to the lack of an efficient algorithm to compute them. Only in the last decade the computation of such vectors has become possible [50][51][52] , and CLVs have been used to investigate e.g. the motion of rigid disk systems 53 , convection 54 , and other atmospheric phenomena [55][56][57][58] . For other theoretical discussions or for reviews on CLVs, see e.g. Refs. [59][60][61][62][63] . Unlike the Lyapunov Exponents (LEs), which are time independent, the correspondent vectors, also known as forward and backward Lyapunov Vectors (LVs), do All this intrinsic information can be summarized using the angle between the CLVs, a scalar field that allows to investigate the spatial structures of the system. The need to pay more attention to the directions between LVs, backward and forward, has also been suggested for the study of turbulence 9 and for the definition of a diagnostic quantity for the study of mixing, the Lyapunov's diffusion 64 . Using three simple examples, we show that the attracting and repelling barriers tend to align along the paths on which the CLVs are orthogonal. The directions of the CLVs along these maxima provide thus information on the attracting or repelling nature of the barriers and can be related to the geometry of the system. Using the CLVs we can define structures, at a given time instant, that are asymptotically the most attractive or repelling. Furthermore, since these structures can be defined for every time instant, it is possible to follow the formation of coherent structures during the evolution of the flow.
In section II we discuss the theory behind the CLVs, and suggest a definition for coherent structures that give asymptotic information. The strategy will be to make use of the scalar quantity defined by the angle between the CLVs to locally identify the structures that asymptotically are maximally attractive or repulsive. In section III we use the CLVs to identify particular patterns in three different systems and we compare the results with FTLEs fields. Finally, in section IV we summarize the conclusions.

II. COVARIANT LYAPUNOV VECTORS
In this section we summarize the theory behind CLVs for two dimensional flows. For a more general and detailed review see for example Ref. 61 . Let the open set D ⊂ R 2 be the domain of the flow, t ∈ R the time and v(x, t) a vector velocity field in D. The dynamical system that describes the motion of a tracer advected by the flow is thus x(x 0 , t 0 ; where x(x 0 , t 0 , t) ∈ D is the trajectory of the tracer starting at the point x 0 at time t 0 . To that maps the initial position x 0 at time t 0 to the position x(x 0 , t 0 , t) at time t. It should be noted that the dependence on the initial condition is very important here, since the vectors will be considered as function of time and of the initial positions. In the following, the contracted form x = x(x 0 , t 0 , t) will be used.
At each point x ∈ D we can identify the tangent space T x D ⊂ R 2 . Infinitesimal perturbations, u(t) ∈ T x D, to a trajectory of this system can be described by the linearized system where J(t) ∈ R 2×2 is the Jacobian matrix composed by the derivatives of the vector field v(x, t) with respect to the component of the vector x. Using the fundamental matrix M(t) ∈ R 2×2 , of (3), that satisfies we define the so called tangent linear propagator F(t 0 , t) maps a vector in x 0 at time t 0 into a vector in x at time t along the same trajectory of the starting system (1), that is According to (5), the propagator is always nonsingular. In terms of the flow map, the tangent linear propagator is Exploiting Oseledec's Theorem 48,65 , it is possible to characterize the system using quantities that are independent on t or t 0 . By virtue of this theorem the far-future operator and the far-past operator are well defined quantities. Note that the product F (t 0 , t)F(t 0 , t) determines the Euclidean norm of the tangent vectors in the forward-time dynamics (a similar role it is played by Operators (8) and (9) probe respectively the future and past dynamics of a certain point, and share the same eigenvalues that, assuming ergodicity, are independent on time and space. Each eigenvalue has multiplicity m i (m 1 + m 2 = 2). Their logarithms correspond to the LEs of the dynamical system (1). If the limits in (8) and (9) are not considered, the resulting eigenvalues are time and space dependent and are called FTLEs.
The two operators (8) and (9) can be evaluated at the same point in space at a given time To overcome these issues, one can build particular spaces, the backward and forward Oseledec subspaces, defined as 48 and In the forward dynamics, the generic vector l + i (t) grows or decays exponentially with an average rate λ i . If the system is evolved forward in time, by means of the tangent linear propagator the evolution of the vector l + 1 (t) will have a non-zero projection inside the space generated by l + 1 (t ), but it will also have a non-zero projection onto the space generated by l + 2 (t ). On the other hand, l + 2 (t) will be transported onto the space generated by l + 2 (t ) and will have zero projection onto the space generated by l + 1 (t ). Repeating similar arguments for the backward Lyapunov Vectors leads to the observation that L − i and L + i are covariant subspaces, Vectors that are covariant with the dynamics and invariant with respect to time reversal will be found at the intersection of the Oseledec subspaces, i.e. at It should be noted that CLVs posses thus the properties of a semi-group. These vectors have an asymptotic grow or decay with an average rate λ i , so that their asymptotic behaviour can be summarized as where τ = |t − t|, which shows their invariance under time reversal. Note also that these vectors do not depend on any particular norm.
Equations (12), (13) and (15) imply a simple relation between CLVs and LVs in a two dimensional system, Since the CLVs highlight particular expansion and contraction directions at each point of the coordinate space, and these directions are not necessarily orthogonal, they can be used to understand the geometrical structure of the tangent space. This geometric information can be summarized by the scalar field of the angle θ(t) between the CLVs. Because the CLVs identify asymptotically the expansion and contraction directions sets of the tangent space associated at every point of the domain, the θ field represents a measure of the hyperbolicity of the system, that is a measure of the orthogonality between these two directions. Note that the orientation of the CLVs is defined as arbitrary. The angle between the CLVs, is thus where {w 1 (t), w 2 (t)} are the first and the second CLVs 51 . It is interesting to point out that when θ(t) = π/2, the CLVs reduce to LVs in two dimensions, and the backward and forward Lyapunov basis coincide. If the computation of the angle is done at every point of the domain, and if we consider the system (1) in which the initial conditions are varied in such a way that all the domain is spanned, we can build a field of the orthogonality between the expansion and contraction directions of the system.
In the following section we will show how CLVs, with their capability of probing the geometric structure of the tangent space, can be used to determine coherent structures that give asymptotic information on the tracer and then, how the CLVs highlights the mixing template of the flow.

A. Hyperbolic Covariant Coherent Structures
Several frameworks are available to study passive scalar mixing. Although all these methods aim to describe the mechanism underlying the chaotic advection, they have significantly different approaches. However, most of them share a fundamental feature called objectivity.
Objectivity is a fundamental requirement to define structures emerging from a flow. In practical applications, objectivity can be used e.g. to move the reference frame into the reference frame of a coherent structure. In particular a structure can be considered objective if it is invariant under a coordinate changes of the form where Q denotes a time-dependent orthogonal matrix and P a time-dependent translation 2 .
The FTLEs, FSLEs and geodesic theory define objective quantities. Sala et all. 62 have shown that, generally, also for a linear transformation of coordinates the new angle nonlinearly depends on the angle of the old reference system. However, for the particular class of transformation we are interested in, (21), the angle θ is an invariant. This means that structures highlighted by θ are objective.
To study the asimptotically most attractive or repelling behaviour of tracer particles near a particular structure identified at a given instant of time, one can consider the lines r(s, t), with s representing the length parameter, defined by where the primes indicates the derivative with respect to s, and characterized by θ = π/2 along their paths. In fact, in this case where, l 1 and l 2 are Lyapunov vectors. In these circumstances the backward and the forward basis are coincident. A sphere of initial conditions around a point of the path will be deformed in an ellipsoid whose axis are aligned with these Lyapunov vectors. If (22) are aligned along a ridge of the hyperbolicity field θ, characterized by θ = π/2, then they are also pointwise the most attractive or repelling lines in terms of the asymptotic behaviour of tracer particles. Consider in fact line (22a) and its tangent CLV w 1 as depicted in Figure 1. Along this curve, for any point P 1 we can chose a point P 2 which is on the normal vector to the curve in P 1 , n. Because of (11), the distance between these two points decreases as Consider now a point P 3 on a nearby curve, and a point P 4 laying on the normal vector to the curve in P 3 . The normal vector n at P 3 can be written as a linear combination of w 1 and w 2 and, considering again (11), one has thus The previous computation can be repeated analogously for w 2 . The same reasoning holds for w 2 .
It must be remarked that (22b) gives asymptotic information, but is defined for a particular instant of time. For every time instant the lines (22), together with the orthogonality condition θ(t) = π/2, describe the so called tensorlines, and can be seen as the asymptotic version of shearless barriers 44 .
Taking into account these properties of the CLVs, we propose a definition to describe these asymptotic coherent patterns:

Definition II.1 (Hyperbolic Covariant Coherent Structures (HCCSs)). At each time t, a
Hyperbolic Covariant Coherent Structure, is an isoline of the hyperbolicity scalar field θ at the level θ = π/2. Its attractive or repelling nature is determined by the CLVs aligned with it.
It is important to emphasize the difference between the information provided by FTLEs CLVs therefore allow the definition of instant structures, which, however, give asymptotic information.
This approach presents some numerical challenges, CLVs are not fast to be computed as the FTLEs, and numerically, in order to identify the isolines, it is necessary to consider an interval of values for the angle. The HCCSs are thus computed here considering isoregion , where δ is 0.087 rads.

III. NUMERICAL EXAMPLES
The algorithm used in this work is the same presented in Ref. 51 , so only a schematic of its structure will be presented here. This algorithm can compute a large number of CLVs converging exponentially fast when invertible dynamics is considered. The basic idea is that, if the backward Lyapunov vectors basis {l − 1 (t), l − 2 (t)} is evolved by means the the operator F, it is always possible to keep it orthonormalized with a QR decomposition and store the corresponding upper triangular projection matrix. It is then possible to exploit the informations contained in these upper triangular matrices to evolve backward in time arbitrary vectors that will converge to the CLVs 51 . This is done in five different phases, which are described in Appendix A.
In the next we investigate the hyperbolicity field θ for three different examples, which include one autonomous Hamiltonian flow and two non-autonomous two dimensional flows.
The HCCSs emerging from the angle between the CLVs are compared with the FTLEs field, while their attractive or repelling nature will be discussed in terms of the CLVs. The so that the dynamical equations of motion are  and repels every close trajectories. We apply the CLVs theory to see if the angle between that vectors is able to detect such kind of structure.  Along the x axis, θ = π/2 remains constant for all the time. Consider in fact the propagator (7) where Along y = 0, both the tangent linear propagator and the CGT (A2) are diagonal. If at time t one has θ(t) = π/2 for y = 0, then where ξ 1 and ξ 2 are the eigenvectors of the CGT, and at time t > t This shows that y = 0 is a repelling HCCSs.

B. Double gyre
The previous analysis is now applied to a with initial condition x(0) = x 0 , y(0) = y 0 . The streamline of the flow at t = 5 are shown in Figure 2b.
We consider two numerical experiments, DV1 and DV2, in which the CLVs are computed in two different time intervals with a non-null intersection. During DV1 we compute the CLVs in the interval t = [5,10] and during DV2 for t = [10,20]. This is important, as will be shown in this section, to highlight that θ is a quantity that depends on a particular time instant and not on the time interval considered for its computation, as long as the interval considered is sufficiently long. Figure 4 shows a comparison between the angle θ and the FTLEs for this system. The first column is relative to DV1, while the second column is relative to DV2.   Figure 5 shows the repelling nature of the HCCSs in the zoomed region, since the contour is aligned with the second CLVs, w 2 , that define the contraction direction. Every particle of the tracer near these HCCSs will tend asymptotically to get away from them.     It is also interesting to consider the evolution of the PDF of the hyperbolicity field, P (θ, t).
If the distribution is peaked around π/2, the expansion and contraction directions are almost everywhere perpendicular between each other. If the distribution is peaked around zero, the expansion and contraction directions are almost tangent everywhere. If the distribution is flat there is not a clear correlation between the expansion and contraction directions. This is the reason why we do not see in this interval a periodicity in the moments of the PDF for the variable θ.

C. Bickley jet
The Bickley jet is an idealized model of a jet perturbed by a Rossby wave [66][67][68] . The velocity field is given in terms of the stream function, ψ(x, y, t), that can be decomposed as the sum of a mean flow ψ 0 (x, y, t) and a perturbation ψ 1 (x, y, t) where and Following the work by Onu et al. 66 , the forcing is chosen as a solution that runs on the chaotic attractor of the Duffing Oscillator From now on the time will be scaled with the quantity L x /U . The streamline of the flow at t = 1.89 are shown in Figure 2c.
As for the double gyre, we show two experiments, BJ1 and BJ2, for which the end of the CLVs computation in the first experiment correspond with the first CLVs computation for the second experiment. The first column of Figure 9 shows the analysis for BJ1 and the second column for BJ2. Figures 9a and 9b  The left panel of Figure 10 shows the superposition of HCCSs computed for t = 3.79 for the experiment BJ2, green lines, and the FTLEs computed in the same experiments ( Figure   9b). This figure remarks the fact that the central jet is not completely a HCCSs. HCCSs are also present at the border of the vortices and in the center part of the upper vortices.
The right panel of Figure 10 show an enlargment into a portion of the central jet containing HCCSs and the CLVs in that part of the domain. In this picture it is possible to appreciate the repelling nature of the HCCSs looking the alignment of the contour with the second CLVs. It is interesting to point out that, looking at the FTLEs field, it is not possible to see ridges in the central jet. In this case it is not possible to find LCSs, if we use definition of LCSs based on the FTLEs field 10,37 , and this once again remark the difference between the HCCSs and LCSs. Figure 11 shows the comparison between the angles at the end of the CLVs computation interval of the experiment BJ1 with the one computed at the beginning of the CLVs calculation interval for the BJ2 experiments. As for the double gyre, the comparison between Figure 11a and Figure 9d shows that regions characterized by θ = π/2, i.e. the HCCSs, have the same structure. The difference between the two fields (Figure 11b), shows that these regions are characterized by smaller errors. In contrast, the regions with low hyperbolicity of Figure 11a appear to be more noisy with respect to the same regions computed for the BJ2 simulation (Figure 9d). In agreement with this, the difference between the two fields shows that the low hyperbolicity regions are characterized by larger errors (Figure 11b). Figure 11. As in Figure 7, but for the BJ2 experiment.  For the two non-autonomous systems we have considered the evolution of the PDFs of the θ field and the evolution of its first four moments. For the Bickley jet the probability distribution of the angle is stationary in time, and so its moments, but for the double gyre it is possible to appreciate a small variation in time for the PDF. The information deriving by the evolution of the θ field, related to the variation of the strength of the hyperbolicity field, could be used to characterize the dynamical mixing of the system.
Finally, future studies will have to address the detection of hyperbolic structures beyond analytical systems, i.e. for two-dimensional turbulent flows. This will be particularly interesting for flows at the transition between balance and lack of balance [69][70][71][72] , where the detection of HCCSs can shed light both on the structure of mixing and on the forward cascade of energy to dissipation, or the lack of thereof. In particular, the evolution of the moments of the PDF of θ could be used to define an index of dynamical mixing for the system under study. Particularly interesting will be the behaviour of the HCCSs in presence of intermittency. We can conjecture that in particular cases, such as e.g. the merging of two vortices, the instantaneous structures underlined by the HCCS can give reliable information of the asymptotic tracer dynamics, that is the dynamics after the merging event. This is however left for future studies.
2. Forward Transition (T2 ): the backward Lyapunov bases are evolved from time t to t . The evolution is done with the help of (1), (6) and the QR decomposition for every evolution step. We indicate with X(t k ) the matrix which columns contain the new bases at time t k , and with R(t k−1 , t k ) the correspondent upper triangular matrix.
During this step both the local Lyapunov bases and upper triangular matrices are stored. The diagonal elements (R(t k−1 , t k )) ii of the upper triangular matrices give information about the local growth rates of the bases vectors at a given time t k , and they are used to compute the FTLEs as a time average where N time step are considered between t and t . If the LEs exist, (A1) will converge to them for a sufficiently long evolution. It is worth to note that sometimes the FTLEs are computed using a different method, that is using the so called Cauchy Green Tensor (CGT) defined as This operator is also known as deformation tensor, whose eigenvalues µ i (t 0 , t), and eigenvectors ξ i (t 0 , t), satisfy From a geometric point of view, a set of initial conditions corresponding to the unit sphere is mapped by the dynamics into an ellipsoid, with principal axis aligned in the direction of the eigenvectors of the CGT and with length determined by the correspondent eigenvalues. The eigenvalues of the CGT determine the FTLEs as where the dependence on the starting position has here been suppressed. This second method for the computation of the FTLEs exhibits some problems, for example, if just one finite local growth rate is taken into account considering a large time interval, (A4) teds to zero and not to the LEs. The first method should be preferred. bases. Using the stored matrices R of the step 3, these matrices are evolved backward in time, until time t , using the following relation: where t n and t n+1 are time step between t and t . This method uses all the information contained in R and not just the diagonal part of the matrix. The D diagonal matrices contain the column norm of C. Using Eq. (A5) it is possible to show 51 that the generic vectors chosen will be aligned with the CLVs. Note that when a trajectory passes close to tangency of an invariant manifold, the matrices C can be ill-defined, and so a little amount of noise on the diagonal element of these matrices, or an average on the diagonal elements of the neighbor matrices is used to correct the problem. bases. In Figure 13 it is shown the convergence for the spatial average of the scalar product between the starting random bases chosen.   Times have been rescaled with L x /U .

Bickley jet
The