Designing a Finite-Time Mixer: Optimizing Stirring for Two-Dimensional Maps

TL;DR

This paper develops methods to design efficient stirrers for optimizing mixing in two-dimensional fluid flows over finite times, using area-preserving maps and spectral analysis to identify near-optimal stirring protocols.

Contribution

It introduces a novel approach to optimize mixing protocols in 2D flows using a sequence of area-preserving maps and spectral control, including a predictive scheme for Harper maps.

Findings

Single vertical shear with horizontal shearing is nearly optimal for Chirikov maps.

A predictive scheme effectively controls Fourier spectrum in Harper maps.

Near-optimal mixing protocols significantly enhance scalar mixing efficiency.

Abstract

Mixing of a passive scalar in a fluid flow results from a two part process in which large gradients are first created by advection and then smoothed by diffusion. We investigate methods of designing efficient stirrers to optimize mixing of a passive scalar in a two-dimensional nonautonomous, incompressible flow over a finite time interval. The flow is modeled by a sequence of area-preserving maps whose parameters change in time, defining a mixing protocol. Stirring efficiency is measured by a negative Sobolev seminorm; its decrease implies creation of fine scale structure. A Perron-Frobenius operator is used to numerically advect the scalar for two examples: compositions of Chirikov standard maps and of Harper maps. In the former case, we find that a protocol corresponding to a single vertical shear composed with horizontal shearing at all other steps is nearly optimal. For the Harper maps, we devise a predictive, one-step scheme to choose appropriate fixed point stabilities and to control the Fourier spectrum evolution to obtain a near optimal protocol.

Figures

Equations

$$F_T=f_T\circ f_{T-1}\circ \dots \circ f_2\circ f_1.$$

$${\rho }^t(z)={\mathcal{P}}_{f_t}[{\rho }^{t-1}](z)={\rho }^{t-1}(f_t^{-1}(z)),t=1,2,\dots ,T.$$

$$s_y(x,y)$$

$$s_x(x,y)$$

$$(\dot{x},\dot{y})=v_x(x,y)=(Y(y),0).$$

$$f(x,y)=s_x\circ s_y(x,y)=(x+Y(y+X(x)),y+X(x)),$$

$$f^{-1}(x,y)=(x-Y(y),y-X(x-Y(y))).$$

$$X(x)$$

$$Y(y)$$

$$X(x)$$

$$Y(y)$$

$$E_s=\underset{D}{\sup }\Vert s(x,y)-(x,y){\Vert }^2.$$

$$\sum _{t=1}^T(a_t^2+b_t^2)=E_H,$$

$$\sum _{t=1}^T{a}_t^2+T=E_C.$$

$$⟨\rho ⟩=\frac{1}{|D|}{\int }_D\rho d\mu ,$$

$$\text{mbox}Var(\rho )=⟨{\rho }^2⟩-{⟨\rho ⟩}^2,$$

$$\Vert \rho {\Vert }_{H^q}^2=\frac{1}{|D|}{\int }_D|(-L^2\Delta {)}^{q/2}\rho {|}^{2}d\mu ,$$

$$\Vert \rho {\Vert }_{H^q}^2=\sum _{m\in {\mathbb{Z}}^2\setminus \{0\}}|m{|}^{2q}|{\hat{\rho }}_m{|}^2,$$

$${\hat{\rho }}_m={\int }_{{\mathbb{T}}^2}\rho (z)e^{-2\pi im\cdot z}{d}^{2}z.$$

$$\Vert \rho {\Vert }^2=-{\int }_D\nabla \rho \cdot {\nabla }^{-1}\rho \le \Vert \rho {\Vert }_{H^1}\Vert \rho {\Vert }_{H^{-1}},$$

$${\rho }_{i,j}^t=P_{ij,i^{\prime }{j}^{\prime }}{\rho }_{i^{\prime }{j}^{\prime }}^{t-1}.$$

$${\hat{\mathcal{P}}}_{s_y}[\rho ]{)}_m$$

$$({\hat{\mathcal{P}}}_{s_x}[\rho ]{)}_m$$

$${\mathcal{F}}_{j,k}[g]\equiv {\int }_0^1{e}^{2\pi i(jx-kg(x))}dx.$$

$${\mathcal{F}}_{j,k}[id]={\delta }_{j,k},$$

$${\mathcal{F}}_{j,k}[-a\sin (2\pi x)]=J_j(-2\pi ka)=J_{-j}(2\pi ka),$$

$$e^{iz\sin \theta }=\sum _{j\in \mathbb{Z}}{J}_j(z)e^{ij\theta }$$

$${\hat{\rho }}_m^t={\hat{\mathcal{P}}}_{s_x}\circ {\hat{\mathcal{P}}}_{s_y}[{\rho }^{t-1}]=\sum _{n\in {\mathbb{Z}}^2}{\hat{\rho }}_n^{t-1}{\mathcal{F}}_{n_1-m_1,n_2}[X]{\mathcal{F}}_{n_2-m_2,m_1}[Y]$$

$${\hat{\rho }}_m^t=\sum _{n\in {\mathbb{Z}}^2}{\hat{\rho }}_n^{t-1}{J}_{m_1-n_1}(2\pi n_{2}a){\delta }_{n_2-m_2,m_1}=\sum _{k\in \mathbb{Z}}{\hat{\rho }}_{k,m_1+m_2}^{t-1}{J}_{m_1-k}\left(2\pi (m_1+m_2)a\right).$$

$${\hat{\rho }}_m^T=\sum _{n\in {\mathbb{Z}}^2}{\hat{\rho }}_n^{t-1}{J}_{m_1-n_1}(2\pi n_{2}a)J_{n_2-m_2}(2\pi m_{1}b).$$

$${\rho }_{i,j}^t={\rho }^t(x_i,y_j)={\rho }^0(f_1^{-1}\circ f_2^{-1}\circ \dots \circ f_t^{-1}(x_i,y_j))$$

$${\hat{\rho }}_m$$

Full text

1 Introduction

The design of an efficient mixer for a passive scalar in an incompressible flow requires effective mechanical “stirring” of the fluid that eventually will lead to homogenization through diffusion. The point is that when stirring gives rise to fine scale filaments, even a small diffusivity can be highly effective in homogenization. As was noted by Spencer and Wiley [SW51]:

When the materials to be mixed are very viscous liquids…neither diffusion nor turbulence can assist very much in mixing. Rather, mixing must be effected by some complex, continuous deformation which serves to disperse the components to the desired degree.

Thus, as a first step in designing a mixer one can ignore diffusion and study only advection by an incompressible flow.

In this paper we are interested in the finite-time case: how can one design the best mixer for a given initial state and a given finite time? For example, the partitioned-pipe mixer [GAPM03] and the rotated-arc mixer (RAM) [SDMC14] are finite-time mixers in which fluid is injected at one end of a cylindrical pipe, flows past a finite number of mixing elements, and is emitted at the other end. Other examples of this type of mixing include “making up rubber formulations on compounding rolls, using an extruder as a mixer, kneading bread dough, and pulling taffy” [SW51]. Our goal is to develop methods for optimizing the geometry of the elements to maximize mixing.

