Mesochronic classification of trajectories in incompressible 3D vector fields over finite times

The mesochronic velocity is the average of the velocity field along trajectories generated by the same velocity field over a time interval of finite duration. In this paper we classify initial conditions of trajectories evolving in incompressible vector fields according to the character of motion of material around the trajectory. In particular, we provide calculations that can be used to determine the number of expanding directions and the presence of rotation from the characteristic polynomial of the Jacobian matrix of mesochronic velocity. In doing so, we show that (a) the mesochronic velocity can be used to characterize dynamical deformation of three-dimensional volumes, (b) the resulting mesochronic analysis is a finite-time extension of the Okubo--Weiss--Chong analysis of incompressible velocity fields, (c) the two-dimensional mesochronic analysis from Mezic et al. \emph{A New Mixing Diagnostic and Gulf Oil Spill Movement}, Science 330, (2010), 486-489, extends to three-dimensional state spaces. Theoretical considerations are further supported by numerical computations performed for a dynamical system arising in fluid mechanics, the unsteady Arnold--Beltrami--Childress (ABC) flow.


1.
Introduction. Chaotic advection is the theory of material transport in fluids based in dynamical systems theory. [3,51,64] It is largely rooted in analysis of geometric structures in flows with a simple time-dependence, time-autonomous or periodic. Since the realization that even flows with a complicated time dependence, e.g., turbulent flows, possess organized Lagrangian structures, it has become increasingly important to detect geometric structures analogous to invariant manifolds in steady flows. In particular, detection of structures using only finite-time information about the flow has been seen as the most practically-useful direction. [2,20,27,31,55,58,59] The need for analysis of geometric structures that organize advection is not purely academic. Transport of material by fluid flows played a crucial role in the fallout from several recent catastrophic events, namely, the volcanic eruption of Eyjafjallajökull (2010), the Deepwater Horizon oil spill (2010), and the nuclear disaster in Fukushima (2011). These events highlighted how important it is to detect organizing geometric structures in (near) real time from data, either measured or generated by detailed simulation models; consequently, such problems have become a very active intersection of dynamical systems and fluid dynamics.
The problem of identifying geometric structures in flows has resulted in several approaches, each focusing on a somewhat-different structure as the objective of its analysis.
The theory of Lagrangian Coherent Structures [32] (LCS) identifies barriers that organize the transport in flows with complex time-dependence. Initially, LCS were closely associated with computation of Finite-Time Lyapunov Exponent fields [59]; more recently, they have been re-formulated using a variational principle [5,28,30], which defines them as certain geodesic lines of the local deformation field induced by the fluid flow. This new definition allows a finer classification of LCS, both in two and three dimensions, based on the type of deformation, e.g., hyperbolic, elliptic, corresponding to different behaviors of fluid parcels in the flow. The recent review by Haller [29] gives a detailed coverage of the current state-of-the-art of techniques centered around LCS.
LCS and associated theories mostly focused on magnitude of non-rotational deformation in the flow. The rotation deformation has been classically studied by Poincaré in topological dynamics and Ruelle [57] in ergodic theory; recently, a Finite-Time Rotation Number [60] has been proposed as a useful quantity in analysis of flows. At the closing of this manuscript we were also notified of the recent work by Farazmand and Haller [15], working along the same lines.
In the effort to study the "stretch-and-fold" mechanism for chaos in finite time, most studies focus on the first-order "stretch". Folding in finite time has not received direct attention; an exception is the study of the Finite-Time Curvature Fields [36][37][38]. At this point, structures observed through all these methods have been connected mostly on phenomenological basis, through comparison of visualizations, showing considerable overlap between observed structures but, also, nonnegligible differences.
Magnitudes of the local material deformation are typically estimated by processing velocity gradients; they cannot be precisely computed in the absence of the detailed data about the velocity field, e.g., when the system is sampled by sparse trajectories only. In sparsely-sampled planar systems, trajectories can be represented by space-time braids -extremely-reduced, symbolic representations of trajectories. The resulting approach, known as braid dynamics [2,6,61] has been successful in providing lower bounds on the amount of topological deformation present in the flow, in limited-data settings. The obtained bounds have been used both in design and analysis of the material advection; unfortunately, there are currently no extensions of braid dynamics to three-dimensional flows.
Instead of looking for barriers to transport, as is the case with the LCS theory, the theory of almost-invariant sets identifies a collection of sets, fixed in space, such that the material placed in them does not leak out. These sets act as routes for the material transport. The approach is based on the Perron-Frobenius transfer operator, which models how the flow moves distributions of points, instead of individual trajectories. The Perron-Frobenius operator is always infinite-dimensional and linear; the identification of almost-invariant sets is then intimately connected with approximating its eigenfunctions. While the Perron-Frobenius operator has been a staple of the ergodic and probability theory since the early 20th century, it was introduced to the applied, non-probabilistic context by Dellnitz and Junge [11,12] as the basis for identification of invariant sets in time-invariant dynamical systems. Since then, the theory has been expanded to include detection of almost-invariant sets of autonomous systems [17,19], and flows with more general time dependencies [18,20].
Spatial invariants of dynamical systems relate to infinite-time averages of functions along Lagrangian trajectories. This connection between the ergodic theory and applied mathematics was initially developed in Ref. [42]. Based on these ideas, Ref. [55] proposed that even finite-time averages of functions can enable detection of geometric structures important for fluid transport, which broadly constitutes the mesochronic, i.e., time-averaged, theory of transport in fluids. The utility of timeaverages has been corroborated on numerical and experimentally-realizable flows with simple time dependence [39-41, 46, 47].
The mesochronic techniques have developed in two directions. One focuses on computations of ergodic invariant sets using long-time averages of a large set of averaged functions. [8,35,42,47] The other, which we follow here, does not aim to compute all invariant sets; rather it uses much shorter averages of the velocity field itself to identify the character of the deformation, i.e., presence or absence of rotation. [45] This is in contrast with mentioned LCS and rotational theories, which describe deformation starting from analysis of the magnitude of the deformation.
Before we dive into calculations, we give a short explanation of the approach. Introductory courses in dynamical systems often discuss the stability of a stagnation point p in a planar systemẋ = f (x) by looking at its linearizationξ = ∇f | p · ξ around a fixed point p whose stability depends on positions of two eigenvalues of the Jacobian ∇f | p . Instead of computing the eigenvalues, their locations can be inferred from the trace and the determinant of ∇f | p . If the flow is incompressible, the trace is zero, so the determinant alone is needed for the full stability analysis. For unsteady systems this analysis may not always hold; however, Ref. [45] showed that even then similar results can be obtained by looking at the Jacobian matrix of the velocity averaged over a Lagrangian trajectory, termed the mesochronic Jacobian matrix. Away from fixed points, the calculation does not compute stability, but rather the spectral class of the Jacobian: hyperbolic (strain) or elliptic (rotation), termed mesochronic classes. Applied to prediction of the oil slick transport in the aftermath of the Deep Water Horizon spill [45], mesochronic analysis showed that regions of mesohyperbolicity correspond to jets which dispersed the slick, while mesoelliptic zones correspond to centers of eddies in which the slick accumulated.

