{'fraction_non_alphanumeric': 0.05675224061977822, 'fraction_numerical': 0.013550053167248974, 'mean_word_length': 4.039191671769749, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Coherent sets for nonautonomous dynamical systems

We describe a mathematical formalism and numerical algorithms for identifying and tracking slowly mixing objects in nonautonomous dynamical systems. In the autonomous setting, such objects are variously known as almost-invariant sets, metastable sets, persistent patterns, or strange eigenmodes, and have proved to be important in a variety of applications. In this current work, we explain how to extend existing autonomous approaches to the nonautonomous setting. We call the new time-dependent slowly mixing objects coherent sets as they represent regions of phase space that disperse very slowly and remain coherent. The new methods are illustrated via detailed examples in both discrete and continuous time.


Introduction
The study of transport and mixing in dynamical systems has received considerable attention in the last two decades; see e.g. [31,46,1,47] for discussions of transport phenomena. In particular, the detection of very slowly mixing objects, known variously as almost-invariant sets, metastable sets, persistent patterns, or strange eigenmodes, has found wide application in fields such as fluid dynamics [34,30,35], ocean dynamics [21, ?], astrodynamics [4], and molecular dynamics [9,41]. A shortcoming of this prior work, based around eigenfunctions of Perron-Frobenius operators (or transfer operators, or evolution operators) is the restriction to autonomous systems or periodically forced systems. In this work, we extend the notions of almost-invariant sets, metastable sets, persistent patterns, and strange eigenmodes to timedependent Lagrangian coherent sets. These coherent sets form a time parameterised family of sets that approximately follow the flow and disperse very slowly; in other words they stay coherent. Coherent sets are the natural nonautonomous analogue to almost-invariant sets.
The standard dynamical systems model of transport assumes that the motion of passive particles are completely determined by either an autonomous or a time-dependent vector field. Traditional approaches to understanding transport are based upon the determination of the location of geometric objects such as invariant manifolds. In the autonomous setting, an invariant manifold of one dimension less than the ambient space will form an impenetrable transport barrier that locally partitions the ambient space. In the periodically-forced setting, primarily in two-dimensional flows, it has been shown that slow mixing in the neighbourhood of invariant manifolds is sometimes controlled by "lobe dynamics" [36,37,46]. In the truly non-autonomous, or aperiodically forced setting, finite-time hyperbolic material lines [24] and surfaces [25] have been proposed as generalisations of invariant manifolds that form barriers to mixing. These material lines and surfaces are known as Lagrangian coherent structures; see also [42] for an alternative definition. The geometric approach can often be used to find co-dimension 1 sets (coherent structures) that form boundaries of coherent sets.
An alternative to the geometric approach is the ergodic theoretic approach, which attempts to locate almost-invariant sets (or metastable sets) directly, rather than inferring their location indirectly from their boundaries. The basic tool is the Perron-Frobenius operator (or transfer operator). Real eigenvalues of this operator close to 1 correspond to eigenmodes that decay at slow (exponential) rates. Almost-invariant sets are heuristically determined from the corresponding eigenfunctions f as sets of the form {f > c} or {f < c} for thresholds c ∈ R. Such an approach arose in the context of smooth autonomous maps and flows on subsets of R d [7,6] about a decade ago. Further theoretical and computational extensions have since been constructed [17,15,16]. A parallel series of work specific to time-symmetric Markov processes and applied to identifying molecular conformations was developed in [40,8,9] and surveyed in [41].
There have been some recent studies of the connections between slow mixing in periodically driven fluid flow and eigenfunctions of Perron-Frobenius operators. Liu and Haller [30] observe via simulation a transient "strange eigenmode" as predicted by classical Floquet theory. Pikovsky and Popovych [34,35] numerically integrated an advection-diffusion equation to simulate the evolution of a passive scalar, observing that it is the sub-dominant eigenfunction of the Perron-Frobenius operator that describes the most persistent deviation from the unique steady state.
The Perron-Frobenius operator based approach has been successful in a variety of application areas, however, as the key mathematical object is an eigenfunction, there is no simple extension of the method to systems that have nonperiodic time dependence 1 Indeed, Liu and Haller [30] state that: "...strange eigenmodes may also be viewed as eigenfunctions of an appropriate Frobenius-Perron operator...This fresh approach offers an alternative view on scalar mixing, but leaves the questions of completeness and general timedependence open." It is this question of general time-dependence that we address in the current work. We extend a standard formalism for random dynamical systems to the level of Perron-Frobenius operators to create a Perron-Frobenius operator framework for general time-dependence. We also state an accompanying numerical algorithm, and demonstrate its effectiveness in identifying strange eigenmodes and coherent sets.
An outline of the paper is as follows. In Section 2 we formalise the notions of nonautonomous systems in both discrete and continuous time. In Section 3 we describe a Galerkin projection method that we will use to produce finite matrix representations of Perron-Frobenius operators. In Section 4 we define the critical constructions for the nonautonomous setting. We show that the nonautonomous analogues of strange eigenmodes are described by the "Oseledets subspaces" or "Lyapunov vectors" corresponding to compositions of the projected Perron-Frobenius operators. In Section 5 we describe in detail a numerical algorithm to practically compute these slowly decaying modes, and demonstrate that in the continuous time setting, these modes vary continuously in time. Our numerical approach is illustrated firstly in the discrete time setting with an aperiodic composition of interval maps, and secondly in the continuous time setting with an aperiodically forced flow on a cylinder. Section 6 provides some further background on almost-invariant sets and coherent sets and Section 7 describes a new heuristic to extract coherent sets from slowly decaying modes in the nonautonomous setting. This heuristic is then illustrated using the examples from Section 5.

Nonautonomous Dynamical Systems
We will treat time dependent dynamical systems on a smooth compact d-dimensional manifold M ⊂ R D , D ≥ d in both discrete and continuous time. In order to keep track of "time" we use a probability space (Ω, H, P), with the passing of time controlled by an ergodic automorphism θ : Ω preserving P (ie. P = P•θ −t for all t ≥ 0). We require this somewhat more complicated description of time for technical reasons: to run the ergodic-theoretic arguments in Theorem 1 that guarantee the existence of the nonautonomous analogues of strange eigenmodes. The requirement that P be an ergodic probability measure rules out obvious choices for Ω and θ: (i) in discrete time, Ω = Z and θ s (t) = t + s, and (ii) in continuous time, Ω = R and θ s (t) = t + s. In both (i) and (ii), there is no ergodic probability measure on Ω preserved by θ. In the next two sections, we will introduce suitable examples of Ω and θ and describe the nonautonomous systems they generate.

Discrete time -Maps
In the discrete time setting, we will think of Ω ⊂ (Z) Z , and θ as a left shift σ on Ω defined by (σω) i = ω i+1 , where ω = (. . . , ω −1 , ω 0 , ω 1 , . . .) ∈ Ω. We assume that σ is ergodic with respect to P. Let T = {T ω 0 } ω 0 ∈Z be a collection of (possibly non-invertible) piecewise differentiable maps on a compact manifold M. For brevity, we will sometimes write T ω in place of T ω 0 . We will define a nonautonomous dynamical system by map compositions of the form For k ≥ 0 (resp. k < 0), Φ(k, ω, x) represents the forward time (resp. backward time) k-fold application of the nonautonomous dynamics to the point x initialised at "time" ω. Whenever T ω is non-invertible, T −1 ω (x) will represent the finite set of all preimages of x. We call Φ a map cocycle. Definition 1. Endow M with the Borel σ-algebra and let µ be a probability measure on This definition of an invariant measure is stricter than is usual for random or nonautonomous dynamical systems (e.g. [2,Definition 1.4.1]). More generally, one may allow sample measures µ = µ ω and insist that µ σ −1 ω • Φ(−1, ω, ·) = µ ω for all ω ∈ Ω.
Example 1 (Aperiodic map cocycle). We construct a map cocycle Φ by the composition of maps T i from a collection T according to sequences of indices ω ∈ Ω. The collection T := {T 1 , T 2 , T 3 , T 4 } consists of expanding maps of the circle S 1 , which we think of as [0, 1] with endpoints identified. The sequence space Ω ⊂ {1, 2, 3, 4} Z is given by Elements of Ω correspond to bi-infinite paths in the graph Figure 1. The shift σ : Ω → Ω is a subshift of finite type. A Borel σ-algebra H is generated by the length-one cylinder sets C i = {ω : ω 0 = i}, i = 1, . . . , 4, and by giving equal measure to these four cylinder sets, we generate a shift-invariant probability measure P.
Let m denote normalised Lebesgue measure on M. To each map T ω we associate a Perron-Frobenius operator P ω : L 1 (M, m) defined by P ω f = y∈T −1 ω x f (y)/| det DT ω (y)|. The operator P ω is a linear operator that acts on integrable functions in analogy to the action of T ω on points. If f ∈ L 1 (M, m) represents a density function for an ensemble of initial conditions, then P ω f represents the density function of the ensemble after the action of T ω has been applied to the ensemble. The map cocycle Φ naturally generates a Perron-Frobenius cocycle P (k) ω = P σ k−1 ω • · · · • P σω • P ω . This composition of k Perron-Frobenius operators capture the action on a function f after k iterations of the non-autonomous system.

Continuous time -Flows
Let F : Ω × M → R d be a sufficiently regular vector field. More precisely, we suppose that F satisfies the conditions of [2, Theorem 2.2.2], which will guarantee the existence of a classical solution of the nonautonomous ODEẋ(t) = F (θ t ω, x(t)), t ∈ R.

Definition 2.
Endow M with the Borel σ-algebra and let µ be a probability measure on M. We call µ an invariant measure if µ • φ(−t, z, ·) = µ for all z ∈ Ξ and t ∈ R. Remark 1. In Definition 2 we are insisting that µ is preserved at all "time instants". As in the discrete time setting, more generally one may allow µ = µ z and insist that µ ξ(−t,z) • φ(−t, z, ·) = µ z . However, as we will soon begin to focus on coherent sets rather than invariant measures, we will restrict the invariant measure to a "time independent" measure for clarity of presentation. This is perfectly reasonable for one of the main applications we have in mind, namely, aperiodically driven fluid flow where µ ≡ Lebesgue, and volume is preserved by the flow at all times.
We may define a family of Perron-Frobenius operators as P

Galerkin projection and matrix cocycles
Following Ulam [45], in the sequel we will consider the finite rank operators π n P (1) ω : L 1 (M, m) → B n and π n P (1) z : L 1 (M, m) → B n , and the matrix representations of the restrictions of π n P (1) ω and π n P (1) z to B n . We denote these matrix representations (under multiplication on the right) by P (ω) and P (z). Extending Lemma 2.3 [29] in a straightforward way to the nonautonomous setting, one has and In particular, these matrices are numerically accessible.
Remark 2. Note we do not concern ourselves at all with the relationship between P (1) ω and π n P (1) ω ; this is a subtle relationship and beyond the scope of this work. See [29, 13, 10, 32, 14, 3, ?] for work in this direction.
The matrices P (ω) and P (z) generate matrix cocycles and

Discretised Oseledets functions and the Multiplicative Ergodic Theorem
In periodically driven flows, Liu and Haller [30] and Pikovsky and Popovych [34], observed that certain tracer patterns persisted for long times before eventually relaxing to the equilibrium tracer distribution. Pikovsky and Popovych [34] recognised these patterns as graphs of eigenfunctions of a Perron-Frobenius operator corresponding to an eigenvalue L < 1. These eigenfunctions decay over time and the closer L is to 1, the slower the decay and the more slowly an initial tracer distribution will relax to equilibrium. We now develop a framework for the considerably more difficult aperiodic setting. Consider some suitable Banach space (F , · ) of real valued functions; F is the function class in which we search for slowly decaying functions. Suppose that the norm is chosen so that for each ω ∈ Ω and k ≥ 0, the operator P (k) ω is Markov; that is, P (k) ω = 1 for all ω and k ≥ 0. For f ∈ F , we calculate the following limit: We refer to λ(ω, f ) ≤ 0 as the Lyapunov exponent of f . If f decays under the action of the Perron-Frobenius operators at a geometric rate of r k , 0 < r < 1, then λ(ω, f ) = log r. The closer r is to 1, the slower the decay. The extreme case of r = 1 (no decay) is exhibited when f is the density of the invariant measure µ that is common to all maps in our nonautonomous dynamical system. We define the Lyapunov spectrum Λ(P, ω) In the aperiodic setting the new mathematical objects that are analogous to strange eigenmodes and persistent patterns will be called Oseledets functions.

Definition 3.
Oseledets functions correspond to f for which (i) λ(ω, f ) is near zero and (ii) the value λ(ω, f ) is an isolated point in the Lyapunov spectrum.
By considering (F , · ) = (B n , · 1 ), the actions of P (k) ω and P (k) z are described by P (k) (ω) and P (k) (z), respectively. We may replace P (k) ω and P (k) z in (12) by P (k) (ω) and P (k) (z), respectively, to obtain a standard setting where the possible values of λ(ω, f ) are the Lyapunov exponents of cocycles of n × n matrices, and and Λ(P, z) := lim exist for P almost-all ω ∈ Ω, and consist of at most n isolated points, λ n < · · · < λ 1 = 0. Of particular interest to us is the function f 2 (ω) (or f 2 (z)) in B n , which represents the function that decays at the slowest possible geometric rate λ 2 .
Remark 3. In certain settings, this matrix cocycle exactly captures all large isolated Lyapunov exponents of the operator cocycle P : (BV, · BV ) . One such setting is a map cocycle formed by composition of piecewise linear expanding maps with a common Markov partition B = {B 1 , . . . , B n }; see [19].
The following example illustrates the concept of Lyapunov spectrum and Oseledets functions in the familiar autonomous setting. For the remainder of this section, we adopt the discrete time notation of σ and ω.
Example 3 ("Autonomous" single map). In [5] individual maps are constructed for which the Perron-Frobenius operator has at least one non-unit isolated eigenvalue when acting on the Banach space (BV, · BV ). A single autonomous map may be regarded as a cocycle over a one-point space Ω = {ω}, and so we may drop the dependence on ω in notation. Keller [27] shows that for a piecewise expanding map T of the interval I, the spectrum of the associated Perron-Frobenius operator P has an essential spectral radius ρ ess (P) equal to the asymptotic local expansion rate sup x∈I lim k→∞ 1/DT k (x) 1/k , and that there are at most countably many spectral points, each isolated, of modulus greater than ρ ess (P). In order to have an isolated spectral point, we construct a map of S 1 which has an almost-invariant set (see Definition 5). The relation between almost-invariant sets and isolated eigenvalues was noted in [7]. Consider the partition Given a ∈ Z 6 , any map T : S 1 → S 1 defined by is Markov with respect to B. Here we take a = (0, 0, 1, 4, 3, 3); see Figure 5. Notice that there is a low transfer of mass between the two intervals [0, 1/2] and [1/2, 1]. Since B is a Markov partition for T , the space of characteristic functions is an invariant subspace of BV for the Perron-Frobenius operator P of T . Thus the action of P ω = P on F = B 6 is represented by the matrix which has non-zero eigenvalues L 1 = 1, The map T is piecewise affine with constant slope 3 and so the logarithm of the local expansion rate is log(1/3). The eigenvalue L 2 ≈ 0.805 of P thus gives rise to an isolated point λ 2 ≈ log 0.805 in the Lyapunov spectrum Λ(P). The corresponding Oseledets function f 2 is given by where w 2 is the eigenvector of P corresponding to the eigenvalue L 2 , see Figure 5. Since |L 3 | ≈ 0.138 < 1/3, this means that log L 2 is the unique isolated Lyapunov exponent in Λ(P). Note that the set {f 2 > 0} corresponds to the set [0, 1/2]. We will discuss this property further in Section 7.
In the nonautonomous setting, we can no longer easily construct Oseledets functions as eigenfunctions of a single operator, or eigenvectors of a single matrix. In fact, the Oseledets functions are themselves (aperiodically) time dependent in the nonautonomous setting. Our model of Oseledets functions for nonautonomous systems is, as the name suggests, built around the Multiplicative Ergodic Theorem, see e.g. [2,Chapter 3,§4]. We now state a strengthened version [19] of the Multiplicative Ergodic Theorem that we require for our current purposes.

Numerical approximation of Oseledets functions
In the autonomous and periodic settings we have seen in Examples 3 and 4 that the subspaces W 2 (ω) = sp{w 2 (ω)} were one-dimensional, and that the vectors w 2 (ω) could be simply determined as eigenvectors of matrices. For truly nonautonomous systems (those that are aperiodically driven), the Oseledets splittings are difficult to compute. In this section we outline a numerical algorithm to approximate the W j (ω) subspaces from Theorem 1. The algorithm is based on the push-forward limit argument developed in the proof of Theorem 1.
To streamline notation, we describe the discrete time and continuous time setting separately.

Discrete time
We first describe a simple and efficient method to construct the matrix P (ω) defined in (8).
1. Partition the state space M into a collection of connected sets {B 1 , . . . , B n } of small diameter.
2. Fix i, j, and k and create a set of Q test points x j,1 , . . . , x j,Q ∈ B j that are uniformly distributed over B j .

Set
We now describe how to use the matrices P (ω) to approximate the subspaces W j (ω). An intuitive description of the ideas behind Algorithm 2 immediately follows the algorithm statement.
1. Construct the Ulam matrices P (M ) (σ −N ω) and P (N ) (σ −N ω) from (19) and (10) for suitable M and N. The number M represents the number of iterates over which one measures the decay, while the number N represents how many iterates the resulting "initial vectors" are pushed forward to better approximate elements of the W j (ω).

Form
as an approximation to the standard limiting matrix (ω) is our numerical approximation to W j (ω).
Here is the idea behind the above algorithm. If we choose M large enough, the eigenspace U (M ) j (σ −N ω) should be close to the limiting (M → ∞) eigenspace U j (σ −N ω). Vectors in the eigenspace U j (σ −N ω) experience stretching at a rate close to L j . Note that the eigenspace U (M ) j (σ −N ω) is the j th singular vector of the matrix P (M ) (σ −N ω), which experiences a "per unit time" average stretching from time −N to −N + M of L j . Choose some arbitrary v ∈ U j (σ −N ω) and write v = ℓ j ′ =j w j ′ with w j ′ ∈ W j ′ (σ −N ω). Pushing forward by P (N ) (σ −N ω) for large enough N will result in P (N ) (σ −N ω)w j dominating P (N ) (σ −N ω)w j ′ for j < j ′ ≤ ℓ. Thus, for large M and N we expect W (M,N ) j (ω) to be close to W j (ω). Remarks 1.

Theorem 1 states that W
2. This method may also be used to calculate the Oseledets subspaces for two-sided linear cocycles, and may be more convenient, especially for large n, than the standard method of intersecting the relevant subspaces of flags of the forward and backward cocycles.
The numerical approximation of the Oseledets subspaces has been considered by a variety of authors in the context of (usually invertible) nonlinear differentiable dynamical systems, where the linear cocycle is generated by Jacobian matrices concatenated along trajectories of the nonlinear system. Froyland et al. [18] approximate the Oseledets subspaces in invertible two-dimensional systems by multiplying a randomly chosen vector by P (N ) (σ −N ω) (pushing forward) or P (−N ) (σ N ω) (pulling back, where P (−N ) (σ N ω) = P −1 (ω) · · · P −1 (σ N −1 ω). Trevisan and Pancotti [?] calculate eigenvectors of Ψ (M ) (ω) for the three-dimensional Lorenz flow, increasing M until numerical convergence of the eigenvectors is observed. Ershov and Potapov [12] use an approach similar to ours, combining eigenvectors of a Ψ (M ) with pushing forward under P (N ) . Ginelli et al. [23] embed the approach of [18] in a QR-decomposition methodology to estimate the Oseledets vectors in higher dimensions. In the numerical experiments that follow, we have found our approach to work very well, with fast convergence in terms of both M and N.

Continuous time
As our practical computations are necessarily over finite time intervals, from now on, when dealing with continuous time systems, we will compute P (k) (z) as π n P (k) z rather than as P (ξ(k − 1, z)) · · · P (ξ(1, z)) · P (z). If the computation of π n P (k) z can be done accurately (this will be discussed further in Section 5.5), then this representation should be closer to P (k) z as there are fewer applications of π n .
We first describe a simple and efficient method to construct the matrix P (ω) defined in (8). 2. Fix i,j, and z and create a set of Q test points x j,1 , . . . , x j,Q ∈ B j that are uniformly distributed over B j .

Set
The flow time M should be chosen long enough so that most test points leave their partition set of origin, otherwise at the resolution given by the partition {B 1 , . . . , B n }, the matrix P (M ) (ξ(−N, z)) matrix will be too close to the n × n identity matrix. If the action of φ separates nearby points, as is the case for chaotic systems, clearly the longer the flow duration M, the greater Q should be in order to maintain a good representation of the images φ(M, ξ(−N, z), B i ) by the test points.
1. Construct the Ulam matrices P (M ) (ξ(−N, z)) and P (N ) (ξ(−N, z)) from (20) for suitable M and N. The number M represents the flow duration over which rate of decay is measured, while the number N represents the duration over which the resulting "initial vectors" are pushed forward to better approximate elements of the W j (z).   (ξ(δ, z)), we represent each as a function and make a comparison in the L 1 norm. We assume that U

Continuity of the Oseledets subspaces in continuous time
(z). We firstly comparef n,ξ(δ,z),M,N +δ andf n,z,M,N .  Proof. The proof runs as a non-autonomous version of the discussion in Remark 7.6.2 [28]. ξ(δ, z), ·) denotes the flow from ξ(δ, z) in reverse time for duration δ. For the moment consider continuous f . Since x → φ(s, z, x) is at least C 1 for each s, z (the derivative of φ wrt to x is continuous with respect to s and x for each fixed z) by [2, Theorem 2.2.2 (iv)] and M is compact, P Proof. This result is more difficult to demonstrate as we need to firstly compare Ψ (M ) (ξ(−N, z)) with Ψ (M ) (ξ (−N + δ, z)). To this end, consider The right hand side converges to zero as δ → 0 by Lemma 1. This result implies that (ξ(δ, z)) are close for sufficiently small δ. Thus f n,ξ(−N,z),M and f n,ξ(−N +δ,z),M are close in L 1 norm. Now we need to push both of these forward by π n P(ξ(−N, z)) (N ) . This will not increase the norm of the difference at all, so f n,ξ(δ,z),M,N − f n,z,M,N 1 will also be small.
As in Examples 3 and 4, the maps T i are piecewise affine with constant slope three, and so ρ(ω) = 1/3. Thus log L 2 and log L 3 may approximate isolated Lyapunov exponents in Λ(P).
We follow Algorithm 2 to approximate the second Oseledets subspace W  (ω) to an Oseledets subspace. In Section 7.2 we will see how to extract coherent sets from these functions.

Oseledets functions in a 2D continuous time nonautonomous system
We consider the following nonautonomous system on M = [0, 2π] × [0, π], t ∈ R + : This equation describes a travelling wave in a stationary frame of reference with rigid boundaries at y = 0 and y = π, where the normal flow vanishes [33,39]. The streamfunction (Hamiltonian) of this system is given by We set c = 0.5, A = 1, and the phase speed to ν = 0.25. The velocity field is 2π-periodic in the x-direction, which allow us to study the flow on a cylinder. The velocity fields in a comoving frame for these parameters are shown in Figure 10. The closed recirculation regions adjacent to the walls (y = 0 and y = π) move in the positive x-direction and are separated from the jet flowing regime by the heteroclinic loops of fixed points, which are given below. This model can be simplified to an autonomous system with a steady streamfunction in the comoving frame by setting X = x − νt and Y = y. The steady streamfunction is then given by S(X, Y, t) = −(c − ν)Y + A sin(X) sin(Y ). Let X s = sin −1 ((c − ν)/A) and Y s = cos −1 ((c − ν)/A). In the comoving frame, the recirculation region at the wall Y = 0 contains an elliptic point q 1 = (π/2, Y s ) and is bounded by the heteroclinic loop of the hyperbolic fixed points p 1 = (X s , 0) and p 2 = (π − X s , 0). Similarly, those elliptic and hyperbolic points at the wall Y = π are q 2 = (3π/2, π − Y s ), p 3 = (π + X s , π), and p 4 = (2π − X s , π), respectively, see Figure 10. One may observe that there is a continuous family of invariant sets in the comoving frame as any fixed level set of the streamfunction bounds an invariant set. In a stationary frame these elliptic and hyperbolic points (and their heteroclinic loops) are just translated in the x−direction. That is, any fixed level set of the time-dependent streamfunction (23) is a (time-dependent) invariant manifold. We note, however, that the recirculation regions are distinguished from the remainder of the cylinder as they are separated from the jet flowing region, which has a different dynamical fate. In the subsequent sections we will perturb this somewhat "degenerate" system to destroy the continuum of invariant sets in the comoving frame and produce a small number of almostinvariant sets (see Definition 5) in the comoving frame, or coherent sets in the stationary frame).

