# The influence of the vorticity-scalar correlation on mixing

Xi–Yuan Yin1, Wesley Agoua1, Tong Wu1,2, Wouter J.T. Bos1 1 LMFA-CNRS, Ecole
Centrale de Lyon, Université de Lyon, Ecully, France
2 Theoretical Physics I, University of Bayreuth, Bayreuth, Germany

###### Abstract

We investigate the role of the correlation between a scalar quantity and the
vorticity in two-dimensional mixing at infinite Péclet number. We assess,
using a diffusivity independent mixing-norm, the dynamics of both Galerkin-
truncated ensembles and freely evolving two-dimensional scalar mixing. Both
statistical mechanics and numerical experiments show how the mixing-rate is
attenuated when vorticity and scalar are initially correlated. Since the
vorticity is shown to be a poorly mixing scalar, the results suggest that, in
general, mixing can be enhanced by minimizing the correlation between
vorticity and passive scalar.

## Introduction.

The present investigation considers the mixing (or stirring) of a passive
scalar. The scalar can represent temperature fluctuations, pollutants or any
other passively advected quantity. The propension of a given flow to mix, is
determined by the dimensionality of space, the properties of the velocity
field and of the scalar, boundary conditions, the type of stirring, and
initial conditions. We discuss these different features now, for the special
case of two dimensions, considering an incompressible fluid, advecting a
passive scalar.

The influence of boundary conditions is clearly very important, as was shown
by Raynal and Gence [1] comparing chaotic mixing in periodic and unbounded
domains. Chaotic mixing in wall bounded domains was considered by Gouillart et
al. [2, 3], where the presence of unmixed fluid near solid walls was shown to
slow-down the mixing. In turbulent flows, the generation of enstrophy near the
walls was shown, on the contrary, to enhance mixing [4]. In the present case
we will focus on the simplest academic case of periodic boundary condition,
avoiding thereby the issues related to the influence of solid walls.

The stirring protocol is another important factor that we will omit from our
analysis. Indeed, to obtain efficient mixing, the stirring protocol can be
optimized under constraints [5, 6]. We will let the flow freely evolve from
initial conditions. The precise type of the initial conditions, in particular
the typical size of the initial scalar compared to the typical velocity length
scales, or domain-size can importantly influence the mixing [7, 8]. Besides,
the study of strange eigenmodes reveal the existence of scalar fields of
slowest decay which attract the long-time evolution of general initial scalar
distributions [9, 10], though the rate of decay of scalar variance is still
sensitive to initial conditions. Again, to keep the description as general as
possible, we will consider random initial velocity and scalar fields. These
fields are constituted by random Fourier modes with a prescribed energy
spectrum for the velocity field, and analogously, a prescribed variance
spectrum for the scalar quantity.

We have at this point simplified and specified the flow, the boundary
conditions, omitted stirring and prescribed random initial conditions. The
problem is then determined, except for one point: to what extent the initial
velocity and scalar fields are correlated? The influence of this initial
correlation, more specifically the vorticity-scalar correlation is the main
subject of this study. We will predict theoretically and assess by numerical
experiments how the initial vorticity-scalar correlation affects the mixing
efficiency. This question has received little attention [11, 12] and only a
very partial understanding of the influence of the scalar-vorticity
correlation on mixing exists. In particular, the question whether such an
initial correlation is beneficial or, on the contrary, detrimental for mixing
has not been answered.

## Governing equations and constants of motion.

We consider the advection of a passive scalar field $\phi({\bm{x}},t)$ in a 2D
periodic box $\Omega$. The scalar is transported by a two-dimensional
incompressible and inviscid flow whose velocity field ${\bm{u}}({\bm{x}},t)$
is described by the Euler equations. The curl of the velocity, or vorticity,
${\omega}(\bm{x},t)=\nabla\times{\bm{u}}({\bm{x}},t)$ evolves according to

$\displaystyle(\partial_{t}+{\bm{u}}\cdot\nabla){\omega}=0$ (1a)
$\displaystyle{\bm{u}}=-\nabla\times\Delta^{-1}{\omega}$ (1b)

where the velocity can be recovered from the vorticity by the Biot-Savart
relation Eq. (1b). The initial condition is
${\omega}({\bm{x}},0)={\omega}_{0}({\bm{x}})$.

The linear advection equation under the velocity field ${\bm{u}}$ then
describes the evolution of the scalar $\phi$:

$(\partial_{t}+{\bm{u}}\cdot\nabla)\phi=0,$ (2)

with $\phi({\bm{x}},0)=\phi_{0}({\bm{x}})$. As we wrote above, we consider a
general case where the initial velocity and scalar field
${\omega}_{0}({\bm{x}}),~{}\phi_{0}({\bm{x}})$ are random. The investigation
of more specific initial conditions is left for further research.

The quadratic invariants, or constants of motion, of a two-dimensional Euler
velocity field are

$\displaystyle
E=\int|\bm{u}|^{2}d{\bm{x}},~{}~{}~{}W=\int{\omega}^{2}d{\bm{x}}.$ (3)

The energy $E$ and enstrophy $W$ characterize the velocity field. One can
define a typical wavenumber

$k_{u}=\sqrt{W/E}.$ (4)

This wavenumber characterizes the typical lengthscale of the velocity field.
In the absence of sources and sinks, it remains constant. The scalar field is
also characterized by two invariants,

$\displaystyle S=\int\phi^{2}d{\bm{x}},~{}~{}~{}Q=\int{\omega}\phi d{\bm{x}}.$
(5)

Whereas the variance $S$ is determined by the scalar field only, the last
quantity $Q$ measures the correlation between the scalar and the vorticity
field.

Two-dimensional turbulence contains an infinite number of additional
invariants, Casimirs of the form $\int f({\omega})d{\bm{x}}$. Equivalently for
the scalar, additional invariants can be defined. The investigation of the
influence of these invariants on mixing will be left for further research. The
choice to retain only the 4 quantities $E,W,S,Q$ is theoretically supported by
the fact that in the Galerkin truncated system that we will consider first,
they are the only surviving constants of motion, conserved, triad-by-triad by
the nonlinear and advection terms [13, 14]. Clearly since we ignore all the
Casimirs, our investigation can be refined in future studies.

The quantity $Q$ is, like the other invariants, conserved in the absence of
sources or sinks, so that its initial correlation persists throughout the
evolution of the system. We can further define the relative correlation,

$h=Q/\sqrt{SW},$ (6)

a quantity which is constrained by the bounds $0\leq h\leq 1$. Our goal is to
measure and understand the influence of $h$ on mixing.

## Mixing norm.

The eventual faith of a scalar in a realistic flow is determined by the
diffusion. However, the mixing is, at sufficiently small values of the
diffusivity, piloted by the stretching and folding of the large-scale scalar
structures. Therefore we want to measure the mixing independently from the
values of diffusion (and viscosity). We use a mixing norm previously
introduced in the study of optimal mixing [15] (see also [16]). This measure
is a function of the inverse gradient $|\nabla^{-1}\phi|$; its value is
representative of the importance of the large scales of $\phi$. We introduce a
characteristic wavenumber,

$k_{\phi}=\sqrt{\frac{\int\phi^{2}d{\bm{x}}}{\int|\nabla^{-1}\phi|^{2}d{\bm{x}}}}.$
(7)

An increase of this characteristic wavenumber indicates that the scalar is
transported towards smaller scales. Typically this implies that the flow is
mixing, since the eventual scalar dissipation in a realistic flow is taking
place at the small scales (or large wavenumbers). We also note that in a
periodic domain, for scalar $\phi$, one can check by Parseval’s identity that
$\int|\nabla^{-1}\phi|^{2}d{\bm{x}}=\sum_{\bm{k}}k^{-2}|\hat{\phi}(\bm{k})|^{2}=\int|-\nabla\times\Delta^{-1}\phi|^{2}d{\bm{x}}$
so that the typical wave number $k_{u}$ is the ratio of the $L^{2}$ to
$H^{-1}$ seminorms for $\omega$, and $k_{\phi}$ is the same ratio for $\phi$.

Since the Euler-equations conserve both energy and enstrophy, $k_{u}$ (defined
in Eq. (4)) does not evolve. However $k_{\phi}$ can evolve and this evolution
will characterize the mixing rate. In the following, we will investigate the
evolution of $k_{\phi}$ for given values of $h$ and $k_{u}$. We will start by
deriving a theoretical prediction.

## Correlated and uncorrelated scalars.

The passive scalar and the vorticity share the same evolution equation
(2),(1). The difference is that the vorticity is related to the velocity by
Eq. (1b). Since the evolution-equation is identical, we can decompose the
scalar into a correlated part which is proportional to the vorticity and an
uncorrelated part $\phi^{\prime}$,