BUDIŠIĆ, SIEGMUND, SON AND MEZIĆ
The main contribution of this paper is the extension of the mesochronic analysis to three-dimensional flows.
There are several connections of the mesochronic analysis with other approaches.
• On a fundamental level, averages of functions along trajectories are intimately related to spectral properties of the Koopman operator [34,43,44], which is adjoint to the mentioned Perron-Frobenius operator. • Greene [25,26] defined the residue criterion in order to predict the order of destruction of Kolmogorov-Arnold-Moser (KAM) tori in perturbed Hamiltonian maps. The computation of the residue is almost identical to that of the mesohyperbolicity indicator for 2D flows. The three-dimensional version of the residue criterion [16] also resembles mesochronic indicators introduced here. • An analysis of oceanic flows based on the Jacobian of the instantaneous, i.e., non-averaged, velocity is well-known as Okubo-Weiss-Chong partition [9,50,63]; this is the limit of the mesochronic theory as the averaging time T → 0. • In the other extreme, as T → ∞, real parts of eigenvalues of the mesochronic Jacobian relate to Lyapunov exponents and rotation numbers). It can be shown under generic conditions, using the polar decomposition of ∇ψ T that the limit of the eigenvalues of the gradient of the flow map are the Lyapunov exponents (see the heuristic discussion in [22], that can be turned into rigorous proofs using ergodic-theoretic techniques in [57], see also [1,Chapt. 3, §9]). The paper is organized as follows. In Section 2 we introduce the precise definitions for dynamical systems we are considering, review the basics of differential geometry needed, and re-state the Okubo-Weiss-Chong analysis in these terms. Section 3 contains the main result, the 3D mesochronic classification, while Section 4 makes connections to Lyapunov, Okubo-Weiss-Chong, and 2D mesochronic analyses. In Section 5 we illustrate the technique by a set of analytic and numerical examples; in particular the steady and unsteady Arnold-Beltrami-Childress flow [13]. Numerical details are given in the Appendices. The paper closes with the discussion in Section 6.

2.
Preliminaries. Consider a time-varying differential equatioṅ For an initial condition (t 0 , x 0 ) ∈ D let t → ϕ(t, t 0 , x 0 ) denote the solution (or trajectory) of the initial value problem (1), with x(t 0 ) = x 0 . Throughout the paper we assume for an arbitrary but fixed initial time t 0 ∈ R, duration T > 0 and open set of initial values usually called Poincaré map if the equation is periodic, i.e., if for the chosen T , We are mainly interested in finite-time dynamics for a fixed duration T > 0 but will also investigate the instantaneous (in the zero-time limit T → 0 + ) and asymptotic (in the infinite-time limit T → +∞) dynamics, assuming the solution ϕ(·, t 0 , x 0 ) exists on [t 0 , ∞).
Observables of (1) are continuous functions F : X → R n , which are evaluated along arbitrary solutions t → x(t) on I. They are used to model physical measurements of a state of the system, e.g., the time trace t → F (t, x(t)) might represent the ocean temperature recorded by a sensor as it is passively carried by ocean currents along the trajectory x(t). A time average or trajectory averageF T : is a function that depends on the initial value x(t 0 ) = x 0 of the trajectory x(t) at time t 0 . Trajectory averages depend on the duration T > 0 and can be analyzed from the instantaneous, asymptotic, or finite-time perspective. The instantaneous case is the most obvious, as lim T →0FT (x 0 ) = F (t 0 , x 0 ), e.g., if F is a component of the velocity field itself. Certain choices of F can still provide valuable information, as we explain in the next paragraph. Asymptotic analysis studies ergodic averages, i.e., limitsF ∞ (x 0 ) := lim T →∞FT (x 0 ), in case they exist. Ergodic theory analyzes limitsF ∞ of observables F which are specified only in general terms, e.g., only by the space of functions from which they are drawn. Even in such general cases valuable information can be recovered, e.g., on time-invariant measures on the state space. [8,47] On the other hand, choosing a particular observable can provide us with more detailed information. Since the components of the velocity field f (t, x) = [f 1 (t, x), f 2 (t, x), f 3 (t, x)] are themselves continuous functions on the time-state space X ⊂ R × R 3 , they are observables. One could argue that they are the most distinguished observables for analysis of dynamical systems, as they directly provide dynamical information about the behavior of the system. We adopt the velocity-as-observable viewpoint and analyze the time average of the velocity field, which was termed mesochronic velocity [45]. We note that the values of the mesochronic velocity appear as quantities of interest in [55] while [45] is the first to look into their gradients.
Note that the spatial derivatives and the averaging over trajectories do not commute, i.e., ∇f T is not equal to the average of the instantaneous Jacobian over trajectories.
The mesochronic velocity in the instantaneous limitf T T →0 + − −−− → f coincides with the velocity field f . The asymptotic limit for T → ∞ exists in many cases, for example, if the dynamical system is autonomous and volume-preserving on a compact domain.
We use the mesochronic velocity to determine the character of the evolution of a material volume (see Section 2.1) by an incompressible dynamical system, which satisfies the Liouville condition Although the limits T → 0 + and T → ∞ have been studied classically, neither theory is applicable to transient behavior. On the other hand, the finitetime analysis of the mesochronic velocity recovers the character of the time-T map ψ : X(t 0 ) → X(t 0 + T ), which captures transients at the time-scale T [45].
2.1. Deformation of a volume cell under a diffeomorphism. We now briefly review basic differential geometry that characterizes deformation under a volumepreserving diffeomorphism ψ using the spectral class of its Jacobian ∇ψ. This theory is later applied to time-T maps ψ T of dynamical systems over finite time intervals.
Let ψ : U → V be a diffeomorphism between two open subsets U ⊂ R 3 , V ⊂ R 3 , with the usual volume measure on R 3 . We are interested in deformation of an infinitesimal volume cell surrounding x ∈ U as x → ψ(x). The central object of our interest is the Jacobian matrix ∇ψ : U → R 3×3 . It is a basic result in differential geometry that volumes of a set S ⊂ U and its image ψ(S) ⊂ V are equal if and only if |det ∇ψ(x)| ≡ 1. We now restrict our attention to orientation-and volumepreserving maps ψ, i.e., maps for which det ∇ψ ≡ 1.
At the coarsest level, we distinguish between a hyperbolic deformation, when the volume cell is deformed along all three spatial dimensions, and the opposite, nonhyperbolic character. Let µ 1 , µ 2 , µ 3 ∈ C denote the eigenvalues of ∇ψ, assuming |µ 1 | ≤ |µ 2 | ≤ |µ 3 |. Different fields of mathematics may interpret presence or absence of hyperbolicity differently, e.g., as the material deformation in continuum mechanics, or stability of the map in dynamical systems and control. Since our analysis could include both domains, we refer to presence/absence of hyperbolicity, and their sub-classification (see below) as the spectral character of the diffeomorphism. Definition 2.2 (Hyperbolicity). The map ψ is hyperbolic at x if no eigenvalues µ 1 , µ 2 , µ 3 of the Jacobian ∇ψ(x) lie on the unit circle in the complex plane, i.e., ∀ i = 1, 2, 3, |µ i | = 1. Otherwise, it is non-hyperbolic at x. In particular, if all eigenvalues lie on the unit circle, i.e., ∀ i = 1, 2, 3, |µ i | = 1, the non-hyperbolic map is elliptic.
Depending on the complexity of eigenvalues, we further distinguish sellar (saddlelike) and helical (spiral/helix-like) character of the deformation under ψ. 1 Definition 2.3 (Sellar and helical deformation). The map ψ is sellar at x if all eigenvalues µ 1 , µ 2 , µ 3 of the Jacobian ∇ψ(x) are real and non-defective (their algebraic and geometric multiplicities match). If, instead, there is a pair of complexconjugate eigenvalues, the map is helical (spiral-like) at x.
Deformation at a sellar point exhibits three distinct spatial axes, meeting at x, whose directions are preserved under ψ. The directions of preserved spatial axes correspond to real-valued eigenvectors of ∇ψ, while the associated (real) eigenvalues µ determine whether the points along the axes are moving away from x, (|µ| > 1), moving towards x, (|µ| < 1), or remain neutral (|µ| = 1).
Deformation at a helical point exhibits only a single preserved spatial axis, around which a volume cell is rotated, resulting in a single real-valued eigenvalue. As
∇ψ(x) is a real matrix, any complex eigenvalues must arise in conjugate pairs, µ,μ. The modulus of the complex pair again determines expansion or contraction of the material in the rotation plane, while the real and imaginary components of the associated eigenvector pair span the rotation plane.
Since we restrict our analysis to volume-preserving maps ψ, existence of an eigenvalue inside the unit circle (contraction) necessarily means that at least one other eigenvalue lies outside of the unit circle (expansion), as |det ∇ψ| = |µ 1 µ 2 µ 3 | ≡ 1. All possible combinations are enumerated in Definition 2.4 in which we use the symbols + and − to denote expansion and contraction directions, respectively. Definition 2.4 (Spectral classes). Let ψ : U → V be a volume-and orientationpreserving diffeomorphism, x ∈ U a point, and ∇ψ(x) the Jacobian of ψ at x with eigenvalues µ i , ordered as |µ 1 | ≤ |µ 2 | ≤ |µ 3 |. The class of the point x is determined according to Table 1 by the number of eigenvalues of ∇ψ(x) that are inside and on the unit circle.
The first four classes are hyperbolic, whereas the remaining cases are non-hyperbolic. Informally, we will refer to signatures [−++] and [−−+] of hyperbolic points as, respectively, flattening and elongating, due to the shape of a volume cell after application of ψ, as sketched in Figure 1.

