Computation and optimal perturbation of finite-time coherent sets for aperiodic flows without trajectory integration

TL;DR

This paper introduces a novel method to compute and optimize finite-time coherent sets in aperiodic flows without trajectory integration, using transfer operators and spectral optimization to control mixing.

Contribution

It develops a space-time generator approach for coherent set computation and provides explicit solutions for perturbation-based optimization of flow coherence.

Findings

Avoids trajectory integration by using space-time generators

Provides explicit solutions for flow perturbation optimization

Demonstrates control of mixing through small velocity field changes

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 time-dependent dynamical system we consider the relevant transfer operators and their infinitesimal generators on an augmented space-time manifold. This space-time generator approach avoids trajectory integration, and creates a convenient linearization of the aperiodic evolution. This linearization can be further exploited to create a simple and effective 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.

Figures

Tables

Table 1. Table 1 : Eigenvalues ( μ k subscript 𝜇 𝑘 \mu_{k} ) of 𝑮 ^ bold-^ 𝑮 \bm{\hat{G}} ordered in ascending magnitude and corresponding approximate singular values ( σ k subscript 𝜎 𝑘 \sigma_{k} ) of 𝒫 0 , τ subscript 𝒫 0 𝜏 \mathcal{P}_{0,\tau} according to ( 31 ).

${\mu }_1$	$0$	${\mu }_4$	$-0.35061$	${\sigma }_1$	$1$	${\sigma }_4$	$0.24599$
${\mu }_2$	$-0.09033$	${\mu }_5$	$-0.44766$	${\sigma }_2$	$0.69674$	${\sigma }_5$	$0.16685$
${\mu }_3$	$-0.34938$	${\mu }_6$	$-0.45702$	${\sigma }_3$	$0.24720$	${\sigma }_6$	$0.16072$

Table 2. Table 2 : Eigenvalues ( μ k subscript 𝜇 𝑘 \mu_{k} ) of 𝑮 ^ bold-^ 𝑮 \bm{\hat{G}} ordered in ascending magnitude and corresponding approximate singular values ( σ k subscript 𝜎 𝑘 \sigma_{k} ) of 𝒫 0 , τ subscript 𝒫 0 𝜏 \mathcal{P}_{0,\tau} according to ( 31 ). The eigenvalues μ 7 , μ 8 , μ 10 subscript 𝜇 7 subscript 𝜇 8 subscript 𝜇 10 \mu_{7},\mu_{8},\mu_{10} correspond to companion modes. They do not yield purely real singular values σ 7 , σ 8 , σ 10 subscript 𝜎 7 subscript 𝜎 8 subscript 𝜎 10 \sigma_{7},\sigma_{8},\sigma_{10} through the exponentiation ( 31 ), because the numerically computed companions ( 35 ) contain a bias induced by discretization; see [ FK17 , Section 7.3] for further details.

${\mu }_1$	$-0.02523$	${\mu }_6$	$-0.29908$	${\sigma }_1$	$0.79690$	${\sigma }_6$	$0.06776$
${\mu }_2$	$-0.21086$	${\mu }_7$	$-0.03534-0.35003i$	${\sigma }_2$	$0.14990$	${\sigma }_7$	$-0.72750+0.00634i$
${\mu }_3$	$-0.25710$	${\mu }_8$	$-0.03534+0.35003i$	${\sigma }_3$	$0.09887$	${\sigma }_8$	$-0.72750-0.00634i$
${\mu }_4$	$-0.25836$	${\mu }_9$	$-0.39208$	${\sigma }_4$	$0.09776$	${\sigma }_9$	$0.02934$
${\mu }_5$	$-0.29905$	${\mu }_{10}$	$-0.21995-0.33451i$	${\sigma }_5$	$0.06778$	${\sigma }_{10}$	$-0.13695+0.01805i$

Equations