Most previous studies of mixing by stirring consider autonomous dynamics and concepts appropriate for the infinite-time limit, i.e., positive Lyapunov exponents or measure-theoretic mixing. Indeed, it was noted by Aref that chaotic advection can lead to effective stirring [Are84]. Examples of the optimal design of a periodic sequence of obstacles in a channel or pipe were studied in [GAPM03, JWLW03, GLGS06, HWH14]. However, it is not clear what role chaos plays in the finite-time, nonautonomous case. Froyland et al do study the finite-time case in a flow with diffusion, but assume the advective part is fixed, only modifying diffusion via local perturbations to enhance mixing [FGTW16].

We view the stirring process as a finite sequence of steps; each corresponding to one element of the mixer. The design corresponds mathematically to the selection of a finite sequence $f_{t}:D\to D$, $t=1,2,\ldots,T$, of maps on a domain $D$. Each map is a fixed-time map of an incompressible fluid flow so that each $f_{t}$ is a volume-preserving diffeomorphism. We will study several families of maps, each a composition of shears, see §2. The mixer acts on a passive scalar, represented by a density $\rho:D\to{\mathbb{R}}$. Our goal is to find a sequence $\{f_{t}\}$ that takes a fixed “unmixed” initial state, $\rho^{0}$, to one that is “optimally mixed” after flowing through the full system, i.e., under the composition

[TABLE]

Ignoring diffusion, sources, and reactions, the evolution of the density is determined by the Perron-Frobenius operator,

[TABLE]

In order to select an “optimal” mixing protocol, we must use a “mix-norm” to define the stirring effectiveness; our choice, following [Thi12], is a negative Sobolev seminorm; see §3.

We are primarily interested in the design of efficient, finite-time mixers in a three-dimensional domain. However, in this paper we treat a simpler problem: we assume that the axial velocity is constant and that $f$ corresponds to the map from a two-dimensional cross-section, through a mixing element, to another such cross-section. Since the geometry of the mixing elements will vary, subsequent maps will differ. Each of the $f_{t}$ will be chosen from a given family that satisfies some geometric constraints. For example the RAM consists of an inner pipe surrounded by a spinning outer pipe, and the mixing elements correspond to windows cut into the inner pipe. Each map represents the effect of the angularly sheared flow induced by the boundary drag. The maps, $f_{t}$, will depend on parameters that represent the window width, length and rotation rate.

In this paper we will study two families, Chirikov and Harper maps, see §2. Each map depends on one or more parameters, say $a_{t}$. Optimization then corresponds to selecting the best sequence from a given family, that is, to vector $\textbf{a}=(a_{1},a_{2},\ldots a_{T})$ of parameters for maps in the family. As a constraint we assume that the total energy expended in the mixing process is bounded. The simple version of this is to assume that some norm, $\|\textbf{a}\|$, is fixed. When the norm is small, any fixed map in one of our families is only weakly chaotic on a small subset of phase space. However, varying the parameters from step-to-step can nevertheless result in effective stirring.

The simplest method to find an optimal protocol is to exhaustively search all possible cases; this corresponds to selecting an optimal parameter vector a that satisfies an energy constraint. Thus we can simulate the optimization process by a random optimization method: simply select the best result from a large number of trials that are equi-distributed on the sphere $\|\textbf{a}\|^{2}=E$, corresponding to a fixed energy.

Our investigation seeks to understand mixing by focusing only on advection created by a sequence of stirring events, accounting for diffusion in the mixing measure itself. Through numerical simulations using simple area-preserving maps and through analytical considerations based on Fourier analysis and key characteristics of the maps, we gain some understanding of how to design efficient stirrers.

In §2 we introduce families of maps, consisting of a composition of shears, to describe a stirring mechanism. In §3 we discuss the use of negative Sobolev seminorms as a measure of mixing. In §4 we explain how the Perron-Frobenius method is used to apply the map to a passive scalar. We present results from numerical experiments with Chirikov’s standard map in §5, and the Harper map in §6.

2 Compositions of Shears

As a simple model, we study a family of two-dimensional maps given by the composition of a pair of shears on a periodic domain. Taking the domain to be the torus with unit period, $D={\mathbb{T}}^{2}=({\mathbb{R}}/{\mathbb{Z}})^{2}$, and $z=(x,y)\in{\mathbb{T}}^{2}$, the vertical and horizontal shears have the form

[TABLE]

where $X,Y:{\mathbb{S}}\to{\mathbb{S}}$ are circle maps. Each such shear is the time-one map of a simple incompressible velocity field; for example, $s_{x}(x,y)$ is the time-one flow map for the ode

[TABLE]

Composing the shears (3) gives the one-element map

[TABLE]

that could be realized by imposing a vector field $v_{y}=(0,X(x))$ for $0<t<1$ followed by the vector field (4) for $1<t<2$, i.e., by two successive shears. When the shears are smooth, (5) is an area-preserving diffeomorphism with the inverse

[TABLE]

This family of maps includes many well-studied cases, including Arnold’s cat map [AA68] with $X(x)=x$ and $Y(y)=y$, Chirikov’s standard map [Chi79, Mei92] for which

[TABLE]

and the Harper map [ACB*+*94, Shi02] for which

[TABLE]

The design of a mixer consists of the choice of a mixing protocol, that is, to a sequence of $T$ maps (1). Formally, $T$ may be infinite, but we will typically take $T={\cal O}(10)$. The goal is to optimize the mixing for some maximal total “energy” of the flow. As a measure of this energy we take a squared norm of the shear

[TABLE]

For example for the Harper maps, $E_{s_{x}}=b^{2}$, $E_{s_{y}}=a^{2}$. Since $s_{x}$ corresponds to the velocity (4), the energy is proportional to the kinetic energy of the fluid. For the sequence (1), the parameters $(a_{t},b_{t})$ of the Harper maps are constrained by

[TABLE]

where $E_{H}$ is some given maximum energy. Note that the energy for a sequence of Chirikov maps includes that required to create the horizontal shear $s_{x}$, thus for (6)

[TABLE]

3 Mixing Measures

Mixing, in a measure theoretic sense, is defined through an infinite time limit [LM94]. A transformation $f:D\to D$ that preserves a measure $\mu$ is mixing if for any two measurable sets $A$ and $B$, the measure of the overlap $\mu(A\cap f^{-t}(B))\to\mu(A)\mu(B)/\mu(D)$ as $t\to\infty$. Equivalently, $f$ is mixing if for any initial density $\rho^{0}\in L^{1}(D)$, $\rho^{t}(x)$, determined by (2), converges weakly as $t\to\infty$ to its mean

[TABLE]

where $|D|=\mu(D)$ is the volume of $D$. Mixing is also implied by a positive Kolmogorov-Sinai entropy, computed—according to Pesin—as the mean of the sum of positive Lyapunov exponents of $f$.

For our purposes, measure-theoretic mixing is not an appropriate yardstick for two reasons. Firstly, mixing implies that every $L^{2}$ initial state, not just the those of physical interest, approaches a uniform state. Moreover, even if the entropy is zero, a physical state may still develop fine-scale structure [MMP05]. Secondly, we are interested in a given state becoming “mixed-up” after a finite number of different transformations, not in the limit $t\to\infty$ for a single map.