2.2.
Instantaneous deformation by a dynamical system. We now apply the classification outlined in Section 2.1 to the time-T map ψ T defined in (2) for a Figure 1. Sketches of deformation depending on the spectral class of a volume-preserving diffeomorphism ψ at the point x ∈ U (see Definition 2.4). Note that the initial and final axes do not need to be parallel in general. Pure shear and reflection (including the identity) are not sketched.
fixed time interval T to recover the Okubo-Weiss-Chong [9,50,63] criterion for classification of velocity fields. Instead of forming the classification of ψ T based on properties of ∇ψ T , we chose instead to base our approach on properties of the mesochronic Jacobian ∇f T , resulting in criteria whose expressions are consistent across all time scales T .
Instantaneous classification refers to the infinitesimally small time interval length T for which the class of x ∈ R 3 under ψ T can be inferred from the spectrum of the Jacobian ∇f (t 0 , x), through the connections with ∇ψ T . Both Jacobians are 3 × 3 matrices; for a general such matrix M , the characteristic polynomial is given by The coefficient tr Cof M is the cofactor trace, also called "second trace" in [16]; intuitively, it is the "variance" of the vector containing eigenvalues. Additionally, it follows from (3) Table 2. Okubo-Weiss-Chong classification based on signs of quantities in the first two columns (see Theorem 2.5).
As mentioned in Definition 2.4, spectral class of the time-T map is determined by locations of eigenvalues of ∇ψ T which, in turn, are roots of the characteristic polynomial P ψ (µ) = det (µId − ∇ψ T (x)).
Expanding ∇ψ(x) into a Taylor series around T = 0 yields for T ≈ 0. Consequently is the characteristic polynomial of the Jacobian matrix ∇f of the instantaneous velocity field. In all cases, Jacobian matrices are evaluated at the time-space point (t 0 , x) which is suppressed in notation. Incompressibility of the flow implies and therefore the spectral character depends only on the determinant d f and the sum of minors m f of the velocity field Jacobian ∇f . We now state the Okubo-Weiss-Chong criteria, using terminology specified in Definition 2.4.
Theorem 2.5 (Okubo-Weiss-Chong [9]). From the Jacobian of the velocity field (velocity gradient tensor) Hyperbolic points can be further classified into four subclasses, depending on signs of the determinant d f and the quantity 27d 2 f + 4m 3 f , as listed in Table 2. The analogous result for planar dynamics is known as the Okubo-Weiss criterion [50,63].
3. Mesochronic classification. The Okubo-Weiss-Chong criterion (Theorem 2.5) classifies the points in the domain of the time-T map ψ T according to the spectral character of the Jacobian ∇f of the velocity field. In general, the spectral character of ∇ψ T for T away from 0 will not be approximated well by the spectral character of ∇f , except in extremely slowly varying flows. In order to capture both deformation at small and large T , we replace ∇f by the Jacobian ∇f T of the mesochronic velocityf T , defined as the Lagrangian average of the velocity f over the duration To see why ∇f T (x) plays an important role in our analysis, we write the solution ϕ(t, t 0 , x) in its integral form which is a consequence of the fundamental theorem of calculus. The same integral appears both in (6) and (7), which we use to write the time-T map as and hence Comparing with (4), ∇ψ T (x) ≈ Id+T ∇f (t 0 , x) which is valid for T ≈ 0, we see that for general time intervals [t 0 , t 0 +T ] the mesochronic velocity Jacobian ∇f T (x) plays the same role as the velocity field Jacobian ∇f (t 0 , x) did for intervals of infinitesimal length.
In analogy to Okubo-Weiss-Chong, we classify x ∈ X(t 0 ) ⊂ R 3 under the action ψ T using only information obtained from the mesochronic velocity and its Jacobian ∇f T (x). Similarly, the mesochronic class of x is the spectral class of ψ T at x, as specified by Definition 2.4.
As explained in Section 2.1, the Liouville incompressibility results in only two, instead of three, "axes" in which we understand the spectral class of ψ T : 1. number of contracting directions, indicated by labels [− + +] and [− − +].
2. presence of a rigid rotation, indicated by helix vs. saddle split. The four mesohyperbolic classes are formed by choosing an option along each of the axes above, with the non-hyperbolic classes acting as boundary cases between them.
For conceptual and computational reasons we will formulate quantities Σ and ∆, each corresponding to one of the "axes" above, such that the signs of their values sort the finite-time dynamics around x into mesochronic classes. The number of contracting directions will be detected by Σ, the presence of rotation by ∆.  Table 3. Mesochronic classification based on signs of Σ and ∆ (see Theorem 3.2). The classes with Σ = 0 are mesohyperbolic. The class for Σ = ±∞ is the neutral (non-mesohyperbolic) saddle.
If Σ is finite, the point x is classified into mesochronic classes according to Table 3. Mesohyperbolic classes are those for which Σ is finite and Σ = 0, i.e., for which both The distinction between the shear and the boundary saddle class when ∆ = 0 depends on eigenvectors of the mesochronic Jacobian and cannot be made purely based on the spectrum.
The proof of the theorem will be the result of two lemmas: • Lemma 3.3 defines Σ and ∆ using coefficients of the characteristic polynomial of the time-T map Jacobian matrix ∇ψ T (x). • Lemma 3.4 establishes relations between characteristic coefficients of ∇ψ T (x) and characteristic coefficients of ∇f T (x).
If Σ is finite, then the point x is classified into mesochronic classes defined by Table 3. Mesohyperbolic classes are those for which Σ is finite and Σ = 0, i.e., t ψ − m ψ = 0 and t ψ + m ψ + 2 = 0. Remark 1. The trace t ψ and the determinant d ψ of the Jacobian matrix ∇ψ T are commonly encountered in matrix analysis due to their simple relationships to eigenvalues. To interpret the cofactor trace m ψ we can expand it as: Under an incompressible flow (det ∇ψ T ≡ 1) the cofactor trace m ψ is the trace of the inverse of the time-T map Jacobian, as discussed after (3): Additionally, we can clarify the meaning of Σ. Rewrite (12) as Then through simple algebraic manipulation it follows that Using m ψ = t ψ −1 and rewriting the expressions for t ψ using eigenvalues µ i of the flow map Jacobian, we obtain that Proof of Lemma 3.3. The mesochronic class of x is determined by eigenvalues of the Jacobian ∇ψ T (x), whose characteristic polynomial is The mesochronic class is determined by relative locations of zeros P ψ (µ * ) = 0 in reference to the unit circle and to the horizontal axis. As mentioned, converting the reference unit circle to a vertical axes simplifies the criteria, as they can then be computed solely from signs of eigenvalues. To this end, we introduce the conformal map Γ(s) = 1+s 1−s which maps the left half-plane in C to the inside of the unit circle in C, while preserving the upper half-plane. It follows that the location of zeros of the composite function P ψ • Γ with respect to axes is the same as the location of zeros of P ψ with respect to horizontal axes and the unit circle. Note that the inverse Γ −1 (µ) = µ−1 µ+1 has a pole at µ * = −1, so no finite zeros in the s-plane will correspond to the zero of P ψ at µ * = −1. For this reason, we will separately treat the case when ∇ψ T has an eigenvalue at −1.
Assuming P ψ (−1) = 0, the composite function is where the second equalities are obtained by incompressibility A point x is non-hyperbolic whenever one of the roots of the characteristic polynomial P ψ (µ) of ∇ψ T (x) is on the unit circle. If P ψ (−1) = 0, this condition implies that a purely imaginary number ia, a ∈ R is a zero of P ψ • Γ. Substituting into the numerator of P ψ • Γ, we obtain the condition n 1 n 2 − n 0 n 3 = 0 for hyperbolicity. This condition translates to a condition on the spectral coefficients of ∇ψ T : t ψ − m ψ = 0. The other case for non-hyperbolicity is P ψ (−1) = 0, which by direct substitution into P ψ translates to t ψ + m ψ + 2 = 0.
When P ψ (−1) = 0, we can proceed to a finer classification of the location of the roots of P ψ . Due to incompressibility, there either have to be two zeros of the polynomial P ψ inside a circle and one outside, or vice-versa. Under a conformal transformation Γ this condition translates into two zeros of P ψ • Γ having matching signs, while the third is of the opposite sign. It follows that the product of zeros s * i of P ψ • Γ is positive when two zeros are negative, corresponding to two directions of contraction under the flow map, while it is negative when there is a single contracting and two expanding directions. By using the zeros s * i of P ψ • Γ to factorize its numerator 3 i=0 n i s i = n 3 3 i=1 (s − s * i ) and equating the zeroth-order coefficients one obtains −n 3 3 i=1 s * i = n 0 . We therefore define an indicator Σ as By this argument, Σ = 0 implies hyperbolicity: Σ > 0 implies two directions of contraction, one of expansion, while Σ < 0 implies one direction of contraction, two of expansion. When Σ = 0 one of the roots of P ψ lies on the unit circle, which we term "neutral", while the other two directions are contracting and expanding. The presence of rotation is indicated by the complexity of zeros of P ψ • Γ, which is determined by the discriminant of the numerator, When D 3D > 0, all zeros are real and distinct, when D 3D < 0, two zeros are complex-conjugates of each other, when D 3D = 0, one real zero is repeated, while the third is distinct. We therefore define our second indicator ∆ to be It follows that ∆ < 0 indicates that all eigenvalues of ∇ψ T (x) are real, implying that there is no rigid rotation; ∆ > 0 indicates that two eigenvalues are complex, so rigid rotation is present. When ∆ = 0 two eigenvalues coincide somewhere along the real line.
Finally, we have to deal with the case when P ψ (−1) = 0 which is not covered by the s-complex plane as defined above, as it implies that the denominator of Σ is zero, i.e., t ψ + m ψ + 2 = 0. If one zero is µ * 2 = −1, then the two others have to be µ * 1 = z, µ * 3 = −1/z, for some z ∈ C to satisfy incompressibility. Therefore the studied point x is non-hyperbolic, with signature [− 0 +]. Notice that if z has an imaginary component, its conjugate isz = −1/z, which cannot be, as zz ≥ 0, therefore one zero at −1 implies that the other two zeros are necessarily real. It follows that P ψ (−1) = 0 indicates a neutral saddle case, even if Σ cannot be evaluated.
We now relate the characteristic polynomials of ∇ψ T and ∇f T .
Lemma 3.4. The characteristic polynomial of the Jacobian matrix ∇ψ T and the characteristic polynomial of the mesochronic Jacobian ∇f T are given by the expres- The coefficients are linked by the expressions t ψ = 3 + T tf where T is the length of the averaging period inf T . Moreover, the incompressibility condition d ψ ≡ 1 imposes the relation on the coefficients of Pf for T = 0.
Proof. The connection between the two Jacobian matrices ∇ψ T and ∇f T is given by relation (9): ∇ψ T = Id + T ∇f T . The characteristic polynomial P ψ of ∇ψ T can be re-written using the characteristic polynomial Pf of ∇f T Using the notation for coefficients of Pf from (14), we can expand the last line to obtain . By comparing coefficients with the general expression for P ψ , the statement of the theorem follows. The incompressibility condition from the statement of the theorem is a consequence of substituting d ψ ≡ 1 for any T > 0.
We now combine the preceding Lemmas to give the proof of Theorem 3.2.
Since tf , mf , and df are related through the incompressibility constraint tf +mf T + df T 2 ≡ 0 derived in (16), we can formulate the condition in two alternative ways: The other mesohyperbolicity condition t ψ − m ψ = 0 translates into either 8df T 3 = 0 or 8T (tf + mf T ) = 0.
Reformulating the expressions for Σ and ∆ in Lemma 3.3 using df , mf and tf through (15) we obtain Since tf , mf , and df are related through the incompressibility constraint (16), we can re-formulate the expressions using either tf or mf To emphasize the connection to the OWC expressions (Theorem 2.5), we will use the df and mf versions of the formulas. The statement of the proof is equivalent to Lemma 3.3 where the criteria are expressed in terms of the spectral coefficients of ∇f T instead of ∇ψ T .
In summary, to identify the mesochronic class of a point x, take the following steps: 1. compute the Jacobian ∇f T (x) (details in Appendices B and C), 2. evaluate ∆ and Σ using (10), 3. use Table 3 to identify the mesochronic class. . As a consequence, the averaged fieldf T tends to f pointwise as T → 0 + and ∇f → ∇f . By continuity of eigenvalues with respect to matrix elements, df → d f , mf → m f , and tf → t f as T → 0 + . As all three of these quantities are finite, the mesochronic incompressibility criterion (16) is trivially satisfied in the limit. Additionally, due to incompressibility of the vector field, it holds that tf  Proof. According to Theorem 2.5, we obtain that d f = 0. Since the maps T → ∇f T , T → df T , and T → mf T are continuous, the instantaneous hyperbolicity will imply mesohyperbolicity for some small T min > 0. The continuity means that there exists T min > 0 on which continuity of the maps T → df (T ) and T → 3t 3 df T +2T 2 mf T −8 implies df T = 0, 3T 3 df T + 2T 2 mf T − 8 = 0, for T ∈ [0, T min ]. By virtue of Theorem 3.2, the point (t 0 , x) is also mesohyperbolic with respect to [0, T ] and the proof is complete.
As a consequence, the signs of the indicators ∆ and Σ (10) have the following limits sgn Σ . By comparing mesochronic classification criteria (Table 3) in the limit T → 0 with the instantaneous OWC criterion (Theorem 2.5), we can conclude that mesochronic classification reduces to OWC classification in the T → 0 limit. Put differently, mesochronic classes generalize OWC classes to time intervals [t 0 , t 0 + T ] with T > 0.