$$x^{\prime}(t)=v(t,x(t))\quad\leadsto\quad\left\{\begin{aligned} \theta^{\prime}(t)&=1,\\ x^{\prime}(t)&=v\big{(}\theta(t),x(t)\big{)}.\end{aligned}\right.$$

$$d{x}_t=v(t,x_t)dt+\epsilon d{w}_t$$

$${\partial }_{t}f(t,x)$$

$$f(0,x)$$

$$\frac{\partial f(t,\cdot )}{\partial n}$$

$$⟨{\mathcal{P}}_{0,t}f,g{⟩}_{L^2(X)}=⟨f,{\mathcal{P}}_{0,t}^{\ast }g{⟩}_{L^2(X)}$$

$$-{\partial }_{t}g(t,x)$$

$$g(\tau ,x)$$

$$\frac{\partial g(t,\cdot )}{\partial n}$$

$$⟨{\nabla }_x(g(t,x)),v(t,x){⟩}_{{\mathbb{R}}^d}$$

$$={\text{div}}_x(g(t,x)v(t,x)).$$

$$-{\partial }_{t}g(t,x)={\text{div}}_x(g(t,x)v(t,x))+\frac{{\epsilon }^2}{2}{\Delta }_{x}g(t,x).$$

$${\partial }_{t}f(t,x)$$

$$f(0,x)$$

$$\frac{\partial f(t,\cdot )}{\partial n}$$

$${\partial }_{t}f(t,x)$$

$$f(\tau ,x)$$

$$\frac{\partial f(t,\cdot )}{\partial n}$$

$$\hat{v}(t,\cdot)=\zeta^{\prime}(t)v(\zeta(t),\cdot)=\left\{\begin{array}[]{ll}v(t,\cdot),&\hbox{$t\in[0,\tau]$;}\\ -v(2\tau-t,\cdot),&\hbox{$t\in(\tau,2\tau]$,}\end{array}\right.$$

$$\zeta(t)=\left\{\begin{array}[]{ll}t,&\hbox{$t\in[0,\tau]$;}\\ 2\tau-t,&\hbox{$t\in(\tau,2\tau]$.}\end{array}\right.$$

$${\partial }_t\hat{f}(t,x)$$

$$\hat{f}(0,x)$$

$$\frac{\partial \hat{f}(t,\cdot )}{\partial n}$$

$$\mathbf{A}=\bigcup _{\theta =0}^{\tau }\{\theta \}\times A_{\theta }\subset [0,\tau ]\times X,$$

$$\mathbf{v}(\mathbf{x})=(1,v(\theta ,x)).$$

$$\hat{A}_{t}=\left\{\begin{array}[]{ll}A_{t},&\hbox{$t\in[0,\tau]$,}\\ A_{2\tau-t},&\hbox{$t\in(\tau,2\tau]$}\end{array}\right.$$

$$\hat{\mathbf{A}}:=\bigcup _{t=0}^{2\tau }\{t\}\times {\hat{A}}_t=\bigcup _{t=0}^{2\tau }\{t\}\times A_{\zeta (t)};$$

$${\int }_0^{\tau }{\int }_{\partial A_t}⟨v(t,x)-b(t,x),n_t(x){⟩}_{{\mathbb{R}}^d}^{+}dS(x)dt,$$

$${\int }_{\partial \mathbf{A}}⟨\mathbf{v}(\mathbf{x}),\mathbf{n}(\mathbf{x}){⟩}_{{\mathbb{R}}^{d+1}}^{+}d\mathbf{S}(\mathbf{x}),$$

$$\hat{b}(t,x)=\left\{\begin{array}[]{ll}\frac{\partial a}{\partial t}(t,r),&\hbox{$t\in[0,\tau]$}\\ -\frac{\partial a}{\partial t}(2\tau-t,r),&\hbox{$t\in(\tau,2\tau]$.}\end{array}\right.$$

$${\int }_0^{2\tau }{\int }_{\partial {\hat{A}}_t}$$

$${\int }_0^{\tau }{\int }_{\partial A_t}|⟨v(t,x)-b(t,x),n_t(x){⟩}_{{\mathbb{R}}^d}|dS(x)dt.$$

$${\int }_{\partial \mathbf{A}}|⟨\mathbf{v}(\mathbf{x}),\mathbf{n}(\mathbf{x}){⟩}_{{\mathbb{R}}^{d+1}}|d\mathbf{S}(\mathbf{x}).$$

$${\int }_0^{2\tau }{\int }_{\partial {\hat{A}}_t}⟨\hat{v}(t,x)-\hat{w}(t,x),{\hat{n}}_t(x){⟩}^{+}dS(x)dt$$

$$={\int }_0^{\tau }{\int }_{\partial A_t}⟨v(t,x)-w(t,x),n_t(x){⟩}^{+}dS(x)dt$$

$$+{\int }_{\tau }^{2\tau }{\int }_{\partial A_{\zeta (t)}}⟨-v(\zeta (t),x)-(-w(\zeta (t),x)),n_{\zeta (t)}(x){⟩}^{+}dS(x)dt$$

$$\text{eqref}outfluxint$$

Full text

1 Introduction

Analysing 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. Non-autonomous time-aperiodic dynamics poses additional difficulties, especially the case of finite time, where asymptotic notions cannot be applied.

The current work has two main contributions.

1. (i)

Extending the work from [FK17], that deals with the time-periodic case, to the situation of general aperiodic finite-time dynamics. We detail a method to compute finite-time coherent sets [Fro13] of aperiodic flows that does not require any time-integration of trajectories. 2. (ii)

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 finite-time coherent sets. This extends optimization results in [FS17], which considered the time-periodic setup of [FK17], to (a) aperiodic dynamics and to (b) infinite-dimensional velocity field space.

1.1 Augmentation

The key construction in [FK17] is the representation of a $\tau$-periodically forced flow on phase space $X\subset\mathbb{R}^{d}$ as an autonomous flow on time-augmented phase space $\tau S^{1}\times X$, where $S^{1}$ denotes the unit circle. On this time-expanded phase space, the time coordinate is simply advanced at a constant rate:

[TABLE]

Finite-time coherent sets on the time interval $[0,\tau]$, as introduced in [FSM10, Fro13], were extracted from singular functions of the transfer operator $\mathcal{P}_{0,\tau}$, 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 $\mathcal{P}_{0,\tau}$ are eigenmodes of $\mathcal{P}_{0,\tau}^{*}\mathcal{P}_{0,\tau}$, where $\mathcal{P}_{0,\tau}^{*}$ denotes the dual of $\mathcal{P}_{0,\tau}$, and that for area-preserving dynamics (corresponding to divergence-free velocity fields $v$) the dual is the transfer operator of the time-reversed dynamics, i.e., the (again, slightly stochastically perturbed) flow governed by the time-reflected velocity field $(t,x)\mapsto-v(\tau-t,x)$. The dynamic interpretation of this operator-based characterization is that (finite-time) coherent are those sets that are to a large extent mapped back to themselves by a noisy forward-backward evolution of the dynamics. This operator-based framework not only gives a qualitative framework for coherence; the singular values of $\mathcal{P}_{0,\tau}$ 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 “forward-time” and “backward-time” velocity fields on time intervals $[0,\tau]$ and $[\tau,2\tau]$—see (9) below—we will construct a system on the augmented space $[0,2\tau]\times X$ that mimics the forward-backward evolution of the dynamics. In this way, the eigenmodes of the (Fokker–Planck-) generator $\bm{\hat{G}}$ of this augmented system yield the eigenmodes of $\mathcal{P}_{0,\tau}^{*}\mathcal{P}_{0,\tau}$; i.e., the singular modes of $\mathcal{P}_{0,\tau}$, and from these singular modes the desired finite-time coherent sets. Again, the corresponding eigenvalues of the generator can be used to give quantitative bounds on coherence/mixing, see Theorem 8. We note that a numerical approach to extract coherent sets from $\mathcal{P}_{0,\tau}$ by solving the Fokker–Planck equation has been described in [DJM16]. In contrast to [DJM16] we do not require time-integration over $[0,\tau]$, which is especially advantageous once we consider the optimization of coherence and mixing.

Connecting the spectral properties of the generator of the augmented-space system with the (finite-time) dynamical properties of the original system is a generalization of the results from [FK17], where it has been done for time-periodic velocity fields on infinite time intervals. Through time-reflection of the finite-time problem we construct a time-periodic one, to which we apply the concepts of [FK17]. We also remove an assumption from [FK17] (the “niceness”) on the so computed sets, thus strengthening the approach. There are several further non-trivial adjustments needed to fit the theory of [FK17] to this time-reflected setting, and the necessary details are covered in sections 2 and 4.

The interplay between the spectra of the dynamics in augmented space and the non-autonomous dynamics in the original space has strong connections to the correspondence between evolution semigroups [How74] and two-parameter evolution families, as elaborated in, e.g., [CL99]; see also [EN00, Section VI.9] for a general introduction. We mention that by a similar construction, spatio-temporal 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 augmented-space dynamics.

1.2 Manipulation of coherence and mixing

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 [MMG*+*07, CAG08, OBPG15]. Other strategies include optimising the diffusion component of the dynamics [FGTW16], the optimal distribution of concentration sources [TP08] and geometric dynamical systems techniques [Bal15]. An interesting theoretical result is that arbitrary mixedness under advection-diffusion 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, finite-time coherent sets are characterized by the singular vectors of the transfer operator, and, equivalently, by the eigenvectors of the generator of the augmented-space 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 $\bm{\hat{G}}$, and can target these eigenvalues if we want to enhance or diminish mixing.

Given a default velocity field $v$ and non-autonomous perturbations $u\in U$ from an admissible space $U$ of divergence-free velocity fields, our approach considers “small” $u$ that change a targeted eigenvalue $\mu$ of $\bm{\hat{G}}$ (thus also the singular value $\sigma$ of $\mathcal{P}_{0,\tau}$) locally optimally. This procedure can be iterated to obtain a larger perturbation in a gradient-method fashion. Optimizing singular values of the transfer operator $\mathcal{P}_{0,\tau}$ directly is difficult as it would necessarily involve the variation of the nonlinear dynamics under the velocity field $u$. Instead, a linearized optimization of the eigenvalue of the generator $\bm{\hat{G}}$ leads to a very simple optimization problem (44), which can be solved via a linear system of the same dimension as $U$. Moreover, the theory holds for infinite-dimensional perturbation spaces $U$ as well, since Fréchet differentiability of the transfer operator and its spectrum with respect to perturbing velocity fields has been established in [KLP19].

1.3 Overview

This work is structured as follows. In section 2 we introduce the $L^{2}$-function based formalism to study advection-diffusion 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 operator-based considerations that follow. Section 4 introduces the new reflected augmented generator needed to handle aperiodic, finite-time driving. We state formal connections between the spectrum of the reflected augmented generator and the reflected transfer operator, and provide spectral-based 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 trajectory-free 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 time-dependent local perturbation of the velocity field. Section 7 specializes the infinite-dimensional 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 Advective-diffusive dynamics

Let $X\subset\mathbb{R}^{d}$ be a bounded and open set with compact and smooth (piecewise $C^{4}$) boundary. We consider the time interval $[0,\tau]$ and the dynamics

[TABLE]

with reflecting boundary conditions for $v\in C^{(1,1)}([0,\tau]\times\overline{X};\mathbb{R}^{d})$111$C^{(1,1)}([0,\tau]\times\overline{X};\mathbb{R}^{d})$ denotes the Banach space of functions $f\,:\,[0,\tau]\times\overline{X}\rightarrow\mathbb{R}^{d}$ that are continuously differentiable in $t$ and cotinuously differentiable in $x$.and $(w_{t})_{t\geq 0}$ being a standard Wiener process in $\mathbb{R}^{d}$. The initial point $x_{0}$ is distributed according to some initial density $f_{0}\in L^{2}(X)$. 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],

[TABLE]

where $\tfrac{\partial}{\partial n}$ is the normal derivative on the boundary. Associated to (2) is an evolution operator $\mathcal{P}_{0,t}:L^{2}(X)\to L^{2}(X)$ that transports a density $f_{0}\in L^{2}(X)$ at time 0 to the solution density of (2) at time $t$. The evolution operator $\mathcal{P}_{0,t}$ is an integral operator with stochastic222Doubly stochastic if the flow is volume-preserving. kernel $k(t,\cdot,\cdot):X\times X\to\mathbb{R}^{+}$ that satisfies [Fro13, Assumptions 1 and 2].

2.1 Construction of a forward-backward process

For simplicity of presentation we assume that the velocity field $v(t,\cdot)$ is divergence free for all $t\in[0,\tau]$. We note that the remaining arguments in this section may be carried through for general velocity fields. Denote by $\langle\cdot,\cdot\rangle_{H}$ the canonical scalar product of a Hilbert space $H$. Following [Fro13] in the volume-preserving setting333The operators $\mathcal{P}_{0,\tau}$ and $\mathcal{P}_{0,\tau}^{*}$ are replaced by normalised versions for nonzero divergence velocity fields; these are denoted by $\mathcal{L}$ and $\mathcal{L}^{*}$ in [Fro13]., coherent sets over the time interval $[0,\tau]$ are extracted from the eigenfunctions of $\mathcal{P}_{0,\tau}^{*}\mathcal{P}_{0,\tau}$ corresponding to large eigenvalues, where $\mathcal{P}_{0,\tau}^{*}$ is the $L^{2}$-adjoint of $\mathcal{P}_{0,\tau}$, defined to be the unique linear operator satisfying

[TABLE]

for all $f,g\in L^{2}(X)$. The eigenvalues of $\mathcal{P}_{0,\tau}^{*}\mathcal{P}_{0,\tau}$ (the singular values of $\mathcal{P}_{0,\tau}$) are known to lie in the interval $[0,1]$ (cf. [Fro13, p.3]). The rationale behind the operator $\mathcal{P}_{0,\tau}^{*}\mathcal{P}_{0,\tau}$ is that $\mathcal{P}_{0,\tau}$ describes evolution in forward time, $\mathcal{P}_{0,\tau}^{*}$ describes evolution under the time-reversed dynamics, and coherent sets are characterized exactly by the property that they are “stable” under a noisy forward-backward evolution of the dynamics.

The adjoint operator $\mathcal{P}_{t,\tau}^{*}$ is the solution operator to the Kolmogorov backward equation [PS08]:

[TABLE]

The operator $\mathcal{P}_{t,\tau}^{\ast}$ maps a density $g_{\tau}$ at time $t=\tau$ backward in time according to (4) to produce a density $g_{t}$ at time $t<\tau$. We may simplify (4) using volume preservation:

[TABLE]

Thus we may write (4) as

[TABLE]

Reversing time in (4) to obtain an initial value problem we get

[TABLE]

using $f(t,x)=g(\tau-t,x)$, $\bar{f}_{0}(x)=g_{\tau}(x)$ and the velocity field $\bar{v}(t,x)=-v(\tau-t,x)$. 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 $[0,2\tau]$ that corresponds to the operator $\mathcal{P}_{0,\tau}^{*}\mathcal{P}_{0,\tau}$. We view $\mathcal{P}_{0,\tau}^{*}$ as evolution on the time interval $[\tau,2\tau]$ and we therefore shift (7) by $\tau$ time units, defining $\tilde{v}(t,x):=\bar{v}(t-\tau,x)=-v(2\tau-t,x)$ to obtain a forward problem on $[\tau,2\tau]$:

[TABLE]

We denote the solution operator of this problem as $\tilde{\mathcal{P}}_{\tau,t}$ ($=\mathcal{P}^{\ast}_{2\tau-t,\tau}$).

Finally, we concatenate the two forward problems (2) and (8) to make a single process over $[0,2\tau]$. We mark objects that live on this extended interval $[0,2\tau]$ with a hat $\hat{\phantom{v}}$. Define the velocity field

[TABLE]

using the reflection map

[TABLE]

The resulting velocity field $\hat{v}$ exhibits discontinuities in $0,\tau$ and $2\tau$ whenever it does not vanish there, but one-sided derivatives exist444The following regularity may not be most general but is meant to give some intuition: $\hat{v}\in C^{(1,1)}(((0,\tau)\cup(\tau,2\tau))\times\overline{X};\mathbb{R}^{d})$ and $\hat{v}\in L^{p}((0,2\tau)\times X;\mathbb{R}^{d})$, for $1\leq p\leq\infty$. In what follows, we will solve the Fokker–Planck equation

[TABLE]

over the interval $t\in[0,2\tau]$; more precisely on $(0,\tau)\cup(\tau,2\tau)$ with $L^{2}$-continuous concatenation at $t=\tau$ (see Prop. 5). Let us summarize the above construction with the following proposition.

Proposition 1**.**

The concatenation $\mathcal{P}_{0,\tau}^{*}\mathcal{P}_{0,\tau}=:\hat{\mathcal{P}}_{0,2\tau}$ with $\tilde{f}_{\tau}=\mathcal{P}_{0,\tau}f_{0}$ comprises initializing (2) at time 0, solving forward using the vector field $v(t,\cdot)$ until time $\tau$, then continuing to evolve (8) for another $\tau$ time units, but now using the reflected and shifted vector field $-v(2\tau-t,\cdot)$ for $t\in[\tau,2\tau]$ corresponding to the reversed dynamics.

3 Cumulative flux from a reflected family of sets

Before proceeding with the operator-based description of finite-time coherence, in this section we analyse 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., space-time) in Proposition 3 is partially to set the stage for the augmented-space operator-based description in section 4. Apart from strengthening the intuition for the forward-backward construction that is used in here, the results of the section 3.2 are not formally needed for the rest of the paper.

3.1 Augmentation and reflection

For a family of sets $\{A_{t}\}_{t\in[0,\tau]}$, $A_{t}\subset X$, we consider the augmented set

[TABLE]

in the augmented state space $\bm{X}:=[0,\tau]\times X$. Let $\bm{n(x)}$ denote the unit outer normal on $\partial\bm{A}$ at $\bm{x}\in\bm{X}$ and $\bm{v}$ the augmented velocity field defined by

[TABLE]

We define a reflected family of sets $\{\hat{A}_{t}\}_{t\in[0,2\tau]}=\{A_{\zeta(t)}\}_{t\in[0,2\tau]}$ in synchrony with the reflected vector field $\hat{v}$:

[TABLE]

and the augmented reflected set

[TABLE]

see Figure 1.

3.2 Outflow flux

We consider a family of $d$-dimensional sets $\{A_{t}\}_{t\in[0,\tau]}$, $A_{t}\subset X\subset\mathbb{R}^{d}$ satisfying the following assumptions (the boundaries are piecewise smooth in space and differentiable in time):

Assumption 2**.**

There exists a co-dimension 1 parameterisation set $R\subset\mathbb{R}^{d-1}$ such that for each $t\in[0,\tau]$ there is a bijective function $a(t,\cdot)\,:\,R\rightarrow\partial A_{t}$ with $a$ being piece-wise smooth, $a\in C^{(1,0)}([0,\tau]\times R;\mathbb{R}^{d})$ 555$C^{1}$ in $t$ and $C^{0}$ in $r$ piecewise. 2. 2.

The mapping $b(t,x):=\frac{\partial a}{\partial t}(t,r)$ is well defined for all $t\in[0,\tau]$ and all $x\in\partial A_{t}$, where $a(t,r)=x$.

The cumulative outflow flux under the vector field $v(t,\cdot),\,t\in[0,\tau]$ from a family of sets $\{A_{t}\}_{t\in[0,\tau]}$ is given by

[TABLE]

where $(\cdot)^{+}$ denotes the positive part, $S(x)$ is the $d-1$ dimensional surface measure and $n_{t}(x)$ is the outer normal unit vector. The result [FK17, Theorem 2] shows that (14) is equal to the instantaneous outflow flux from the set $\bm{A}$, defined by

[TABLE]

with $\bm{S}(\bm{x})$ denoting the $d$-dimensional surface measure.

We extend this result to the reflected velocity field $\hat{v}(t,\cdot),t\in[0,2\tau]$, generated by a general aperiodic velocity field.

In particular, recalling $\zeta$ from (10), for every $t\in[0,2\tau]$ the boundary $\partial\hat{A}_{t}$ has a parametrization $\hat{a}(t,\cdot)=a(\zeta(t),\cdot)\,:\,R\rightarrow\partial\hat{A}_{t}$ and

[TABLE]

At $t=\tau$ the right- and the left-sided partial derivatives of $a$ with respect to $t$ exists but they may not be equal. The family of normal vectors is mirrored in time (see Figure 1): $\hat{n}_{t}(x)=n_{\zeta(t)}(x)$ for $t\in[0,2\tau]$.

Proposition 3**.**

The cumulative outflow flux from the family of sets $\hat{A}_{t},t\in[0,2\tau]$ under the vector field $\hat{v}(t,\cdot),t\in[0,2\tau]$, is equal to the cumulative absolute flux in and out of the family of sets $A_{t}$, $t\in[0,\tau]$, under the vector field $v(t,\cdot),t\in[0,\tau]$; that is,

[TABLE]

Furthermore, the time-integrated flux (17) is equal to the instantaneous absolute flux in augmented space:

[TABLE]

Proof.

Let us first prove the first equality. Therefore we do not need objects of the augmented setting yet. We split the integral

[TABLE]

and see that the first integral is already a part of what is needed. So we only need to treat the second integral (19). We use Fubini’s theorem and substitute using $\zeta$ with $g_{\zeta(t)}(r)$ being the Gram determinant. (It is important to note that we need and use substitution in one dimension, the time dimension, because we need the sign we get from substitution, which we would not get in higher dimensions.)

[TABLE]

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 so-called spectral mapping theorem for our reflected augmented process (Proposition 5) and derive a bound for the finite-time coherence of a family of sets $\{A_{t}\}_{t\in[0,\tau]}$ in terms of the second eigenvalue of a generator (the infinitesimal operator) of our augmented reflected advection-diffusion process (Theorem 8). To do this we build on discrete-time theory from [Fro13] with the periodic continuous-time theory from [FK17].

4.1 The evolution operator for the reflected process

It is well known that $\mathcal{P}_{s,s+t}$, $t>0$, is a compact, integral preserving, real and positive operator on $L^{2}(X)$ while $t\mapsto\mathcal{P}_{s,s+t}f$ is continuous as a mapping from $[0,\infty)$ to $L^{2}(X)$ for any fixed $f\in L^{2}(X)$666See 19 for compactness and continuity, and [LM13] or [KKS16] for the other properties.. Furthermore $\smash{\hat{\mathcal{P}}_{0,2\tau}=\mathcal{P}^{\ast}_{0,\tau}\mathcal{P}_{0,\tau}}$ is a self-adjoint operator on $L^{2}(X)$ with simple largest eigenvalue $\lambda_{1}(\hat{\mathcal{P}}_{0,2\tau})=1$. Following [Fro13] one has that the second eigenvalue $\smash{\lambda_{2}(\hat{\mathcal{P}}_{0,2\tau})}$ satisfies

[TABLE]

where in the volume-preserving777In the nonzero divergence case, $\mu_{0}$ is a reference measure describing the initial mass distribution of the (possibly compressible) fluid being evolved, and $\nu_{\tau}$ is the forward evolution of $\mu_{0}$ under a normalised version of $\mathcal{P}_{0,\tau}$ denoted by $\mathcal{L}$ in [Fro13]. setting $\mu_{0}$ and $\nu_{\tau}$ are both simply the Lebesgue measure. We now consider the problem (11) introduced in section 2 as a time-periodic problem on $2\tau S^{1}\times X$ (we extend $\hat{v}$ periodically). Following the considerations of section 2 the evolution operator $\hat{\mathcal{P}}_{s,s+t}$ starting from time $s$, w.l.o.g. $s\in[0,2\tau]$, flowing for time $t\geq 0$ to $s+t=k\tau+r$, $k=\lfloor\frac{s+t}{\tau}\rfloor\in\mathbb{N}\cup\{0\}$ and $(s+t)\mod\tau=r\in[0,\tau)$, is given by

[TABLE]

The situation when $t$ is exactly $2\tau$ is of particular importance:

[TABLE]

Note that $\hat{\mathcal{P}}_{s,s+2\tau}$ is self-adjoint for $s=k\tau$, $k\in\mathbb{N}\cup\{0\}$.

4.2 The time-augmented generator and evolution family

We now turn to the augmented reflected system

[TABLE]

in $\bm{\hat{X}}:=2\tau S^{1}\times X$, and note that $(\hat{w}_{t})_{t\geq 0}$ is a standard Wiener process in $\mathbb{R}^{d}$; in particular it is not constructed by time reflection. We define augmented versions of $\hat{x}_{t},\hat{v}$, $\varepsilon$, and $\hat{w}_{t}$, denoting them with bold symbols:

[TABLE]

and $\bm{\hat{w}}_{t}$ is a $d+1$ dimensional standard Wiener process. The augmented system is

[TABLE]

Considering the time-periodic version of problem (11) with $\hat{v}$ on $\bm{\hat{X}}$ we formulate a Fokker–Planck equation in augmented space. To avoid confusion with the (new) state variable $\theta$, we write dependence on time $t$ as a subscript of the augmented function, i.e., $\bm{f}_{t}\,:\,\bm{\hat{X}}\to\mathbb{R}$ for all $t\geq 0$. The augmented Fokker–Planck equation is

[TABLE]

Note there is no diffusion in the $\theta$ direction, as per the definition of $\bm{\varepsilon}$.

We will now consider the augmented Fokker–Planck equation as a linear differential equation in the space $L^{2}(\bm{\hat{X}})$. Its right-hand side is given by the so-called (augmented) infinitesimal generator $\bm{\hat{G}}\,:\,\mathcal{D}(\bm{\hat{G}})\subset L^{2}(\bm{\hat{X}})\rightarrow L^{2}(\bm{\hat{X}})$, with domain $\mathcal{D}(\bm{\hat{G}})$ defined as the subspace of $L^{2}(\bm{\hat{X}})$ on which the generator is well-defined in terms of semigroup theory [Paz83, EN00]. The augmented Fokker–Planck equation (24) in augmented space then reads as

[TABLE]

We may also write (24) and (25) in terms of the non-autonomous (“unaugmented”) dynamics:

[TABLE]

where $\hat{G}({\theta})$ is the right-hand side operator of (11) at time $t=\theta\in 2\tau S^{1}$, i.e., the time-$\theta$ differential operator of the Fokker–Planck equation (on $[0,2\tau]$).

Remark 4**.**

1. (a)

The periodicity in $\theta$ and the boundary conditions from (25) are, as is common in semigroup theory, encoded in the domain of $\smash{\bm{\hat{G}}}$. This domain $\smash{\mathcal{D}(\bm{\hat{G}})}$ enforces continuity conditions in $\theta$ at $0,\tau,2\tau$. 2. (b)

For the purposes of the current work we will not require the well-posedness of the problem (25). We will only need the operator $\bm{\hat{G}}$ defining the right-hand side of this equation, and its relation to the transfer operator family $\hat{\mathcal{P}}_{s,s+t}$, $s<t$. This will be the focus of section 4.3. For additional theory for these augmented problems, objects, and related results we refer the reader to **[CL99]**. 3. (c)

Any non-constant solution $\bm{f}\,:\,(t,\theta,x)\mapsto\bm{f}_{t}(\theta,x)$ to (25) has three input variables, $t$, $\theta$ and $x$, and may have different regularity properties in each variable. In the following we will focus on eigenfunctions of $\bm{\hat{G}}$ that are of course constant in $t$.

We now comment on the crucial connection between solutions of (25) and $\hat{\mathcal{P}}_{s,s+t}$, neglecting the issue of solvability. Note that the stochastic augmented differential equation (23) allows for evolving the non-autonomous equation (1) from any initial time $s$ by setting $\theta_{0}=s$. In an analogous manner, the augmented Fokker–Planck equation (25) with initial condition $\bm{f}_{0}$ evolves every initial condition $\bm{f}_{0}(s,\cdot)$, $s\in[0,2\tau)$,—i.e., a configuration of initial conditions—by the non-autonomous reflected Fokker–Planck equation (11). More precisely, the following holds for the evolution of (25):

[TABLE]

In the terminology of semigroup theory [CL99, EN00] the solution operators of (25), here formally denoted by $\smash{(e^{t\bm{\hat{G}}})_{t\geq 0}}$, form an evolution semigroup (or Howland semigroup), and it is given exactly by (27). Informally, the action of $\smash{e^{t\bm{\hat{G}}}}$ in the context of $\hat{\mathcal{P}}_{\theta,\theta+t}$ can be described as follows. On the left-hand side of (27), $\smash{e^{t\bm{\hat{G}}}}$ takes the initial configuration $\bm{f}_{0}$ (on all of $\bm{\hat{X}}$) and evolves the entire configuration for the time duration $t$, to obtain $\bm{f}_{t}$. The result is then evaluated at the $\theta+t\mod 2\tau$ fiber. The $\theta+t\mod 2\tau$ fiber of $\bm{f}_{t}$ corresponds to the $\theta$ fiber of $\bm{f}_{0}$ evolved for time $t$ due to the constant drift in the $\theta$ variable: $-\partial_{\theta}$, see (26). That equation (27) indeed gives the solutions to (25), is a conseqence of [CL99, Theorem 6.20] adapted to the concatenation of the forward and backward evolutions described by the reflected system (11). As we will not require a result of this generality, we omit the details. For our purposes it will be sufficient to consider the special case, where $\bm{f}_{0}=\bm{f}$ is an eigenfunction of $\bm{\hat{G}}$. This is done next.

4.3 Eigenfunctions of the time-augmented generator

Analogously to [FK17, Lemma 22] the following result holds. It can be obtained from (27) by noting that every eigenpair $(\mu,\bm{f})$ of $\bm{\hat{G}}$ gives a solution to (25) by $\bm{f}_{t}=e^{\mu t}\bm{f}$. However, to highlight the intuitive connection between the augmented generator and the non-autonomous problem, we prove the following proposition by other means.

Proposition 5**.**

Let $\bm{f}$ be an eigenfunction of $\bm{\hat{G}}$ corresponding to the eigenvalue $\mu\in\mathbb{C}$. One has then

[TABLE]

for all $s\in 2\tau S^{1}$ and $t\geq 0$.

Proof.

We will prove (28) following ideas from [FK17]. First we apply Theorem 19 piecewise on $[0,\tau]$ and $[\tau,2\tau]$ concerning well-posedness and regularity. Therefore we consider the original problem (2) and the reflected, shifted, time-reversed problem (11). Now theorem 19 guarantees for any initial condition $f_{0}\in L^{p}(X)$, $p\in(1,\infty)$, the unique existence of a function $f$ with the regularity

[TABLE]

and the properties

[TABLE]

Further $f(t)\in\mathcal{D}(\hat{G}(t))$ holds for all $t\in(0,2\tau]$. Now we can proceed as in [FK17, Lemma 22].

Let $\mu\in\mathbb{C}$ and $\bm{f}\in\mathcal{D}(\bm{\hat{G}})$ with $\bm{\hat{G}}\bm{f}=\mu\bm{f}$. According to the construction above (26) we know

[TABLE]

and this implies, in accordance with (25), for all $\theta\in 2\tau S^{1}\backslash\{0,\tau\}$

[TABLE]

Now $\hat{\mathcal{P}}_{\theta,\theta+t}$ is the evolution operator to the evolution equation

[TABLE]

Therefore the function $e^{-\mu t}\hat{\mathcal{P}}_{\theta,\theta+t}\bm{f}(\theta)$ solves (29) uniquely and 19 guarantees continuity in $\theta$. Therefore we can connect the eigenfunctions of the augmented reflected generator $\bm{\hat{G}}$ with the evolution given by $\hat{\mathcal{P}}_{s,s,+t}$ for all $s,t$ as stated in the claim. ∎

Let $\mu$ be an eigenvalue of $\bm{\hat{G}}$ with an eigenfunction $\bm{f}$. Inserting $s=0$ and $t=2\tau$ into Proposition 5 yields

[TABLE]

This is a spectral mapping theorem type of result, as it connects the eigenvalues and eigenfunctions of the evolution operator $\hat{\mathcal{P}}_{s,s+2\tau}$ with those of an associated (infinitesimal) generator $\bm{\hat{G}}$. We refer the reader to standard literature on classical results for operator semigroups [Paz83, EN00].

Recalling that $\hat{\mathcal{P}}_{0,2\tau}=\mathcal{P}_{0,\tau}^{\ast}\mathcal{P}_{0,\tau}$ is a compact self-adjoint positive operator, it must be that $e^{\mu 2\tau}=\sigma^{2}$ for some $0<\sigma\in\mathbb{R}$. This implies

[TABLE]

from which it follows that $\Im(\mu)=\frac{k\pi}{\tau}$ for some $k\in\mathbb{Z}$.

Remark 6**.**

Theorem 19 guarantees that for initial conditions $f_{s}\in\mathcal{D}(G(s))$ the solution $\hat{\mathcal{P}}_{s,s+t}f_{s}=f(t)$ is in $t$ a continuous mapping to the domain of the generator, $\mathcal{D}(G(s+t))=\mathcal{D}_{p}$. Theorem 20 further gives for each eigenfunction $\bm{f}$ that $\bm{f}:\theta\mapsto\bm{f}(\theta,\cdot)\in C(2\tau S^{1};\mathcal{D}_{p})$ and that $\bm{f}\in C(\bm{\hat{X}})$. This regularity is utilized in the proof of Theorem 8 below.

4.4 Coherent families of sets

In the specific case where the velocity field $v$ is periodic in time, [FK17] shows that the families of sets

[TABLE]

have an escape rate (see [FK17, Definition 8]) of at most $\text{Re}(\mu_{2})$, where $\mu_{2}$ is the first nontrivial eigenvalue of $\bm{\hat{G}}$ corresponding to the eigenfunction $\bm{f}$. 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 $v$ we will quantify the coherence of families $\{A^{\pm}_{\theta}\}_{\theta\in[0,\tau]}$ and provide a construction of highly coherent families with associated rigorous coherence bound.

Definition 7** (Coherence ratio).**

Let $\{A_{t}\}_{t\in[0,\tau]}$ be a family of measurable sets. Denote by $\mathbb{P}_{m}$ the law of the process $\{x_{t}\}_{t\in[0,\tau]}$ generated by the SDE (1) initialised with $x_{0}\sim m$, where $m$ denotes normalised Lebesgue measure on $X$. For $m(A_{0})\neq 0$ we define the coherence ratio of the family $\{A_{t}\}_{t\in[0,\tau]}$ as

[TABLE]

It was shown in [FK17, Appendix A.6] that for a family of sets with sufficient regularity (called “sufficient niceness” therein) the quantity (32) is well defined. Here we will alleviate this requirement entirely by showing regularity of the augmented eigenfunctions $\bm{f}$, and showing that this is sufficient to prove the desired results.

Theorem 8 makes a link between the coherence of a particular family of sets defined by zero super/sublevel sets of an eigenfunction of $\bm{\hat{G}}$ and the corresponding eigenvalue $\mu$. It shows that the probability of a trajectory remaining in a family of sets constructed from the positive and negative parts of eigenfunctions of $\bm{\hat{G}}$ decays no faster than the rate given by the corresponding eigenvalues. This result extends to aperiodically driven continuous-time systems, similar results for autonomous systems in discrete time [FS10] and continuous time [FJK13], and periodically driven dynamics in continuous time [FK17].

Theorem 8**.**

Let $\bm{\hat{G}}\bm{f}=\mu\bm{f}$ with $\mu<0$. If $\bm{f}$ is scaled such that $\|\bm{f}(\tau,\cdot)\|_{L^{1}}=2$, then it holds for the family $\{A_{t}^{\pm}\}_{t\in[0,\tau]}$ of sets with $A_{t}^{\pm}=\{\pm\bm{f}(t,\cdot)\geq 0\}$ that

[TABLE]

where $|A|$ denotes the non-normalised Lebesgue measure of the set $A$. In particular, eigenfunctions at eigenvalues $\mu\approx 0$ yield families of sets with high coherence ratio.

Proof.

See Appendix B. We note that by Corollary 21 $\bm{f}$ is continuous on $\bm{\hat{X}}$. ∎

Intuitively, the left hand side of (33) quantifies the likelihood of escape from the family of sets, and the right hand side of (33) is a (scaled) measure of mixing; the bound says that the likelihood of escape is less than the mixing incurred over the same time duration. The bound is not intended to be sharp; we remark that one could optimise the level set cutoff to improve the ratio $\rho_{m}$, as has been done in previous work on coherent sets [FP09].

In Section 5.3, we consider the dominant 6 eigenvectors of $\bm{\hat{G}}$ and apply sparse eigenbasis approximation (SEBA) [FRS18] to find a sparsity-inducing rotation of this eigendata and separate individual slow escape / slow mixing subdomains. The following proposition generalises Theorem 8 so that it may apply to vectors formed from linear combinations of eigenvectors.

Proposition 9**.**

Let $\bm{\hat{G}}\bm{f}_{i}=\mu\bm{f}_{i}$, $i=1,\ldots,k\in\mathbb{N}$, with $\mu_{k}\leq\ldots\leq\mu_{1}\leq 0$. For $\bm{f}=\sum_{i=1}^{k}\alpha_{i}\bm{f}_{i}$ with $\alpha_{i}\in\mathbb{R}$, the statement of Theorem 8 remains true with $\mu=\mu_{k}$ if

[TABLE]

The proof Proposition 9 is deferred to Appendix B.

In the computations performed in the next sections we will numerically approximate $\bm{\hat{G}}$, compute its upper spectrum and associated eigenfunctions, and plot super/sublevel sets of the eigenfunctions. In Section 5.3 we will additionally apply SEBA, plot the sparse basis functions, and one of the superlevel sets.

5 Computational aspects

5.1 Numerical discretization

We use the “Ulam’s discretization 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 discretization for the generator yields a matrix that may be interpreted as a rate matrix of a finite-state, continuous-time Markov chain, with the states corresponding to a partition of $\tau S^{1}\times X$ into hypercubes (hyperrectangles) in $\mathbb{R}^{d+1}$. 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 $1/h$, where $\tau S^{1}$ is discretized into intervals of length $h$. The entries in the remaining $d$ 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 $G_{n}^{\mathrm{drift}}$ in [FK17, Section 7.2]. The entries of the rate matrix corresponding to diffusive dynamics (in the $d$ space coordinates only, there is no diffusion in the time coordinate) are computed from a finite-difference approximation of the Laplace operator; see the expression for $G_{n}^{\mathrm{diff}}$ in [FK17, Section 7.2]. We then set $G_{n}:=G_{n}^{\mathrm{drift}}+G_{n}^{\mathrm{diff}}$. The matrix $G_{n}$ can also be interpreted as a rate matrix for a finite-state Markov chain; it has an eigenvalue [math] and its spectrum is confined to the left half of the complex plane.

The reflected velocity field $\hat{v}$ from (9) may be substituted for the velocity field $v$ used in [FK17] and the methodology of [FK17] employed; this is the approach taken in the numerical experiments below.

Remark 10**.**

The computation of $G_{n}^{\mathrm{drift}}$ uses only the outward-pointing 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 (15)—, discarding the inward-pointing parts. Because of the reflected structure of $\hat{v}$, 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 outward-pointing direction). The outward-pointing components would be used on the time interval $[0,\tau]$, while the inward-pointing components would be used on the time interval $(\tau,2\tau)$, where they are outward-pointing 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. However, the assembly of the generator matrices is relatively fast anyway, and we have not tried to optimize our implementation of Ulam’s method for the generator in this reflected setting.

5.2 Example: Periodically driven double gyre

We consider the periodically driven double gyre system [SLM05]:

[TABLE]

on the time interval $[0,\tau]=[0,4]$. The forcing is $f(t,x)=\gamma\sin(2\pi\Omega t)x^{2}+(1-2\gamma\sin(2\pi t))x$ and the parameters are $A=0.25$, $\Omega=2\pi$, and $\gamma=0.25$, implying the forcing period $1$, on the spatial domain $X=[0,2]\times[0,1]$. This system has been a standard example of coherent sets [FPG14]. The purpose of this section is to show that our method reliably computes the singular functions and values of $\mathcal{P}_{0,\tau}$; further analysis is deferred to later sections. In particular, we will revisit this example in context of optimal manipulation of these coherent sets in sections 7.3 and 7.4.

The augmented reflected generator approach with a resolution of $40\times(100\times 50)$ and noise intensity $\varepsilon=0.1$ gives us the non-trivial dominant eigenvectors of $\bm{\hat{G}}$ at time $t=0$ shown in Figure 2.

Ordered by ascending magnitude, the $3$rd, $4$th and $6$th eigenvalues (Table 1) and eigenvectors (not shown) correspond to features also detected in [FK17], where they were connected to complex non-companion eigenvalues—the concept of companion eigenvalues will be introduced around (35) below. These features become less coherent, i.e., their respective real eigenvalues decrease compared with the others, as the length $\tau$ of the time interval increases.

5.3 Example: Bickley jet

We now apply the reflected augmented generator approach to a perturbed Bickley Jet [RBBV*+*07]. 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 $v=(-\frac{\partial\Psi}{\partial y},\frac{\partial\Psi}{\partial x})$ is induced by the stream function

[TABLE]

The constants are chosen according to [RBBV*+*07]. The length unit is $Mm$ ($1\,Mm=10^{6}m$) and the time unit is days. For the amplitudes $A_{n}$ and the speed of the Rossby waves $c_{n}$ we choose

[TABLE]

with $r_{e}$ being the Earth’s radius. Further we choose

[TABLE]

The state space is periodic in $x$ direction and is given by $X=\pi r_{e}S^{1}\times[-3,3]$ (in accordance with other literature), and the time interval will be chosen as $[0,\tau]=[0,9]$. For good numerical tractability, we resolve our reflected space-time manifold with a spatially somewhat coarse $108\times(120\times 36)$ grid, that is uniform in space ($120\times 36$) leads to square boxes needed for isotropic diffusion) and sufficiently finely resolved in time ($108$). We choose $\varepsilon=0.1$.

The system described above is equipped with homogeneous Dirichlet boundary conditions instead of homogeneous Neumann conditions on $\partial X$. This leads to a slightly different spectral structure of the generator, which now generates a semigroup of sub-Markovian operators. Its leading eigenvalue is strictly less than zero. We expect Theorem 8 to hold in the case of Dirichlet boundary conditions. One possible theoretical justification would require “close” the open system that is represented by the homogeneous Dirichlet boundary conditions by introducing a virtual “external” state, then apply the Neumann theory to that system. The details would lead beyond the scope of this work, and will be discussed elsewhere.

We highlight that the computations we are about to perform here are different to those performed in [FK17, Section 7.6] in at least two respects. Firstly, the Bickley jet under investigation is aperiodically driven, in contrast to the periodically driven Bicklet jet in [FK17, Section 7.6] (we use slightly different parameters in the velocity field). Secondly, we wish to find functions that decay the least under finite-time evolution, in contrast to the problem considered in [FK17], which sought functions that decayed at the slowest time-asymptotic ($t\to\infty$) rate under periodic driving. In particular, even for a periodically driven Bickley jet, the finite-time question considered in the present paper is different to the infinite-time question addressed in [FK17]. Thus, even though we chose a flow interval of length 9 as in [FK17], the problem in consideration is different.

Analogously to [FK17, Section 7] our time-augmentation 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 Ulam-discretization we refer to [FK17, Section 7.3]. We will use and verify the relations derived there. Therefore we calculate the eigenvalues and vectors of $\bm{\hat{G}}$ with the smallest magnitude instead of largest real part using eigs(G,10,’SM’) in Matlab.

Table 2 shows a gap after the first and sixth eigenvalue. Let us first discuss the leading $6$ eigenvectors. Figure 3 shows the eigenvectors corresponding to the dominant (i.e., smallest real part) $6$ eigenvalues. The first eigenvector, the quasistationary (or conditionally invariant) distribution highlights (in blue, see Figure 3(a)) parts of the domain that get pushed out of the region $X$ of consideration. The red regions are those parts of phase space that remain longest in $X$ 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 $\bm{\hat{G}}$ contain genuinely new dynamical information, we checked that the complex eigenvalues of $\bm{\hat{G}}$ from $\mu_{7}$ to $\mu_{50}$ are all companion eigenvalues equal to (recall from section 5.1 that $h$ is the temporal grid spacing)

[TABLE]

for an eigenvalue $\mu$ of $\bm{\hat{G}}$ and a $k\in\mathbb{Z}$, as derived in [FK17, Section 7.3]. Under the assumption that the eigenvector $\bm{w}$ is sufficiently smooth in time and time is sufficiently resolved (i.e., $w_{t}\approx w_{t-h}$), each eigenpair $(\mu,\bm{w})$ of the discretized generator $\bm{\hat{G}}$ has an approximate companion pair $(\mu-\mu^{(k)},\bm{w}\psi_{k})$, where $\bm{w}\psi_{k}$ is understood as pointwise multiplication and $\psi_{k}(t)=\omega^{kt}$, which only varies in time but not in space. For additional verification we can check the correlation of the companion eigenvectors

[TABLE]

as in [FK17, Section 7.5]. For instance, by looking at $\mu_{1}$ in Table 2 and noting that $\mu^{(\pm 1)}=0.01015\pm 0.34887i$ ($h=18/108$), we find that $\mu_{7},\mu_{8}\approx\mu_{1}-\mu^{(\pm 1)}$ are candidates for companion eigenvalues for $\mu_{1}$. The small difference in the numerical values of $\mu_{1}-\mu_{7,8}$ and the shift $\mu^{(\pm 1)}$, around $2.3\cdot 10^{-3}$ in magnitude, is due to the first eigenvector not being constant in time, i.e., merely $w_{t}\approx w_{t+h}$. Nonetheless, the complex eigenvalues $\mu_{7}$ and $\mu_{8}$ are companion eigenvalues to $\mu_{1}$. This is supported by the correlation for the corresponding eigenvectors

[TABLE]

while the correlation with other eigenfunctions $n\in\{2,\ldots,10\}\backslash\{7,8\}$ satisfies $|c_{n}^{1}(\pm 1)|\leq 0.00323$. The construction above (20) implies that every singular value is real. Our numerical calculations strongly suggest that within the first $50$ eigenvalues every complex eigenvalue is a companion to a real eigenvalue with smaller magnitude. The correlations using (36) yield results similar to those stated in the special case above.

5.4 Vortex isolation by sparse eigenbasis approximation

The space-time signatures of six coherent vortices in the Bickley flow are captured in the leading six singular vectors shown in Figure 3 (note only the initial time slice is displayed). In order to isolate these six vortices in space-time, we apply an orthogonal rotation and some sparsification to the six-dimensional subspace of $\mathbb{R}^{108\times(120\times 36)}$ spanned by the leading six (space-time) eigenvectors shown in Figure 3. The orthogonal rotation is chosen so as to construct an approximating basis of six sparse vectors. To find such a sparse approximating basis, we applied the SEBA (Sparse EigenBasis Approximation) algorithm (see [FRS18, Algorithm 3.1]). The six sparse basis vectors $\bm{\varphi}_{k}$, $k=1,\ldots,6$, produced by this algorithm are shown in Figure 4, and each of these vectors strongly isolates a single vortex. We emphasise that we input the full space-time vectors to the SEBA algorithm, but in Figure 4 display only the initial time slice.

In Figure 5 we seed particles inside the calculated vortical features (in the super-level set $\{\bm{\varphi}_{k}(0,\cdot)>0.4\}$ if $\bm{\varphi}_{k}$ is scaled to have maximum-norm 1) and evolve them forward in time to visualize the coherence. In addition to the deterministic evolution we also visualize a stochastic evolution using $\varepsilon=0.1$ as in our calculations above. Both simulations use a fourth-order Runge–Kutta (-Maruyama) scheme with step size $\frac{9}{4\cdot 108}=\frac{1}{48}$. Figures 5(b) and (c) demonstrate the coherence of the single vortex in Figure 5(a).

We wish to apply Proposition 9 to further demonstrate that the positive parts of the $\bm{\varphi}_{k}(0,\cdot)$ represent coherent sets. Because the $\bm{\varphi}_{k}(0,\cdot)$, $k=1,\ldots,6$ do not exactly span the leading six-dimensional eigenspace of $\bm{\hat{G}}$, Proposition 9 does not directly apply. Nevertheless, using the linear combinations of eigenfunctions $\bm{f}_{i}$ that result in the SEBA-features $\bm{\varphi}_{k}(0,\cdot)$ we found that the hypotheses of Proposition 9 were satisfied with the following modification: wherever the contributions (34) were negative, the corresponding $\alpha_{i}$ in (34) were set to zero. Any $\alpha_{i}$ that needed to be treated in this way was very close to zero. We do not depict these slightly modified linear combinations as they are still very close to the SEBA-features in Figure 4.

6 Optimization

Having developed an efficient means of computing singular vectors of $\mathcal{P}_{0,\tau}$ as eigenfunctions of the augmented generator $\bm{\hat{G}}$, we turn our attention to manipulating these eigenfunctions. These manipulations will be used to control the mixing properties of aperiodic flows. Theorem 8 provides a construction of a family of coherent sets $\{A_{t}^{\pm}\}_{t\in[0,\tau]}$ from eigenfunctions of $\bm{\hat{G}}$, with a coherence guarantee controlled by the corresponding eigenvalues. Our goal now is to either enhance or diminish the coherence of a family of sets related to an eigenvalue $\mu_{k}$ by small time-dependent perturbations $u(t,x)$ of the velocity field $v(t,x)$. This will be achieved by optimally manipulating $v(t,x)$ to increase or decrease the real part of the second eigenvalue $\mu_{2}$ of $\bm{\hat{G}}$.

For $m\geq 2$ let $H^{m}((0,\tau)\times X,\mathbb{R}^{d})$ denote the Sobolev space of vector fields on $(0,\tau)\times X$ whose weak derivatives of order up to $m$ are square integrable888In order to use results from [KLP19] we will later make an additional assumption on $m$.. As we have previously assumed that $\mathrm{div}_{x}v=0$ to simplify our presentation, for consistency we consider the subspace $\mathcal{D}_{0}$ of $H^{m}((0,\tau)\times X,\mathbb{R}^{d})$ consisting of spatially divergence-free vector fields that satisfy homogeneous Neumann boundary conditions: $\frac{\partial u}{\partial n}=0$ on $\partial X$. We consider small perturbations $u$ lying in a bounded, closed and strictly convex subset $\mathcal{C}\subset\mathcal{D}_{0}$.

We adopt the approach of [FS17], who select a perturbation $u$ so as to maximise the derivative of the real part of $\mu_{2}$ (or a group of leading eigenvalues) with respect to the perturbation. The perturbation was made to the Ulam-discretized generator of the vector field, and optimized using linear programming; the perturbed velocity field could then be inferred from the optimized generator. In [FS17] the velocity field was assumed to be time-periodic, however, the same optimization approach could be applied to aperiodic velocity fields using the time-reflected velocity field introduced in this work; one would need only to additionally impose the relevant time-reflection constraints in the optimization problem for the Ulam-discretized generator. In the present work, we consider aperiodic velocity fields and in contrast to [FS17] we perturb the velocity field directly and solve the resulting optimization problem by Lagrange multipliers. This potentially allows for greater flexibility in the discretization scheme.

To theoretically justify our approach in the infinite-dimensional setting we need some results from [KLP19] regarding the regularity of the spectrum of $\mathcal{P}_{0,\tau}$ with respect to perturbations of the velocity field. Transferring these results to the regularity of the spectrum of $\bm{\hat{G}}$ with respect to velocity field perturbations, we derive a first variation of $\mu_{k}$ with respect to $u$ and detail the steps below in section 6.1. In section 6.2 we specify the constraints and we then proceed with discussing necessary and sufficient (section 6.3) conditions for the optimization. Finally we summarize the result of the construction of our optimization.

6.1 Objective functional and its smoothness

We choose an eigenvalue $\mu_{k}$ of $\bm{\hat{G}}$ corresponding to a feature we want to enhance or diminish. We want to alter the real part of $\mu_{k}$ with a perturbation $u$ of the original velocity field $v$ as much as possible within our constraints. This is because by Theorem 8 the real part of $\mu_{k}$ is a measure for the coherence of the family of features highlighted by the corresponding eigenvector. Our chosen objective functional should be a good measure of the change of $\mu_{k}$ with respect to $u$. The response of an eigenvalue or a singular value with respect to a perturbation in this infinite-dimensional setting is in general complicated. Therefore we approximate it locally via linearization; that is by computing a first variation or first-order Taylor expansion. In what follows, we assume that $\mu_{k}$ is real; the obvious modifications can be made if $\mu_{k}$ is complex by considering the real parts.

Our domain $X$, the drift $v$, the perturbation $u$, for $m\geq 1+\frac{d}{2}$ such that the space $H^{m}((0,\tau)\times X;\mathbb{R}^{d})$ is conitnuously imbedded in $C^{(1,1)}([0,\tau]\times\overline{X};\mathbb{R}^{d})$ [AF03, Thm. 4.12], and the noise $\varepsilon I_{d\times d}$ are smooth enough to apply the results of [KLP19] to $\mathcal{P}_{0,\tau}$. Assuming that the singular value $\sigma_{k}$ is simple and isolated, [KLP19, Theorem 5.1] and the paragraph following it guarantee Fréchet differentiability of $\sigma_{k}$ and the corresponding singular function with respect to $u$. Using (20) and the expression $\hat{\mathcal{P}}_{0,2\tau}=\mathcal{P}^{\ast}_{0,\tau}\mathcal{P}_{0,\tau}$ we can relate the singular values and functions of $\mathcal{P}_{0,\tau}$ and the eigenvalues and eigenfunctions of $\smash{\hat{\mathcal{P}}_{0,2\tau}}$. In particular, the eigenvalues of $\smash{\hat{\mathcal{P}}_{0,2\tau}}$ are Fréchet differentiable with respect to $\hat{u}$. The spectral mapping property999See [CL99, Chapters 2,3 and 6] for analogous results in the context of evolution semigroups, [FK17, Lemma 22] for periodically forced systems, and [Paz83, EN00] for spectral mapping results for one-parameter semigroups. of Proposition 5 asserts for corresponding eigenvalues $\mu(\bm{\hat{G}})$ of $\bm{\hat{G}}$, eigenvalues $\lambda(\hat{\mathcal{P}}_{0,2\tau})$ of $\hat{\mathcal{P}}_{0,2\tau}$ and singular values $\sigma(\mathcal{P}_{0,\tau})$ of $\mathcal{P}_{0,\tau}$ that

[TABLE]

which extends the differentiability to the spectrum of $\bm{\hat{G}}$ and, in particular, $\mu_{k}$, hence $\mu_{k}=\tfrac{1}{2\tau}\ln\lambda_{k}=\tfrac{1}{\tau}\ln\sigma_{k}$.

Having establishing the Fréchet differentiability of $\mu_{k}$ with respect to $u$, we now calculate the first variation of $\mu_{k}$ with respect to $\hat{u}$; in other words, we compute the Gâteaux derivative of $\mu_{k}$ at $\hat{v}$ in the direction induced by $u$, which exists and coincides with the Fréchet derivative of $\mu_{k}$ at $\hat{v}$ applied to the direction induced by ${u}$. Consider $u$ as above and some small $\delta>0$. We insert the reflected perturbed velocity field $\hat{v}+\delta\hat{u}$ into (26):

[TABLE]

where $\bm{\hat{v}}=(1,\hat{v})$, $\bm{\hat{u}}=(0,\hat{u})$, $\hat{E}(t)(\hat{f})=-\hat{f}\text{div}_{x}(\hat{u})-\hat{u}\nabla_{x}\hat{f}=-\hat{u}\nabla_{x}\hat{f}$, using $\text{div}_{x}(\hat{u})=0$ and the perturbation generator $\hat{\bm{E}}=-\text{div}_{(\theta,x)}\bm{\hat{u}}-\bm{\hat{u}}\nabla_{(\theta,x)}$. Let $\hat{\bm{g}}^{k}(\delta)$ and $\hat{\bm{f}}^{k}(\delta)$ be the left and right eigenfunctions, respectively, of $\bm{\hat{G}}+\delta\bm{\hat{E}}$ for the eigenvalue $\mu_{k}(\delta)$; that is,

[TABLE]

normalising so that $\langle\hat{\bm{f}}^{k}(\delta),\hat{\bm{f}}^{k}(\delta)\rangle_{L^{2}}=\langle\hat{\bm{g}}^{k}(\delta),\hat{\bm{f}}^{k}(\delta)\rangle_{L^{2}}=1$. For $\delta=0$, we use the shorthand $\hat{\bm{f}}^{k}=\hat{\bm{f}}^{k}(0)$ and $\hat{\bm{g}}^{k}=\hat{\bm{g}}^{k}(0)$. To estimate the effect of the perturbation $u$ on $\mu_{k}$ we linearise $\mu_{k}(\delta)$ at $\delta=0$. We have

[TABLE]

using the eigenproperties of $\hat{\bm{f}}^{k}$ and $\hat{\bm{g}}^{k}$, and the normalisations above; see also [FS17, Section 4.3]. Now,

[TABLE]

where $c:\mathcal{C}\to\mathbb{R}$ is a linear function of $u$. If $\mu_{k}$ is complex, then one considers the real part of the functional $c$.

Lemma 11**.**

The objective functional $c\,:\,H^{m}((0,\tau)\times X)\rightarrow\mathbb{R}$, with $m\geq 1$, is continuous, Fréchet differentiable and the Fréchet derivative is Lipschitz continous.

Proof.

Using (38), the following estimate shows that $c$ is continuous.

[TABLE]

Here $\|\cdot\|_{\infty}$ denotes the canonical $L^{\infty}((0,2\tau)\times X)$ norm. The Fréchet differentiability of $c$ is straightforward because $c$ is linear. ∎

We will prove further relevant properties of $c$ in section 6.3.

6.2 Constraints

As mentioned above we consider perturbations $u\in\mathcal{C}$, a bounded, closed and strictly convex subset of $\mathcal{D}_{0}\subset H^{m}((0,\tau)\times X;\mathbb{R}^{d})$. For our objective functional $c$ to be a valid estimate of the change in $\mu_{k}$ due to the perturbation $u$, we restrict $u$ to be small using $R>0$. We consider a ball or ellipsoid in the form of the following energy constraint. For multi-indices $\alpha$ and a weight vector $\omega=(\omega_{\alpha})_{|\alpha|\leq m}$ with $0<\omega_{\alpha}\in\mathbb{R}$ for all $|\alpha|\leq m$ we require

[TABLE]

where $D^{\alpha}=(\partial_{1}^{\alpha_{1}}\cdots\partial_{d}^{\alpha_{d}})$. The functional $h$ is continuously Fréchet differentiable in $u$ from $H^{m}$ to $\mathbb{R}$ since $B_{\omega}$ is a bounded positive definite bilinear form. We consider $\mathcal{C}=U\cap\{h\leq 0\}$ where $U$ is a proper subspace of $H^{m}((0,\tau)\times X;\mathbb{R}^{d})$ and might only have a relative interior.101010For a rigorous definition of “relative interior” and related concepts see [Roc70, II.6]. For a (not entirely correct but) short intuition: The relative interior of a set $A$ can be thought of as the interior with respect to the trace topology of the smallest subspace that contains $A$. $U$ describes the set of admissible perturbations to the original vector field $v$. By construction in (39) there are constants $\gamma>0$ ($B_{\omega}$ is bounded) and $\beta>0$ ($B_{\omega}$ is positive definite) such that

[TABLE]

Now, Theorem 22 ensures the existence of a linear, bounded, injective and self-adjoint operator $J_{B}(\omega)$ with

[TABLE]

We will use $J_{B}(\omega)$ to derive an explicit formula for the optimal solution.

6.3 Optimality conditions

We have an optimization problem with a continuous linear objective $c$ and a closed, bounded, strictly convex feasible set $\mathcal{C}$, defined by a single constraint $\{h\leq 0\}$ on a subspace $U$. Lagrange multipliers provide a convenient and explicit solution to this problem, e.g. [Lue97], however, first we establish a general existence and uniqueness result.

6.3.1 Unique optimum

Lemma 12**.**

Let $\mathcal{C}$ be a closed, bounded, and strictly convex subset of $H^{m}$ containing the zero element in its (relative) interior and $c\,:\,\mathcal{C}\rightarrow\mathbb{R}$ be a bounded linear functional that does not uniformly vanish on $\mathcal{C}$. Then the optimization problem $\min_{u\in\mathcal{C}}c(u)$ has a unique solution $u^{\ast}\in\mathcal{C}$.

Proof.

Continuity of $c$ and boundedness of $\mathcal{C}$ imply $\inf_{u\in\mathcal{C}}c(u)=\alpha>-\infty$. Let $u_{k}\in\mathcal{C}$ be such that $\lim_{k\rightarrow\infty}c(u_{k})=\alpha$. This sequence is bounded, and so there is a weakly convergent subsequence $u_{n_{k}}\rightharpoonup u^{\ast}$. The set $\mathcal{C}$ is closed and convex and therefore also weakly closed, which implies $u^{\ast}\in\mathcal{C}$. By the definition of weak convergence $c(u^{\ast})=\lim_{k\rightarrow\infty}c(u_{n_{k}})=\alpha$ follows. Therefore $u^{\ast}$ is a solution for our optimization problem.

In order to demonstrate uniqueness assume we have two solutions $u_{1}\neq u_{2}$ with $c(u_{1})=c(u_{2})=\alpha$. Strict convexity of $\mathcal{C}$ implies that $u_{3}:=u_{1}/2+u_{2}/2\in\text{int}(\mathcal{C})$, where the relative interior is meant. Linearity of $c$ implies $c(u_{3})=\alpha$. Let $r>0$ be such that an open ball of radius $r$ centred at $u_{3}$ is contained in $\text{int}(\mathcal{C})$.

Because $c$ does not vanish on $\mathcal{C}$ and the zero vector is in the relative interior of $\mathcal{C}$, there exists a $v\in\mathcal{C}$ such that $c(v)<0$. By linearity of $c$ we have $c(u_{3}+(r/2)v)<\alpha$, contradicting optimality of $u_{3}$ and establishing uniqueness of the optimum. ∎

6.3.2 Necessary conditions

The following property will be used for the necessary conditions in the theory of Lagrange multipliers in Lemma 14.

Lemma 13**.**

The unique optimal solution $u^{\ast}$ is a regular point for $h:U\to\mathbb{R}$.

In order to distinguish between the point at which the derivative is taken and the input of the resulting mapping, we will use square brackets for the reference point of the derivative and round brackets for the input of the resulting mapping.111111So $h^{\prime}[u](v)$ denotes the derivative of $h$ at $u$, which again is a linear mapping, applied to $v$.

Proof.

Following the definition on [Lue97, p. 240], the point $u^{\ast}$ is a regular point for $U\to\mathbb{R}$ if the derivative of $h$ at $u^{\ast}$, denoted by $h^{\prime}[u^{\ast}]:U\rightarrow\mathbb{R}$ is surjective. Since $h(u)=B_{\omega}(u,u)-R^{2}$, and the functional $h^{\prime}[u]$ acts as $v\mapsto h^{\prime}[u](v)=\langle J_{B}(\omega)v,u\rangle_{H^{m}}+\langle J_{B}(\omega)u,v\rangle_{H^{m}}$, $h^{\prime}[u]$ is obviously surjective onto $\mathbb{R}$ for all $u\neq 0$. ∎

The uniqueness argument in the proof of Lemma 12 shows that we may replace our constraint $h(u)\leq 0$ with $h(u)=0$. We now use [Lue97, Theorem 1, p. 243], stated below.

Lemma 14**.**

If the continuously Fréchet differentiable functional $c$ has a local extremum under the constraint $h(u)=0$ at the regular point $u^{\ast}$, then there exists an element $z\in\mathbb{R}$ such that the Lagrangian functional $L(u)=c(u)+zh(u)$ is stationary at $u^{\ast}$; that is, $c^{\prime}[u^{\ast}]+zh^{\prime}[u^{\ast}]=0$.

We thus obtain the two necessary conditions:

[TABLE]

6.3.3 Sufficient conditions

We now prove that the necessary conditions (40) and (41) are in fact also sufficient. Because our objective is linear and our constraint is of inner product form, we take a direct approach to developing sufficient conditions, avoiding more complicated general theory.

Proposition 15**.**

Let $\mathcal{C}=U\cap\{h\leq 0\}$. There are exactly two elements of $\mathcal{C}$ that satisfy (40) and (41). One is the unique minimizer (with $z>0$) and the other is the unique maximizer (with $z<0$).

Proof.

Lemma 12 guarantees the existence of at least two extrema (one minimum and one maximum), and therefore at least two distinct elements $u,w\in\mathcal{C}$ satisfying (40) and (41). We show that these are the only such elements. There exist $z_{u},z_{w}\in\mathbb{R}$ such that

[TABLE]

Subtracting these two equations and using $c^{\prime}[u](\cdot)=c^{\prime}[w](\cdot)=c(\cdot)$, we obtain the functional equation $h^{\prime}[w]=(z_{u}/z_{w})h^{\prime}[u]$. Thus, the linear functional $h^{\prime}[w]$ is a scalar multiple of the linear functional $h^{\prime}[u]$. Since $h^{\prime}[u](\cdot)=2\langle\cdot,J_{B}(\omega)u\rangle_{H^{m}}$, by the Riesz representation theorem, we have $w=(z_{u}/z_{w})u$. However, the necessary condition $h(u)=h(w)=0$ implies that $h^{\prime}[u](u)=2\langle u,u\rangle_{m,\omega}=2\|u\|_{m,\omega}^{2}=2R^{2}$ and similarly that $\|w\|_{m,\omega}^{2}=R^{2}$. Thus, either $z_{u}=z_{w}$ or $z_{u}=-z_{w}$. If $z_{u}=z_{w}$, then we have $u=w$, while if $z_{u}=-z_{w}$, then $u=-w$. Thus, the only possibility for distinct $u$ and $w$ is that $u=-w$, and therefore that there are at most two functions satisfying the necessary conditions. Finally, without loss, assume that $u$ is a minimum. Since $c^{\prime}[u](u)+z_{u}h^{\prime}[u](u)=c(u)+z_{u}R^{2}=0$ and $c(u)<0$ if $u$ is a minimum, we must have $z_{u}>0$. Therefore $z_{w}<0$, implying that $c(w)>0$ and that $w$ is a maximum. ∎

Using the injective operator $J_{B}(\omega)$ we can solve the necessary and sufficient conditions (40) and (41) for the optimal solution $u^{\ast}$, leading to (44) below. First we can transform (40) into an equation in $H^{m}$ using the Riesz representation theorem. We now know that $u^{\ast}$ exists and fulfills the following equation.

[TABLE]

Here $c_{R}\in H^{m}$ is the Riesz representation of the functional $c$ on $H^{m}$. Thus $c_{R}$ is in the range of $J_{B}(\omega)$ and we can apply $J_{B}(\omega)^{-1}$ to $c_{R}$. Now we can solve (42) for $J_{B}(\omega)u^{\ast}$ and for $u^{\ast}$ because Proposition 15 guarantees $z\neq 0$, giving

[TABLE]

Inserting these expressions into (41) leads to

[TABLE]

Solving (43) for $z>0$ (minimizer) and using this $z$ leads to the following explicit expressions

[TABLE]

7 Optimization of $\mu_{k}$ numerically

In this section we apply the results derived in section 6 to some examples.

7.1 Discrete optimization problem

Before discretising the objective functional we want to construct $\mathcal{C}_{N}$, a finite dimensional version of the constraint set $\mathcal{C}$. Therefore we choose finitely many basis elements $\{\varphi_{\ell}\}_{\ell=1,\ldots,N}$ spanning the admissible subspace $U_{N}:=\text{span}\{\varphi_{\ell}\}_{\ell=1,\ldots,N}$ for our perturbations. We then intersect $U_{N}$ with $\mathcal{C}$ and represent elements by their coefficient vectors in $\mathbb{R}^{N}$ with respect to the chosen basis. Hence we define $\mathbb{R}^{N}\supset\mathcal{C}_{N}\simeq\mathcal{C}\cap U_{N}\subset H^{m}$.

Coefficient vectors will be denoted with a bar $\bar{\phantom{u}}$ in the following. We will omit the $\hat{\phantom{X}}$ in the following calculations, but note that it can be done analogously using augmented reflected objects. The energy neighborhood constraint (39) can be expressed as a quadratic constraint in the coefficient vector $\bar{u}\in\mathbb{R}^{N}$,

[TABLE]

This constraint describes a strictly convex set (ball or ellipsoid) in $\mathbb{R}^{N}$. Regarding the objective functional

[TABLE]

we have to account for the two possibly different bases for discretization: (i) the discretization of $u$, and (ii) the discretization of $\bm{\hat{G}}$, which can also involve a test and an ansatz basis. Let us denote the basis functions for the discretization of $\bm{f}^{k}$ by $\{\chi_{j}\}_{j=1,\ldots,m}$ and $\bm{g}^{k}$ by $\{\xi_{i}\}_{i=1,\ldots,n}$ for now. Then for $\bm{f}^{k}=\sum_{j=1}^{m}\bar{f}_{j}^{k}\chi_{j}$ and $\bm{g}^{k}=\sum_{i=1}^{n}\bar{g}^{k}_{i}\xi_{i}$ we have

[TABLE]

Using this last equation we can calculate a finite-dimensional representation of $\bm{E}$ acting from $\text{span}\{\xi_{j}\}_{j}$ to $\text{span}\{\chi_{i}\}_{i}$.

Due to linearity we may decompose in $\ell$ and separately compute $\bm{E}_{\ell}$ in a similar way to the numerical approximation of $\bm{G}$, outlined in Section 5.1. We use Ulam’s method to discretize the generators, taking $m=n$ and $(\xi_{j})_{j}=(\chi_{i})_{i}$ to be indicator functions of space-time boxes. The cost vector can be constructed by

[TABLE]

The discretized optimization problem then has a linear objective and a single quadratic constraint

[TABLE]

Since the energy constraint is induced by a scalar product, the matrix $B_{\omega}$ is invertible by typical arguments for Galerkin discretization (i.e., $B_{\omega}$ is symmetric positive definite). Thus the optimal solution can be obtained with Lagrange multipliers, identical to the analysis leading to (44), and is given by

[TABLE]

In practice we choose the “admissible energy” $R$ sufficiently small to ensure the approximate validity of our linearised objective functional. We can then iterate the optimization process as a gradient ascent/decent method to invest more cumulative energy in the perturbation. Each step consists of constructing $\bar{c}$ and solving the equations (45). The construction of $\bar{c}$ requires the calculation of $\bm{G}$, $\bm{g}^{k}$, $\bm{f}^{k}$, and $\bm{E}_{\ell}$ (or the respective augmented reflected objects); the latter are fixed through all optimization steps and do not need to be updated. We apply this procedure in the following examples.

Remark 16**.**

All of our finitely many perturbation ansatz functions $(\varphi_{\ell})_{\ell=1,\ldots,N}$, introduced in the next section 7.2, are $C^{\infty}$ thus $\|\cdot\|_{H^{m}}$ and $\|\cdot\|_{L^{2}}$ are equivalent on $\text{span}\{\varphi_{l}\}_{l}$. Therefore for numerical convenience, in the numerical examples we use $m=0$ to calculate $B_{\omega}$, although the functional $c$ is only strictly well defined for $m\geq 1$.

7.2 Perturbing fields

We construct a suitable basis for velocity field perturbation as follows. Our spatial domain will be a rectangle. In order to impose zero velocity normal to the boundary and divergence-freeness, we construct the spatial components of the basis vectors $\{\varphi_{\ell}\}_{\ell=1,\ldots,N}$ with a possible constant movement in $x$ or $y$ direction from smooth stream functions $\Psi_{kl}$ and then multiply these components with time-dependent scalar (amplitude) functions $\phi_{r}$. For a rectangular domain $[a_{x},b_{x}]\times[a_{y},b_{y}]$ we take the streamfunctions

[TABLE]

$k=1,\ldots,K$, $l=1,\ldots,L$, which are slightly modified Fourier modes that induce a velocity field $\psi_{kl}:=(-\frac{\partial\Psi_{kl}}{\partial y},\frac{\partial\Psi_{kl}}{\partial x})$ with $k$ horizontal gyres and $l$ vertical gyres that satisfy the homogeneous Neumann boundary conditions in space and are divergence free in space. These fields may travel in $x$ direction with speed $c_{x}$ or in $y$ direction with speed $c_{y}$.

We use $L^{2}((0,\tau)\times X)$-normalized versions of the functions

[TABLE]

where $\phi_{r}$ is a scalar (temporal) modulation of the amplitude of the spatial Fourier modes:

[TABLE]

We omit using $r=1$ for the $\sin$-modulation, since it would have both positive and negative values, meaning a sign change of the perturbing velocity field during the evolution. Such perturbing fields proved to be less efficient in early numerical experiments. Thus, the increasing time-linear modulation is assigned $r=-1$ to avoid confusion. In summary, time-modulation of the perturbing fields is described by $\phi_{r}(t)$, $r\in\{-1,0,2\}$, and we have $3KL=N$ basis functions $\{\varphi_{\ell}\}_{\ell=1,\ldots,N}$ in total.

7.3 Increasing coherence: forced double gyre flow

Extending the example from section 5.2, our goal in this section is to increase coherence of the left-right separation captured by the $2$nd eigenvalue of $\bm{\hat{G}}$ (equivalently, the left $2$nd singular vector of the transfer operator $\mathcal{P}_{0,\tau}$), Figure 2 (a).

The original velocity field has a total energy (space-time $L^{2}$ norm) of $\approx 1.6$. We will use an optimization budget of $R=0.05$per step to increase the $2$nd eigenvalue $\mu_{2}$ of $\bm{\hat{G}}$ which encodes the left-right separation of the domain. We iterate our optimization procedure $8$ times to invest a total energy of $0.4$, which is $25\%$ of the energy of the original velocity field. This seems like a moderate investment, but we note that while in general it is easy to destroy coherence by almost any perturbation, to increase it, the dynamics and the perturbation need to work together—making it harder to increase coherence than to decrease it. Following the formula (45) to optimize coherence (i.e, minimize mixing) we use $-\bar{u}$, where we recall that $\bar{u}$ denotes the coefficient vector representing the optimal perturbation. Of course the other eigenvalues of the generator also change as the velocity field is perturbed. In each iterative step we check a posteriori that the second eigenvalue did indeed increase, which we consider an indicator for the validity of our objective functional. Our perturbation library consists of the functions from (46) for $k=1,\ldots,5$, $l=1,2,3$ with $c_{x}=0=c_{y}$ and the time modulation $r=-1,0,2$, hence $N=45=5\cdot 3\cdot 3$. After the 8 iteration steps we arrive at the solution vector $\bar{u}_{(8)}\in\mathbb{R}^{N}$ and the effective change of the second eigenvalue singular value:

[TABLE]

Each step of the iteration increases the eigenvalue roughly by $0.003$. The result is visualized in Figure 6.

We seeded $200,000$ particles on the right side of the line $x=1$ and evolved them forward in time with the noisy flow, using $\varepsilon=0.1$ (Runge–Kutta–Maruyama with time step size $h=\frac{1}{100}$) for $\tau=4$ time units using the original and the optimized drift. For the original velocity field roughly $15\%$ end up on the left side of the domain, where as for the optimally increased left-right separation only about $10\%$ of the particles end up on the left side (results not shown). We repeated the same noisy evolution procedure for $\varepsilon=0.01$; then the $9\%$ of particles changing sides originally were reduced to $5\%$ for the coherence-improved velocity field. This is shown in Figure 6 (b)–(d), where the seeded particles are shown at initial time and colored according to whether they end up left or right of the line $x=1$ after this noisy evolution, at time $t=\tau$. Note that the time direction is unimportant and we could have colored at the initial time $t=0$ and evolved forward.

The optimization of the left-right coherence in the periodically forced double gyre has also been considered in  [FS17, Section 6.3]; see, in particular Figure 18 therein, where the regular regions of the flow were dramatically increased.

7.4 Decreasing coherence: forced double gyre flow

We now wish to diminish the coherence of various coherent features. Firstly, the left-right separation discussed in Section 7.3. Using the same optimization protocols as in Section 7.3, but switching the sign of the objective, we produce a velocity field that should increase mixing across the left-right separatrix. This is indeed indicated in Figure 6 (d). These results are consistent with the results of [FS17], where the “lobes” of stable and unstable manifold intersections greatly increased [FS17, Figures 11 and 14].

We now turn our attention to the two central vortices, which are encoded in the $5$th eigenmode (and after one iteration step of the optimization in the $6$th eigenmode). We use the same perturbation basis and energy criterion for the iteration as above ($8$ iterative optimization steps). After a first optimization step of the iteration described at the end of section 7.1 the gyre feature (initially $5$th eigenvalue) is pushed to the $6$th eigenvalue spot, i.e., there is an interchange position with the mentioned features in terms of ranking with respect to their coherence. Thus, we have to keep track of the ranking of the eigenvalues during our iterative optimization procedure, which we do here manually between each step. This could be done in an automated fashion similarly to section 5.3, by computing correlations between eigenvectors of successive iterates.

We obtain the following change in the eigen- and (corresponding) singular values:

[TABLE]

Next we seed particles in the vortices induced by the level sets of the $5$th eigenvector shown in Figure 2 and evolve them with and without perturbation. The results are visualized in Figure 7.

We note that increasing mixing in the double gyre flow has been considered in [FS17, Section 6.2]; in particular our Figure 7 could be compared with Figures 13, 16, and 17 therein.

7.5 Targeted manipulation of distinguished coherent features

In certain situations it may be of interest to manipulate the mixing of parts of phase space that do not arise as eigenmodes. For example, these parts of phase space may be individual (or combinations of) SEBA vectors as described in section 5.4, or they may be related to the phase space geometry. We describe a procedure to accomplish this.

We recall Proposition 5 and Theorem 8, where we showed

[TABLE]

and that the sign structure of $\bm{f}$ (recall $\bm{f}$ necessarily has zero mean) indicates a family of finite-time coherent sets. Further, the coherence ratio of the family is bounded by an expression involving $\mu$, indicating more coherence the closer $\mu$ is to [math]. By normalising $\bm{f}$ so that $\|\bm{f}\|=1$ we see smaller $\|\bm{\hat{G}}\bm{f}\|$ corresponds to more strongly coherent features encoded in $\bm{f}$.

Now let $\bm{\varphi}\in\mathcal{D}(\bm{\hat{G}})$ be a general normalised, zero-mean space-time feature; that is, $\bm{\varphi}$ is not an eigenfunction. We might think of $\bm{\varphi}$ being mean-removed SEBA vector or a mollified (such that it is in the domain $\mathcal{D}(\bm{\hat{G}})$ of the generator) version of $\bm{\varphi}=\mathds{1}_{\bm{C}}-|\bm{C}|\mathds{1}_{\bm{X}}$, the mean-centred121212We note that the removal of the mean from $\bm{\varphi}$ makes no difference for the optimization of the objective function in (48) below, since $\bm{\hat{G}}\mathds{1}_{\bm{X}}=0$, however we keep this for the intuitive connection with “eigenfeatures” and Theorem 8. indicator function of a possibly coherent family $\bm{C}\subset\bm{X}$ of sets in augmented-space representation, where $|\bm{C}|$ denotes the augmented-space Lebesgue measure of $\bm{C}$. As we would in general like $\bm{\varphi}$ to represent a finite-time coherent set, we should restrict our attention to features satisfying $\bm{\varphi}(t,\cdot)\approx\bm{\varphi}(2\tau-t,\cdot)$.

Analogously to the case of an eigenfunction $\bm{f}$, to quantify the coherence of a feature $\bm{\varphi}$ that is not necessarily an eigenfunction, we employ the heuristic of measuring $\|\bm{\hat{G}}\bm{\varphi}\|$. The rationale for this is as follows. If a family of sets encoded by the eigenvector $\bm{f}$ is completely coherent (in the absence of diffusion), then the temporal change (“movement”) of the sets at any time $\theta$, namely $\partial_{\theta}\bm{f}(\theta)$, would be identical to how the dynamics transports the mass located in the set, i.e., $\partial_{\theta}\bm{f}(\theta)=\hat{G}(\theta)\bm{f}(\theta,\cdot)$ by the Fokker–Planck equation (11). Thus, if the coherence of the feature $\bm{\varphi}$ is strong, one has $\partial_{\theta}\bm{\varphi}(\theta,\cdot)-\hat{G}(\theta)\bm{\varphi}(\theta,\cdot)\approx 0$ for all $\theta$, leading to $\|\bm{\hat{G}}\bm{\varphi}\|\approx 0$. Section 3 gives a geometric view on the very same situation: in (14) and (17), if the boundary of a time-dependent set moves with a velocity $b(t,x)$ that is approximately equal to the velocity field $v(t,x)$ driving the dynamics, then the outflow from this family of sets will be small—and this can analogously be quantified by the space-time flux (18).

Thus, to destroy a coherent feature encoded in $\bm{\varphi}$ we could maximize $\|\bm{\hat{G}}\bm{\varphi}\|^{2}$ with respect to the perturbing fields $u$. Again, as this is a nonlinear problem, we approach it by local optimization, and aim to maximize the objective function given by the local linear change,

[TABLE]

subject to constraints on the perturbation $u$. Conversely, if we wish to enhance a coherent feature $\bm{\varphi}$ we should minimize $\|\bm{\hat{G}}\bm{\varphi}\|^{2}$. If we would simultaneously like to destroy coherence of a feature $\bm{\varphi}_{1}$ and enhance the coherence of other features encoded in $\bm{\varphi}_{2}$, then we would maximize

[TABLE]

with weights $\alpha_{1},\alpha_{2}>0$.

7.6 Non-eigenfeature optimization: traveling wave

We consider a traveling wave example [Pie91, SW06, FLS10] given by

[TABLE]

on the domain $[0,8]\times(2\pi S^{1}\times[0,\pi])$. Here, two rotating gyres move in $x$ direction with speed $\nu=0.25$ and are superimposed with a constant drift $c_{\text{drift}}=1$. To account for this constant speed and periodicity in $x$-direction we choose $k=2,4,\ldots,20$, $l=1,\ldots,5$ and $c_{x}=\nu$ for our perturbation dictionary; see Section 7.2 and equation (46). We remark that Balasuriya [Bal15] has investigated a similar dynamical system and considered single “one at a time” (as opposed to general linear combinations of) perturbations drawn from a family similar to ours. In [Bal15], flux out of a small “gate” connecting a stable and unstable manifold is taken as a measure of mixing. In contrast, we measure mixing through the $L^{2}$-norm decay of an initial concentration field, and we find the unique perturbation in a convex subset of a 150-dimensional subspace that maximizes the change in mixing rate.

We use the resolution $80\times(80\times 40)$ and $\varepsilon=0.1$. We now take a feature that is not described by an eigenfunction and aim to increase its coherence. The feature we choose is the mean-centered version of the following time-constant and horizontally constant profile

[TABLE]

which is shown in the coloring of Figure 8 (a). We iteratively update the perturbation by solving problem (48) in each step. We iterate for 35 steps with an energy budget of $R=0.1$ ($\sim 1\%$ of the original energy of $v$) per step. Figure 8 (b) shows the final time slice of the original (deterministic) evolution of the particles in Figure 8 (a), while Figure 8 (c) shows the final time slice of the optimized evolution.

We next investigate which perturbing basis functions are favored by the optimization. First, we note that the basis function induced by $\varphi_{0,2,1}$ is equal to the original velocity field up to the constant horizontal drift $c_{\text{drift}}$. Thus, it is conceivable that this basis function is heavily used in order to partly cancel the original velocity field and slow the flow down.

The amplitudes of all of the streamfunctions in the optimized solution are shown in Figure 9, ordered according to amplitude and grouped according to the spatial mode. The streamfunctions of the spatial modes with the largest amplitudes in the optimal solution (first to fourth) are shown in Figure 10. Altogether, these modes combined to perturb the traveling double gyre towards a laminar horizontal flow (Figure 11 (c)).

We additionally extended the above computation to $100$ steps to investigate the asymptotic behavior under large perturbations. Around step $95$ the increase in coherence diminishes rapidly, as the value of the objective function $\|\bm{\hat{G}}(v+u)\bm{\varphi}\|^{2}$ approaches a plateau. The corresponding velocity field approaches a purely laminar flow (results not shown here), as the optimal perturbation cancels the rotational part of the original velocity field.

8 Conclusions

Froyland and Koltai [FK17] introduced a time-augmented construction to enable the efficient numerical construction of the infinitesimal generator of a periodically driven flow; [FK17] built on earlier work for steady flows [FJK13]. In the first part of this work, we extended these results to a finite-time flows with general aperiodic driving, by a novel “reflected” process. Proposition 3 provided a formula for the accumulated outflow from a general time-dependent family of sets. Proposition 5 derived a spectral mapping theorem for our reflected process, connecting the spectrum of the time-augmented generator with the spectrum of the corresponding reflected evolution operator. Using the sign structure of the time-augmented eigenfunctions of the augmented generator we built a family of coherent sets and Theorem 8 lower bounded the coherence ratio of this family in terms of the corresponding eigenvalue.

In the second part of this work, we built on the optimization techniques of [FS17] to optimally enhance or destroy coherent features encoded in eigenfunctions. We directly manipulated the underlying drift field, subject to energy constraints, and proved that the “small perturbation” problem has a unique optimum (Proposition 15). Using Lagrange multipliers, we derived an explicit formula for the infinitesimal drift field perturbation; this would be very difficult to achieve without using a generator framework. In Section 7 we implemented a multi-step gradient descent method, utilising the efficient time-augmented generator framework, and the explicit Lagrange multiplier solution, to optimize over relatively large energy budgets.

An advantage of our optimization approach is that the basis of perturbing velocity fields is fixed; thus their generators need only be computed once. On the other hand, the generators are very large, sparse matrices, of which we need to compute the spectrum with the smallest real part. Hybrid spectral discretization techniques as in [FK17] or multilevel solvers can be a remedy to this. Finally, the approach could be extended to the non-zero divergence case, and to open systems by considering homogeneous Dirichlet boundary conditions in the Fokker–Planck equation (2).

Acknowledgments

We thank Andreas Denner for discussions leading to the inception of this work, and Nicolas Perkowski for suggestions concerning Theorem 8. GF thanks the DFG Priority Programme 1881 “Turbulent Superstructures” for generous travel support, and the Department of Mathematics at the Free University Berlin for hospitality. GF is partially supported by an ARC Discovery Project. DP180101223. MS is supported by the DFG Priority Programme 1881 “Turbulent Superstructures” and by the DFG CRC 1294. PK is partially supported by the DFG CRC 1114, project A01.

Appendix A Non-autonomous Cauchy Problems

In this appendix we summarize important results from the theory of non-autonomous abstract Cauchy problems (NACPs) of the form:

[TABLE]

over some finite time interval $[s,T]$ in $L^{p}(X)$, for $1<p<\infty$. The linear operators $G(t)$ are so-called infinitesimal generators, defined below. We state the relevant assumptions for this section, which all have been made and used through out the paper:

• •

$X\subset\mathbb{R}^{d}$ is an open bounded domain with piece-wise $C^{4}$ boundary.

• •

$v\in C^{(1,1)}([0,\tau]\times\overline{X};\mathbb{R}^{d})$.

To see that our Fokker–Planck equations (2) and (11) fit in the setting of non-autonomous abstract Cauchy problems we will make use of the following notation:

[TABLE]

Equivalently, $C((0,\tau);L^{p}(X))$ will denote the space of functions mapping $(0,\tau)$ continuously to $L^{p}(X)$. This means that instead of considering the PDE pointwise in every $(t,x)$ the NACP considers $f$ as a mapping from $[0,\tau]$ into a the Banach space $L^{p}(X)$ this leads to a differential equation (50) in $L^{p}(X)$.

Lemma 18 below establishes that our family of operators $(G(t))_{t\in[0,\tau]}$ considered over $L^{p}(X)$, $1<p<\infty$, defined by

[TABLE]

for $\varepsilon>0$, satisfies the assumptions made in the references [Tan96] and [Lun95] using the domain \smash{\mathcal{D}(G(t))=\mathcal{D}_{p}:=\big{\{}f\in W^{2,p}(X)\,\big{|}\,\frac{\partial f}{\partial n}=0\text{ on }\partial X\big{\}}} for all $t\in[0,\tau]$. Theorem 19 is concerned with the existence of unique solutions to (50) and some regularity properties of the solution. Theorem 20 states further regularity results and Corollary 21 is concerned with the regularity of eigenfunctions of $\bm{\hat{G}}$. The latter is needed in the proof of Theorem 8 in Appendix B.

Remark 17**.**

The results in **[Tan96]** are formulated for the time interval $[0,T]$ but of course also hold for arbitrary time intervals $[a,b]$. 2. 2.

We formulate all of the following on $[0,\tau]$. The reflected problem on $[\tau,2\tau]$ can be done analogously. 3. 3.

The case $p=1$ can be treated as well with the theory of Tanabe **[Tan96]**, cf. **[FK17, Appendix]**, but is more complex and will be omitted here.

Lemma 18**.**

We consider the family $(G(t))_{t\in[0,\tau]}$ of unbounded linear operators over $L^{p}(X)$ with the domain $\mathcal{D}(G(t))=\mathcal{D}_{p}$, defined by (51). The following conditions are fulfilled:

1. (1)

The spatial domain $X$ is a bounded open set of $\mathbb{R}^{d}$ of class $C^{4}$ **[Tan96, p.279]** (locally) (and globally uniformly regular of class $C^{2}$ **[Tan96, sec.5.2]**. 2. (2)

The coefficients of the differential operator $G(\cdot)$ and the boundary operator $B(t,x)(\xi):=\langle n(x),\xi\rangle$, with $n(x)$ being the unit outer normal at $x\in\partial X$ and $\xi\in\mathbb{R}^{d}$, are Hölder-continuous in $t$ **[Tan96, sec. 6.13]**. 3. (3)

The coefficients of $G(\cdot)$ are bounded and uniformly continuous in $x$ on $\overline{X}$ **[Tan96, sec. 6.13]**. 4. (4)

The differential operator $G(\cdot)$ and the boundary operator $B(t,x)(\xi):=\langle n(x),\xi\rangle$, satisfy the complementing conditions **[Tan96, p.131]** for all $t\in[0,\tau]$. 5. (5)

The family $(G(t))_{t\in[0,\tau]}$ is uniformly strongly elliptic **[Tan96, Def.5.4]** and $G(t)$ satisfies the root condition **[Tan96, p.130]** for all $t\in[0,\tau]$. 6. (6)

The conditions **[Tan96, (P1), p.221]**, **[Tan96, (P2), p.222]** and **[Tan96, (P4), p.256]** are fulfilled.

Before we continue with the proof, let us remark, that most of the proof has already been done in [Tan96] and some conditions ((3),(4) and (5)) become trivial in our case because the highest order part (called principal part in [Tan96]) of $G(t)$ is the time-independent Laplace operator $\Delta$.

Proof.

One important aspect of this lemma is to establish that our choice $\mathcal{D}_{p}\subset W^{2,p}(X)$ for the domain $\mathcal{D}(G(t))$ is appropriate.

Most of the statements above are proven at some point in [Tan96].

1. (1)

The regularity required is compatible with our general assumption on the regularity of $X$. The condition for uniformly $C^{2}$ in [Tan96] is irrelevant for our context as it is only important for unbounded domains. 2. (2)

$B(t,x)$ is constant in $t$ and therefore Hölder continuous. Our assumptions guarantee that $v$ and $\partial_{x}v$ are Hölder continuous in $t$. 3. (3)

This is satisfied by the assumptions we make on $v$. 4. (4)

This can be verified by straightforward calculations. The important thing to note is that the principal part (called $L^{0}$ in [Tan96]) of $G(t)$ is time-independently $\frac{\varepsilon^{2}}{2}\Delta=:L^{0}$, which implies the roots $r_{1}=i,r_{2}=-i$ of $L^{0}(x,\xi+rn(x))$ for $\xi$ perpendicular to $n(x)$ with unit norm, and that $B(t,x)(\xi+rn(x))=r$ for $\xi$ being perpendicular to $n(x)$. 5. (5)

It is well known that the operator family $(G(t))_{t\in[0,\tau]}$, defined by (51), is uniformly elliptic in our setting. Now [Tan96, Thm. 5.4] states that every strongly elliptic operator satisfies the Root Condition. 6. (6)

[Tan96, sec. 6.13] shows that the assumptions (P1), (P2) and (P4) are satisfied for advection-diffusion type operators as in our setting.

∎

Theorem 19**.**

The NACP (50) with the operators defined by (51) and the assumptions stated in this paper has a unique solution

[TABLE]

given by $f(t)=\mathcal{P}_{s,t}f(s)$, and that $f(t)\in\mathcal{D}(G(t))$ for $t\in(0,\tau]$. Further, the family $(\mathcal{P}_{s,t})_{t\geq s}$ is a family of linear, bounded and even compact operators (for $t>s$) on $L^{p}(X)$.

Proof.

Lemma 18 ensures that we can apply the results from [Tan96] that we will use in the following. The results from [Tan96, sec. 6.13] that use more abstract results from Acquistapace and Terreni ([AT86] and [AT87]) give the existence of a unique solution $f$ to (50) [Tan96, Thm. 6.6] and a corresponding two parameter family of solution operators (called fundamental solution in [Tan96]) $(\mathcal{P}_{s,t})$.

• •

The function $f$ is a classical solution [Tan96, Def.6.1], i.e.

[TABLE]

and $f(t)\in\mathcal{D}(G(t))=\mathcal{D}_{p}$ for $t>0$ follows directly from [Tan96, Thm. 6.6].

• •

For all $t>s$ we the solution can be expressed as $f(t)=\mathcal{P}_{s,t}f(s)$ with the two parameter solution family $(\mathcal{P}_{s,t})_{t\geq s}$. This follows from [Tan96, Thm. 6.5].

• •

For every $t>s$ the operator $\mathcal{P}_{s,t}$ is compact on $L^{p}(X)$. The result [Tan96, Thm. 6.6], i.e., $\mathcal{P}_{s,t}f(s)=f(t)\in\mathcal{D}_{p}$ for $f(s)\in L^{p}(X)$), gives that for $t>s$ $\mathcal{P}_{s,t}\,:\,L^{p}(X)\rightarrow\mathcal{D}_{p}\subset W^{2,p}(X)$ is a bounded linear operator from $L^{p}(X)$ to $W^{2,p}(X)$. Now the Rellich–Kondrachov embedding theorem [AF03, Thm. 6.3] states that $W^{2,p}(X)$ is compactly embedded in $L^{p}(X)$. Thus $\mathcal{P}_{s,t}$ is a compact operator from $L^{p}(X)$ to $L^{p}(X)$, because it maps bounded sets in $L^{p}$ to bounded sets in $W^{2,p}(X)$ which are relatively compact sets in $L^{p}(X)$.

∎

Due to the common domain of all $G(t)$ we get a better regularity for the the solution if the initial condition is more regular. More precisely, we use the result [Lun95, Cor. 6.1.9 (iv)] to prove the following.

Theorem 20**.**

1. (i)

It holds that $t\mapsto f(t)=\mathcal{P}_{0,t}f_{0}\in C([0,\tau];\mathcal{D}_{p}))\cap C^{1}([0,\tau];X)$ if and only if $f_{0}\in\mathcal{D}_{p}$ and $G(0)f_{0}\in L^{p}(X)$. 2. (ii)

