Chiral Separation by Flows: The Role of Flow Symmetry and Dimensionality

Separation of enantiomers by flows is a promising chiral resolution method using cost-effective microfluidics. Notwithstanding a number of experimental and numerical studies, a fundamental understanding still remains elusive, and an important question as to whether it is possible to specify common physical properties of flows that induce separation has not been addressed. Here, we study the separation of rigid chiral objects of an arbitrary shape induced by a linear flow field at low Reynolds numbers. Based on a symmetry property under parity inversion, we show that the rate-of-strain field is essential to drift the objects in opposite directions according to chirality. From eigenmode analysis, we also derive an analytic expression for the separation conditions which shows that the flow field should be quasi-two-dimensional for the precise and efficient resolutions of microscopic enantiomers. We demonstrate this prediction by Langevin dynamics simulations with hydrodynamic interactions fully implemented. Finally, we discuss the practical feasibility of the linear flow analysis, considering separations by a vortex flow or an extensional flow under a confining potential.

Scientific RepoRts | 6:35144 | DOI: 10.1038/srep35144 (1) tt tr rt rr te re t r where ≡ U U r ( ) denotes the ambient fluid velocity at  r , and Ω = ∇ × U r ( ) 1 2 is half the fluid vorticity. Here the ambient flow is approximated as a linear Stokes flow = + U r U r J ( ) (0) with a constant velocity gradient tensor = ∇ U J 21 . The rate-of-strain field E is the symmetric part of J, given as E = (J + J T )/2 with the transpose operator T. The mobility tensors, μ's, are related to the strengths of thermal noises, ξ  's, through the fluctuation dissipation theorem ξ ξ µ 〈 〉= k T 2 i m j n B ij mn , with m, n = t, r and k B T being the thermal energy at temperature T. Here the superscripts, t and r, stand for 'translational' and 'rotational' , respectively; μ's couple two degrees of freedom corresponding to its superscripts. The third-rank resistance tensors, ζ's, determine hydrodynamic friction force and couple the translational and rotational motions of the object with the rate-of-strain field (for which we assign 'e' as a superscript symbol, and see also index notation of Carrasco et al. 23 ).
The tensor product in Eq. (1) is defined as (ζ : E) i ≡ ∑ j, k ζ ijk E jk . In calculating μ and ζ, the no-slip condition on object surface is imposed, and they depend on the position vectors of the surface elements relative to O. The mobility and resistance tensors are thus independent of the choice of the origin, and they are functions only of the geometry of the particle such as object orientation φ and handedness α. Therefore, the chirality-dependent drift can be induced only by the rate-of-strain field E. We write it in shorthand notation, E tt te tr re with chirality index α = R or L and refer to it as drift velocity.
Symmetry under parity inversion. Using a symmetry property under parity-inversion operation, we first show that the chiral separation (objects move in opposite directions according to handedness indeed) occurs by the parity-even rate-of-strain field. In particular, this is proven even with considering the rotational motions of object orientations. Using parity-inversion () operation about the point O, a linear flow is decomposed into two parts; one is a parity-odd flow whose direction is changed by  operation, and the other is a parity-even flow that remains invariant under  operation. In Fig. 1, translation motions of a chiral molecule and its chiral partner under parity-odd and parity-even flows are illustrated, respectively. Shown on the left panel of Fig. 1(a) is a left-handed (L) object in a parity-odd flow, which is supposed to drift to the left by the flow. Applying  operation on this system, the flow direction is inverted, and the object handedness as well as the drift direction are reversed (middle panel). If the flow field becomes again inverted (right panel), the original flow field is recovered and the right-handed (R) object then moves to the left because of the linearity of relations between the object velocity and the flow fields. One clearly sees that the drift motions of the chiral pair are identical under the parity-odd flow, leading to no chirality resolutions. On the contrary, if  operation is applied on an object of L in a parity-even flow ( Fig. 1(b)), the object handedness is converted into R and the drift direction is changed as well, while the flow itself remains invariant. This enables us to have the opposite motions of the chiral pair under the same ambient flow and thus demonstrates that the parity-even flow field is essential for chiral separation.
According to Eq. (1), the translational velocity  v of an object includes two flow components U and E. It is obvious that U is parity-odd while E is parity-even as J remains invariant under  operation. Therefore, from the analysis above, we find that the drift velocity of a left-handed object is opposite to its chiral partner under a given rate-of-strain field: where ϕ′ denotes parity-inverted orientation of φ . It is to be noted that the angular motion remains invariant under  operation, and thus the probability distributions of orientation are symmetric, , at any time t (see the Supplemental Information for details). Accordingly, the orientation-averaged drift velocity of the left-handed object is given by where the negative sign of the average drift velocity implies a possibility of chiral separation. In a similar way, the separation direction can be specified if a linear flow field has a mirror-reflection symmetry. Shear flow, γ = Û r yx ( ) with shear rate γ  , belongs to the case, if choosing the xy-plane as a mirror-symmetry plane. Regard a relation, where the double prime indicates quantities transformed by mirror-reflection () operation. Here, α v E, and α ⊥ v E, are, respectively, the parallel and perpendicular components of a drift velocity α  v E to a chosen mirror-reflection plane. For a mirror-symmetric flow field (E″ = E), taking average over object orientations for , one has a relation, Combining with Eq. (4), one finds that 〈 〉 α Φ v E, = 0, and the separation of objects with opposite chirality can occur along the direction perpendicular to the mirror-symmetry plane, as indeed in the case of the shear flow. We again note that Eqs (4) and (6)   v E to be well separated on macroscopic scales, at least, larger than experimental resolutions of microfluidic devices. Below we argue that this is achieved with quasi-two-dimensional flow fields.
Eigenmode analysis. It is more intuitive to examine the equation for  v for the case of diagonalizable J (for a non-diagonalizable case the analysis can be performed in a similar way, leading to the qualitatively same conclusions; see the Supplemental Information for details). The equation for  v in Eq. (1) can be written in the diagonalizing basis of J, and one of its components along the direction of the q-th eigenvector of J reads as, where λ q is the q-th eigenvalue of J, and quantities with the subscript λ q symbolize the q-th component of transformed vectors = λ X X S ( ) q q with (SJS −1 ) q,q′ = λ q δ q,q′ and X = r, v E , ξ t . For notational simplicity, we shall drop the eigenvalue index q and chirality index α hereafter. Without loss of generality, we have set = U(0) 0. The formal solution of Eq. (7) is given by Here we define the displacements resulting from the drift motion, , and from the thermal noise . Taking average over the thermal noises and initial positions, we obtain the average displacement with the average of initial positions located at origin. Mean separation distance between particles having opposite chirality is given by , and separation precision can be quantified by  The chiral separation is achieved if the quantifiers satisfy following relations: The first condition requires that the mean separation distance t ( )  should be much larger than the dispersion Σ (t) for efficient separation. The second condition states that the separation distance t ( )  ought to be much larger than the initial width of distribution σ 0 which roughly amounts to experimental resolution such as the distance between recovery outlets of devices for the separated particles.
It is clear that an exact analysis of the separation conditions (10) requires to obtain the precise form of v E,λ (t) which is highly nonlinear and has no known solution to the best of our knowledge. Notwithstanding the difficulty, one can still put forward a reasonable analysis: Suppose that there exists a finite maximum drift velocity v m of v E,λ (t), viz. v E,λ (t) ≤ v m , which is reasonable because the drift velocity as given in Eq. (2) is determined by the product of the bounded quantities. It is also assumed that up to the leading order, the maximum drift velocity v m is proportional to the magnitude of flow field given as λ ∑ J i i 3 2 , and for a fixed |J|, only weakly depends on an individual λ i . This assumption may be supported by an observation that as increasing the flow gradient tensor J by a factor of a certain constant, the resulting drift velocity will be increased proportionally. Lacking in mathematical rigour, this ansatz is effective to extract an essential flow property to induce chiral separation, as evidenced m t On the other hand, the positional dispersion can be written as σ σ σ λ where σ(t) determines the dispersion when v E,λ (t) is purely deterministic to annihilate fluctuations of X E,λ (t).
Obviously, an inequality, Σ (t) ≥ σ(t), follows, which together with Eq. (10) then constitutes a necessary condition for chiral separation, 0 which is given in a greatly simplified form to allow an analytic approach and helps to extract essential flow factor inducing chiral separation. It should also be mentioned that in high Péclet number regime the drift velocity as a function of object orientation remains roughly constant, as suggested by Marcos et al. 11 . For the case, v m can be interpreted as the constant drift velocity and hence, conditions in Eq. (14) are equivalent to Eq. (10). In the following, we proceed our analysis of the separation conditions, (14) with d(t) and σ(t) determined by Eq. (12) and Eq. (13), respectively.
Short time behavior. We examine the behavior of the separation precision for t ≪ t λ ≡ |λ| −1 , where t λ is the saturation time scale at which the separation precision with a non-zero λ approaches a constant value. In Fig. 2, we illustrate the characteristic behaviors of d/σ depending on λ. For λ = 0, d/σ grows unboundedly with time t, while for λ ≠ 0, it approaches a constant value after a time t λ . That is, in the case of λ ≠ 0, the feasibility of d/σ ≫ 1 depends on the system parameters such as initial dispersion and diffusion constant, even at large t. For t ≪ t λ , d is monotonically increasing function of time, and irrespective of the sign of the eigenvalue λ, d/σ is simplified as   . From the separation conditions, Eq. (14), and the equation of the separation precision, Eq. (15), one can readily define the separation time scale t s at which both d/σ and d/σ 0 become of the order of unity (see Fig. 2 s D and therefore, an efficient and precise chiral resolution (i.e., d/σ ≫ 1 and d/σ 0 ≫ 1) occurs when This constitutes the condition of macroscopic chiral separation for short times, i.e., t ≪ t λ .
Long time behavior. In the long time regime of t ≫ t λ , if the eigenvalue is not zero (λ ≠ 0), d/σ is saturated to a certain value as This can also be intuitively understood from Fig. 2, i.e., only when t λ ≫ t s , the saturated value of d/σ can be much greater than one. Then, the condition of . On the other hand, when λ < 0, the first condition of Eq. (14) leads to t λ ≫ t D , and the second condition does to t λ ≫ t σ .
Remarkably, we find that in both of short-time and long-time regime, the macroscopic chiral separations are dictated by a single criterion, irrespectively of sign of λ: s s D which means that the eigenvalue of the velocity gradient tensor should be significantly smaller than the inverse of the separation time scale. We note that the eigenvalue index q has been omitted for notation simplicity. Therefore, Eq. (20) should be interpreted as conditions required to be satisfied by respective λ q in order to obtain the separation along the corresponding coordinate λ r q . In other words, the chiral separation may occur if at least one of the eigenvalues satisfies the condition.
Dimensional analysis and estimates of parameters. Let us now express the separation condition, Eq. (20), in terms of physical parameters such as the linear size of the object  and magnitude of the flow velocity gradient V defined through J (for example, see Eqs (22) and (23) of the next section). One may define V also as ∑ = V J i j ij , 1 3 2 , where the proportionality constant is of the order of unity, and it is irrelevant to the present analysis. We introduce dimensionless parameters where δ is initial dispersion σ 0 in units of , and it seems reasonable to assume δ ≫ 1 in most cases of practical interest of small particles. ε is the amplitude of λ relative to the magnitude of the flow velocity gradient, which can characterize the dimensionality of the flow field; if ε = 0, the corresponding flow becomes a two-dimensional flow. c measures the chiral-dependent drift velocity relative to a variation of ambient flow velocity over .
. As a result, the separation condition of Eq. (20) can be written in dimensionless form as . Among the parameters, c affects both of t σ and t D , and its estimation is important for numerical evaluations of Eq. (21). The exact value of c, of course, relies on system details such as object shape and chirality, and it should be a tremendous task to obtain the general expression of c analytically. However, as we show in the Supplemental Information, a possible upper bound of c can be envisioned, which turns out to be of the order of 10 −2 . In addition, we numerically evaluate c for different objects and flow patterns considered in simulations, which is indeed found to be small as consistent with the proposed upper bound (see the next section of Simulation results). Combining these facts, we can take a conservative bound of |ε| for an efficient separation as  ε − ⪅ (10 ) 3 for various flow strengths and/or object sizes, even though a rather unrealistically narrow initial distribution (δ ~ 1, i.e., σ  0 ) is assumed. Considering the measurement and control accuracy of current microfluidic devices (to our knowledge, of the order of 0.1% at best), one might view this range of ε to be synonymous for a singular flow. It implies that the chiral separation is possible only by quasi-two-dimensional flows described by a velocity gradient tensor J with a vanishingly small eigenvalue.
Before going further, we explicitly mention the validity range of our theory, especially, in terms of the relevant object size . Our theoretic formulation assumes the low Reynolds and high Péclet number conditions, which give the range of appropriate object size as η η ρ  v V and ρ is the fluid density. For a water at room temperature, µ µ . 0 01 m 1 m for V ~ 10 6 /s, and µ µ . 0 1 m 100 m for V ~ 10 2 /s. The present theory is based on the assumption of linear flow field which is hard to be realized in most cases over extended length and time scales. For example, the presence of hydrodynamic boundaries arising from confining walls of microfluidic devices leads to nonlinear flows, which might yield nontrivial effects on the separation. However, analysis on general nonlinear flow fields is beyond the scope of the present study, and the relevance of our linear flow analysis will be discussed later in more detail in section of Possible applications.
Simulation results. In order to demonstrate the arguments, we perform Langevin dynamics simulations by integrating Eq. (1) with explicitly taking into account full hydrodynamic interactions among object elements at the level of the Rotne-Prager tensor 23 . As typical examples of enantiomers, helix and tetrahedral structures are concerned. We decompose the structures into arrays of N closely packed beads with radius a. The decomposition allows us to calculate their grand mobility tensors which in turn yield mobility tensors for a rigid body motion of Eq. (1), according to the conversion equations, Eq. (19) to (22) in ref. 23 (see the Supplemental Information for details). In simulations, we rescale all lengths by the bead radius a, giving linear size  of the helix and tetrahedral structure as 225a and 254a, respectively. Characteristic time scale is τ = 6πηNa 3 /k B T. The distribution of initial positions of objects is assumed to be a Gaussian centered at the origin with an initial width σ 0 .
We consider two different types of flow fields, one with a diagonalizable velocity gradient tensor J A and the other with a non-diagonalizable tensor J B , parameterized by a dimensionless variable ε. First, flow-A is described by the velocity gradient tensor, A A which has distinct eigenvalues, λ 3 /V A = ε and λ ε ε 2 2 . If ε = 0, the eigenvalue λ 3 vanishes and the corresponding J A describes a two-dimensional flow. Flow-B is represented by where a finite ε describes a deviation from the shear flow (ε = 0), and eigenvalues are degenerated, λ 3 /V B = 2ε and λ 1,2 /V B = − ε. Both flows are incompressible, i.e., J is traceless, and have a reflection symmetry about xy-plane. According to our symmetry argument, the chirality-dependent drift is expected to occur along the z-direction which is the eigenmode direction of λ 3 . Dimensionless flow velocities, for i = A, B, are set to 30 in order to realize high Pe; Pe ~ 3.2 × 10 4 for helix and Pe ~ 3.9 × 10 4 for tetrahedral. The magnitude of flow velocity gradient, V, can be defined in a basis independent way as = ∑ V EE i j ij ij 0 , . However, the current expression of V equals V 0 up to a prefactor of the order of one and does not lead to any qualitative difference in results. For a wide range of flow velocity V and object size , c is found to be small as − ⪅ c 10 3 for the helix and − ⪅ c 10 2 for the tetrahedral, in accord with the estimate on the possible maximum value of c (see the Supplemental Information). The previous criterion of Eq. (21) predicts that for separation to occur, the upper bound of |ε| is given as ε − ⪅ O(10 ) 3 even when a very narrow initial distribution of molecular size is assumed (δ ~ 1). We test this prediction through numerical simulations, as varying the value of ε of the flows considered above. Figure 3 explicitly shows the flow fields at different flow parameters, ε = 0 and ± 0.01, for the flow-B (a-c), and the time evolutions of corresponding probability distributions of R or L helices, along the z-direction, obtained from the simulations (d-f). Note that, in (b-c), the arrows indicating difference fields are magnified by 50 times for clear visibility. As consistent with Eq. (21), a small but finite value of ε leads to qualitatively different behaviors, despite the apparent similarity to the separable shear flow with ε = 0. For either small positive or negative ε, the probability distributions continue to substantially overlap even at t ≥ t λ , while for ε = 0, the two distributions are well discriminated and d increases with time, enabling complete separation. Now we propose to quantify the degree of the separation by the Jensen-Shannon divergence (JSD): where p α (z) represents the probability distribution to find an object with chirality α = L, R at a projected position z along a chosen axis of observation. Unlike d/σ, JSD is bounded as 0 ≤ JSD ≤ ln2, and its value depends on overlapping area between p L (z) and p R (z). If p L (z) = p R (z) (perfect overlap), JSD vanishes. If p L (z) has no overlapping region with p R (z) (complete separation), JSD reaches its maximum value ln 2. Any values of JSD less than ln 2 signal that finite overlap between p L (z) and p R (z) exists, and resulting separation is inaccurate. JSD therefore provides the well-defined scale of separation precision, taking into account shape details of p α (z). Fig. 4 is JSD as a function of ε for different flow patterns and chiral objects. It follows from Eqs (15) and (18) that when λ = 0, d/σ > 1 for t > t s , while for λ ≠ 0, it is saturated to a constant value after t λ . The occurrence of macroscopic separation can therefore be determined by measuring JSD at time t larger than t s (for λ = 0) and t λ (for λ ≠ 0) which are estimated as follows. For the high Pe regimes considered in this work, t s is given by t σ ~ δ/cV. In simulations, the upper bound of δ is (1)  , the lower bound of c is − , and the smallest value of |ε| other than zero is 0.01 in Fig. 4. Thus, JSD's are measured at  τ t (10 ) 3 for ε ≤ 0, τ (10 )  for 0 < ε ≤ 0.2, and  τ ( ) for ε > 0.2, which are long enough for JSD to reach its stationary value. For numerical evaluations of JSD, the probability distributions p(z) are discretized into histograms with the bin size of the order of the object size. The complete separation with the maximum JSD of ln 2 is achieved in all considered cases only when ε (and thus at least one of the eigenvalues) is vanishingly small. A finite value of JSD results for small  , it is more pronounced that JSD yields very small values unless ε is vanishingly small. As shown in Fig. 4, the simulation results clearly demonstrate our claim that the complete separation indeed occurs by quasi-two-dimensional flows satisfying the condition, Eq. (21).

Depicted in
Possible applications. Finally, we discuss the relevance of our linear flow analysis in terms of possible practical applications. In most cases, it is challenging to realize linear flows persisting for extended length and time scales in microfluidic devices. A simple shear flow is a well-known exception. We exhibit here two other examples where the linearity of flows is easily assured and at the same time, the present analysis can be useful for practical purposes. In particular, we show that for chiral separations, it is enough to have a linear flow only in a localized region if a confining potential is applied together.
Vortex flow. According to our theoretic formulations, the chiral separation occurs when both of the following conditions have to be satisfied: First, one of the eigenvalues of velocity gradient tensor should be much smaller than the inverse of the separation time scale (the eigenmode analysis). Secondly, there should exist a non-vanishing rate-of-strain field that induces a finite drift velocity (the parity-inversion argument). Now we present a salient example explicitly showing the indispensable role of the rate-of-strain field: consider a pure rotational flow with a perfect circular streamline, defined by the velocity gradient tensor, This flow satisfies the first condition but obviously not the second condition because it has a vanishing rate-of-strain field, E = 0. Hence it cannot separate any kind of chiral pairs. On the contrary, a vortex flow deformed by a finite rate-of-strain field, e.g., with the following velocity gradient tensor, fulfills the both conditions and could, therefore, lead to the separations, as indeed confirmed by the Langevin dynamics simulations (Fig. 5). The vortex flow can be an obvious solution for chiral separations, if not perfectly circular, and is another example of linear flows that persist for an extended period of time in microfluidic setups. As shown here, the vortex flows with circular streamlines have very distinct separation powers, depending on whether or not they have a finite rate-of-strain field. Our prediction on the separation power of vortex flows can be tested by a microfluidic four-roll mill device suggested by Lee et al. which can produce the entire spectrum of flow types, from purely rotational flow to purely extensional flow, by varying flow rate ratio 24 . Considering the dimensions of the device, the sub-micron helical objects ( µ  ⪅ m), for example, should be observed to be separated for elongated vortex flow but not for purely rotational flow. Also using various chiral objects of different shapes in this setup can be a feasible way to confirm our theory predicting that the flow properties rather than object shapes are essential to chiral separation.
Confining potential. At low Reynolds numbers, the resistance formalism linearly relates the hydrodynamic net force and torque exerted on a rigid object to the flow parameters 21 . Due to the linearity of the relations, originating from the linearity of the Stokes equations, the external forces such as confining force can be separately added into the equations of motion (see the Supplemental Information). Consequently, even in the presence of external forces our analysis can be performed in a similar way and the separation criterion remains intact in general.
In order to numerically verify this, we perform Langevin dynamics simulations with Flow-A and Flow-B previously defined, considering helical and tetrahedral objects under an external potential, φ(x, y) = k(x 2 + y 2 ). Here, we set ε = 0, and then Flow-A and B represent a two-dimensional extensional flow and a shear flow, respectively. The external potential plays a role as a two-dimensional confinement, constraining particle positions near to the origin in the xy plane. We note that for both flows, the separation occurs along the z-axis. The right panel of Fig. 6 shows the simulation results of temporal evolution of JSD in the presence of the potential while the left panel exhibits the results without the potential. As deduced, one can see that the separation behaviors remains the same qualitatively. This suggests that for practical applications of our analysis, a linear flow field does not necessarily have to persist over extended length and time scales. The two-dimensional confinement potential makes the separation (e.g., along the z-axis) take place only in a limited space of the xy-plane where the linearity of flows can be rather easily assumed for a time longer than the separation time scale.

Discussion
In practical applications of harnessing microfluidic devices for separations, one of the central questions will be to determine which flow has separation capability. Despite increasing attention to microfluidic chiral resolutions, most of the previous works are restricted to be considering an object of a specific shape in a given flow pattern. The complicated mathematical structure of the equations of motion present difficulties in understanding chiral separation phenomena even for a specific flow, and a comprehensive picture of common mechanisms and general conditions for chiral separations in terms of flow properties has been lacking. We have tackled this problem for an arbitrary linear flow using simple ideas, namely, considering symmetry properties and adopting the eigenmode analysis of flow fields. These enable us to draw an intuitive physical picture of the underlying mechanism of chiral separation dynamics. According to our results, the common features of separable flows are summarized as i) flows with a finite strain-rate tensor and ii) quasi-two-dimensional flows with small eigenvalues obeying Eq. (20). The typical examples satisfying both conditions are shear flow, vortex flow, and two-dimensional extensional flow, all of which are indeed demonstrated here to cause separations via the Langevin dynamics simulations. The present study thus provides a theoretical understanding of why two-dimensional flows such as shear and vortex flow are efficient to induce chiral separations. Our results provide simple criteria that would allow us to categorize and decide what kind of flows have a separation power or not. This is the prediction that cannot be easily made, without complicated numerical calculations or extensive simulations, from the theoretical studies known hitherto. It is also important to note that the separation criteria only concern the properties of flows, not of objects, so that they are applicable for objects of different shapes. In the presence of confining potential φ(x, y) = k(x 2 + y 2 ) with k = 8k B T/a 2 , JSD shows the very similar behaviors and the complete separation is also achieved for all cases. For helix and tetrahedral, the same objects as in Fig. 4 are considered, the number of ensembles is 10 5 , and the initial distributions are given as Gaussian centered at the origin with dispersion of 100a.