Asymptotic limits.
For autonomous dynamical systems defined over compact domains, asymptotic rates of deformation are defined by Lyapunov exponents. The finite-time analogs, termed Finite-Time Lyapunov Exponents (FTLE) [27,59] are defined using the polar decomposition of the Jacobian matrix of the time-T map ψ T . Let ∇ψ T = R |∇ψ T | be the polar decomposition of ∇ψ T where |∇ψ T | 2 = ∇ψ T ∇ψ T . Eigenvalues of |∇ψ T | are non-negative, singular values of ∇ψ T , so we can define Finite-Time Lyapunov Exponents σ i to be the exponential growth/decay rates of the singular values, i.e., we represent the singular values of ∇ψ T as e σiT .
It is a well-known fact in matrix analysis that when a matrix is normal, i.e., unitarily diagonalizable, its singular values are equal to absolute values of its eigenvalues. In the language of this paper, this means that when ∇ψ T is normal, positions of the eigenvalues µ i of ∇ψ T in reference to the unit circle are determined by the signs of the Finite-Time Lyapunov Exponents σ i . It follows that Σ > 0 ( Table 3) implies that both two eigenvalues of ∇ψ T are outside of the unit circle and that two Finite-Time Lyapunov Exponents are positive, when ∇ψ T is normal.
Unfortunately, ∇ψ T is not generally normal for any finite T , i.e., its eigenvectors are not orthonormal. However, Ref. [22] shows, although non-rigorously, that for smooth ergodic systems which have distinct Lyapunov exponents both real parts of eigenvalues and singular values of ∇ψ T can be written as e σiT as T → ∞. As the sign of Σ indicates the number of eigenvalues of ∇ψ T outside the unit circle, which, in turn, is determined by real parts of logarithms of those eigenvalues, we conclude that, when T → ∞, the sign of Σ indicates whether two or one Lyapunov exponents are positive, assuming that the conjecture in Ref. [22] holds.