$\phi(\bm{x},t)=(h\sqrt{S/W}){\omega}(\bm{x},t)+\phi^{\prime}(\bm{x},t),$ (8)

with the following relations

$\int\phi^{\prime}{\omega}d{\bm{x}}=0\textrm{~{}~{}~{}and~{}~{}~{}}S^{\prime}\equiv\int\phi^{\prime
2}d{\bm{x}}=S(1-h^{2}).$ (9)

The quantities $S$ and $S^{\prime}$ are independently conserved. The
decomposition of the scalar (8) allows to simplify the problem and to
formulate predictions for the statistical mechanics of the scalar using known
results only. Indeed, the evolution of $\phi^{\prime}$ to a statistical
equilibrium state is governed by one invariant, $S^{\prime}$. The evolution of
the correlated part will be fully determined by the evolution of ${\omega}$,
since $h$ is constant.

(a) |
---|---
(b) |
Figure 1: Results and predictions for equilibrium ensembles of the Galerkin-
truncated system. The wavenumber range is confined to $k\in[1,42]$. (a)DNS
results for the Equilibrium spectra of the scalar variance $E_{\phi}(k)$.
Theory predicts a spectrum of $h^{2}(S/W)k^{2}E(k)+c_{\phi}k$ from eq. (13).
For clarity of the image, the $h=0.38$ and $h=0.77$ data are shifted up by a
half and full decade respectively. (b) Influence of the relative correlation
$h$ on the mixing scale $k_{\phi}$ predicted from eq. (14). Results are shown
for different values of the characteristic velocity scale $k_{u}$.

## Theoretical predictions for the Galerkin truncated system.

We consider that the velocity and scalar excitation are projected on a finite
number of Fourier modes in two-dimensional Fourier-space. The smallest and
largest wavenumber are $k_{0},k_{1}$ respectively. The evolution-equations are
also restricted to keep the excitation confined on these modes. The scalar
variance is then obtained by

$\int_{k_{0}}^{k_{1}}E_{\phi}(k)dk=S.$ (10)

And a similar expression relates $E_{\phi^{\prime}}(k)$ and $S^{\prime}$. In
such a set-up, the equilibrium solution of the uncorrelated passive scalar is
given by the equipartition solution

$E_{\phi^{\prime}}(k,t=\infty)=c_{\phi}k,$ (11)

where $c_{\phi}=2S^{\prime}/(k_{1}^{2}-k_{0}^{2})$ (see for instance [17]).
This simple expression describes the equilibrium spectrum associated with
$\phi^{\prime}$. The correlated part is fully determined by the vorticity
field and its dynamics is determined by both the kinetic energy and the
enstrophy. Indeed, the energy spectrum determines both quantities according to

$\int_{k_{0}}^{k_{1}}E(k)dk=E,~{}~{}~{}\int_{k_{0}}^{k_{1}}k^{2}E(k)dk=W.$
(12)

For the scalar spectrum we have therefore

$E_{\phi}(k,t)=E_{\phi^{\prime}}(k,t)+h^{2}\frac{S}{W}k^{2}E(k,t).$ (13)

To verify this relation, we carry out numerical simulations of the system.
Standard Cartesian pseudospectral solvers are the ideal tools to verify the
present predictions of equilibrium statistical mechanics, since these methods
are based on a Galerkin truncation [18] and since they are known to correctly
conserve the invariants of the system. We use the two-dimensional code Ghost
[19], based on a vorticity-streamfunction formulation and second order time-
integration. We start from a Gaussian-shaped initial spectrum
$E(k,0)\sim\exp(-b(k-k_{L})^{2})$, with random Fourier-phases and $k_{L}=4$
and $b=0.125$, setting the total energy at 1 and enstrophy at approximately
28, i.e. $k_{u}\approx 5.3$.

We plot in Fig. 1(a) the spectra $E_{\phi}(k,t=t_{\infty})$ and the estimate
Eq. (13), with $E_{\phi^{\prime}}(k)=c_{\phi}k$. For each value of $h(t=0)$ a
different statistical realizations of the initial conditions for $E(k,0)$ is
used. The time $t_{\infty}$ is determined by running the simulations for
$t\approx 8000t_{*}$, with $t_{*}=W^{-1/2}$, a time sufficient for the
statistics to converge to a statistical equilibrium. It is observed that the
predicted spectrum is well approached by the theoretical estimate. We can then
use this theoretical estimate of the scalar spectrum to compute analytically
the mixing scale $k_{\phi}$.

