Computation and optimal perturbation of finitetime coherent sets for aperiodic flows without trajectory integration
Abstract
Understanding the macroscopic behavior of dynamical systems is an important tool to unravel transport mechanisms in complex flows. A decomposition of the state space into coherent sets is a popular way to reveal this essential macroscopic evolution. To compute coherent sets from an aperiodic timedependent dynamical system we consider the relevant transfer operators and their infinitesimal generators on an augmented spacetime manifold. This spacetime approach avoids trajectory integration, and creates a convenient linearization of the aperiodic evolution. This linearization can be further exploited to create a simple spectral optimization methodology for diminishing or enhancing coherence. We obtain explicit solutions for these optimization problems using Lagrange multipliers and illustrate this technique by increasing and decreasing mixing of spatial regions through small velocity field perturbations.
1 Introduction
Analyzing complicated flows through their transport and mixing behavior has been and still is attracting a great amount of attention [MMP84, RKLW90, Wig92, HP98, Are02, JW02, Wig05, SLM05, FP09, Thi12, FPG14, KK16, KR18, HKK18], both from geometric and probabilistic points of view. Nonautonomous timeaperiodic dynamics poses additional difficulties, especially the case of finite time, where asymptotic notions cannot be applied.
The current work has two main contributions.

A technique to find a small perturbation of the underlying aperiodic vector field in a prescribed ball in a space or subspace of vector fields, which optimally enhances or destroys the existing finitetime coherent sets. This extends optimization results in [FS17], which considered the timeperiodic setup of [FK17], to aperiodic dynamics.
The key construction in [FK17] is the representation of a periodically forced flow on phase space as an autonomous flow on timeaugmented phase space , where denotes the unit circle. On this timeexpanded phase space, the time coordinate is simply advanced at a constant rate:
Finitetime coherent sets on the time interval , as introduced in [FSM10, Fro13], were extracted from singular functions of the transfer operator , where the transfer operator is the linear operator describing the evolution of distributions under the dynamics subject to a small random perturbation. The crucial observation is that singular modes of are eigenmodes of , where denotes the dual of , and that for areapreserving dynamics (corresponding to divergencefree velocity fields ) the dual is the transfer operator of the timereversed dynamics, i.e., the (again, slightly stochastically perturbed) flow governed by the timereflected velocity field . The dynamic interpretation of this operatorbased characterization is that (finitetime) coherent are those sets that are to a large extent mapped back to themselves by a noisy forwardbackward evolution of the dynamics. This operatorbased framework not only gives a qualitative framework for coherence; the singular values of provide quantitative bounds for coherence [Fro13, FPG14]—namely the closer the singular value is to one, the less mixing occurs between the coherent set and its exterior under the noisy dynamics.
By concatenating the “forwardtime” and “backwardtime” velocity fields on time intervals and —see (9) below—we will construct a system on the augmented space that mimics the forwardbackward evolution of the dynamics. In this way, the eigenmodes of the (Fokker–Planck) generator of this augmented system yield the eigenmodes of ; i.e., the singular modes of , and from these singular modes the desired finitetime coherent sets. Again, the corresponding eigenvalues of the generator can be used to give quantitative bounds on coherence/mixing, see Theorem 7. We note that a numerical approach to extract coherent sets from by solving the Fokker–Planck equation has been described in [DJM16]. In contrast to [DJM16] we do not require timeintegration over , which is especially advantageous once we consider the optimization of coherence and mixing.
Connecting the spectral properties of the generator of the augmentedspace system with the (finitetime) dynamical properties of the original system is a generalization of the results from [FK17], where it has been done for timeperiodic velocity fields on infinite time intervals. The interplay between the spectra of the dynamics in augmented space and the nonautonomous dynamics in the original space has strong connections to the correspondence between evolution semigroups [How74] and twoparameter semigroups, as elaborated in, e.g., [CL99]; see also [EN00, Section VI.9] for a general introduction. We mention that by a similar construction, spatiotemporal dynamical patterns were extracted in [GD17] by considering the generator of the Koopman operator (the adjoint of the transfer operators considered here) associated with the augmentedspace dynamics.
There are several different ways to measure mixing and mixedness under (stochastically perturbed) dynamics, such as considering dispersion statistics or the change of variation in a concentration field; see e.g. [Pro99, LH04, TDG04, Thi08]. Multiscale norms of mixing measure how “oscillatory” a concentration field is [MMP05]; see [Thi12] for a review. The most widely used approaches to the problem of mixing optimization search for switching protocols between some fixed velocity fields in order to optimize some topological [BAS00] or other mixing measure [MMG07, CAG08, OBPG15]. Other strategies include the optimal distribution of concentration sources [TP08] and geometric dynamical systems techniques [Bal15]. An interesting theoretical result is that arbitrary mixedness under advectiondiffusion can be achieved in finite time solely by sufficiently increasing the strength of the (otherwise fixed) advective flow [CKRZ08]. If there are no restrictions to the choice of the velocity field, one can choose the one that is optimally mixing the actual concentration at every time instance [LTD11]. We also note that a related problem to mixing enhancement arises in statistical mechanics [LNP13] where the convergence toward the stationary distribution should be accelerated, e.g., to increase the efficiency of sampling.
Instead of focusing on one fixed concentration field, we will bound the mixing characteristics of a flow in terms of the objects that most inhibit mixing: coherent sets. As we mentioned above, finitetime coherent sets are characterized by the singular vectors of the transfer operator, and, equivalently, by the eigenvectors of the generator of the augmentedspace process, while the corresponding eigenvalue delivers an upper bound on transport between the coherent set and its exterior. Thus, we can quantitatively access the mixing behavior of a flow on finite time through the spectrum of the augmented generator , and can target these eigenvalues if we want to enhance or diminish mixing.
Given a default velocity field and nonautonomous perturbations from an admissible space of divergencefree velocity fields, our approach considers “small” that change a targeted eigenvalue of (thus also the singular value of ) locally optimally. This procedure can be iterated to obtain a larger perturbation in a gradientmethod fashion. Optimizing singular values of the transfer operator directly is difficult as it would necessarily involve the variation of the nonlinear dynamics under the velocity field . Instead, a linearized optimization of the eigenvalue of the generator leads to a very simple optimization problem (42), which can be solved via a linear system of the same dimension as . Moreover, the theory holds for infinitedimensional perturbation spaces as well, since Fréchet differentiability of the transfer operator and its spectrum with respect to perturbing velocity fields has been established in [KLP18]. It should be noted that additive perturbations in the sense might be nonphysical or technically not realizable, however that is the topic of flow control [BMT01, KB07], and certainly beyond the scope of the current work.
This work is structured as follows. In section 2 we introduce the function based formalism to study advectiondiffusion systems by the Fokker–Planck equation and its evolution operator, the transfer operator, in forward and backward time. In section 3 we consider purely advective transport between a family of sets and its exterior in terms of fluxes through the boundary of the family; this is a geometric analogue of the operatorbased considerations that follow. Section 4 introduces the new reflected augmented generator needed to handle aperiodic, finitetime driving. We state formal connections between the spectrum of the reflected augmented generator and the reflected transfer operator, and provide spectralbased bounds on the maximal possible coherence of sets in phase space under the aperiodic dynamics. Section 5 contains a numerical demonstration of the efficacy of our trajectoryfree approach. Section 6 describes the formal setup of the optimization problem designed to manipulate the position of the dominant spectral values of the reflected augmented generator, culminating in an explicit expression for the optimal timedependent local perturbation of the velocity field. Section 7 specializes the infinitedimensional results of section 6 to the numerical setting via discretization and includes a variety of examples of coherence reduction and enhancement. We conclude in section 8.
2 Advectivediffusive dynamics
Let be compact with smooth boundary (at least and piecewise ). We consider the time interval and the dynamics
(1) 
with , for some , and being a standard Wiener process in . The initial point is distributed according to some initial density . The evolution of the density of the governing equation (1) is given by the Fokker–Planck equation or Kolmogorov forward equation [LM13, Section 11.6].
(2)  
(3)  
where is the normal derivative on the boundary. Associated to (2) is an evolution operator that transports a density at time 0 to the solution density of (2) at time . The evolution operator is an integral operator with stochastic^{1}^{1}1Doubly stochastic if the flow is volumepreserving. kernel that satisfies [Fro13, Assumptions 1 and 2].
2.1 Construction of a forwardbackward process
For simplicity of presentation we assume that the velocity field is divergence free for all . We note that the remaining arguments in this section may be carried through for general velocity fields. Denote by the canonical scalar product of a Hilbert space . Following [Fro13] in the volumepreserving setting^{2}^{2}2The operators and are replaced by normalised versions for nonzero divergence velocity fields; these are denoted by and in [Fro13]., coherent sets over the time interval are extracted from the eigenfunctions of corresponding to large eigenvalues, where is the adjoint of , defined to be the unique linear operator satisfying
for all . The eigenvalues of (the singular values of ) are known to lie in the interval (cf. [Fro13, p.3]). The rationale behind the operator is that describes evolution in forward time, describes evolution under the timereversed dynamics, and coherent sets are characterized exactly by the property that they are “stable” under a noisy forwardbackward evolution of the dynamics.
The adjoint operator is the solution operator to the Kolmogorov backward equation [PS08]:
(4)  
The operator maps a density at time backward in time according to (4) to produce a density at time . We may simplify (4) using volume preservation:
(5)  
thus we may write (4) as
(6) 
Reversing time in (4) to obtain an initial value problem we get
(7)  
using , and the velocity field . Comparing (2) and (7) we see that the natural evolution of the adjoint problem (the Kolmogorov backward equation) corresponds to the forward problem (the Kolmogorov forward equation) of the time reversed dynamics.
We wish to construct the process over the time interval that corresponds to the operator . We view as evolution on the time interval and we therefore shift (7) by time units, defining to obtain a forward problem on :
(8)  
We denote the solution operator of this problem as ().
Finally, we concatenate the two forward problems (2) and (8) to make a single process over . We mark objects that live on this extended interval with a hat . Define the velocity field
(9) 
using the reflection map
(10) 
The resulting velocity field exhibits discontinuities in and whenever it does not vanish there, but onesided derivatives exist. In what follows, we will solve the Fokker–Planck equation
(11)  
over the interval ; more precisely on with continuous concatenation at . Let us summarize the above construction with the following proposition.
3 Cumulative flux from a reflected family of sets
Before proceeding with the operatorbased description of finitetime coherence, in this section we analyze the reflected dynamics by its flux through the boundary of a moving (possibly coherent) set. Our intention behind connecting this to a flux in augmented space (i.e., spacetime) in Proposition 3 is partially to set the stage for the augmentedspace operatorbased description in section 4.
We consider a family of dimensional sets , satisfying the following assumptions (the boundaries are piecewise smooth in space and differentiable in time):
Assumption 2.