4.3.
Recovering the 2D mesochronic deformation criterion. The supplement to Ref. [45] presented a derivation of the criteria for mesohyperbolicity for 2D (planar) differential equations. Planar differential equations can be trivially embedded into a 3D state space by adding a third state with trivial (zero) dynamics. We use such an embedding to demonstrate how the 3D mesochronic criteria (Theorem 3.2) specialize to the 2D criterion.
Since df ≡ 0 as noted above, it follows that the flow is not 3D-mesohyperbolic, as it is expected, as the third coordinate is always preserved due to the construction of the 3D flow.
The indicators Σ and ∆ evaluate to . Therefore, the sign of the ∆ indicator is determined by the 2D mesohyperbolicity criterion. If ∆ > 0, D 2D < 0 and the two eigenvalues are non-real and lie on the unit circle due to incompressibility constraints. If ∆ < 0, D 2D > 0 and the eigenvalues are real, one is contracting and the other expanding.
For planar dynamics, D 2D = 0 implies that the dynamics are either a pure shear, when the geometric multiplicity of the eigenvalue is 1, or trivial, i.e., ψ T (x) ≡ x.
Incompressibility in the planar case is more restrictive than in 3D, yielding only two structurally stable cases: mesohyperbolic λ i ∈ R and mesoelliptic λ i ∈ C, which intersect at the pure reflection/shear case λ 1 = λ 2 = 1. The derivation of the mesohyperbolicity criterion for the 2D case therefore relied only on detection of real vs. complex eigenvalues, which is the reason why D 2D in (18) is taken as the sole 2D mesohyperbolicity criterion, without resorting to a more complicated calculation. Table 4 we compute explicitly the values of the indicators Σ and ∆ for a simple class of linear time-invariant systems whose ∇ψ T is constant, and given by the polar decomposition:

Linear time-invariant velocity fields. In
In this parametrization, we can independently manipulate rates of strain λ 1,2,3 as well as the rate of rotation ω present in the system. As two of the rates have the same sign, unless one of the rates is zero, we choose to order the directions by setting sgn λ 1 = sgn λ 2 , which means that the third direction is of the opposite sign λ 3 = −λ 1 − λ 2 , due to incompressibility. While these systems do not represent a broad range of dynamical systems, we have a good understanding of their dynamics so it is instructive to see how their properties are reflected in the mesochronic classification.
First, when all rates λ i are non-zero, all points are mesohyperbolic as Σ is constant and non-zero; presence or absence of rotation determines whether a point is a (mesohyperbolic) saddle (ω = 0) or a helix (ω = 0). The signature [− − +] or [− + +] of the saddle is then determined by the sign of the pair λ 1,2 . When the two rates match exactly, λ 1 = λ 2 , it implies ∆ = 0. Since the quantities ∆ and Σ are functions of the spectrum, they alone are not enough to detect whether associated directions align (shear) or not (saddle). In the case of the systems derived, we know that those two directions correspond to independent eigenvectors, which means that the point is a saddle.
If one of the rates is equal to 0, it always implies that Σ = 0, which is classified as one of the neutral 3D mesochronic classes. In that case, ∆ corresponds to the 2D mesohyperbolicity indicator D 2D , as described in Section 4.3. Again, the presence ω = 0 or the absence of rotation ω = 0 is reflected on the sign of ∆, where ∆ > 0 corresponds to the former, and ∆ < 0 to the latter case. The magnitudes of ∆ and Σ grow exponentially in most of the cases; however, notice that in the presence of rotation, a periodic function multiplies the exponentiallygrowing magnitude, resulting in ∆(T ) = 0 periodically. This means that there is a potential for resonance, i.e., if T is a multiple of the period of oscillation, dynamics momentarily appears to be on the boundary behavior between [−++] and [−−+] saddle mesohyperbolicity, or even a pure reflection when ∆ = Σ = 0. This choice of T is, of course, highly unlikely without a prior knowledge of ω.
In summary, analysis of simple linear systems shows that mesochronic classes correctly reflect our intuition about presence of stretching and rotation in linear, time-invariant flows.

5.2.
Arnold-Beltrami-Childress Flow. The Arnold-Beltrami-Childress (ABC) flow [13] is a kinematic model of an incompressible fluid flow evolving in a threedimensional periodic domain. Even though the system of ODEs specifying the ABC flow is simple, it exhibits a variety of different behaviors and has been used as a test-bed for various computational algorithms [7,8,19,27].
The ABC flow evolves on a 3-torus in periodized state variables (x, y, z) ∈ [0, 2π] 3 ∼ = T 3 . Dynamics depend on parameters A, B, C, D ∈ R and are specified by differential equationsẋ = A(t) sin z + C cos ẏ where the time-varying parameter A(t) is given by If D = 0, the equations are autonomous; if, additionally, any other parameter is 0, the system is integrable. [13] The linearization along a solution p(t) = (x(t), y(t), z(t)) of (19) is given bẏ The determinant and the sum of minors are given by the expressions These expressions can be used to evaluate the OWC criterion for the instantaneous hyperbolicity according to Theorem 2.5.
Thus, for the initial condition p(0) = (x 0 , y 0 , z 0 ), z(t) = z 0 + (C sin y 0 + B cos x 0 )t for all t ∈ [0, T ]. Furthermore, if we write σ := x + y, δ := x − y, then the ABC flow can be rewritten as a decoupled second order system which is a direct product of two pendulum-like equations. Solutions of the first two components can be written in terms of integrals of Jacobi elliptic functions, and it follows that the system is integrable. A similar argument follows in case when either A, B, or C are zero, in addition to D = 0. Lemma 5.1. All points (x, y, z) in the state space of the system (19) with A(t) ≡ 0 are non-mesohyperbolic over any interval [t 0 , t 0 + T ].
Proof. When A(t) ≡ 0, the matrix defining the linear system of equations (20) is block-diagonal, with blocks and 0 on the diagonal. The fundamental matrix of the system is, therefore, also block diagonal, with value 1 on the diagonal corresponding to the exponential of the block 0 in (20). Since 1 is then in the spectrum of the time-T map Jacobian, (x, y, z) is non-mesohyperbolic for all T .
Remark 2. If A(t) ≡ 0, the mesochronic class of a point (x, y, z) does not depend on the value of z since the Jacobian matrix in (20) does not depend on the z-coordinate of the solution around which we linearized.