As was noted by Danckwerts, in the presence of diffusion, the decrease of variance

[TABLE]

is one measure of mixing [Dan52]. Indeed, diffusion will cause $\mbox{Var}({\rho})$ to decrease at a rate proportional to $\|\nabla\rho\|^{2}$ [Thi12]. Thus to be effective, a stirrer should increase gradients so that diffusion is activated. This “intensity of segregation” has been used as to optimize mixing in diffusive flows [KGA*+*01, GAPM03, GLGS06]. However, in the absence of diffusivity and sources, the variance (and indeed each $L^{p}$-norm) of the density is constant.

One class of mixing measures include the mix-norm of Mathew [MMP05] and equivalent Sobolev seminorms [Thi12]. For functions that have zero mean, the squared $H^{q}$-seminorm is defined by

[TABLE]

where $L$ is a normalizing length scale. When $q/2$ is not a positive integer, the Laplacian $\Delta$ is to be interpreted by its Fourier symbol. For the case $D={\mathbb{T}}^{2}$, we set $L=1/(2\pi)$ and (13) reduces to

[TABLE]

where the Fourier coefficient of $\rho$ for mode $m=(m_{1},m_{2})\in{\mathbb{Z}}^{2}$ is

[TABLE]

An effective stirring mechanism will increase $\|\rho\|^{2}_{H^{1}}=\|\nabla\rho\|^{2}$; however, this norm is difficult to compute effectively since it relies on accurate evaluation of high Fourier modes. However, following [Thi12], note that the Cauchy-Schwarz inequality implies that

[TABLE]

assuming $\left<\rho\right>=0$. Since the $L^{2}$ norm is invariant under the Perron-Frobenius dynamics (2), if $\|\rho\|_{H^{-1}}\to 0$, then $\|\rho\|_{H^{1}}\to\infty$. Thus a measure of the effectiveness of the stirring is the degree to which the $H^{-1}$ seminorm decreases. Since high Fourier modes are scaled by $|m|^{-1}$, the $H^{-1}$ norm will not be as sensitive to fine-scale structure, and thus is more amenable to computation. It is possible to show that any negative Sobolev seminorm has similar properties [LTD11], and that if $\|{\rho^{t}}\|_{-1}\to 0$ as $t\to\infty$ for all $\rho^{0}\in L^{2}$, then $f$ is mixing [MMP05].

Below, we use the $H^{-1}$ seminorm as our primary measure of mixing.

4 Computing the Perron-Frobenius Operator

To compute the image $\rho^{t}(x,y)$ and the seminorm (13) numerically, one must represent the density in some finite basis set. One possible technique is the “mapping” [SW51] or “Ulam” method [Ula60, Fro05, FP09, FS10]: choose a grid of $N\times N$ points $(x_{i},y_{j})=(\tfrac{i}{N},\tfrac{j}{N})$ on the torus, with $\rho_{ij}=\rho(x_{i},y_{j})$ and approximate the operator ${\cal P}$ (2) by an $N^{2}\times N^{2}$ stochastic matrix $P$,

[TABLE]

Here $P_{ij,i^{\prime}j^{\prime}}$ is the fraction of the $i^{\prime}j^{\prime}$ cell that is mapped onto the $ij$ cell by $f$. This can be estimated—for a fixed map $f$—by a one step iteration of a large number of trial points in each cell. Such a discretization introduces an effective diffusion that will decrease the variance (12) [KGA*+*01], and explicit diffusion can also be included in the process [GGA12]. However, this method is not practical for our purposes, since the stochastic matrix would need to be recomputed for each new map in the sequence (1).

A second discrete method is to use a Fourier basis. This seems natural since the Sobolev norm that we seek to minimize is most easily evaluated in Fourier space (14). Each shear (3) has a relatively simple representation in Fourier space. Using (15) with $m=(m_{1},m_{2})$, we obtain

[TABLE]

where for each $(j,k)\in{\mathbb{Z}}^{2}$, we define a functional ${\cal F}_{j,k}:C^{0}({\mathbb{S}})\to{\mathbb{C}}$ that gives a Fourier coefficient of the exponential of a function:

[TABLE]

For a linear shear $id(x)=x$ this gives

[TABLE]

the Kronecker-delta, while for a sinusoidal shear such as $X(x)$ in (6) and (7) we have

[TABLE]

using the standard identity

[TABLE]

for the Bessel function. Note that in general ${\cal F}_{j,0}[g]=\delta_{j,0}$, which implies, in particular, that the $(0,0)$-Fourier mode is preserved—as expected for any volume-preserving map.

Composing the two results (16) gives the general formula for updating Fourier coefficients for the map (5):

[TABLE]

since $f^{-1}=s_{y}^{-1}\circ s_{x}^{-1}$. For example for the Chirikov map, (19) becomes, using (18),

[TABLE]

For the Harper map, (19) becomes

[TABLE]

To implement this method numerically, one must truncate the Fourier basis, keeping say, the $N^{2}=(2M+1)^{2}$ modes with $m_{i}\in[-M,M]$. Since the sums in (19) naturally generate modes that are not in the basis set, these must be discarded. Conversely, modes that are not in the truncated basis also may contribute to those in the truncation, and these give errors in the numerical scheme.

A final technique for implementing (2) is to simply evaluate an analytically defined function $\rho^{0}$ at a point $f^{-1}(x,y)$. Given a grid of points $(x_{i},y_{j})$ we can define

[TABLE]

whenever the initial state $\rho^{0}$ is known everywhere. Numerically this is efficient if one is interested only in the final mixed state, $\rho^{T}$, but less so if one would also like intermediate states $t<T$ since $\rho^{t}_{ij}$ cannot be used to compute the next image. Nevertheless, to the extent that one can accurately compute the composition of the maps—subject to the growth of errors due to sensitive dependence on initial conditions—the resulting values of $\rho^{t}$ are accurate when $\rho^{0}$ is Lipchitz, though of course we only know the values on the selected grid.

We choose a regular grid of $N^{2}$ points, so that the Sobolev norm, (14) can be easily evaluated using a fast Fourier transform. Using discrete lattice with $N^{2}$ points, we have $z_{k}=\frac{k}{N}$ for $k\in[0,N)^{2}$, and $\rho_{k}=\rho(z_{k})$, with

[TABLE]

5 Chirikov Map

In this section we discuss various techniques to find optimal mixing protocols for a sequence of Chirikov maps (6). For most of the results in this section, we use an initial Gaussian profile,

[TABLE]

This distribution has standard deviation $\sigma=\sqrt{\frac{1}{2\varepsilon}}\approx 0.2236$. We will also assume that

[TABLE]

so that the total energy (10) is $E_{C}=12.03$. Unless otherwise specified, we use a $301\times 301$ grid and choose the maximum Fourier mode $M=150$ for the computations.

5.1 Single Step Stirring

If the initial state consists of a single, nonconstant Fourier mode, $\rho^{0}(x,y)=\hat{\rho}_{k,l}^{0}e^{2\pi i(kx+ly)}$, then, using (20), the squared seminorm (14) with $q=-1$ at $t=1$ is