Using Parseval’s relation between the mixing norm and the scalar spectrum
$\int|\nabla^{-1}\phi|^{2}d{\bm{x}}=\int_{k_{0}}^{k_{1}}k^{-2}E_{\phi}(k)dk$,
we find after some elementary manipulation for the mixing scale,

$k_{\phi}=\left[h^{2}k_{u}^{-2}+(1-h^{2})k_{\phi^{\prime}}^{-2}\right]^{-1/2}.$
(14)

where
$k_{\phi^{\prime}}=\left((k_{1}^{2}-k_{0}^{2})/(2\ln(k_{1}/k_{0}))\right)^{1/2}$.
The two limits of these expression are, for $k_{\phi}(h=1)=k_{u}$, where the
mixing scale is fully determined by the typical velocity length-scale. And for
the uncorrelated case we have $k_{\phi}(h=0)=k_{\phi^{\prime}}$.

Expression (14) is plotted for fixed resolution (or $k_{\phi^{\prime}}$) in
Fig. 1(b). We see there that the influence of finite $h$ is for most cases
deleterious, i.e, $k_{\phi}$ is a decreasing function of $h$. Only for large
initial $k_{u}$, this effect is not observed. Nevertheless, this latter case
corresponds to initial conditions where the initial scalar field is
concentrated at scales much larger than the excited velocity scales. In
realistic flows this will lead to very slow mixing. The reason that
statistical mechanics predict enhanced mixing for this case is that only final
states are predicted, irrespective of the time it takes to reach such an
equilibrium state. This constitutes the limit of equilibrium statistical
mechanics to predict realistic flows and we will therefore in the following
use numerical simulations to understand the dynamics.

## Simulations of scalar mixing by the Euler-equations.

Figure 2: Results from CMM computations. Figure (a) shows the final time
scalar spectra for various correlation values. Figure (b) shows the mixing
enhancement $M(\phi,t_{*})$ as function of the scalar-vorticity correlation
when transported by various flows for a given amount of time. The average
enhancement over a sample of 30 scalars is shown in full line, the shaded area
covers one standard deviation. Figures (c) and (d) show the initial and final
time scalar fields for the $h=0$ and $h=1$ cases respectively.

What we can retain from the results of equilibrium statistical mechanics is
that it indicates a tendency of less efficient mixing when the vorticity-
scalar correlation is strong. The worst mixing is attained when the velocity
scale is largest, $k_{u}/k_{0}\approx 1$. We note that this is a relevant
case, since two-dimensional flows are often dominated by space-filling
structures.

We will assess now if, as suggested by statistical mechanics, mixing is
reduced when $h\neq 0$, dynamically, when the system is not Galerkin
truncated. Simulating the evolution of vorticity and the mixing of a scalar
without the influence of viscosity and diffusion is numerically challenging.
Nevertheless, recent advances in numerical analysis [20] have led to the
development of a numerical method allowing to compute, in a precise and
efficient way, exactly this evolution. The method we therefore apply is the
characteristic mapping method [21], a method which, instead of advancing the
Eulerian fields such as vorticity and velocity, evolves the inverse flow-map.
An advantage of this method is that once the flow map is determined for a
given flow, it allows to determine directly the scalar field $\phi(\bm{x},t)$
for all possible initial conditions $\phi_{0}(\bm{x})$. Furthermore, the
solution of this advection is not frequency truncated and solves the non-
diffusive advection. The Casimirs and invariants of the system are naturally
conserved so that we can compare the statistical predictions for the truncated
system with the long-time solutions of the non-truncated one.

We start from an initial velocity field given by an enstrophy spectrum
concentrated around $k_{u}$ with random phases and normalized to unit kinetic
energy. The initial scalar spectrum is defined by the bump function with
$k_{\phi}\approx 3.3$. The evolution of the energy spectrum and the scalar
spectrum are shown in Fig. 2(a) for a range of correlation values. We see that
as correlation decreases, more scalar energy is moved towards the small scales
indicating a better mixing.

We quantify the extent of the mixing enhancement using the mixing norm of
$\phi$. Define the enhancement