5.2.2.
Steady non-integrable case. The structure of the invariant sets in the state space of the ABC flow for parameters A = √ 3, B = √ 2, C = 1, D = 0 is well studied analytically [13] and numerically [8]. The state space contains six interwoven vortices with the space between them filled by chaotic dynamics (Figure 2). We place a grid of 400 × 400 initial conditions on the (x, y) face of the periodicity cell and calculate tf , mf , df for time intervals of different lengths. Other details about numerics are given in Appendix C.  To give a sense of time scales involved in the system, Figure 3 shows several trajectories (pathlines) within a single vortex, simulated for various durations T . Trajectories inside vortices take approximately T = 3 to cross one periodicity cell. The two other panels in Figure 3 show that the vortex rotates around its axis while the inner layers move at slightly faster speeds than its outer layers.
To detect non-mesohyperbolic behavior we need to numerically evaluate condi- written here using the traces of the flow map and the relation m ψ = t ψ −1 to estimate their growth more clearly. These conditions are difficult to reliably compute in the face of numerical errors that will almost-certainly result in non-zero quantities. Numerically, we determine when the criteria are satisfied using a numerical tolerance determined by estimating the growth rate of t ψ and t ψ −1 with increasing   length T of time intervals. We infer the behavior of these quantities in linear timeinvariant systems in order to establish the correct baseline -the criterion has to reflect correctly our knowledge of linear systems if we are to even consider using it for classification of nonlinear systems. In what follows, we derive an empirical estimate of mesohyperbolicity, using a quantity termed numerical hyperbolicity. Our considerations will further rely on T → ∞ arguments as our intent is to estimate the behavior of the flow beyond the time interval in which we are sampling it.
In both cases, as T → ∞, To account for exponential growth, we set the numerical tolerance of mesohyperbolicity based on logarithms of expressions (21) (In all expressions we omit dependence on state variables and time interval for shortness). Non-zero values of either h 1 or h 2 are signs of mesohyperbolic behavior; conversely, we need either one of them to be small to declare non-mesohyperbolicity. In nonlinear flows, we do not expect that h 1,2 will be entirely independent of the value of T . Nevertheless, in ergodic regions, [66] we expect convergence in mean as T → ∞.
Even the rate of convergence to the mean is not uniform: in regular ergodic regions, e.g., vortices of the ABC flow, the expected decay is O(T −1 ); in strongly mixing regions, conjectured to be embedded within the chaotic region, the expected decay is O(T −1/2 ), i.e., similar to the Central Limit Theorem for i.i.d. random variables. Initial conditions that are neither regular nor strongly mixing may potentially have an even slower decay of variance [54], T −α for any 0 < α < 1/2. Values of h at those points would then still grow as h ∼ T 1/2−α . The volume of weakly mixing zones is small in systems containing Kolmogorov-Arnold-Moser-type dynamics, [8,53,62] and therefore we do not expect those values to occur as major features in the histogram of h.
A good quantitative criterion for deciding whether a point is mesohyperbolic should estimate whether the smaller of the two h 1,2 is "sufficiently" close to zero. Under a conjecture that some variant of the Central Limit Theorem (CLT) holds for the estimate of max i |λ i | we can test whether the deviation of our estimator min{h 1 , h 2 } from the hypothesis of non-mesohyperbolicity max i |λ i | = 0 is normally distributed, i.e., whether numerical mesohyperbolicity h, defined by h = |min{h 1 , h 2 }| √ T , where, as before is small, h < ε, for some (small) constant ε. When this is indeed the case, then we are reasonably confident that with longer T the estimated eigenvalues of the flow map would indeed have a unit modulus and we empirically declare that the point is not mesohyperbolic, as classified by Theorem 3.2.
Is there a fixed value of ε that could help us decide whether the hypothesis of mesohyperbolicity for a point holds? To answer this question standard CLT results from probability theory, e.g., Lindeberg-Feller [14,Thm. 3.4.5], would employ higher order moments of the distribution of samples of the random processes h 1,2 . Such a constant ε would then turn our empirical criterion into a statistical test, where h would be the z-score for testing the hypothesis of whether a point under consideration is mesohyperbolic or not. Unfortunately, we cannot assume to know how the higher moments behave, despite existing work on CLTs in the context of dynamical systems [23,24,56,65,66].
In lieu of rigorous results, we choose to proceed empirically and set the cutoff value ε based on distributions of the numerical mesohyperbolicity h. Figure 4a shows histograms of numerical mesohyperbolicity h for a range of values of T , conforming well to expectations. As T increases, the distribution of h changes from a fairly flat distribution (T = 1) to a bimodal distribution with well-separated  peaks. Figure 5 shows that each mode of distribution of h corresponds, respectively, to vortices and to the large chaotic region between them as T → ∞. Based on these results, we declare numerical non-mesohyperbolicity using the cutoff parameter value h < ε, with ε = 10 −0.3 . Therefore, a solution ϕ(·, 0, p) with p = (x, y, 0) is instantaneously hyperbolic everywhere except along the lines As T ≈ 0, the mesochronic partition is virtually identical to the Okubo-Weiss-Chong partition.
The mesochronic partition for T ≈ 0 is illustrated in Figure 6, and due to short integration period T , matches exactly the Okubo-Weiss-Chong partition. Notice that the conventional intuition about vortices being "elliptic" structures cannot be inferred from short integration times, as for T ≈ 0 almost the entire space is mesohyperbolic (non-mesohyperbolic lines (23) are difficult to sample numerically). Nevertheless, neutral mesohelical regions roughly coincide with locations of vortices, while the chaotic region between vortices contains a mixture of all four classes of mesohyperbolicity. Comparing with Figure 2, we see that the boundaries of mesochronic classes do not align with boundaries of invariant sets.
Increasing T results in the sequence of images shown in Figure 7, where the mesohyperbolic regions are shown in the left column, and the non-mesohyperbolic regions in the right column. As the averaging period is increased to T = 1, partitions deform, but remain largely uncorrelated with invariant features. As we increase T beyond 1, non-hyperbolic behavior significantly re-appears along the interface between different mesohyperbolic classes. Parts of boundaries of mesohyperbolic zones start to align with invariant vortices. Since level sets of any function averaged for a sufficiently long time will partition the state space into invariant sets [8,47], this is expected. Notice that the non-mesohyperbolic regions appear almost exclusively inside invariant vortices. As T is increased even further, the non-mesohyperbolic zones grow inside the invariant vortices and eventually completely occupy them. In the chaotic zone, we see disappearance of mesohelical dynamics, with only saddle mesohyperbolicity remaining, which matches the asymptotic analysis in Section 4.2.