The regularity $f\in C([0,\tau];W^{m,p}(X))$ implies Hölder continuity131313$C^{\alpha}(\overline{X})$ denotes the space of conitinous functions that are Hölder continuous on $\overline{X}$ with exponent $\alpha\in(0,1)$.* in space for every time slice, i.e., $f(t,\cdot)\in C^{\alpha}(\overline{X})$ for all $p$ such that $0<\alpha\leq m-\frac{d}{p}$, with uniformly bounded Hölder-norm in $t$.* 3. (iii)

Further, $f\in C([0,\tau];W^{2,p}(X))$ implies $f\in C([0,\tau]\times X)$.

Proof.

1. (i)

This result is basically [Lun95, Cor. 6.1.9 (iv)] with $\mathcal{D}_{p}=\mathcal{D}(G(t))$ being the common domain of the family of unbounded operators $(G(t))_{t\in[0,\tau]}$ in the Banach space $L^{p}(X)$ ($\mathcal{D}_{p}$ is dense in $L^{p}(X)$), and the resulting two parameter evolution family is $(\mathcal{P}_{s,t})_{t\geq s}$. 2. (ii)

For every $t\in[0,\tau]$ we have $f(t,\cdot)\in W^{m,p}(X)$. From the appropriate Sobolev embedding theorem [AF03, Thm. 4.12] (also cf. Morrey) follows that for domains as smooth as our $X\subset\mathbb{R}^{n}$ the space $W^{m,p}(X)$ is continuously embedded in $C^{\alpha}(\overline{X})$ for $0<\alpha\leq m-\frac{d}{p}$, i.e., there exists a $K>0$ such that