$M(\phi,t)=\frac{k_{\phi(t)}}{k_{\phi(0)}}=\sqrt{\frac{\int|\nabla^{-1}\phi(0)|^{2}d{\bm{x}}}{\int|\nabla^{-1}\phi(t)|^{2}d{\bm{x}}}},$
(15)

where the second equality is obtained from the fact that
$\int\phi^{2}d{\bm{x}}$ is a conserved quantity. The enhancement number
$M(\phi,t)$ measures the relative increase of the typical wavenumber in time
due to transport by the fluid. The essential difference with the results of
equilibrium statistical mechanics is that the enhancement measure $M(\phi,t)$
takes into account the initial conditions, whereas equilibrium statistical
mechanics only gives information on the final state, irrespective of initial
conditions. We note that in the $h=1$ fully correlated case,
$\phi=(Q/W)\omega$, in which case $k_{\phi}=k_{u}=\sqrt{W/E}$ which is a
constant of motion and hence $M(\phi,t)$ must remain 1.

To test the influence of the scalar-vorticity correlation on the mixing
enhancement, we perform the scalar advection using three flows obtained from
the Euler equations where the characteristic length scale $k_{u}$ is
approximately $1.7,3.3$ and $9.4$ respectively. We generate for each flow a
sample of 30 initial conditions $\phi^{\prime}_{0}$ of random phase and bump
function spectrum such that $\phi^{\prime}_{0}$ are fully decorrelated from
the vorticity $\omega_{0}$ and $k_{\phi^{\prime}_{0}}\approx 3.3$. Each
$\phi^{\prime}$ correspond to the $h=0$ case, and together with $\omega$, we
generate a family of initial conditions $\phi$ ranging from $h=0$ to $h=1$ and
whose initial characteristic length $k_{\phi_{0}}$ range from
$k_{\phi^{\prime}_{0}}$ to $k_{u}$. This effectively tests for the mixing
enhancement generated by a flow which moves at larger, same, and smaller
spatial scales than the advected quantity. Fig. 2 (b) shows the results from
our tests. The time $t$ at which the enhancement is assessed is fixed at
$t\approx 20t_{*}$ with $t_{*}=W^{-1/2}$ as to normalize the average turnover
over the three flows. The same tests can be performed by instead scaling $t$
with $E^{-1/2}$ but have similar results due to the scaling symmetry of the
Euler equations. In all cases, we observe that the mixing enhancement is
monotonously decreasing with respect to the correlation.

## Conclusion.

The fact that initial velocity conditions are important for the efficiency of
mixing is not surprising. However, the fact that for a given random scalar
field and a given random velocity field, fixing $E,W,S$, the mixing can as
dramatically be affected by the value of $Q$ as observed in Fig. 2 is
something which might not have been expected from the outset. Indeed, the
mixing scale for $h=1$ remains constant, while for the fully uncorrelated case
the mixing enhancement, for all our tests, is typically a factor 2. The
underlying reason, which we now understand, is that, if the vorticity is fully
correlated with the scalar, the scalar cascade is entirely locked to the
enstrophy cascade. It seems from the test cases considered here, that due to
its link with velocity field, the vorticity is indeed the least mixed scalar
field in 2 dimensional flows.

In the definition of the mixing-scale and the invariants of the system, we
have not explicitly assumed that the system is turbulent. We do therefore
think that these ideas can be applied to non-turbulent systems as well. A
practical implication of the present insights is that, in optimization
problems (e.g. [22, 23, 5]), the constraint of minimizing the correlation
between scalar and vorticity, might be a powerful way to approach the optimal
mixing regime in both turbulent and non-turbulent flows.

Acknowledgments. All DNS simulations were carried out using the facilities of
the PMCS2I (École Centrale de Lyon). For the purpose of Open Access, a CC-BY
public copyright licence has been applied by the authors to the present
document and will be applied to all subsequent versions up to the Author
Accepted Manuscript arising from this submission.