[TABLE]

This results in a function that asymptotically decreases as the shear amplitude, $a_{1}$, grows; however, it does not decrease monotonically because of the oscillations of the Bessel functions, see Fig. 1(a). Indeed, if the total energy (8) is small enough, then the optimally mixed protocol, after one step, may be that with $a_{1}=0$.

When a full spectrum initial condition is used in (20), the oscillations seen in Fig. 1(a) are mostly averaged away. Figure 1(b) shows the seminorm for Cauchy (black curve) and Gaussian (red curve) initial spectra. Note that the resulting curves are not symmetric about $a_{1}=0$; for these real initial spectra, a positive shear is more effective than a negative shear. Of course, even for $a_{1}=0$, the density is still sheared since $s_{x}$ is not the identity for (6). For large shear, the Sobolev norm decreases, but—though it is not seen in the figure—the decrease saturates. As $a_{1}\to\infty$, $\|{\rho^{1}}\|_{-1}\to 0.758$, for the Gaussian initial spectrum, and $\to 0.809$ for the Cauchy distribution. It is not possible to “completely” mix after one step, even with arbitrarily large energy.

5.2 Constant Amplitude Stirring

When the amplitude is constant, $\textbf{a}=(a_{o},a_{o},\ldots,a_{o})$, and satisfies the energy constraint (24), then $a_{0}=0.05$. If we think of this as the autonomous standard map, then its phase portrait is almost completely regular; its only features are elliptic islands around the fixed point at $(x,y)=(0,0)$ and the period-two orbit at $(\tfrac{1}{2},\tfrac{1}{2})$, see Fig. 2. The initial density (23) is centered over this period-two orbit, and evolves as shown in Fig. 3. The computations were done using the grid method (22) on an $N=2M+1=301$ squared grid. It is clear from the figure that the stirring is dominated by the horizontal shear $s_{x}$ with $Y(y)=y$ from (6), though some stretching from the saddle fixed point at $(\tfrac{1}{2},0)$ can also be seen.

Figure 4 shows the decay of the Sobolev norm (14) with time. Here the computations were done with a varying number of grid points $N$, as shown in the figure. Note that as $N$ is increased, the computations quickly converge: the value of $\|\rho^{12}\|_{-1}$ varies by less than $10^{-3}$ for $N\geq 51$. This indicates that our standard choice $N=301$ gives sufficient accuracy for $T=12$. The decrease in the mix-norm saturates with time, just as the single step case of §5.1 saturated with increasing amplitude. The point is that repeatedly reapplying the same shear amplitudes may not be the best mixing strategy. This constant case will serve as the benchmark of comparison for our optimization methods.

5.3 Varying Amplitudes: Random Optimization Method

Choosing parameters at random from some distribution can serve as a crude optimization scheme: simply select the best protocol—the solution with smallest seminorm (14)—from the trials. For the standard map (6), there is one parameter, $a_{t}$, for each step of the map, so the parameter space will have $n=T$ dimensions. For convenience we think of the stirring parameters as a vector, letting

[TABLE]

denote a given stirring protocol. The energy (10) then becomes

[TABLE]

To satisfy the energy constraint, we must sample amplitude vectors from the surface of a ball in $\mathbb{R}^{n}$. To ensure uniform sampling on ${\mathbb{S}}^{n-1}$, we randomly select $n$-dimensional vectors from a standard iid normal distribution, scaling each realization to ensure it has the correct energy.

Figure 5(a) shows the distribution of $\|{\rho^{T}}\|_{-1}$ for some two million randomly chosen amplitude vectors with the fixed energy (24) using the Gaussian initial profile (23). Almost all of the trials, 94.8%, perform better than the constant protocol, and the best case achieves the norm $\|{\rho^{T}}\|_{-1}=0.0545$, which is $3.35$ times smaller than the constant protocol. We will refer to this vector as the “random optimal protocol” (RO) even though it is certainly not truly optimal.

The evolution of $\|{\rho^{t}}\|_{-1}$ is shown in Fig. 5(b) for the constant protocol, the RO protocol, and the worst performing solution of the random trials, which we call “RW”. Note that though the seminorm for the RW case decreases on the first step, it remains essentially constant after that. Aside from the first step of the map, the optimal trial significantly out-performs the constant case.

The twenty most successful amplitude vectors are shown in Fig. 6. For these vectors, the largest magnitudes of $a_{t}$ occur primarily in steps $t=2$ or $3$; the first step and all following steps have smaller amplitudes. This suggests a simple scheme: concentrate the energy in the vertical shear in a small number of key steps of the map. We investigate this next.

5.4 Sparse Stirring

Since the most successful protocols found at random have energy concentrated in a small number of map steps, we can modify the optimization scheme to assume sparsity in a. Fixing $d\leq T$ to be the number of nonzero $a_{t}$, we first randomly select $d$ values from the set $\{1,2,\ldots,T\}$ to correspond to the nonzero $a_{t}$, and then perform the random optimization as described above, but with $n=d$. For $d=1$, there are simply $2\times{12}=24$ discrete choices since exactly one of the $a_{t}=\pm 0.03$. For $d>1$, in addition to the discrete selection of $\binom{12}{d}$ shearing times, we must choose sparse vectors on ${\mathbb{S}}^{d-1}$ to fix the energy. Consequently the number of trials should increase with $d$. The best amplitude vectors that we found for each case are shown in Fig. 7; the optimal protocol for $d=1$ corresponds to shearing on the second step with $a_{2}=-\sqrt{0.03}\approx-0.173$. The performance of the optimal protocols as $d$ varies is shown in Fig. 8; note that $\|{\rho^{t}}\|_{-1}$ does not not depend significantly on $d$. Moreover, the optimal sparse protocols achieve even a slightly smaller mix-norm than the full set of random trials. This shows that the latter technique did not find the true optimum, even with more than two million trials.

If we think of the composite map (1) as a single, area-preserving map, for $d=1$ and $a_{2}\neq 0$ it has the form

[TABLE]

when $T\geq 2$. Defining new variables $(\xi,\eta)=(x+y,Ty)$, this reduces to a single standard map in $(\xi,\eta)$ with effective amplitude $Ta_{2}$. Several phase portraits of the map (26) are shown in Fig. 9. As $T$ increases the elliptic regions decrease in size, becoming invisible for $T\geq 8$. Contrast this with the regularity of the constant amplitude case shown in Fig. 2, which has the same total energy. Indeed when $T=12$ the dynamics are nearly uniformly chaotic since, $12|a_{2}|\approx 2.08\gg 0.972/(2\pi)\approx 0.155$, which is the threshold for global chaos for the standard map [Mei92]. This is one reason why the sparse case is such an effective stirrer.

For comparison, Fig. 10 shows the sparse amplitude vectors, a, that maximize the final mix-norm. Note that these have $a_{2}>0$, the opposite of the optimal sparse case. These protocols vary more with $d$ than the optimal protocols did, but the difference in the final mix-norms is small. For the worst $d=1$ case $\|{\rho^{T}}\|_{-1}=0.2041$, which is less than that of the full RW case shown in Fig. 5. A sparse amplitude vector is not as effective in inhibiting mixing as one that is more nearly constant.