[TABLE]

for all $g\in W^{2,p}(X)$. This with the regularity $f\in C([0,\tau];W^{m,p}(X))$ further implies that the constant for the continuous embedding for $(f(t))_{t\in[0,\tau]}$ into $C^{\alpha}$ can be taken as independent of $t$. 3. (iii)

We consider a function $f\in C([0,\tau];W^{2,p}(X))$. Now (ii) gives $f(t,\cdot)\in C^{\alpha}(X))$ for every $t$ with $0<\alpha\leq 2-\frac{d}{p}$. Together with (i) this implies

[TABLE]

This immediately implies $f\in C([0,\tau];C(\overline{X}))$ which is equivalent to $f\in C([0,\tau]\times\overline{X})$. Let us briefly show the direction of the last statement that we need: For $f\in C([0,\tau];C(\overline{X}))$ we have

[TABLE]

∎

Corollary 21**.**

Let $\bm{f}$ be an eigenfunction of $\bm{\hat{G}}$ considered on $L^{p}$, $p>\frac{d}{2}$, then $\bm{f}\in C([0,2\tau]\times X)$.

Proof.

For an eigenfunction $\bm{f}$ follows

[TABLE]

from theorem 5. Further theorem 19 implies $e^{-\mu\tau}\hat{\mathcal{P}}_{0,2\tau}\bm{f}(0,\cdot)=\bm{f}(0,\cdot)\in\mathcal{D}(\hat{G}(2\tau))$. Now theorem 20 (i) and (iii) with $\mathcal{D}(\hat{G}(2\tau))=\mathcal{D}_{p}$ gives the claim. ∎