## References

  * [1] Florence Raynal and Jean-Noël Gence.  Energy saving in chaotic laminar mixing.  International journal of heat and mass transfer, 40(14):3267–3273, 1997.
  * [2] E. Gouillart, N. Kuncio, O. Dauchot, B. Dubrulle, S. Roux, and J.-L. Thiffeault.  Walls inhibit chaotic mixing.  Phys. Rev. Lett., 99(11):114501, 2007.
  * [3] E Gouillart, O Dauchot, B Dubrulle, S Roux, and JL Thiffeault.  Slow decay of concentration variance due to no-slip walls in chaotic mixing.  Physical review. E, Statistical, Nonlinear, and Soft Matter Physics, 78(2 Pt 2):026211–026211, 2008.
  * [4] Benjamin Kadoch, Wouter J. T. Bos, and Kai Schneider.  Efficiency of laminar and turbulent mixing in wall-bounded flows.  Phys. Rev. E, 101:043104, 2020.
  * [5] Dimitry PG Foures, Colm P Caulfield, and Peter J Schmid.  Optimal mixing in two-dimensional plane poiseuille flow at finite péclet number.  Journal of Fluid Mechanics, 748:241–277, 2014.
  * [6] Maximilian F Eggl and Peter J Schmid.  Mixing enhancement in binary fluids using optimised stirring strategies.  Journal of Fluid Mechanics, 899:A24, 2020.
  * [7] Alexander A. Schekochihin, Peter H. Haynes, and Steven C. Cowley.  Diffusion of passive scalar in a finite-scale random flow.  Physical Review E, 70(4):046304, October 2004.
  * [8] Peter H Haynes and Jacques Vanneste.  What controls the decay of passive scalars in smooth flows?  Physics of Fluids, 17(9), 2005.
  * [9] Weijiu Liu and George Haller.  Strange eigenmodes and decay of variance in the mixing of diffusive tracers.  Physica D: Nonlinear Phenomena, 188(1):1–39, January 2004.
  * [10] Jai Sukhatme and Raymond T Pierrehumbert.  Decay of passive scalars under the action of single scale smooth velocity fields in bounded two-dimensional domains: From non-self-similar probability distribution functions to self-similar eigenmodes.  Physical Review E, 66(5):056302, 2002.
  * [11] Greg Holloway and Stefan S Kristmannsson.  Stirring and transport of tracer fields by geostrophic turbulence.  Journal of Fluid Mechanics, 141:27–50, 1984.
  * [12] M. Lesieur and J.R. Herring.  Diffusion of a passive scalar in two-dimensional turbulence.  J. Fluid Mech., 161:77, 1985.
  * [13] David Fyfe, David Montgomery, and Glenn Joyce.  Dissipative, forced turbulence in two-dimensional magnetohydrodynamics.  Journal of Plasma Physics, 17(3):369–398, 1977.
  * [14] Robert H Kraichnan and David Montgomery.  Two-dimensional turbulence.  Reports on Progress in Physics, 43(5):547, 1980.
  * [15] George Mathew, Igor Mezić, and Linda Petzold.  A multiscale measure for mixing.  Physica D: Nonlinear Phenomena, 211(1-2):23–46, 2005.
  * [16] Jean-Luc Thiffeault.  Using multiscale norms to quantify mixing and transport.  Nonlinearity, 25(2):R1, 2012.
  * [17] C.E. Leith.  Diffusion approximation for turbulent scalar fields.  Phys. Fluids, 11(7):1612, 1968.
  * [18] S. A. Orszag.  Numerical methods for the simulation of turbulence.  Physics of Fluids, 12:II–250, 1969.
  * [19] Pablo D Mininni, Duane Rosenberg, Raghu Reddy, and Annick Pouquet.  A hybrid mpi–openmp scheme for scalable parallel pseudospectral computations for fluid turbulence.  Parallel computing, 37(6-7):316–326, 2011.
  * [20] Olivier Mercier, Xi-Yuan Yin, and Jean-Christophe Nave.  The characteristic mapping method for the linear advection of arbitrary sets.  SIAM Journal on Scientific Computing, 42(3):A1663–A1685, 2020.
  * [21] Xi-Yuan Yin, Olivier Mercier, Badal Yadav, Kai Schneider, and Jean-Christophe Nave.  A characteristic mapping method for the two-dimensional incompressible euler equations.  Journal of Computational Physics, 424:109781, 2021.
  * [22] Evelyn Lunasin, Zhi Lin, Alexei Novikov, Anna Mazzucato, and Charles R Doering.  Optimal mixing and optimal stirring for fixed energy, fixed power, or fixed palenstrophy flows.  Journal of mathematical physics, 53(11), 2012.
  * [23] Zhi Lin, Jean-Luc Thiffeault, and Charles R Doering.  Optimal stirring strategies for passive scalar mixing.  Journal of Fluid Mechanics, 675:465–476, 2011.