To further understand stirring in the $d=1$ sparse case, we can observe the effects of the horizontal shear in Fourier space. When $a=0$, (20) gives

[TABLE]

This implies that each Fourier coefficient has its $y$-wavenumber shifted by $m_{1}$ and, in particular, that the coefficients $\hat{\rho}_{0,m_{2}}^{t}$ remain unchanged. The effect of applying the map with $\textbf{a}=0$ is that the Fourier amplitudes with $m_{1}>0$ are shifted up in $y$-wavenumber; thus, for a smooth initial condition, the amplitudes of the low modes become asymptotically small. Of course the $m_{1}=0$ amplitudes are unaffected, so an optimal protocol must use the nonzero $a_{t}$ to minimize the magnitude of these coefficients. Conversely, a poor protocol must accomplish the opposite.

Let us now consider the effect of a single nonzero $a_{t}$ on the amplitudes with zero $y$ wave number. Since the horizontal shear has no effect on these modes, we may write the updated coefficients in terms of the initial vertically shearing step. Using (20) the coefficients $\hat{\rho}_{0,m_{2}}$ are given by

[TABLE]

After one step with $a_{0}=0$, the largest Fourier coefficients from the initial condition (23) with $m_{1}=0$ are $\hat{\rho}_{0,\pm 1}\approx-0.1204$. For a nonzero $a_{2}$, the $0,1$ mode becomes

[TABLE]

where we truncate at the second mode number. For the optimal case, $a_{2}=-0.1732$, (27) gives

[TABLE]

For this negative $a_{2}$, the Bessel functions $J_{-1}$ and $J_{0}$ have opposite signs. This makes the largest terms nearly cancel, causing the value $\hat{\rho}_{0,1}^{2}$ to decrease by a factor of $4$. This can be seen in Fig. 11, which shows the progression of the density in both physical space (top row) and Fourier space (bottom row). The larger Fourier amplitudes in the figure are shown as darker squares, and moving from the $t=1$ column, to the $t=2$ column, shows that the $(0,\pm 1)$ mode amplitudes decrease significantly. The following steps of shear then push the other large coefficients, the $(\pm 1,\mp 2)$ modes in the figure, to larger $y$-wavenumbers, so that the magnitude of the low Fourier modes diminishes significantly after just the four time steps shown.

By contrast, the $d=1$ RW protocol has the opposite effect. When $a_{2}=+0.1732$, (27) gives

[TABLE]

The sign change causes this mode amplitude to grow in magnitude, resulting in a clumping of density that can not be undone by the remaining steps of purely horizontal shear. This is also illustrated in Fig. 12, which shows the first four steps of the sparse RW protocol in physical and Fourier space.

5.5 Other Initial Conditions

Given this intuition for optimizing mixing for the Gaussian initial density (23), we now determine if the sparse protocol is still effective for other initial conditions. We study initial densities constructed by randomly superimposing periodic bump functions of the form

[TABLE]

so that the center of the bump can range over the torus. Adding a set of sinusoidal perturbations, we construct the density by defining

[TABLE]

We randomly choose the centers $(u,v)\in[0,1)^{2}$, the widths $\varepsilon\in[0.1,5]$, and $K$ and $L$ from $1$ to $5$. The wavenumbers $n_{l}$ ranged over $\{-10,10\}^{2}$ and the phases $\phi_{l}$ over $[0,2\pi)$. The amplitudes $C_{l}$ of the sine terms are kept smaller than $20\%$ of the maximum bump amplitude (to create a noisy effect), and each off the function (29) is scaled so that $\|{\rho^{0}}\|_{-1}=1$.

The question becomes: is a $d=1$ sparse protocol, in which we apply all the energy at a single step, still nearly optimal for varying initial conditions? If this is the case, than the optimal $d=2$ case, in which we divide the energy between two steps should resemble a $d=1$ protocol, in that the energy at one step should be much larger than that at the second. For each of $1200$ initial conditions, we found the best $d=2$ protocol from $10^{4}$ random trials with the energy (24). Figure 13(a) shows a histogram of step times $t_{max}$ at which the larger of the two nonzero amplitudes, $|a_{t_{max}}|$, was applied. The majority of these optimal protocols have $t_{max}\leq 2$: as for (23) the largest amplitude is used early in the twelve applications of the map. Figure 13(b) shows the fraction of protocols for which $|a_{t_{max}}|^{2}>rE$; that is, the fraction of trials for which the largest amplitude step contains more than a fraction $r$ of the available energy. Nearly three-quarters of protocols apply more than 70% of the available energy in one step. This supports the conclusion that the scheme applying a large shear early in the evolution is nearly optimal for the majority of initial conditions.

6 Harper Map

The Harper map (7) differs from the standard map in that the horizontal shear is nonlinear with the additional amplitude parameter $b$. The mixing protocol is now determined by two parameter vectors a and b so that the dimensionality of the optimization problem is doubled, allowing for more flexibility. Moreover, since $b_{t}$ can vary, the energy (9) in the horizontal shear is no longer fixed as it was for (10). In this section, we search for an optimal mixing protocol by means of the random optimization method and also develop other methods that require less computation.

For most of the results in this section, we again use the Gaussian density profile (23), the time constraint $T=12$, and fix

[TABLE]

Unless otherwise specified, we use a $301\times 301$ grid, setting $M=150$.

6.1 Fixed Points

Whenever $ab\neq 0$, the Harper map has exactly four fixed points, $(0,0)$, $(0,\tfrac{1}{2})$, $(\tfrac{1}{2},0)$ and $(\tfrac{1}{2},\tfrac{1}{2})$. For small $a$ and $b$ the map can be approximated by the system

[TABLE]

with Hamiltonian

[TABLE]

Contours of this Hamiltonian are shown in Fig. 14 for various parameters. For the flow, the stability of the fixed points depends only on the relative signs of $a$ and $b$. When $a$ and $b$ are of the same sign, the points $(0,0)$ and $(\tfrac{1}{2},\tfrac{1}{2})$ are hyperbolic and $(0,\tfrac{1}{2})$ and $(\tfrac{1}{2},0)$ are elliptic. When $a$ and $b$ are of opposite sign, the stabilities are are reversed. Note that $(\tfrac{1}{2},\tfrac{1}{2})$, the center of the domain, is elliptic when $ab>0$ and hyperbolic when $ab<0$. When $|a|>|b|$ there is more vertical shearing; conversely, when $|b|>|a|$, the horizontal shear is dominant.

The dynamics of the map (5) with (7)—for fixed parameter values—is similar to the Hamiltonian flow when $|a|$, $|b|\ll 1$, except for chaotic layers near the separatrices; two examples are shown in Fig. 15. For the map, the elliptic fixed points undergo a period doubling bifurcation at $|ab|=\pi^{-2}$, becoming hyperbolic with reflection. However, for the total energy levels that we will choose, none of the shear amplitudes will be large enough for this to occur.

6.2 Constant Shear Amplitude