Unsteady case.
We now keep the parameters A = √ 3, B = √ 2, C = 1 as before, but set D = 1, which results in the unsteady variation in the coefficient A(t) = √ 3 + t sin t. The unsteady ABC flow has not received as much analytic attention as the steady case; however, it was used as a demonstration of numerical techniques for computation of the flow map Jacobian ∇ψ T in [7]. Figure 8 shows the same sets of initial conditions used in demonstrating the steady flow ( now advected by the unsteady flow for comparison. Notice that the initial advection patterns are similar, until A(t) starts significantly deviating from its constant term A. Figure 4b shows the histogram of numerical mesohyperbolicity while Figure 9 shows the spatial distributions of (non-)mesohyperbolicity classes, determined using the same numerical criterion as in the steady case (22). For short intervals T = 0.5, 1.0, mesochronic classification of the flow is similar to the steady case. This is expected as the magnitude of A(t) is dominated by the steady component.   split between two behaviors that was observed in the steady flow ( Figure 7) is not present here. Remnants of the axial vortex in the left and two "eyes" of vortices in the right sides of images are visible both in mesochronic and FTLE partitions. Figure 10 illustrates transport of initial conditions sampled from several regions in the state space of the unsteady ABC flow. The first two rows show advection from patches chosen as subsets of regions that are mesohyperbolic for integration times T = 0.5, 1, 5. The last row shows advection from a patch that straddles several mesochronic regions at T = 5.
Advection up to t = 5 demonstrates that the initial conditions selected from single mesochronic regions (central column, top two frames) do not disperse much, i.e., stay coherent. Initial conditions from a mixture of regions show considerably more dispersion.
Advection up to t = 7, shown in the third column, demonstrates how the patches of initial conditions evolve beyond the interval T = 5 that was used to generate mesochronic classes. All material patches at this point show considerable growth; arguably, the patch in the last row again shows the largest dispersion. We conclude this section with Figure 11 showing images of orbits, i.e., pathlines, of material points advected by the unsteady ABC flow through space over 5 time units. Two mesohelical patches used to initialize points are the same as in Figure 10; the two additional mesosellar patches were initialized similarly. The two mesosellar sets are sampled from a narrower region, so it is not surprising that they remain more tightly packed than mesohelical sets. We believe that the major distinction between the top and the bottom row is in the internal order of orbits within each bundle. Mesohelical sets seem to preserve the internal order, i.e., despite the torsion of the bulk, one can still observe an ordered rainbow of colors toward the end of the orbits. On the other hand, both mesosellar sets show more internal mixing of colors. While this is potentially a circumstantial observation here, it may be an interesting point to explore in the future, due to the association of hyperbolic saddles with mixing and bulk rotation with lack of mixing.
6. Discussion. Mesochronic analysis is a local analysis in that it classifies a single point based on the properties of the local deformation gradients along trajectory emanating from it. Nevertheless, the hope is that sets of initial conditions selected using the local mesochronic criteria will collectively stay coherent on a macroscopic level and, potentially, deform as a bulk in the way suggested by the local analysis. The contributions of this paper are in extension of two-dimensional theory of [45] to three-dimensional flows.
We have extended the concepts of mesohyperbolicity (local strain over finite times) and mesoellipticity (local rotation over finite times) to three dimensions, which allows for co-existence of the rotation and the strain. In turn, we differentiate between the mesosellar behavior, involving three distinct directions of strain, and the mesohelical behavior, involving a plane in which material simultaneously rotates  . Orbits (pathlines) of 100 points initialized from two mesohelical patches shown in Figure 10 and two additional mesosellar patches, advected for time t ∈ [0, 5]. Color is added to illustrate internal ordering of the material. and, possibly, uniformly strains. Quantities Σ and ∆ that we defined simplify classification of domain points, side-stepping explicit calculation of eigenvalues of the flow. The eigenvectors of ∇ψ T also have a role in the foregoing analysis. A complex conjugate pair determines a plane whose normal is the vector of finite-time rotation. A real eigenvector indicates the direction of finite-time stretching under the map.
Rotation of the material has a much richer presence in 3D than in 2D classification. In two dimensions, any non-uniformity in strain, accompanied by rotation or not, manifests itself as a hyperbolic deformation. In other words, only rigid-body rotations are highlighted as non-hyperbolic in two dimensions, in addition to shear. In three dimensions, rotation in a plane can be accompanied by strain in the normal direction: the type of that strain distinguishes between three mesohelical classes: [− + +], [− − +] and neutral.
Consequently, the invariant vortices in the steady ABC flow initially show a significant amount of hyperbolicity in them, indicating that layers within them move at different speeds in the axial direction ( Figure 3); however, over the longer timescales they appear neutrally mesohelical, which corresponds to interpretation of them as rigid rotating structures. This suggests that the separation between different layers of vortices is asymptotically sub-exponential.
In the unsteady ABC flow discussed in Section 5.2.3, material is significantly less coherent as the unsteadyness destroys long-term invariant structures. Nevertheless, for interval lengths T = 1, 5, while the unsteady term A(t) = A + Dt sin(t) is of the same order of magnitude as the steady A, the distributions of mesohyperbolic areas show loose correspondence between steady and unsteady ABC flows. However, as T → 10, the regions that turned neutrally mesohelical in the steady case are instead replaced by a mixture of hyperbolic mesosellar regions, indicative of chaotic mixing.
We have demonstrated that visualization of mesochronic classes corresponds with well-known behavior of the fluid-like ABC flow. In particular, it is interesting to see how well the mesohelical regions correspond with the vortex zones in which Kolmogorov-Arnold-Moser-type structures are known to exist. At this point, the theory is immediately applicable to kinematic analysis of the geometric structure in chaotic advection in fluids; however, numerical algorithms used in this paper serve the proof-of-concept purpose, and more accurate and efficient methods for computation of the flow Jacobian [7] can and should be readily used, if applicable.
Note that the mesochronic classification is invariant to Galilean transformations but not to rotating frame transformations [29], it therefore discovers transport properties for three-dimensional incompressible flows which are observed in the given frame of reference the system is in. It is interesting to consider the relationship between the notion of exponential dichotomy [10] and mesohyperbolicity discussed here. Mesohyperbolicity (or -ellipticity, -helicity) are notions that are valid on finite time intervals. They are best used as tools to study the bifurcations of local dynamics of a trajectory in time, when changing time intervals are selected, i.e. mesohyperbolicity is a function of two parameters, the beginning and final time. It appears plausible that, adapting the techniques used in [52], it can be shown that if a trajectory is not mesohyperbolic for any interval [t 1 , t 2 ] inside the interval of interest [T 1 , T 2 ], then there is no (appropriately defined [52]) finite-time exponential dichotomy on that interval. Thus, the concepts defined here have the potential to parametrize finite-time stability properties.
The true test of any method for detection of geometric structures is its usefulness to the applied communities, e.g., physical oceanography, and flow engineers. We therefore plan to apply the technique to the unsteady testbed flows in [4], to transitory systems [48] and to more physically-relevant models in order to further verify practical use of the mesochronic classification, beyond confirmations obtained for the 2D case. Furthermore, it remains to be understood if the highlighted quantities, e.g., the trace, the cofactor trace, and determinant of the mesochronic Jacobian, are also useful in dynamic analysis of turbulent fluid flows.