Appendix B Proof of Theorem 8 and Proposition 9

Proof of Theorem 8.

The proof strongly follows the lines of those of [FK17, Theorems 16 and 19], which, in turn, borrows ideas from [FS10, FS13]. It consists of two main steps. First, we consider the events

[TABLE]

for a dense sequence of times, and show that $\mathcal{E}_{n}\downarrow\mathcal{E}:=\bigcap_{r\in[0,\tau]}\{\omega\,|\,x_{r}(\omega)\in A_{r}^{+}\}$. Second, we use this approximation to bound the retention probability in the family, hence the coherence ratio.

Step 1. Let $(r_{i})_{i\in\mathbb{N}}$ be a dense sequence in $[0,\tau]$ such that $r_{1}=\tau$; this latter condition is needed such that the decomposition in (53) below is always possible. The events $\mathcal{E}_{n}$ are clearly measurable. Since the paths $t\mapsto x_{t}(\omega)$ of the process (1) are continuous, so is $t\mapsto\bm{f}(t,x_{t}(\omega))=:F_{t}(\omega)$. By corollary 21 we have that $\bm{f}\in C([0,2\tau]\times\overline{X})$. To see $\mathcal{E}_{n}\downarrow\mathcal{E}$, note that by the continuity of $t\mapsto F_{t}$ it holds