If the shears have constant amplitudes $|a_{t}|=|b_{t}|=\gamma_{o}$, then the energy constraint (9) implies that $\gamma_{o}=\sqrt{\frac{E_{H}}{2T}}$. It would seem to be counterproductive to have successive signs of the horizontal (or vertical) shears oscillate, so we will consider here two characteristic cases, $a_{t}=b_{t}=\gamma_{o}$, shown in Fig. 16, and $a_{t}=-b_{t}=\gamma_{o}$, shown in Fig. 17. The contours of the Hamiltonian (30) are also shown in these figures to indicate the approximate dynamics of the next map to be applied—indeed as Fig. 15(a) shows, these contours reasonably approximate the dynamics since the amplitude, $\gamma_{o}=0.1443$, is small when $E_{H}=0.5$.

Since the initial density is centered at $(\tfrac{1}{2},\tfrac{1}{2})$, which is elliptic when $a$ and $b$ have the same sign, and the standard deviation of the Gaussian profile (23) is $0.2236$, which is not too large, the dominant dynamics seen in Fig. 16 is a sheared rotation about the center which results in very little mixing. By contrast, when $a$ and $b$ have opposite signs, the point $(\tfrac{1}{2},\tfrac{1}{2})$ is a saddle, so the density in Fig. 17 is stretched along the unstable manifold until it reaches the saddle at $(0,0)\equiv(1,1)$, where it divides and then moves back towards $(\tfrac{1}{2},\tfrac{1}{2})$ along its stable manifolds.

Figure 18 shows the decay of the mix-norm for these two cases and various values of the maximum mode number, $M$. As with the standard map, the decrease in the mix-norm slows as $t$ increases. As the grid size increases, the values of $\|{\rho^{T}}\|_{-1}$ converge—little difference is seen when $M>25$, well below our standard resolution.

Among the four possible constant amplitude cases, $(a_{t}=\pm 0.1443,b_{t}=\pm 0.1443)$ for $E_{H}=0.5$, the optimal mixing protocol for the centered Gaussian (23) corresponds to the opposite sign case depicted in Fig. 17, achieving $\|{\rho^{T}}\|_{-1}=0.0636$. Of course, since the initial condition is symmetric about the fixed point, flipping both signs results in the same Sobolev norm. This raises the question: will the opposing amplitude stirring scheme be optimal for different energy budgets? The mix-norm as function of constant amplitudes $(a_{t},b_{t})=(a,b)$ is shown in Fig. 19. It is interesting that though the optimal protocols for both large and small energy have $a\approx-b$, there is an intermediate regime where this scheme is not optimal.

6.3 Random Optimal Protocol

Computing the optimal mixing protocol for the Harper map, requires choosing the $n=2T$ dimensional vector $(\textbf{a},\textbf{b})$ on the energy sphere

[TABLE]

that has the smallest $\|{\rho^{T}}\|_{-1}$. Following the same procedure as §5.3, Fig. 20 shows a histogram of the mix-norm for our standard Gaussian initial density. Note that the two constant amplitude protocols (orange and red dots in the figure) correspond to near extremes of the range of $\|{\rho^{T}}\|_{-1}$.

The evolution of the density for the best protocol from $1.7\times 10^{6}$ trials is shown in Fig. 21; again we call this the RO protocol. Note that for this protocol $a_{1}$ and $b_{1}$ have opposite signs, causing the initial density to be stretched along the unstable manifold of the point $(\tfrac{1}{2},\tfrac{1}{2})$; this is as expected from our earlier discussion of the constant case. It is interesting that the density evolution in Fig. 21 exhibits not only the stretching seen in the constant amplitude hyperbolic case of Fig. 17, but also considerable folding, the second hallmark of two-dimensional chaos.

It is also important that the amplitudes of the RO protocol do not have equal magnitudes. Indeed, for several of the steps ($t=5,8,10,$ and $12$) one of the amplitudes is at least three times smaller than the other, so that the dynamics is nearly a horizontal or vertical shear (Recall that the Hamiltonian contours for these steps are overlaid on the previous density.). In §6.5, we exploit this by examining stirring protocols that allow only one of the shear amplitudes to be nonzero for each step.

In addition, the energy of the shears in the RO protocol varies considerably with time. Defining the one-step energy by

[TABLE]

we see in Fig. 22, most of the energy is expended in steps $t=1$, $2$, $10$ and $11$.

The evolution of $\|{\rho^{t}}\|_{-1}$ for the RO protocol is shown in Fig. 23(a); this figure also includes the two constant amplitude cases from §6.2. As seen in Table 1, at $t=12$, the RO protocol just outperforms the hyperbolic constant case, showing that it is possible to do better than just stretching the density along the manifolds of the central fixed point. It is interesting, however, that the constant amplitude case has lower mix-norm when $4<t<9$. During this period, as can be seen in Fig. 21, the RO protocol is setting up improved stirring for the last three steps by folding the density contours; this is accomplished using disparate amplitudes of the two shears.

As the energy varies, so does the achievable minimum mix-norm. The optimal protocol as a function of energy gives a Pareto frontier, shown in Fig. 23(b). This figure was constructed using $6.5\times 10^{5}$ trials for the initial condition (23) spread over 20 values of $E_{H}\in[0,2]$. Also shown, for each energy, are the mix-norms for the worst trial and for the constant amplitude, hyperbolic protocol. Note that the latter protocol appears to be optimal when $E_{H}<0.3$; however, for larger energies the RO protocol beats the constant case.

6.4 Stepwise Minimization

The decay of the mix-norm for the RO protocol shown in Fig. 23(a) is not monotone. Though $\|{\rho}\|_{-1}$ decreases sharply during the first three steps, it appears to increase slightly during steps four and five, finally decreasing for the remaining steps. This is odd, since one might expect the optimal twelve-step protocol to minimize the mix-norm at each step as much as possible. Instead of considering the time-history as a whole, could it be as effective to choose $a_{t}$ and $b_{t}$ independently to maximize the decrease in the mix-norm from step $t-1$ to step $t$?

For simplicity, we will look for a protocol that has same stepwise values of the energy as the RO protocol. Thus, if the amplitudes $(a_{t},b_{t})$ for the RO protocol at time $t$ have energy (31), then we will search on this circle for the optimal one-step protocol: for each time $t$, we apply (2) to find

[TABLE]

Recall that the energy $E^{RO}_{t}$ was shown in Fig. 22. Since (32) requires succession of a one-dimensional searches at each time, it is much quicker than the full random optimization. However, the restriction that to a known energy distribution means that we will not be looking for more general protocols. We call this new stepwise minimization protocol “SW”.

Figure 24 illustrates the stepwise minimization by comparing the RO protocol of §6.3 to a SW protocol at three times. For this illustration, we used the density $\rho^{t-1}$ obtained by the RO protocol as the initial state for the SW minimization, shown in the second column. So to obtain the stepwise minimization at $t$, we perform $t-1$ steps of the RO protocol and then use (32). The amplitudes $(a_{1},b_{1})$ for the SW protocol are essentially the same as those for RO–this can be seen in the third column which shows $\|{\rho^{t}}\|_{-1}$ as a function of $(a,b)$. The RO amplitudes, the green dot in this panel, and the selected SW amplitudes nearly overlap.

