Eulerian indicators for predicting and optimizing mixing quality

Many different methods exist for measuring, quantifying and predicting the quality of a mixture resulting from a fluid mixing device. In general, these rely on Lagrangian information from the system, as they involve the evolution of individual trajectories. In this paper, we propose tools that are formulated using only Eulerian information, and yet can predict the Lagrangian mixing properties of a mixing device. The proposed Eulerian indicators are motivated by the recent linked twist map (LTM) formulation of fluid mixing. We describe their features and application in the context of four different fluid flows: the blinking vortex flow and the partitioned pipe mixer (two paradigm flows from the early studies of chaotic mixing), a model of a channel flow containing ‘separatrices’ in the cross-sectional flow (which is a common flow feature in a large number of micromixers), and a flow generated by pulsed source–sink pairs (which has been used as a model for mixing in DNA microarrays).


Introduction
Mixing of fluids is concerned, fundamentally, with the relative displacement of fluid particles. Essentially all definitions of mixing, as well as measures of mixing efficiency, utilize some sort of mathematical operation to quantify the rearrangement of arbitrary 'blobs' of fluid in space and time as a result of their evolution in the flow. Consequently, mixing is inherently a Lagrangian phenomenon whose description is often computationally intensive as it requires the simulation of large numbers of fluid particle trajectories for, possibly, long times. As a consequence, the prediction of the mixing properties of a flow tends to be problematic.
In terms of quantifying mixing, we note that there are many different measures in the literature that can be used to quantify the extent to which two substances are mixed. For example, these include scale and intensity of segregation (Danckwerts 1953a(Danckwerts , 1953b) (see also the discussion in Krasnopolskaya et al (1999)) and topological braiding exponents (Thiffeault 2005) and the MixNorm (Mathew et al 2005) (see also the closely related method in Stone and Stone (2005)). The differences and relative merits of many of these have been discussed in a recent review paper (Funakoshi 2008). However, the purpose of this paper is not to develop another diagnostic to measure the quality of mixing. Rather, we seek to develop Eulerian indicators that will serve to predict the quality of mixing of a flow without the need to first simulate a mixing experiment. We will do this by examining the kinematic characteristics of a flow that promote mixing and classifying them according to whether they are 'Eulerian' or 'Lagrangian' in nature. The particular aspect of this classification that is relevant to our goals is that Eulerian characteristics can be determined solely from characteristics of the flow field, while an understanding of Lagrangian characteristics requires the simulation of the advection 3 of 'blobs' of fluid particles. Armed with this particular classification of characteristics, the goal is then to determine for which type of flows a knowledge of certain Eulerian characteristics suffices to predict that the flow will mix well. It is significant that, for flows with this property, we can predict that they will mix well without simulating mixing experiments for the flow. It is often easier and more efficient to modify mixing devices by changing the Eulerian properties of the resulting flows, and in this way our Eulerian approach should lead to methods to create 'optimal' mixing devices more efficiently. We will refer to these particular Eulerian quantities as Eulerian indicators, since they are Eulerian in nature and they indicate that these particular Eulerian features of the flow play a significant role in influencing mixing.
We remark that Yannacopoulos et al (1998) and King et al (2001) have previously introduced an Eulerian diagnostic to describe regions of a three-dimensional (3D) flow where the particle trajectories are not integrable. The idea originates from the work described in Mezic and Wiggins (1994) and is based on a diagnostic that measures deviations of the flow from a symmetric case, for which the flow is integrable.
Our Eulerian indicators are motivated by the recent linked twist map (LTM) formulation of fluid kinematics described in Wiggins (2004a, 2004b) and Sturman et al (2006). The 'mapping approach' following from dynamical systems theory is now well established as a framework for studying mixing in flows purely as a result of advection (or, at least, where diffusion is weak); see Ottino (1989a), Wiggins (1992) Wiggins and Ottino (2004). In this context, the motion of fluid particles (or particles following fluid pathlines) is considered as a map, which is typically referred to as a Poincaré map. Practically, this means that rather than observe a fluid particle tracing out a continuous curve in the domain, we instead 'stroboscopically' view the particle at fixed, discrete, intervals of time. The collection of discrete points in the flow domain obtained in this way is referred to as a Poincaré section. We will see that Poincaré sections can give one a useful picture of the mixing properties of a flow (but not the complete picture, as can be seen by comparing our Poincaré sections with the variance of the concentration). For example, if it is possible to choose a point at an arbitrary location in the flow domain such that under iteration by the Poincaré map it appears to fill the domain, then the flow would seem to 'mix well' (i.e. it would appear to be ergodic). In fact, the mathematical notion of mixing is stronger than ergodic. Roughly, mixing means that under evolution by the flow nearby particles separate and 'forget' that they started near each other. Ergodicity merely implies that an arbitrary particle comes arbitrarily close to every point in the domain (eventually), so in an ergodic flow neighbouring points need not separate in time. It is important to realize that chaotic advection alone is not adequate for a flow to mix perfectly unless the chaos exists throughout the entire flow domain. It is well known that chaos can co-exist with regular regions, or 'islands'. The significance of LTMs in this context is that they are one of the few (non-uniformly hyperbolic) dynamical systems for which it is possible to prove rigorously the existence of the mathematical property of strong mixing (and even the Bernoulli property, in some sense the strongest type of chaos) on regions of the flow having positive area (in fact, more mathematically, full Lebesgue measure (Burton and Easton 1980, Devaney 1980, Przytycki 1983, Springham 2008, Sturman et al 2006, Wojtkowski 1980. The significance to fluid mechanics stems from the fact that the proofs rely on geometrical and kinematical features of the maps in a way that can lead to practical design principles for 'optimal' mixing devices (Hertzsch et al 2007, Ottino andWiggins 2004a). For example, LTMs embody the fluid kinematical situation of 'monotonic and sufficiently strong shear rates' and 'streamline crossing', for which there is a significant amount of experimental evidence indicating that these are important kinematical mechanisms for creating 'chaotic mixing' (Ottino 1990, Sturman et al 2006. In order to make this more quantitative we first need to define the notion of an LTM in a manner that is sufficient for our needs.

The mathematical definition of a linked twist map on the plane
We consider LTMs defined on a subset of the plane. This is appropriate for the specific applications that we have in mind since the flows will be defined in a region with a solid boundary. LTMs can also be defined on domains with double-periodic boundary conditions, i.e. the torus. This is appropriate for the so-called egg-beater flows discussed in Franjione Ottino (1992), Ottino (1989b) and Sturman et al (2006), and also for mixers based on source-sink pairs with reinjection, described in McQuain et al (2004), Stremler and Cola (2006), Cola et al (2006) and Raynal et al (2004), analysed in Hertzsch et al (2007) and discussed in section 3.3 of the present paper. Whether the LTM is defined on a (subset of) the plane or a torus is significant since the proofs of strong mixing properties differ in the details in each case (but not in the overall strategy and approach). This is discussed in Sturman et al (2006) and Springham (2008).
Let A 1 denote an annulus whose inner (r 1 ) and outer (r 2 ) boundaries represent streamlines centred at a point in the plane denoted by O 1 , and let A 2 denote an annulus centred at a point in the plane denoted by O 2 with boundaries ρ 1 and ρ 2 , as shown in figure 1. As in this figure, the boundaries r 1 , r 2 and ρ 1 , ρ 2 must be chosen so that the annuli A 1 and A 2 intersect in two distinct regions.
For the map F, the radial coordinate r is constant, i.e. it is unchanged under iteration by the map, and the angular coordinate θ increases by f (r ) on each iteration, which is an amount that depends on the radial coordinate. From the point of view of fluid mechanics the annulus A 1 corresponds to a region of closed streamlines. A similar description holds for the map G on the annulus A 2 . We assume that the functions f and g increase monotonically from f (r 1 ) = g(ρ 1 ) = 0 to f (r 2 ) = 2π m and g(ρ 2 ) = 2π n, where m and n are integers, respectively. This is the so-called twist condition that we make more precise below. With these assumptions the maps F and G are each examples of twist maps.
A linked twist map is then defined to be the composition H = G • F, defined on the union A 1 ∪ A 2 . Note that to compose the maps F and G formally, we require a change of variable transformation. As we have already mentioned, LTMs can be shown to have global strong mixing properties under certain conditions. These conditions were first given in Wojtkowkski (1980), and are discussed in detail in Sturman et al (2006) in the context of fluid mixing. New results in this context have recently been obtained by Springham (2008). The significance of these conditions is that they are inherent in the geometrical and analytical definition of the LTM and, most significantly for us, they can be classified in terms of Eulerian and Lagrangian characteristics of the map. We summarize these characteristics heuristically and without detail as follows: Eulerian characteristic 1: monotonic and sufficiently strong twists/shear rate The absolute value of the derivatives df /dr and dg/dρ must be bounded from below by a sufficiently large constant. In the LTM context, we refer to this property by saying that the twist is sufficiently strong and monotonic. In the fluid mechanical context, this property implies that the shear rate is sufficiently strong and non-vanishing. Eulerian characteristic 2: transversely intersecting streamlines/'streamline crossing' The boundaries r 1,2 and ρ 1,2 must be chosen so that the annuli A 1,2 intersect in two distinct regions. This has the effect of ensuring that the intersection of each circle r = constant in A 1 with each circle ρ = constant in A 2 is sufficiently transverse. Lagrangian characteristic: 'good' recurrence to 'saddle-like' overlap regions The annuli A 1 and A 2 intersect in two distinct regions. It can be proven that almost all (that is, a set of full Lebesgue measures) points in A 1 and A 2 visit these regions upon iteration by the LTM infinitely often. Moreover, as a direct result of the two Eulerian characteristics above, points visiting these regions experience 'saddle-point-like' or 'hyperbolic' behaviour, and since almost all trajectories visit these regions infinitely often, it can be shown that almost all trajectories are saddle-like or hyperbolic in stability type. There are numerous technical issues here. For example, we did not comment on the sign of the derivatives df /dr and dg/dρ, which would give rise to a discussion of co-rotating and counter-rotating cases. These are discussed in some detail in Sturman et al (2006) and Springham (2008) and lie at the heart of what we mean by 'good' recurrence. We do not want the stretching that occurs when points visit these regions to be 'undone' (e.g. stretching in the 'opposite' direction) on later visits. Admittedly, this discussion around good recurrence is somewhat vague-but the issues here are clearly Lagrangian in nature and lie at the heart of the ergodic theory analysis required for the Lagrangian aspect of this problem. In the ergodic theory of dynamical systems, this principle is made precise by the method of invariant cones (Alekseev 1969).

6
The main point here that is relevant to this paper is that the first two Eulerian characteristics, coupled with certain geometrical features of the intersecting annuli, imply that this Lagrangian property holds for almost all points in A 1 ∪ A 2 , and from this it follows that the entire domain A 1 ∪ A 2 exhibits strong mixing properties.
In Wiggins and Ottino (2004) and Sturman et al (2006), it is shown that the kinematics of a large class of flows can be formulated as some type of LTM. The composition of the two twist maps corresponds to considering a flow that is defined by two alternating or blinking flows. The alternating or blinking action can correspond to either blinking in time (for example, the blinking vortex flow) or blinking in space (for example, flow in a spatially periodic duct).
The blinking vortex flow (Aref 1984), an example of a flow modelled by blinking between two different flows in time, is a paradigm in the study of chaotic advection, and the fluid particle motion in this flow is described precisely by an LTM on the plane (Sturman et al 2006, Wiggins 1999, Wiggins and Ottino 2004. Each vortex generates a family of closed streamlines, where the time for a fluid particle to traverse a closed streamline decreases monotonically with the distance from the vortex (this is the origin of the phrase 'twist map'), and families of streamlines of the two vortices cross transversely when superimposed as described above.
In 3D, the mixing of fluid can be studied with 2D mappings for channel, or duct, flows, where the flow is steady and the cross-sectional flow varies periodically in space. The cross-sectional flow varies between the beginning and end of a spatially periodic element of the duct and the map takes particles from the cross section of the flow at the beginning to the cross section of the flow at the end of a period along the flow trajectories. The classic example of this situation is the partitioned pipe mixer, first studied in Khakhar et al (1987). The partitioned pipe mixer was considered to be a blinking flow in which particles were advected under the cross-sectional velocity fields at the beginning and end of a spatially periodic element of the duct for a length of time needed for the particle to travel the length of the periodic element. A general construction of such maps is carried out in detail in Wiggins and Ottino (2004) and Sturman et al (2006).
In general, these two types of blinking flows (i.e. roughly, 'blinking in time' or 'blinking in space') can be thought of in terms of generalized LTMs. Physically, we could think of LTMs as providing a mathematical framework for rigorously analysing flows consisting of 'transversely oriented shears', which are encompassed by Eulerian characteristics 1 and 2 described above. The Lagrangian characteristic ensures that almost all points visit the regions of transversely oriented shears infinitely, often in a way that almost all trajectories will be saddle-like in stability.
In recent years, there have been a wide range of mixing devices created based on both the blinking in time and blinking in space approaches. The highly viscous, 'slow flow' regime is ideal for the blinking approach to modelling flows. One way in which flows of this type are created is through the motion of boundaries of cavities containing a fluid. The literature for chaotic advection in this type of boundary-driven 2D cavity flow is large; see e.g. Chien et al (1986), Jana et al (1994aJana et al ( , 1994b, Leong and Ottino (1989) and Stremler and Chen (2007). This type of flow regime fits well with applications in microfluidics, and examples of blinking 2D electro-osmotic flows driven by electromagnetic fields are given in Bau (2002, 2005).
More examples of spatially periodic steady 3D flows for which 2D Poincaré maps can be constructed describing the motion of fluid particles from the beginning of one spatial period to the next spatial period are given in Khakhar et al (1987), Kusch and Ottino (1992), Franjione 7 and Ottino (1991), Jana et al (1994a), Stroock et al (2002), Metcalfe et al (2001Metcalfe et al ( , 2003Metcalfe et al ( , 2006, Kim et al (2004), Kang et al (2007) and Lopez and Graham (2008).
While the mathematical framework of LTMs has proven ideal for analysing and proving the existence of strong mixing properties on sets of nonzero area in the flow, many of the examples from the fluids literature given above contain separatrices in the two flow patterns that are superimposed, and these types of maps are outside the area where rigorous mathematical results on strong mixing have been proven. In order to gain intuition on such flows, we introduce in section 3.2 a simple kinematic model of a blinking flow where the two flow patterns contain separatrices in a way that captures the same kinematical properties of some of the mixers described above (e.g. Bau (2002, 2005), Stremler and Chen (2007), Stroock et al (2002)).
This paper is organized as follows. In section 2, we discuss the elements of mixing that are considered by our Eulerian indicators through their definitions in sections 2.1-2.3. In section 2.4, we consider a fundamentally Lagrangian property of mixing that is not captured by our Eulerian indicators and discuss its significance. In section 3, we consider four examples: the blinking vortex flow (section 3.1), the 'separatrix model' (section 3.2), the pulsed source-sink pair flow (section 3.3) and the partitioned pipe mixer (section 3.4). We consider each flow for four different operating conditions (obtained by selecting the parameter values appropriate for each flow) and compare the behaviour of the Eulerian indicators with Lagrangian measures of mixing given by the variance of the concentration and Poincaré maps. In section 4, we conclude with final remarks and suggestions for future study.

The characteristics of mixing and Eulerian indicators
In this section, we introduce our Eulerian indicators that we use to quantify the Eulerian characteristics described above for specific flows. We will also discuss some of the issues surrounding the Lagrangian characteristic also introduced above.

An Eulerian indicator for quantifying Eulerian characteristic 1: monotonic and sufficiently strong twists/shear rate
A measure of the amount of shear, or strength of twist, at a point in the flow can be obtained by the following simple numerical scheme. For a flow defined by the streamfunction i , i = 1, 2, at each point (x, y) in the flow domain, the flow has direction (u (i) x (x, y), u (i) y (x, y)). We define a pair of points (x i,1 , y i,1 ) and (x i,2 , y i,2 ) a distance δ i 1 from (x, y), each on opposite sides of the line defined by the vector (u (i) x (x, y), u (i) y (x, y)), such that the line segment with endpoints (x i,1 , y i,1 ) and (x i,2 , y i,2 ) is perpendicular to the vector (u (i) x (x, y), u (i) y (x, y)). Define, using the usual Euclidean length, the lengths of vectors l i,1 = |(u (i) x (x i,1 , y i,1 ), )|, and then the strength of twist at each point can be given by β i = |l i,1 − l i,2 |/2δ i . This construction is illustrated in figure 2. Since for a blinking flow we have two streamlines passing through each (x, y), we can define the product of shearsβ(x, y) = β 1 β 2 . To produce a single value for the flow, we averageβ over the flow domain to give an Eulerian indicator for Eulerian characteristic 1: (1)

An Eulerian indicator for quantifying Eulerian characteristic 2: transversely intersecting streamlines/streamline crossing
At each point (x, y) of the flow domain, we can associate an angle between streamlines as follows. The angle between tangents to these streamlines can be easily computed, since the tangent to the streamline defined by the streamfunction i through (x, y) is given by the vector v i = (u (i) x (x, y), u (i) y (x, y)). Then the angle between such vectors is given simply by An important point to note is that we do not take into account the direction of flow along the streamline-we are simply interested in the transversality of the intersection. For this reason, we replace α(x, y) witĥ so that 0 α π/2. See figure 3. Finally, to produce a single value representing the degree of transversality in a flow, we average this quantity over all points in the flow domain to produce our Eulerian indicator for Eulerian characteristic 2: 9 Figure 3. Computation ofα. We disregard the direction of streamlines (dotted lines) in this computation, so that in this figure,α(x 1 , y 1 ) =α(x 2 , y 2 ).

Combining the Eulerian indicators for monotonic and sufficiently strong twists/shear rate and transversely intersecting streamlines/streamline crossing
The exponential stretching associated with chaotic mixing stems from shearing in two complementary directions (although as we discuss in the following section, shearing in two different directions is not enough to guarantee chaotic mixing). Thus mixing might be expected to be most efficient when the strongest shears act in the most transverse directions (that is, closest to perpendicular). With this in mind, we defineγ (x, y) =α(x, y)β(x, y), the product of the first two Eulerian indicators-the measures of transversality of the streamlines with the strength of the twist. We define our third proposed Eulerian indicatorγ by again averaging over the whole domain of the flow:

Issues associated with the Lagrangian characteristic of 'good' recurrence to saddle-like overlap regions
Here we discuss the issue associated with the Lagrangian characteristic of 'good' recurrence to saddle-like overlap regions. We re-emphasize that the focus of this paper is on providing design principles for optimizing mixing based on Eulerian quantities. However, there are situations where this goal cannot be realized, and this Lagrangian characteristic plays the role of 'spoiler' in those situations, which we now describe with a concrete example. We will consider two different LTMs, each defined on the entire torus, denoted by T 2 . Coordinates on the torus are realized by taking the unit square, S = {(x, y) : 0 x 1, 0 y 1}, and identifying the horizontal and vertical boundaries, i.e. x = 0 = 1 and y = 0 = 1 (more mathematically, (x, y) ∼ (x + 1, y) ∼ (x, y + 1)), i.e. we take doubly periodic boundary conditions on the unit square. We explain what we mean by the LTM being defined on the entire torus in the context of our first example. We consider the map defined on T 2 . This map acts as a shear, from left to right, where the strength of the shear increases monotonically from the bottom of S to the top of S, as indicated in figure 4(a). General LTMs need not be defined for the entire range of variables on the torus in the direction of increasing shear (see Sturman et al 2006). However, we will consider F 1 to be defined on the entire torus, i.e. 0 x 1 and 0 y 1.  The action of the maps F 1 and G 1 whose composition, , is the Arnold cat map.
A significant advantage in considering LTMs defined on the entire torus is that recurrence is trivial in that case (every iterate of every orbit remains in the torus), and that focus on what we mean by 'good' recurrence without having to provide technical details about recurrence properties of certain sets of points. To explain this more fully, the notion of recurrence is fairly intuitive. We choose a region of the domain of the dynamical system, in this case the entire torus (i.e. the entire domain). Points are said to be recurrent with respect to this domain if, under time evolution, they repeatedly return to this domain (forever). In this case, they are trivially recurrent since any point in the torus remains in the torus under iteration by H 1 .
Next we consider the map also defined on the entire torus, as in the discussion of F 1 above. This map acts as a shear from bottom to top, where the strength of the shear increases monotonically from the left of S to the right of S, as indicated in figure 4(b). Now we compose F 1 and G 1 to obtain an LTM, H 1 , defined on the entire torus: which is the famous Arnold cat map (Arnold and Avez 1968). This is significant because it has been proven that the Arnold cat map is ergodic, (strongly) mixing, it is (very) chaotic in the sense that it has the Bernoulli property and it mixes (very) quickly in the sense that correlations decay super-exponentially in time. These properties have been proved in many different studies in the literature, and it is often difficult to track down the original references; see e.g. Sturman et al (2006). Now we contrast the strongly mixing Arnold cat map with a different (but markedly similar) LTM defined on the entire torus. We first consider the twist map F 2 : which is identical to F 1 , and the twist map G 2 :  The action of the maps F 2 and G 2 whose composition, H 2 (x, y) ≡ G 2 • F 2 (x, y), is such that every point is a periodic point of period 6.
with each defined on the entire torus, as in the discussion above. In contrast to G 1 , G 2 acts as a shear, from top to bottom, where the strength of the shear increases monotonically from the left of S to the right of S, as indicated in figure 5(b).
We construct an LTM defined on the entire torus by the composition of F 2 and G 2 : While H 2 may superficially appear very similar to H 1 (the only difference is the relative direction of the shears-co-rotating versus counter-rotating in fluid mixing terminology), the dynamics of H 2 is very different. One can actually compute the trajectory of any point explicitly (as well as its stability type). For an arbitrary point (x, y), iteration under H 2 gives the following: Hence, every point on the torus is a fixed point of the sixth iterate of H 2 , or a period-six point for H 2 . Linearization of the sixth iterate about any of its fixed points shows that the eigenvalues of the matrix associated with the linearization are identically one. Hence, there is no hyperbolic behaviour. The LTM H 2 does not mix at all. The obvious question is: Why does H 1 mix so well, and H 2 does not mix at all? The 'streamlines' of the two maps whose composition gives the LTM are orthogonal for both H 1 and H 2 . Therefore we have the best possible 'streamline crossing'. Both LTMs are constructed from maps that give 'transverse shears', where the magnitude of the shear rates for the horizontal twist maps (i.e. F i , i = 1, 2) and the vertical twist maps (i.e. G i , i = 1, 2) is identical for both LTMs, H i , i = 1, 2. It is therefore easy to verify that the Eulerian indicatorsβ,ᾱ andγ are identical for H 1 and H 2 . Hence, our Eulerian indicators clearly are not predictors of good or bad mixing for this case. Here we need to examine more closely what we mean by the Lagrangian characteristic of 'good recurrence to saddle-like overlap regions'. Since the LTMs are defined on the entire torus, the overlap regions (for the domains of F i and G i , i = 1, 2) are the entire torus.
The rather simple nature of the maps H 1 and H 2 allows us to obtain more insight into this process. Each of these maps is linear, and therefore the matrix associated with the linearization This figure illustrates the difference in behaviour between maps H 1 , which is mixing, and H 2 , which is periodic. Note that the Eulerian properties of the two maps are equivalent, in the sense that computingᾱ,β andγ for each gives identical results. The difference in mixing properties comes entirely from Lagrangian behaviour.
is the same for any point about which the map is linearized. We can compute the eigenvalues associated with the linearization, and they are 1 ± √ 2 for H 1 and 1 2 ± i √ 3 2 for H 2 . Hence, for H 1 , each point experiences saddle-like stretching and contraction in directions that are constant in space, and for H 2 , each point experiences rotation (with no stretching or contraction). This effect is illustrated in figure 6. Figure 6(a) depicts a short trajectory for H 1 , shown as black circles connected with dotted lines. We have endowed the initial point (bottom left) with a triangular wedge representing a set of direction vectors. (Strictly, these direction vectors belong to tangent space, but for the two-torus the tangent space at each point can be viewed as R 2 .) We operate on the wedge of directions with the Jacobian of H 1 at each iterate of H 1 itself. It can be seen that along the orbit of the initial hyperbolic point, the wedge narrows and lengthens (its area remains invariant as the map is area-preserving), until eventually all vectors in the initial wedge lie in the same direction (that of the unstable manifold of the initial point). Conversely, the same construction for H 2 is shown in figure 6(b). A period-6 point is endowed with a wedge of direction vectors, and at each iteration the wedge is rotated, until it finally returns to its original configuration as the trajectory returns to its initial point. Since H 1 and H 2 are linear, similar arguments apply to (almost) every point in T 2 . The marked difference in dynamical behaviours between H 1 and H 2 stems entirely from the Lagrangian features of the systems.
The linearity of H 1 and H 2 makes it easy to determine their stretching and contraction properties by examining their eigenvalues and eigenvectors. In the more general case where the LTMs are nonlinear there is a more general technique that generalizes the linear notion 13 of eigenvalues and eigenvectors. This is the notion of invariant cone fields (Alekseev 1969), which allows one to characterize the linear behaviour in the case where the Jacobian varies from point to point, i.e. we can study the behaviour of the Jacobian along a trajectory-which is clearly a Lagrangian property. The construction in figure 6(a) forms the skeleton of an invariant cone argument, since each subsequent iterate of the original wedge, or cone, is mapped into the previous one, while the lengths of vectors in the cone are increased.

Examples
In this section, we illustrate the use of the Eulerian indicators on four contrasting models. What these flows have in common is that they can all be considered blinking flows, and as such are examples of linked twist maps (Sturman et al 2006). Where they differ is in the mechanism by which the flow patterns are altered. We give some details of the systems themselves in the following sections. Although there are many variable parameters for each system, and much research has been undertaken to describe various properties under the alteration of various parameters, we fix all parameters as constant, except that in each case we retain a single parameter that serves to alter the streamline structure and thus the Eulerian properties.
For each model we select four parameter values producing different types of mixing behaviour. For each such parameter, we show the superimposed streamline patterns (which serves to give a visual impression of the extent of 'streamline crossing' in the flow) and two plots indicative of the Lagrangian mixing properties. The first of these is a Poincaré map: for carefully chosen initial conditions, we plot many iterates (typically 10 5 ), and the resulting plot makes distinct the decomposition of the domain into chaotic seas and unmixed islands. The second illustrates the effect of the mixing process on two initially segregated sets of initial points. For each model we seed the process by dividing the domain into two halves, left and right. The left half forms a set L containing 10 6 points coloured red, whereas the right half forms a set R containing 10 6 points coloured green. The plot we show is the images of L and R after ten iterates of the Poincaré map. We note here that although the mixing properties of the system itself will be dependent on the initial sets chosen, the combination of the Poincaré map and the mixing plot gives a reasonable summary of the mixing properties of the system. For each parameter value, we also show three contour plots, ofα,β andγ . These are intended to illustrate how the Eulerian properties of the blinking flows change as the parameter is altered. In each of the contour plots, darker colours indicate greater values ofα,β andγ . Scales are chosen to be consistent within each row of figures and to best display the variations present. In the computation ofβ, we typically take the parameter δ i = 0.01 (i = 1, 2), although the calculation is robust to changes to δ i .
For further comparison of the mixing properties of each system with the Eulerian indicators, we provide figures illustrating the behaviour of each Eulerian indicator. There are many different ways to quantify the mixing performance of a mixing device (see, for example, Funakoshi (2008) for a review). We opt, for simplicity, to compute the (normalized) variance of concentration V. Specifically, we take the tenth iterate of the sets L and R defined above and compute the variance of the concentration of the resulting mixture. This is done by a simple numerical scheme: place a grid over the domain (typically with around 10 4 gridpoints), and estimate the concentration at each gridpoint by counting the number of red and green points inside a circle of (fixed, small) radius r . The concentration c is then given by c = no. of red points no. of red points + no. of green points at each gridpoint. Since initially the domain was half red and half green, a perfect mixture would give c = 0.5 everywhere. Perfectly segregrated it would lead to c = 0 or 1. We then take the variance of this value over all the gridpoints and rescale by multiplying by 4 (the product of the average concentrations of red and green initially). This quantity then corresponds exactly to intensity of segregation (see, for example, Danckwerts (1953a)). This should give a variance of 0 for a perfect mixture (every concentration is 0.5, equal to the mean), and 1 for perfect segregation. Of course, the value of the variance in this scheme depends on the choice of r , that is, V = V(r ). If r is too small, there are a few red and green points in each circle, and the value for concentration is unlikely to be accurate; if r is too large, we sample from too large an area, and differences in concentration are not seen. We have typically chosen r = 0.1, which gives roughly 100 red and green points in each circle while still being small compared to the domain. Numerically, the results for variance are reasonably robust to small changes in r . We also note that the value of V is dependent on the choices of initially segregated sets. We have opted for a fairly generic choice (dividing the domain into left and right halves), but have checked that other choices typically seem to give similar results. However, it should be noted that armed with the information in a Poincaré plot, it is always possible to choose initial sets to coincide with elliptic islands, and therefore skew the resulting value for V.
We compute the variance of concentration for a range of parameter values for each system, and compare this directly with the values of the Eulerian indicatorsᾱ,β andγ in turn. Note that we do not expect any quantitative consensus, merely a broad qualitative agreement.

The blinking vortex flow
The well-known time-dependent flow known as the blinking vortex flow and introduced in Aref (1984) is the most obvious mixing device to take the form of a linked twist map. To specify this flow we consider a circular container of radius R centred at the origin with a point vortex of strength placed, without loss of generality, at (a, 0) in Cartesian coordinates. This has complex potential given by the equation (Aref 1984) and hence the equations of motion are given bẏ The essence of the blinking vortex is to alternate the operation of two such vortices with different locations for the vortices. Here we take a pair of vortices centred at (±a, 0). We discretize this in time by considering the map T = T −a • T a , where T a (respectively T −a ) is the time-1 map of the vortex flow defined above with the vortex located at (a, 0) (respectively (−a, 0)). Throughout this section we fix the radius at R = 2, and the strength of the vortices at = 1.0 (hence, co-rotating vortices). The remaining parameter a is the only one we vary. The streamlines corresponding to both T a and T −a are circles, and chaotic advection is possible because the superposition of these streamlines features streamline crossing. Figures 7(a)-(d) show this streamline crossing for the four values a = 0.4, 0.8, 1.2 and 1.6. Figures 7(e)-(h) show Poincaré plots corresponding to these streamline configurations, whereas figures 7(i)-(l) illustrate the mixing of the two initially segregated sets L and R.
We quantify the streamline configuration and illustrate the Eulerian quantities defined in section 2. Firstly, for the same parameter values used previously, we plot the angleα between streamlines at each point in the circular domain in figures 7(m)-(p). In each figure, darker colours correspond to angles closer to π/2. Note that each of these figures has a horizontal line of tangentiality along y = 0-this is forced by the symmetry of the patterns. These figures suggest that the degree of transversality increases as a increases (as the vortices move to the boundaries of the domain). This appears to contradict figures 7(i)-(l), which clearly show the mixing quality decreasing as a increases. However, as discussed in the previous section, it is important for efficient mixing that transverse streamlines should be linked with strong shear rates. Secondly, we show in figures 7(q)-(t) the productβ = β 1 β 2 of the shears due to each vortex. Since the shear rate due to each vortex increases towards the centre of that vortex, we see the product of the two shear rates decrease as the vortices are moved away from each other. Finally, figures 7(u)-(x) show the value ofγ =αβ, the product of transversality and shear rate. As before, darker colours indicate greater values ofγ .
We now compare the performance of the Eulerian indicators with a standard diagnostic for mixing quality, the variance of the concentration V as discussed above. This comparison is shown in figure 8. In each figure, we plot V as a function of the parameter a on the left-hand axis in red. Each Eulerian indicator is plotted on the right-hand axis in red. Note that the axis scales are different-we do not expect quantitative agreement between the curves. Note also that the right-hand axes are inverted-larger values of the Eulerian indicators correspond to more efficient mixing, and in turn to smaller values of V.
In figure 8(a), we compareᾱ, the average transversality, to V. Clearly there is little agreement; in fact, there is a negative correlation. However,β andγ broadly coincide with the curve for V. For this ubiquitous and fundamental flow, two Eulerian indicators predict mixing behaviour well.

A kinematic model of a blinking flow containing a separatrix
In this section, we construct a kinematic model for a linked twist map with a separatrix that models many of the flows described in section 1. The idea is simple. In a rectangular domain, we first construct a flow consisting of two circulation cells divided by a vertical streamline that acts as a separatrix (flow 1). We then construct a second flow, on the same domain, of the same type as the first, but with the separatrix in a different location (flow 2). Fluid particles are advected by 'blinking' between flow 1 and 2 alternately. More precisely, a fluid particle is advected for time t 1 under flow 1, then for time t 2 under flow 2 and so on.
We first consider the 'standard' streamfunction for a cellular flow, 1 (x, y) = cos(π x/2)cos(π y/2). The velocity field defined by this streamfunction will have 'free-slip' boundary conditions on the boundary of the unit cell, S = [−1, 1] × [−1, 1]. In order to model Mixing properties and kinematic behaviour of the blinking vortex for four contrasting parameter values. Columns 1-4 are for a = 0.4, 0.8, 1.2 and 1.6, respectively. The first row shows the two superimposed streamline patterns for the blinking flow. The second shows a Poincaré section for initial conditions chosen to reveal the structure of the dynamics. The third row shows the images after ten iterations of the map T of two initially segregated sets of initial conditions. The fourth, fifth and sixth rows show contour plots of the anglê α between streamlines, the product of shearsβ and the product of transversality and shearsγ , respectively. Darker colours indicate greater values plotted. the 3D, spatially periodic steady channel flows whose cross-sectional flows have a separatrix, we need to consider a flow having the no-slip condition on the boundary. Towards this end, we modify the standard cellular flow streamfunction as follows: 2 (x, y) = − cos(cos(π x/2) cos(π y/2)) = − cos( 1 (x, y)), where it is still considered on the unit cell S = [−1, 1] × [−1, 1]. The velocity field corresponding to this streamfunction is given by where u x (x, y) = − π 2 cos π x 2 sin π y 2 sin cos π x 2 cos π y u y (x, y) = π 2 sin π x 2 cos π y 2 sin cos π x 2 cos π y 2 .
This velocity field is zero on the boundaries of S, and has a stagnation point at (x, y) = (0, 0). Using this cellular flow, we construct a flow in the rectangular domain Q = [−1, 2] × [−1, 1] having the desired properties. It consists of two cellular flow patterns, divided by a separatrix consisting of a vertical streamline at x = s. The velocity field corresponding to the desired flow is given by where the parameters are such that the cellular flows defined by L and R have stagnation points of centre-type linearized stability at x = b 1 , y = 0 and x = b 2 , y = 0, respectively, in a domain of half-width c 1 and c 2 , respectively. The domain being defined as Q ≡ [−1, 2] × [−1, 1] requires that these parameters must satisfy b 1 + 1 = c 1 and 2 − b 2 = c 2 . Requiring that the cellular flows are divided by a separatrix at x = s implies that s = 2 − 2c 2 and s = 2c 1 − 1. Note that particles on the streamlines defined by L move in a counter-clockwise sense, whereas particles on the streamlines defined by R move in a clockwise sense. The separatrix x = s is a line of stagnation points. The location of the separatrix is controlled by the choice of the parameters b i , c i , i = 1, 2, but b i and c i are related as given above. Looking at this another way, the values of b i and c i depend only on the position of the separatrix, so to reflect this we introduce a single parameter such that c 1 = (3/4) − ( /2) and c 2 = (3/4) + ( /2). Thus when = 0, the separatrix lies in the middle of the rectangle Q, at x = 0.5. Increasing has the effect of moving the separatrix to the left, until it meets the boundary x = −1 at = 1.5, whereas decreasing from zero has the effect of moving the separatrix to the right, until it meets the wall x = +2 at = −1.5. Figures 9(a) and (b) illustrated two cases. To explicitly denote this dependence of the location of the separatrix on the parameter , we write L(x, y) = L (x, y) and R(x, y) = R (x, y).
With this flow it is easy to construct a blinking flow by considering a flow that alternates between two different flows corresponding to two different choices of , with each flow acting for times t 1 and t 2 , respectively. For simplicity we will take t 1 = t 2 = 1 since our goal is to illustrate the general principles underlying 'good mixing' in blinking flows with separatrices, but in optimizing the mixing properties of specific flows it may be beneficial to make specific choices of t 1 and t 2 .
More formally, we consider the time-1 map of the flow, given by The blinking flow is then given by the composition of two such maps, T 1 , 2 = T − 2 • T 1 (x, y). We choose four representative values of 1,2 and show the streamline crossing structure in figures 10(a)-(d). These show the separatrices moving from the centre of Q to the edges of the Figure 10. Mixing properties and kinematic behaviour of the separatrix model for four contrasting parameter values. Columns 1-4 are for 1,2 = 0.1, 0.5, 0.8 and 1.2, respectively. The first row shows the two superimposed streamline patterns for the blinking flow. The second shows a Poincaré section for initial conditions chosen to reveal the structure of the dynamics. The third row shows the images after ten iterations of T 1 , 2 of two initially segregated sets of initial conditions. The fourth, fifth and sixth rows show contour plots ofα,β andγ , respectively. Darker colours indicate greater values plotted.
boundary-specifically we take 1,2 = 0.1, 0.5, 0.8 and 1.2. Figures 10(e)-(h) show Poincaré maps for these parameter values, whereas figures 10(i)-(l) show the plots of mixing of initially segregated sets. For 1,2 = 0 and 1,2 = 1.5, the blinking system is integrable-there is no mixing and no streamline crossing. It is clear that near these extremes (in the first and fourth columns of figure 10) mixing is poor, and that many streamline crossings are near-tangential. The latter point is made clear in figures 10(m)-(p), which show the angle between streamlinesα for The Eulerian indicators themselves involve averaging each of these measures. In figure 11, we show how these indicators compare to Lagrangian diagnostics commonly used to quantify the extent of mixing-the variance of concentration and Poincaré maps. Figure 11(a) shows the result of comparing the Eulerian indicatorᾱ with the Lagrangian measure of the variance of the concentration of the mixture. The initial condition for the mixture is the same as that for the simulations of figures 10(i)-(l) (except that we use 10 4 , rather than 10 6 initial conditions for computational speed), and the variance is computed after ten iterations. At 1,2 = 0 and 1.5 the flow is integrable, so the variance should be one. However, numerical artefacts associated with the existence of a boundary and a finite grid size introduce a slight error. Nevertheless, both the Eulerian indicator and the Lagrangian variance indicate that the most effective configuration for mixing is for 1,2 ≈ 0.6. The agreement between these two curves is clearly very good, which indicates that our Eulerian indicator does indeed predict the value of 1,2 for which this class of flows will result in the best mixing. However, it is clear that while the Eulerian indicatorᾱ varies smoothly with 1,2 , the variance curve has more structure, which is not reflected by the Eulerian indicator. Figure 11(b) shows the result of comparing the Eulerian indicatorβ with the Lagrangian measure of the variance of the concentration of the mixture as in figure 11(a). Here there is little agreement between the red and the green lines, indicating thatβ alone is a poor indicator in this case. However, figure 11(c) compares the variance with the Eulerian indicatorγ , and now the combination of transversality and shear mimics the variance curve fairly closely. Recall that in each graph of figure 11, the right-hand axis is inverted, so that large values of the Eulerian indicator (strong shear and/or transversality) correspond with small values of variance (good mixing).

Pulsed source-sink system
A contrasting mixer design, also in the form of a blinking flow, is commonly used in hybridization of DNA molecules. In a circular boundary of radius R consider a point source and point sink of equal strength, so that fluid particles travel along curved streamlines from the source to the sink. On arriving at the sink a particle is reinjected into the flow from the source according to a prescribed reinjection procedure. If the reinjection procedure is chosen so that a particle leaves the source along the same streamline as it entered the sink, then no chaotic motion occurs-there is no crossing of streamlines. A mechanism for creating crossing of streamlines is to introduce a second source-sink pair into the system. In this situation, the system is rendered time-dependent by switching one source-sink pair on for a time t 1 while the second source-sink pair is switched off, and then switching the first source-sink pair off and the second pair on for a time t 2 . This produces two different streamline patterns, which are effectively superimposed as the source-sink pairs are alternately switched on and off following this procedure. Many variants on this idea have been proposed and studied; see for example McQuain et al (2004), Stremler and Cola (2006), Cola et al (2006), Raynal et al (2004) and Evans et al (1997). The ergodic mixing properties of this class of mixer were investigated as linked twist maps in Hertzsch et al (2007), in which an instantaneous reinjection is assumed, and the system modelled as a map on deformation of the torus.
In Hertzsch et al (2007) the authors discuss how the choice of pumping times t 1 and t 2 affects the quality of mixing. Varying these parameters has no effect on the Eulerian features of the system, so we fix t 1 = t 2 = 0.1 and select a single parameter that affects the geometry of the streamlines. We assume a circular boundary of radius R = 2, and consider source-sink pairs S 1 and S 2 given by where S + (x, y) represents a source at (x, y), and S − (x, y), represents a sink at (x, y). We vary the parameter v ∈ [0, −1.5] to create different streamline patterns. The top row of figure (2007), but is due to the fact that in this paper we choose t 1 = t 2 = 0.1, and therefore do not consider the effect of the variation of the pumping times. Our aim in this paper is not necessarily to select parameters to optimize mixing, but rather to show comparisons that illustrate effectively the use of the Eulerian indicators. Note also that mixing in this type of system is not only due to the stretching and folding of chaotic advection. The reinjection mechanism also creates the possibility of mixing via 'chopping and shuffling' (Ashwin et al 2002). Figures 12(m)-(x) show contour plots ofα,β andγ for the same parameter values. As before, we compare the performance of the Eulerian indicators with the variance of the concentration of two initially separated blobs. Here we initially seed the system with 10 4 red initial conditions in the left-hand half of the circular domain, and 10 4 green initial conditions in the right-hand half and compute the variance of concentration V after 20 applications of both source-sink pairs. The variance V changes as v is varied and is shown as the red curves in figure 13, plotted on the left-hand axis. As before, we show the comparison with, in figure 13(a), the indicatorᾱ, in figure 13(b), the indicatorβ and, in figure 13(c), the indicatorγ . In this situation, the transversality measureᾱ alone models the variance poorly, while including the In figure 14, we show the streamlines, Poincaré sections, mixing behaviour 4 and Eulerian indicators for the parameter values ψ = 0.25, 1.0, π/2 and 2.2. It is apparent from the Poincaré sections and mixing plots that the symmetric configuration ψ = π/2 is not necessarily the most efficient for this mixer, a conclusion also reached in Mizuno and Funakoshi (2002). Figure 15 shows the performance of the Eulerian indicatorsᾱ,β andγ .

Conclusions and outlook
In this paper, we introduced three Eulerian indicators that indicate the extent of (averaged) 'streamline crossing' (ᾱ), (averaged) 'shear rate' (β) and the (averaged) productγ . These were applied to four different flows-the blinking vortex flow, the 'separatrix model', the pulsed source-sink pair flow and the partitioned pipe mixer. All parameters associated with each flow were fixed, except for one, and we considered the nature of the mixing properties of the individual flows as the chosen parameter was varied. The variance of the concentration and Poincaré maps were taken as Lagrangian measures for quantifying mixing. We found that for each flow the Eulerian indicatorγ broadly agreed with the Lagrangian measures of mixing in predicting the value of the parameter for which the best mixing occurred. For two flows, the 'separatrix model' and the partitioned pipe mixer,ᾱ also served to predict the parameter value for which the best mixing occurred, butβ was not a useful predictor of good mixing for these flows. Contrastingly, for the other two flows, the blinking vortex flow and the pulsed source-sink pair flow,β also served to predict the parameter value for which the best mixing occurred, but α was not a useful predictor of good mixing for these flows. This shows that in flows that can be modelled as 'blinking flows' (i.e. by some mechanism the flow is made to alternate between two different flow patterns) it is 'transversely oriented shear' that generates good mixing.
We end with a few remarks and suggestions for future study.
1. A consideration of the manner of comparison between the Lagrangian measures of mixing (i.e. variance of the concentration and Poincaré maps) shows why Eulerian indicators offer the promise of constructing optimal mixers. For each flow it is relatively easy to compute the value of the Eulerian indicators for a large range of parameter values. The Lagrangian measures of mixing require computation of large numbers of trajectories and are computationally time consuming. 2. We want to emphasize that a justification of the Eulerian indicator that we develop as an indicator of 'optimal mixing' does not follow from any of the known LTM theorems. Rather, it is merely suggested by the kinematical features of the linked twist maps for which there are rigorous theorems. We hope that our work serves as motivation for the development of new theorems along the lines of the original LTM theorems. 3. As discussed in section 2.4, the issue of 'co-rotation' and 'counter-rotation' can have a profound effect on mixing properties. A detailed discussion of the subtleties involved in this distinction for LTMs is given in chapter 2 of Sturman et al (2006). The separatrix model discussed in section 3.2 simultaneously exhibited both co-and counter-rotation. Mixing properties and kinematic behaviour of the partitioned pipe mixer for four contrasting parameter values. Columns 1-4 are for ψ = 0.25, 1.1, π/2 and 2.2, respectively. The first row shows the two superimposed streamline patterns for the blinking flow. The second shows Poincaré sections for certain initial conditions. The third shows the tenth iterate of initially segregated sets. The fourth, fifth and sixth rows show contour plots ofα,β andγ , respectively. Darker colours indicate greater values. Also, the known LTM theorems all require 'monotonicity' with respect to the speed along streamlines. In the separatrix model, there are regions of the flow where monotonicity breaks down. Nevertheless, the known LTM theorems should provide a framework for developing a rigorous mathematical theory for this new type of LTM. 4. The discussion of the co-rotating and counter-rotating cases in chapter 2 of Sturman et al (2006) brings to light the issue of the importance of the domain on which the flow is defined. Nevertheless, for the type of flows considered in this paper having solid boundaries (e.g. a cavity or duct), our Eulerian indicator provides an excellent predictor of the equality of mixing based on the geometrical characteristics of the overlapping flow patterns. Consequently, it provides a tool for constructing optimal mixers of the type used in the microfluidic applications listed in section 1. We remark that optimization of mixing in settings where the LTM theorems is applicable is considered in Sturman et al (2006), chapter 9. 5. While we have only been considering blinking flows in this note, there are several ways in which our Eulerian indicator could be applied to continuous, spatially periodic flows in 3D. In particular, we could compute the Eulerian indicator for the overlapping flow patterns of the cross-sectional flows at the beginning and end of a periodic element. Or we could easily compute an 'average transversality' for a system whose streamline pattern changes continuously and periodically in space with time. 6. It would be interesting to check the performance of the Eulerian indicator in the presence of diffusion since the consideration of the Eulerian indicators does not take into account diffusion. 7. The example flows that we considered could be considered as 'blinking' between two different flow patterns. However, the Eulerian indicator approach could easily be extended to flows that blink between n flow patterns (n > 2). Some examples of LTMs for compositions of three maps can be found in Sturman et al (2006).