[TABLE]

For a finite measure $\pi$ let $\mathbb{P}_{\pi}:=\mathbb{P}(\cdot\,|\,x_{0}\sim\pi)$ denote the law of the process (1) with initial distributions $\pi$. If $\pi$ is not a probability measure, then we denote $\mathbb{P}_{\pi}=\pi(\Omega)\mathbb{P}_{\pi/\pi(\Omega)}$, where $\Omega$ is the entire event space. If, additionally, $\pi$ is a signed measure, then $\mathbb{P}_{\nu}:=\mathbb{P}_{\pi^{+}}-\mathbb{P}_{\pi^{-}}$, where $\pi^{+}$ and $\pi^{-}$ denote the positive and negative parts of the signed measure $\pi$, respectively, in the sense of the Hahn decomposition, i.e. $\pi=\pi^{+}-\pi^{-}$. Now, the $\sigma$-additivity of $\mathbb{P}$ yields also $\mathbb{P}_{\pi}(\mathcal{E}_{n})\to\mathbb{P}_{\pi}(\mathcal{E})$ as $n\to\infty$.

Step 2. By Proposition 5 we have that $\bm{f}(0,\cdot)$ is an eigenfunction of the integral preserving operator $\mathcal{P}_{0,\tau}^{*}\mathcal{P}_{0,\tau}$ at the eigenvalue $e^{2\mu\tau}<1$. Thus we have $\int_{X}\bm{f}(0,\cdot)\,dm=0$, and with it $\int_{X}\bm{f}(t,\cdot)\,dm=e^{-\mu t}\int_{X}\mathcal{P}_{0,t}\bm{f}(0,\cdot)\,dm=0$ again by the integral-preserving property. We define the signed measure $\nu$ via $d\nu(x)=\bm{f}(0,x)\,dm(x)$. With the given scaling of $\bm{f}$ we have that $\int_{A_{\tau}^{\pm}}\bm{f}(\tau,\cdot)\,dm=\pm 1$. By Proposition 5 we have that $\mathcal{P}_{0,t}\bm{f}(0,\cdot)=e^{\mu t}\bm{f}(t,\cdot)$, thus for every $n\in\mathbb{N}$ we have