This is also true at $t=3$, though there is a larger difference, presumably because none of the selected random trials landed on the best point. However the difference at $t=2$ is nontrivial: the sign of $b_{2}$ flips between the RO and SW cases. Note that the central point $(\tfrac{1}{2},\tfrac{1}{2})$ is elliptic for the RO amplitudes, whereas it is a saddle for the SW amplitudes. This supports our hypothesis that the RO protocol, by looking ahead, does not minimize the mix-norm at every step: for this step it has emphasized folding instead of stretching.

Fully implementing (32) (without reinitializing to the RO solution at every step) for the initial density (23) gives $\|{\rho^{t}}\|_{-1}$ shown in Fig. 25(a). This figure also shows the mix-norm of the RO protocol from Fig. 23 for comparison. Though the overall performance of the two protocols is about the same at $t=12$, the SW protocol reaching $\|{\rho^{T}}\|_{-1}=0.0401$, the RO protocol does slightly outperform the stepwise minimizer. It seems that an optimal mixing scheme cannot simply consider each step of the map independently, but must instead forgo a one-step improvement to setup the density for future steps of the map.

6.5 Horizontal or Vertical Shearing

For many of the time steps of the RO protocol, as illustrated in Fig. 21, one of the shear amplitudes is much smaller than the other, giving either primarily horizontal or vertical shearing. In this section we investigate protocols in which there is only one nonzero shear at each step. To simplify the problem even further, we assume that the energy per step is constant. Defining the amplitude $\gamma_{o}=\sqrt{E_{H}/T}$; there are now just four options for the parameter pair at a given step: $(a,b)\in\{(\gamma_{o},0),(-\gamma_{o},0),(0,\gamma_{o}),(0,-\gamma_{o})\}$.

To find an optimal stirring protocol with the new constraints is easier than a full RO since there are only four choices for parameters at every step. In fact except for $t=1$, only three of these can decrease the mix-norm. For example if $(a_{1},b_{1})=(\gamma_{o},0)$, then the choice $(a_{2},b_{2})=(-\gamma_{o},0)$ simply undoes the previous step, so it cannot be an optimal choice. Therefore, only the three options $(\gamma_{o},0),(0,\gamma_{o}),(0,-\gamma_{o})$ need be considered at $t=2$. In this way, there are $4\times 3^{T-1}$ distinct shear protocols to consider, i.e., $\sim 7\times 10^{5}$ when $T=12$. The scheme found from randomly sampling this discrete space is compared to the RO protocol in Fig. 25(b). The optimal shear protocol (labeled “Shear RO” in the figure) just beats the full RO protocol, as shown in Table 1.

An even simpler method is to do a stepwise minimization, as in §6.4, now imposing the single-shear restriction. In this case there are only $4+(T-1)3$ choices, i.e., $40$ when $T=12$. The optimal result in this case gives the mix-norm curve labelled CS SW (Constant Shear SW) in Fig. 25(b). This simplified protocol achieves $\|{\rho^{T}}\|_{-1}=0.0693$, which is not comparable to the better performing protocols. Nevertheless, it still performs better than 99.6% of the randomly sampled amplitude vectors shown in Fig. 20. Note that, for this specific initial density profile, the hyperbolic constant scheme is slightly better (see Table 1); however, its effectiveness depends strongly on the localization of (23) near the saddle point.

To gain some insight into how to choose the direction of shear based on the current density distribution, we can consider what the optimal one-step parameter choice would be when the density can be described by a single sinusoid. Consider the initial condition

[TABLE]

for wave vector $n$, and phase shift $\phi$. The resulting Fourier expansion will have only two nonzero terms, both of modulus one. Using this initial condition and fixing the amplitude $\gamma_{o}=0.2041$, we compute the mix-norm $\|{\rho^{1}}\|_{-1}$ as a function of mode numbers for $200$ randomly chosen phase shifts. Figure 26 shows a density plot of the number of times choosing $a\neq 0$ results in better one-step mixing than choosing $b\neq 0$. According to Fig. 26, we should choose $a\neq 0$ whenever $|n_{1}|<|n_{2}|$ and $b\neq 0$ if $|n_{1}|>|n_{2}|$. The lines $|n_{1}|=|n_{2}|$ are ambiguous.

To understand why this is true, consider the simple case of (33) with $\phi=0$, so that the Fourier expansion will contain the only nonzero terms $\hat{\rho}_{n}^{0}=\hat{\rho}_{-n}^{0}=1$. Recalling (21), the Fourier coefficients at $t=1$ will be

[TABLE]

When $b=0$ and $a$ is nonzero, the only terms that survive are those that satisfy $m_{2}=\pm n_{2}$ so that

[TABLE]

for $k\in{\mathbb{Z}}$. On the other hand, if $a=0$ and $b$ is nonzero, we are left with the terms in which $m_{1}=\pm n_{1}$, giving

[TABLE]

Since the $H^{-1}$ norm (14) has the wavenumber in the denominator, the implication is that we would then want to choose $b=0$ when $|n_{1}|<|n_{2}|$, since the resulting terms $\hat{\rho}_{k,\pm n_{2}}$ will have higher combined wave number than the terms $\hat{\rho}_{\pm n_{1},k}$. The Bessel function terms will also tend to be smaller, at least when the shears have comparable amplitude and the mode numbers are large. This supports the more general results of Fig. 26.

To apply this insight to more general density profiles which have a full Fourier spectrum, we define the partial sums

[TABLE]

The shear choice is then determined by the relative sizes of these partial norms: if $S_{<}$ is smaller we choose $(|a|,|b|)=(0,\gamma_{o})$; otherwise, if $S_{>}$ is smaller we choose $(|a|,|b|)=(\gamma_{o},0)$.

To make this a complete protocol, we must also select the signs of the shears. Note that successive steps $(a_{t},0)$ to $(0,b_{t+1})$ are equivalent to one step with amplitudes $(a_{t},b_{t+1})$, so that the choice of sign is equivalent to the selection of the location of elliptic and saddle points of the map. The opposite switch, from $(0,b_{t})$ to $(a_{t+1},0)$ does not give a single Harper map of the form (7), but is instead equivalent to its inverse with the signs of both parameters switched. Thus the locations of the elliptic and saddle cases are still the same as those shown in Fig. 14.

If the density is concentrated in a specific location, then it seems reasonable that a sign choice that makes the nearest fixed point a saddle will beat the sign choice that makes it elliptic. Indeed, for the Gaussian centered at $(\frac{1}{2},\frac{1}{2})$, it was important to choose opposite signs to allow this fixed point to be hyperbolic, since making it elliptic would result in virtually no mixing.