There exists a codimension 1 parameterisation set such that for each there is a bijective function with being piecewise smooth, ( in and in ) piecewise.

is well defined for all and all , where .
The cumulative outflow flux under the vector field from the family of sets is given by
(12) 
where denotes the positive part, is the dimensional surface measure and is the outer normal unit vector.
In augmented state space we consider the augmented set
(13) 
and define the instantaneous outflow flux from the set by
(14) 
with denoting the dimensional surface measure and denoting the augmented velocity field defined by
We extend this result to the reflected velocity field generated by a general aperiodic velocity field. We define a reflected family of sets in synchrony with the reflected vector field :
(15) 
In particular, recalling from (10), for every the boundary has a parametrization and
(16) 
At the right and the leftsided partial derivatives of with respect to exists but they may not be equal. The family of normal vectors is mirrored in time (see Figure 1): for .
Proposition 3.
The cumulative outflow flux from the family of sets under the vector field , is equal to the cumulative absolute flux in and out of the family of sets , , under the vector field ; that is,
(17) 
Furthermore, the timeintegrated flux (17) is equal to the instantaneous absolute flux in augmented space:
(18) 
Proof.
Let us first prove the first equality. Therefore we do not yet need objects of the augmented setting. We split the integral
(19) 
and see that the first integral looks promising. So we only need to treat the second integral (19). We use Fubini’s theorem and substitute using with being the Gram determinant. (Transformation theorem would not work, we need the sign we get from substitution.)
Combining the two calculations above we get the desired result. The second equality involving the objects of the augmented setting follows analogously to [FK17, Theorem 2]. ∎
4 Coherent families of sets and the generator on augmented phase space
In this section we create a spectral mapping theorem for our reflected augmented process (Proposition 4) and derive a bound for the finitetime coherence of a family of sets in terms of the second eigenvalue of a generator (the infinitesimal operator) of our augmented reflected advectiondiffusion process (Theorem 7). To do this we build on discretetime theory from [Fro13] with the periodic continuoustime theory from [FK17].
4.1 The evolution operator for the reflected process
It is well known that is a compact, integral preserving, real and positive operator^{3}^{3}3See [FK17, Lemma 30, at the end of the proof] for compactness, and [LM13] or [KKS16] for the other properties.. Furthermore is a selfadjoint operator with simple largest eigenvalue . Following [Fro13] one has the second eigenvalue satisfies
(20) 
where in the volumepreserving^{4}^{4}4In the nonzero divergence case, is a reference measure describing the initial mass distribution of the (possibly compressible) fluid being evolved, and is the forward evolution of under a normalised version of denoted by in [Fro13]. setting and are both simply the Lebesgue measure. We now consider the problem (11) introduced in section 2 as a timeperiodic problem on (we extend periodically). Following the considerations of section 2 the evolution operator from time to time is given by
(21) 
The situation when is exactly is of particular importance:
(22) 
Note that is selfadjoint for , .
4.2 The timeaugmented generator
We now turn to the augmented reflected system
(23) 
on , and note that is a standard Wiener process on ; in particular it is not constructed by time reflection. We define augmented versions of , , and , denoting them with bold symbols:
and is a dimensional standard Wiener process. The augmented system is
Considering the timeperiodic problem with on we formulate a Fokker–Planck equation in augmented space:
(24) 
Note there is no diffusion in the direction, as per the definition of . Let denote the augmented generator associated with the augmented Fokker–Planck equation (24) in augmented space:
(25)  
with continuous concatenation conditions for the other parts of the boundary, . We may also write (24) and (25) in terms of the nonautonomous (“unaugmented”) dynamics:
(26) 
where is the righthand side operator of (11) at time , i.e., the generator of the reflected Fokker–Planck equation. Analogously to [FK17, Lemma 22] the following result holds.
Proposition 4.
Let be an eigenfunction of corresponding to the eigenvalue . One has then
(27) 
for all and .
Proof.
We will prove (27) using some results and ideas from [FK17]. First we apply [FK17, Lemma 30] piecewise on and concerning wellposedness and regularity. Therefore we consider the original problem (2) and the reflected, shifted, timereversed problem (11). Now [FK17, Lemma 30] guaranteees the unique existence of a function with
and the regularity
Now we can proceed as in [FK17, Lemma 22]. Let and with . According to the construction above (26) we know
and this implies for all
(28) 
Now is the evolution operator to the evolution equation
Therefore the function solves (28) uniquely. This gives the claim. ∎
Let be an eigenvalue of with an eigenfunction . Inserting and into Proposition 4 yields
(29) 
Recalling that is a compact selfadjoint positive operator, it must be that for some . This implies
(30) 
from which it follows that for some .
4.3 Coherent families of sets
In the specific case where the velocity field is periodic in time, [FK17] shows that the families of sets
have an escape rate (see [FK17, Definition 8]) of at most , where is the first nontrivial eigenvalue of corresponding to the eigenfunction . Because we consider the dynamics on a finite time interval, this notion of escape rate is replaced by the concept of a coherence ratio [Fro13].
In the general setting of aperiodic we will quantify the coherence of families and provide a construction of highly coherent families with associated rigorous coherence bound. To ensure the existence of our coherence quantifier, we place an assumption on the family of sets^{5}^{5}5Such a family of sets was called “sufficiently nice” in [FK17]. Here we simply call it “nice”..
Definition 5 (Nice family of sets).
A nice family of sets is such that