[TABLE]

The last equality follows from the decomposition of the event $\{x_{r_{1}}\in A_{\tau}^{+}\}=\{x_{\tau}\in A_{\tau}^{+}\}$ into disjoint events $\{x_{r_{i}}\in A_{r_{i}}^{+}\text{ for all }r_{i}>r_{j},\text{ but }x_{r_{j}}\notin A_{r_{j}}^{+}\}$, $j=2,\ldots,n$, and $\{x_{r_{i}}\in A_{r_{i}}^{+}\text{ for all }i=1,\ldots,n\}$. One can see [FJK13] that $p_{j}\leq 0$, because the set of initial conditions $x_{r_{j}}\notin A_{r_{j}}^{+}$ is contained in the non-positive support of $\nu$. It follows that

[TABLE]

Thus, by step 1,

[TABLE]

The same bound can be obtained for the family $\{A_{t}^{-}\}_{t\in[0,\tau]}$. Noting that $\mathbb{P}_{\nu^{\pm}}\leq\|\bm{f}(0,\cdot)\|_{L^{\infty}}\mathbb{P}_{|X|m}$, the claim follows, since $|A_{0}^{\pm}|=|X|\,m(A_{0}^{\pm})$. ∎

Proof of Proposition 9.

The proof is entirely analogous to that of Theorem 8, except that the system of equations containing (53) is altered. The deviating part is the system of inequalities, which follows by the assumptions of the proposition:

[TABLE]

X ∎

Appendix C Bilinear form

The following theorem is also mentioned in [McI68, Appendix] but only partially proven there. We prove the part that is important for this work.

Theorem 22**.**

For a bounded bilinear form $B$ acting on a Hilbert-Space $(H,\langle\cdot,\cdot\rangle_{H},\|\cdot\|_{H})$,

[TABLE]

there exists a linear, bounded operator $T\,:\,H\rightarrow H$ such that

[TABLE]

If $B$ is symmetric, then $T$ is self-adjoint. Further if $B$ is additionally positive definite

[TABLE]

then $T$ is also injective and has a continuous inverse.

Proof.

First let us fix $v\in H$ and consider the linear functional $\ell_{v}\,:\,u\mapsto B(u,v)=\ell_{v}(u)$. By the Riesz representation theorem there exists and element $z_{v}\in H$ such that $\ell_{v}(u)=\langle u,z_{v}\rangle_{H}$ for all $u\in H$. It remains to show that the mapping $T\,:\,v\mapsto z_{v}=T(v)$ is linear and bounded.

Consider $v+\lambda w\in H$ for $\lambda\in\mathbb{R}$ and $v,w\in H$, then

[TABLE]

holds for all $u\in H$. Thus $T$ is linear. The boundedness follows trivially from

[TABLE]

Let us assume that $B$ is symmetric, then

[TABLE]

shows that $T$ is self-adjoint.

Now assume that $B$ is additionally positive definite. The estimate

[TABLE]

immediately implies injectivity. Thus $T^{-1}$ is defined on the range of $T$ ($\text{ran}(T)$) and it is continuous.

[TABLE]

∎