To determine the relationship between the location of density peaks and the optimal sign choice, we can do a simple one-step experiment with various localized initial densities—we use the initial state (28), varying the center $(u,v)$ of the profile. Figure 27 shows the best sign choice as a function of $(u,v)$. The color at a given point indicates the sign choice that most often resulted in the best one-step mixing under randomly chosen values of the width parameter $\varepsilon$. As expected, if the density is localized near $(\frac{1}{2},\frac{1}{2})$, we should choose $a$ and $b$ of opposite sign. This generalizes: if the density is localized near any fixed point, the relative signs should be chosen to make that fixed point hyperbolic. The regions in which we choose a given sign combination form nearly perfect squares around each fixed point. To understand why this is, consider the approximate eigenvalues and eigenvectors given by the Hamiltonian (30). The fixed point $(\frac{1}{2},\frac{1}{2})$ has eigenvalue-eigenvector pairs

[TABLE]

Now assume that we are using the constant shearing scheme such that $|a|=|b|=\gamma_{o}$. Also assume that the density is concentrated near $(\frac{1}{2},\frac{1}{2})$ so that we at least know we want $a=-b$ to make the point hyperbolic. This causes the eigen-pairs to reduce to

[TABLE]

The vector $(\text{sign}(b),1)$ corresponds to the stable manifold and $(-\text{sign}(b),1)$ to the unstable. The closer the concentration of density is to a hyperbolic fixed point, the more it gets stretched along the unstable manifold of that point. Ideally the density should be localized to straddle the stable manifold so that it gets pushed towards the fixed point and then spread out along the two branches of the unstable manifold. This implies that if the concentration is located in the $1^{st}$ or $3^{rd}$ quadrants, relative to $(\frac{1}{2},\frac{1}{2})$, we should choose $b>0$ to make the local stable manifold have positive slope. By contrast, if the concentration lies in the $2^{nd}$ or $4^{th}$ quadrants we choose $b<0$ to make the local stable manifold have negative slope. Performing a similar analysis for each of the fixed points implies the results obtained experimentally in Fig. 27.

From this we can extract an algorithm to determine effective sign choices. We denote the $16$ regions in Fig. 27 by $R_{1},R_{2},\ldots,R_{16}$, and define the indicator functions

[TABLE]

For each region we compute a partial mix-norm

[TABLE]

Under the hypothesis that region with the largest $S_{i}$ dominates the need for mixing, we pick the signs of $(a,b)$ for this region according to Fig. 27.

A combination of the Fourier and sign algorithms gives a full protocol choice for the multistep shear problem. We attempt two different schemes. The first allows for the signs of $(a,b)$ to change between steps of the map if the location of the density has changed. We refer to this as the “CS Full” method. The second scheme chooses signs for $(a,b)$ based solely on the initial condition maintaining this sign choice throughout the protocol. We refer to this as the “CS FS” (Fixed Signed) method. the CS Full method has the advantage of tracking the density and switching the fixed points accordingly. However, the CS FS method allows for more consistency in the direction of stirring.

The CS Full algorithm performs better when we limit the instances in which a sign change can occur. Switching the signs, and thereby changing the direction of stirring, often causes a momentary increase in the mix-norm. Performing a sign change when it is not strictly necessary can negatively affect the results. However, forgoing a sign change can cause the mix-norm to level off. Recall that the RO solution outperformed the Hyperbolic Constant scheme because it sacrificed some steps to cause folding which will allow for more effective stretching later. Experimentally we find that CS Full performs the best when we only switch signs if the difference between the two largest $S_{i}$ values is more than 10% of the total mix-norm.

Figure 28 shows the performance of the algorithm over $10^{5}$ random choices of initial conditions (29) normalized such that $\|{\rho^{0}}\|_{-1}=1$. The minimum, mean, and standard deviation of $\|{\rho^{T}}\|_{-1}$ for each distribution depicted in Fig. 28 can be found in Table 2. The scheme with sign switching beats out the one with fixed signs slightly in terms of mean final mix-norm. However, the SW minimizer outperforms both analytic solvers. The standard deviation of mix-norms using the sign switching methods is smaller than that for fixed signs. Indeed, the distribution of mix-norms for the fixed sign method has a right tail that extends nearly to $\|{\rho^{T}}\|_{-1}=1$, suggesting that for some initial conditions, the solver performs extremely poorly and results in virtually no mixing. Using sign switching helps to reduce the size of this tail.

7 Conclusions

We investigated the design of mixers by studying the optimization of stirring behavior of a two-dimensional device modeled by a finite number of nonlinear shears, under a constraint of fixed energy. The minus-one Sobolev seminorm, one example of a “mix-norm”, provides an cost-effective way to measure mixing efficiency without explicitly computing diffusion. We computed the results of stirring by evaluating the Perron-Frobenius operator on a uniform grid, assuming that the initial density profile is known analytically. This allowed us to avoid the aliasing errors and numerically induced diffusivity of other methods. One of the most important insights is that effective mixing is highly dependent on the behavior of the lower order Fourier modes: in order to construct a good stirring scheme, one must target the low-order modes with the highest amplitudes.

For the first model, a composition of Chirikov standard maps, concentrating the nonlinear shear energy in a small number of steps was shown to give an efficient mixing scheme. For this map, the horizontal shear is assumed fixed, and using the energy for vertical shearing early in the protocol takes maximal advantage of the shift in Fourier coefficients to higher mode numbers.

For the Harper map model both vertical and horizontal shears can be controlled, allowing for more diverse stirring mechanisms and more interesting behaviors. A simplified scheme assumed that each step corresponds to either a horizontal or a vertical shear with constant energy. Within this reduced parameter space, we devised an algorithm that combines insight on the action of the shears in Fourier space with knowledge of the character of the fixed points, to produce a near-optimal stirring protocol. The method was shown to remain effective when applied to a variety of initial density profiles.

Table 1 compares the optimal stirring protocols produced by the different methods when when applied to the initial condition $(\ref{eq:InitialGaussian})$, and Table 2 compares the overall performance of solvers for varying initial conditions (29). Our results suggest that both consistency and change—contradictory requirements—are necessary to produce an optimal stirring protocol. While the former leads to stretching near saddle points, the latter can result in enhanced folding. Both were required to obtain the optimal protocol for the Harper map. With only stretching, the mix-norm tends to level off; but with too much folding, stirring can be reversed, causing the mix-norm to grow. An optimal protocol will tend to use hyperbolic behavior if a smooth portion of the density is concentrated near a fixed point—if a nearby fixed point is elliptic, mixing efficiency can be decreased even though there will be (weak) shearing.

Though it is clear that both stretching and folding are required to design an effective mixer, it is less clear how to best combine these processes to optimize mixing. Observing stretching and folding in a wider variety of mixing problems may lead to more general results about how one should balance the two mechanisms to induce mixing.

The models that we studied here have their fixed points at fixed locations. Another stirring mechanism to investigate would be one that allows these locations to vary. For instance, one could introduce variable phases into the shears to control the locations of the fixed points at every step. We hope to use results gained from studying these simple, theoretical models to aid investigations of more physically applicable, three-dimensional models such as the partitioned pipe and rotated arc mixers [GAPM03, SDMC14].

It would be worthwhile to compare the $H^{-1}$ norm to a model that implements diffusion. The effectiveness of our measure depends on the assumption that diffusion is small, i.e., that the Péclet number is large. In the future, we plan to include diffusion to get a sense of the errors induced by using the mix-norm for a purely advective model when the Péclet number is finite.