for every the set is closed; and

the function —where denotes the Hausdorff distance—is either left or rightcontinuous at for every .
Definition 6 (Coherence ratio).
Let be a family of measurable sets. Denote by the law of the process generated by the SDE (1) initialised with , where denotes a possibly signed Borel measure. For we define the coherence ratio of the family as
(31) 
It was shown in [FK17, Appendix A.6] that for a nice family of sets the quantity (31) is well defined. If no initial distribution is specified, denotes the coherence ratio with respect to the uniform distribution (normalised Lebesgue measure on ).
The following theorem shows that the probability of a trajectory remaining in a family of sets constructed from the positive and negative parts of eigenfunctions of decays no faster than the rate given by the corresponding eigenvalues.
Theorem 7.
Let with , and assume that the family of sets with is nice. If is scaled such that , then
(32) 
Proof.
By Proposition 4 we have that is an eigenfunction of the integral preserving operator at the eigenvalue . Thus we have , and with it again by the integralpreserving property. With the given scaling of we have that , and let . Note that by Proposition 4 we have that , thus [FK17, Equation (A.10)] gives
Noting that , the claim follows. ∎
5 Computational aspects
5.1 Numerical discretisation
We use the “Ulam’s discretisation for the generator” approach developed in [FJK13] for autonomous flows and extended in [FK17, Sections 7.2 and 7.3] for nonautonomous flows. In brief, referring to the above papers for further details, the Ulam discretisation for the generator yields a matrix that may be interpreted as a rate matrix of a finitestate, continuoustime Markov chain, with the states corresponding to a partition of into hypercubes (hyperrectangles) in . The entries which correspond to rates between boxes adjacent in the temporal coordinate direction (in which the time evolution is a rigid rotation of constant velocity 1) are given by , where is discretised into intervals of length . The entries in the remaining space directions are computed from the rate of flux out of the hypercube faces by numerical integration of the component of the velocity field normal to (and pointing out of) the face; see the expression for in [FK17, Section 7.2]. The entries of the rate matrix corresponding to diffusive dynamics (in the space coordinates only, there is no diffusion in the time coordinate) are computed from a finitedifference approximation of the Laplace operator; see the expression for in [FK17, Section 7.2]. We then set . The matrix can also be interpreted as a rate matrix for a finitestate Markov chain; it has an eigenvalue and its spectrum is confined to the left half of the complex plane.
The reflected velocity field from (9) may be substituted for the velocity field used in [FK17] and the methodology of [FK17] employed; this is the approach taken in the numerical experiments below.
Remark 8.
The computation of uses only the outwardpointing velocity field values on the faces of the partition elements—similarly to how the outward flux is defined through the positive part of the inner product in (14)—, discarding the inwardpointing parts. Because of the reflected structure of , a slightly more efficient implementation would be to store the evaluations of the velocity field normal to hypercube faces in both directions (not only in the outwardpointing direction). The outwardpointing components would be used on the time interval , while the inwardpointing components would be used on the time interval , where they are outwardpointing because of the sign flip in (9)—similarly as it happens in the proof of Proposition 3. This would reduce by half the computational effort in evaluating the velocity field components normal to the hypercube faces. As we report below, the assembly of the generator matrices is relatively fast anyway, and we have not tried to optimise our implementation of Ulam’s method for the generator in this reflected setting.
5.2 Example: Bickley jet
We now apply the reflected augmented generator approach to a perturbed Bickley Jet [RBBV07]. The following model describes an idealized zonal jet in a band around a fixed latitude, assuming incompressibility, on which two traveling Rossby waves are superimposed. The velocity field is induced by the stream function
The constants are chosen according to [RBBV07]. The length unit is () and the time unit is days. For the amplitudes and the speed of the Rossby waves we choose
with being the Earth’s radius. Further we choose
The state space is periodic in direction and is given by . For good numerical tractability, we resolve our reflected spacetime manifold with a spatially somewhat coarse grid, that is uniform in space () leads to square boxes needed for isotropic diffusion) and sufficiently finely resolved in time (). We choose .
The system described above is equipped with outflow boundary conditions instead of homogeneous Neumann conditions on . This leads to a slightly different spectral structure of the generator, which now generates a semigroup of substochastic operators. Its leading eigenvalue is strictly less than zero.
This Bickley Jet differs from the one investigated in [FK17, Section 7.6], as it is not periodic in time. Further, we expect the eigenvectors computed in [FK17, Section 7.6] to induce sharper and smoother separations, because they utilize the periodicity and describe asymptotic properties. The singular vectors computed here give rougher results, because they only have finitetime information available and only describe finitetime dynamics.
Analogously to [FK17, Section 7] our timeaugmentation produces companion eigenvalues. Companion eigenmodes denote eigenmodes that are “higher order harmonics” of existing eigenmodes differing only in temporal modulation, and encoding the same coherence information; see below. For more details on the companion eigenvalues for the Ulamdiscretisation we refer to [FK17, Section 7.3]. We will use and verify the relations derived there. Therefore we calculate the eigenvalues and vectors of with the smallest magnitude instead of largest real part using eigs(G,10,’SM’).
Table 1 shows a gap after the first and sixth eigenvalue. Let us first discuss the leading eigenvectors. Figure 2 shows the eigenvectors corresponding to the dominant (i.e., smallest real part) eigenvalues. The first eigenvector, the quasistationary (or conditionally invariant) distribution highlights (in blue, see Figure 2(a)) parts of the domain that get pushed out of the region of consideration. The red regions are those parts of phase space that remain longest in under the diffusive dynamics (1). This effect is due to the outflow conditions. Note that this example is in this sense explorative, as our theory in the previous sections was only considering reflecting and not outflow boundary conditions.
The second eigenvector indicates an upper/lower separation. The other four eigenvectors show combinations of coherent vortices.
To investigate which elements of the spectrum of contain genuinely new dynamical information, we checked that the complex eigenvalues of from to are all companion eigenvalues equal to (recall from section 5.1 that is the temporal grid spacing)
(33) 
for an eigenvalue of