A coherent family: Mixing case
We modify the traveling wave model in the previous section to allow mixing in the jet flowing region. We add a perturbation to the system in the following way: x = c − A(z(t)) sin(x − νz(t)) cos(y) + εG(g(x, y,z(t))) sin(z(t)/2) y = A(z(t))) cos(x − νz(t))) sin(y).
The following definition generalises invariant sets to almost-invariant sets. In the continuous time case we define: Definition 5. Let µ be preserved by the autonomous flow φ. We will say that a set A ⊂ M is ρ 0 -almost-invariant over the interval [0, τ ] if for all s ∈ [0, τ ], 2. A is connected.
If A ⊂ M is almost-invariant over the interval [0, τ ], then for each s ∈ [0, τ ], the probability (according to µ) of a trajectory leaving A at some time in [0, s], and not returning to A at time s is relatively small. In the discrete time setting, τ = 1, and the obvious changes are made in Definition 5. By convention we ask that A is connected; if A is not connected, we consider each connected component to be an almost-invariant set for suitable ρ 0 .
We now begin to discuss the nonautonomous setting. The notion of an invariant set is extended to an invariant family. Definition 6.

Continuous time:
We will call a family of sets {A ξ(t,z) }, Motivated by a model of fluid flow, we imagine coherent sets as a family of connected sets with the property that the set A ω is approximately mapped onto A σ k ω by k iterations of the cocycle from "time" ω; that is, Φ(k, ω, A ω ) ≈ A σ k ω . The definition of coherent sets combines the properties of almost-invariant sets and an invariant family. As we now have a family of sets we require one more property beyond those of Definition 5, in addition to modifying the almost-invariance property. In the continuous time case we define: Definition 7. Let µ be preserved by a flow φ and 0 ≤ ρ 0 ≤ 1. Fix a z ∈ Ξ. We will say that a family {A ξ(t,z) } t≥0 , A ξ(t,z) ⊂ M, t ≥ 0 is a family of ρ 0 -coherent sets over the interval [0, τ ] if: for all s ∈ [0, τ ] and t ≥ 0, In discrete time, we replace (26) with τ necessarily becomes 1, and we make the obvious changes to the other items in Definition 7. We remark that by selecting some A ⊂ M of positive µ measure and defining A ξ(t,z) := φ(t, z, A), t ≥ 0, the family {A ξ(t,z) } t≥0 , is a family of 1-coherent sets. Such a family is not of much dynamical interest, as there is nothing distinguishing this family from one constructed with another connected subset A ′ ⊂ M. We are not interested in these constructions of coherent sets, and in practice the numerical algorithm we present in the next section is unlikely to find such sets for chaotic systems.
We may extend Algorithm 5 in order to find a sequence of coherent sets {A σ i ω } K i=0 . Since we require the measure of a sequence of coherent sets to be constant, we seek to maximize the mean value of ρ µ over a given time range as we vary the measure of the sets.
Algorithm 6 (To determine a sequence of maximally coherent sets over a range of times ω, . . . , σ K ω).
It is interesting to compare the locations of the coherent sets with their corresponding maps in the mapping cocycle, see Figure 14. As for the previously constructed family Figure 14: Graphs of T σ k ω showing A σ k ω × A σ k+1 ω for k = 0, . . . , 5.
R s(ω,k) ([0.5, 1]) k∈N , the coherent sets alternate between two positions separated by a rotation of approximately 0.25. However, the mean value of ρ µ is greater for the sequence A σ k ω constructed from Algorithm 6 since in each case the coherent set matches up well with local maxima and minima of the preceding map.
In Figure 16 we note that the grey set A + z on the left at time z = (0, 1, 1.5) flows approximately to the light grey set A + ξ(10,z) on the right at time ξ (10, z). Similarly for the black sets A − z and A − ξ (10,z) . This carrying of the time z coherent sets to the time ξ(10, z) coherent sets by the aperiodic flow is only approximate, as ρ µ (Â + z (c * ),Â + ξ(10,z) (η(c * ))) = 0.9605 and ρ µ (Â − z (c * ),Â − ξ(10,z) (η(c * ))) = 0.9599. Thus, we expect a loss of about 10% under the advection of the flow. Figure 16 also zooms onto A + z and A − z to demonstrate this loss of mass. To make this loss even more apparent, we continue to flow forward for 50 time units. Figure 17 shows that the (black) coherent sets A + z and A − z do indeed disperse over time, however at a much slower rate than the arbitrarily chosen (grey) sets. The coherent sets A + z , A − z are just single elements of a time parameterised family {A + ξ(t,z) , A − ξ(t,z) } t≥0 of coherent sets that at any given initial time describe those sets that will disperse most slowly over a duration of 10 time units.

Final Remarks
We have formulated a new mathematical and algorithmic approach for identifying and tracking coherent sets in nonautonomous systems. Our new approach generalises existing suc- (grey) and A − ξ(10,z) (black) are identified by thresholding the Oseledets functions. Overlays of φ(10, z, A + z ) (grey) and φ(10, z, A − z ) (black) are also shown. [Bottom left] Zooms of A + (z) and A − (z) . [Bottom right] Overlays of φ(10, z, A ± z ) (grey/black dots) on A ± ξ(10,z) (white), displaying the loss of mass over 10 time units duration from z = (0, 1, 1.5). Figure 17: Trajectories of the perturbed system (24) for ε = 1. The large light grey (A + ξ(t,z) ) and black (A − ξ(t,z) ) blobs are the coherent sets identified by our approach. The other (medium grey) blobs are chosen nearby the coherent sets to show strong mixing away from the coherent regions.
cessful transfer operator methodologies that have been used in the autonomous setting. Our constructions address the question raised by [30] of how to study strange eigenmodes and persistent patterns observed in forced fluid flows in the general time-dependent situation. Future work will include applying these techniques to detect and track mobile coherent regions in oceanic and atmospheric flows, extending significantly the flow times studied in [21, ?] and [?].