Simultaneous positioning and orientation of a single nano-object by flow control: theory and simulations

In this paper, we theoretically describe a method to simultaneously control both the position and orientation of single nano-objects in fluids by precisely controlling the flow around them. We develop and simulate a control law that uses electro-osmotic flow (EOF) actuation to translate and rotate rigid nano-objects in two spatial dimensions. Using EOF to control nano-objects offers advantages as compared to other approaches: a wide class of objects can be manipulated (no magnetic or electric dipole moments are needed), the object can be controlled over a long range (>100 μm) with sub-micrometer accuracy, and control may be achieved with simple polydimethylsiloxane (PDMS) devices. We demonstrate the theory and numerical solutions that will enable deterministic control of the position and orientation of a nano-object in solution, which can be used, for example, to integrate nanostructures in circuits and orient sensors to probe living cells.


Introduction
We theoretically describe a technique for simultaneously positioning and orienting single nanoobjects in a fluid in two spatial dimensions by manipulating the flow around them. We address two object control goals with this technique. The first is the ability to move the object across large distances (tens of micrometers). The second is the ability to accurately control the position and orientation of objects of a variety of shapes and material properties-for example, semiconductors [1], conductors [2] and dielectrics [3].
Previously, we showed simulations demonstrating position control of spherical objects [4] in a microfluidic device using electro-osmotic flow (EOF) control. The ideas developed in that study enabled an experimental demonstration of position control of micrometer-sized spherical objects to sub-micrometer accuracy [5] and subsequently the control of nanoscopic particles (single quantum dots) to nanometer precision [6]. In this paper, we explain and simulate a technique that shows how the translational and shear components of the flow field in the device can be manipulated, to trap an object at a desired position and orientation or manipulate both its position and orientation along a desired trajectory.
Existing approaches for simultaneously controlling the position and orientation of a single (or a few) object(s) can be classified into actuation strategies that include optical, magnetic and electrical techniques. If objects have a higher refractive index than their surrounding medium, a laser beam can be used to attract the object into the region of highest light intensity. Variations of this basic principle have been used to transfer translational and angular momentum from the laser beam to the object either by making use of special optical properties of the object or by manipulating the wave front of the incident light [7]- [13].
Magnetic fields in combination with fluidic forces have been used to control magnetic objects [14] or with a combination of optical forcing for translation and electromagnets for rotation [15]. Apart from the requirement that the object be magnetic, the objects may require specially designed shapes to enable fluidic actuation [16]. Alternatively, magnetic features are lithographically patterned on objects to allow manipulation of their orientation by magnetic fields [17]. The electrodes actuate a flow that translates and rotates the object from its current to its desired position and orientation. The right panel shows the feedback loop that achieves flow control in the device. At every instant, a camera captures an image of the object and an image processing algorithm computes the position and orientation of the object and transmits that information to the controller. The controller uses this information to actuate a flow in the device (by creating an electric field that moves the fluid in the device) that translates and rotates the object to the desired position and orientation.
If the object has a significant dipole moment, dielectrophoresis (a type of electric actuation) can be used to position and orient nanowires [18]- [23] and biological cells [24]. In this technique, a high electric field gradient interacts with the object's dipole moment to translate and rotate the object. Since the electric field gradients have to be high, the object is controlled near an electrode in a device where the gradient is steep.
The position and orientation control that we present below uses electrical actuation to modulate the flow around an object. Our technique depends on controlling viscous drag, a force that applies to every object, and can hence be used to control a general class of objectthe object does not have to be charged or magnetic or have any other special properties. Our simulations will show that the object can be controlled across a large (≈100 µm) region. The approach uses feedback control of shear flow in a microfluidic device to translate and rotate the object. The fluidic shear force acting on the object (which is assumed to be inertia-less) suspended in the fluid rotates the object to any desired orientation (left panel of figure 1), while the translational component of the fluidic drag moves it to any desired location in a two-dimensional (2D) control region. The position and orientation of the object, which are randomly perturbed due to Brownian motion, are measured at regular intervals. From this visual measurement, a feedback control loop determines and applies the fluid flow that will translate and rotate the object from where it is towards where it should be (see right panel of figure 1).
The flow is actuated electro-osmotically. Here, the electric field moves the fluid, which moves the object by viscous forces [25,26] (which is different from electrophoretic or dielectrophoretic actuation where the electric field creates a force directly on a charged or polarizable object). By using multiple (here eight) electrodes in concert (see figure 2), it is possible to create complex electric field patterns in space, and in our thin planar devices those patterns are faithfully transmitted to the fluid flow in the control region. Feedback control of  those patterns in time has allowed us to control the position of single objects [6] and is here being extended to also allow control of their orientation. The rest of the paper is organized as follows. Section 2 describes electro-osmotic flow. Section 3 describes the effect of the fluid on the motion of the object. Section 4 discusses the control algorithm that is used to manipulate the flow around the object. We show numerical examples demonstrating object control in section 5 and end with a discussion of additional considerations towards experiments in sections 6 and 7.

Electro-osmotic flow
At the interface between solid and an electrolytic fluid, the surface energy of the solid surface is reduced by the adsorption of ions from the fluid onto the oppositely charged ions at the surface of the solid [25,26]. This results in a charge imbalance in a thin (<100 nm) fluid layer, termed the diffuse (or Debye) layer, adjacent to the wall-fluid interface of the device.
A potential difference applied at the electrodes creates a planar electric field E(x,ŷ) in the plane of the device. That electric field moves the ions in the Debye layer, which in turn drags the rest of the fluid in the device due to viscosity. The flow is laminar (Reynolds number Re <10 −4 ), has a steady-state velocity profile U (x,ŷ) of a plug flow (that is constant alongẑ apart from the variation in the thin Debye layer as shown in figure 3) that is linearly proportional to the applied electric field and has the value [25] [25]): The negatively charged surface of the device is shielded by positively charged ions from the electrolyte solution. The ions in the thin diffuse (Debye) layer near the device-fluid interface move under the influence of the electric field and drag the rest of the fluid by viscous forces [25]. The resulting flow profile is uniform alongẑ (except for the variation in the thin diffuse layer, not drawn to scale in the figure) with the flow velocity U proportional to the applied electric field E.
where µ and are the dynamic viscosity and permittivity of the fluid, respectively, and ζ is the potential difference acrossẑ between the edge of the Debye layer and the device-fluid interface.
Since the electric field is irrotational (curl-free) the fluid velocity in the device, which follows the local electric field, is also irrotational. This means that EOF can only impart a translational velocity to a spherical object-one cannot rotate a sphere using EOF. However, as described in the next section, EOF can translate as well as rotate a non-spherical object such as an ellipsoid or a nanorod.

Translational and rotational velocity of an ellipsoid in electro-osmotic flow (EOF)
Similar to the inertia-free spherical particles 6 used previously to demonstrate positional control [4,5], an inertia-free ellipsoid will instantaneously translate along the stream lines of any flow that is set up in the device. If the flow in the device was rotational, one could create a vortex flow and the rotational velocity of a spherical object would be proportional to the vorticity 6 A simple calculation using Newton's laws [4] shows that a spherical particle of diameter 10 µm with a density equal to that of water takes 0.04 ms to reach the ambient translational flow velocity. This interval is small compared to the expected 10 ms between successive control updates, which allows us to assume that the nanorods (major axis length of 10 µm and minor axis length of 200 nm) have negligible inertia and their velocity instantaneously conforms to the local flow field. If the moment of inertia of the sphere is denoted as I s and the rotational drag coefficient as s , then a similar calculation for rotation shows that a time t r = I s ·log(β) s = ρa 2 log(β) 15 µ is needed for the sphere to reach within 1 β of the ambient rotational flow velocity of the fluid (ρ and µ are the density and viscosity of the fluid, respectively, and a is the sphere radius). For the above sphere in water, t r ≈ 0.01 ms for β = 1000.  The total torque acting on the body is the sum of infinitesimal torques due to the shear force acting throughout the boundary of the ellipsoid. Due to the unequal axes lengths of the ellipsoid, a saddle flow (illustrated above) can be shown to rotate the ellipsoid clockwise [29]. The body-fixed frame of reference is x-y, and the device-fixed frame isx-ŷ.
of the flow. However, for an irrotational flow, one can nevertheless rotate an ellipsoid due to the interaction of the individual shear components of the flow and the reduced symmetry of the ellipsoid (as compared to a sphere) as shown below. After computing the flow velocity, exact expressions for the force and torque acting on an ellipsoid [29] can be obtained by summing the infinitesimal shear force components (and the resultant torques) acting on every point of the ellipsoid's boundary. For the saddle flow shown in figure 4, the net fluidic torque will rotate the ellipsoid clockwise. The slow, viscous, incompressible and isothermal flow that is electro-osmotically set up in the device is well described by Stokes flow [29], which neglects the momentum of the fluid in the Navier-Stokes equations of fluid dynamics. The ellipsoid is controlled in a plane parallel to the floor of the device and is assumed to lie far enough (>200 nm) from the floor and ceiling of the device to neglect increased drag due to wall effects [30]. This claim is supported by recent calculations [31], which can be applied, for example, to compute the drag correction terms for a cylindrical particle of diameter 200 nm and length 2 µm, with the nearest point on its surface located at a distance d = 200 nm from a plane wall, and its axis tilted at an angle φ to the plane of the wall. Detailed simulations [31], performed for such cylinders (with aspect ratios of 10), show that the correction term for the (x, y) translational drag coefficients and the rotational drag coefficient about the z-axis are each less than 15% as compared to the unbounded fluid drag coefficients for φ as large as 45 • . The correction terms decrease as the aspect ratio of the cylinder increases because a proportionately lower area of the rod is close to the wall [31]. Hence, we disregard the wall-correction term and use the unbounded fluid drag coefficients in what follows.
The flow setup in the device is perturbed by the presence of the ellipsoid [32]- [34]. The linear nature of Stokes flow can be exploited to obtain the force and torque acting on any body (not just an ellipsoid) that is immersed and free to move in the flow. Denote the surface of the ellipsoid by x 2 a 2 1 + y 2 a 2 2 + z 2 a 2 7 semi-major axis length and a 2 is the semi-minor axis length of the ellipsoid). The unperturbed flow field u( r ) (i.e. the flow that would be observed if no object were present) in the device is approximated to be a superposition of uniform and pure-shear flow fields, i.e. spatial variations of u( r ) of O(a 2 1 ) or higher are neglected. Denote the translational velocity of the ellipsoid along the device-fixed axesx,ŷ (see figure 4) by Ux and Uŷ, respectively. As shown in figure 4, θ is the angle between the body and device fixed frames of reference. Denote the angular velocity of the ellipsoid about thê z-axis by ωẑ. Denote byû andv the components of the unperturbed uniform flow along thê xandŷ-axes, respectively, evaluated at the position that is occupied by the center of the ellipsoid, while the terms ∂û ∂x , ∂û ∂ŷ , ∂v ∂x and ∂v ∂ŷ represent the shear components of the flow. The constant e = a 2 a 1 is the ratio of the minor to major axis lengths of the ellipsoid. For the flow u( r ) described above, ignoring parasitic pressure flows and using the assumption of negligible inertia of the ellipsoid, it can be shown [29] that the translational velocity of the ellipsoid matches the uniform flow field component Ux =û,

Uŷ =v
(2) and the rotational velocity of the ellipsoid is given by Equations (2) and (3) are derived by integrating the shear and pressure distributions on the surface of the inertia-less ellipsoid after solving the quasi-static Stokes equations [27] (see section A of the supplementary information for more details, available from stacks.iop.org/NJP/13/013027/mmedia). Even if higher order spatial flow variations were considered (in addition to uniform and pure shear flow terms), successive correction terms in equations (2) and (3) would differ by the operator , where x i are the axis coordinates [34]. In the proposed device, these higher-order terms are of the order O(a 2 i /r 2 dev ) and higher compared to the linear terms included in equations (2) and (3); here r dev is the distance from the center of the control region to the midpoint of any straight edge at the boundary in figure 5. Since a 1 r dev is 0.1 or smaller in the proposed device, we ignore higher spatial order flow variations.
Equation (3) can be simplified by making use of two relations between the four shear components. By continuity, the divergence of the flow field is zero, i.e. ∂û ∂x + ∂v ∂ŷ = 0. We also have that in EOF, since the fluid velocity is proportional to the electric field (as explained in the previous section), the flow field is curl-free-i.e. the component of the the vorticity about theẑ-axis ( ∂v ∂x − ∂û ∂ŷ ) (the first term of equation (3)) is identically zero. This simplifies equation (3) to Thus, if we apply an electric field that creates the unperturbed flow field (û,v) in the device, then this flow will instantaneously turn the ellipsoid with the rotational velocity ωẑ in equation (4). . The leftmost edge is maintained at 1 V and all the others are maintained at 0 V. The electric field E(x,ŷ) (which is responsible for the object's translation) is shown with white arrowheads. The color plot shows the ellipse rotation that will be created. For an ellipsoid oriented at θ = 0 • , the shear component ∂ ∂ŷ (Ex (x,ŷ)) decides the direction of the ellipsoid's rotation (see equation (4)). Thus, in the reddish hued region, the ellipsoid will turn counterclockwise (seeing into the paper), and in the bluish hued region it will turn clockwise. We plot sign( ∂ ∂ŷ (Ex (x,ŷ))) log 10 (| ∂ ∂ŷ (Ex (x,ŷ))|) to show both the sign and magnitude of the rotation creating term, which varies between ≈ ±10 9 V m −2 .
Since the flow velocity componentsû andv are linearly dependent on the electric field, the translational and rotational velocities of the object are also linearly dependent on the electric field-a fact that is used in the control algorithm. For any arbitrary orthotropic object (an object with three mutually perpendicular planes of symmetry) like the ellipsoidal rod (or other bodies like right elliptical cylinders and rectangular parallelepipeds), the exact same analysis that has been used in this section can be applied to obtain the translational and rotational velocities; the only difference will be a different shape-dependent constant, instead of 1−e 2 1+e 2 , in equation (4) [28]. For non-orthotropic objects there will be an additional dependence of the rotational velocity of the object on the fluid velocity (and not just the fluid shear as in equation (4)). This is because the force at a particular point on the boundary of the body can give rise to a torque that is not balanced by an opposing torque (due to a lack of sufficient symmetry in the body), thus causing the body to rotate [29]. There will also be an additional dependence of the translational velocity of the object on the fluid shear (and not just the fluid velocity as in equation (2)). As explained later, this kind of coupling between the translational and rotational motions of the object can be included by experimenting with an appropriate gain parameter in the control algorithm.

Translational and rotational Brownian motion of an ellipsoid
Thermal equilibrium between the fluid and any object suspended in it is maintained by the random collisions between the object and the surrounding fluid molecules. The object translates and rotates in a time interval dt along random directions by an amount that is, on average, proportional to √ dt. Expressions for the translational and rotational diffusion coefficients [35] of an ellipsoid are given in section B of the supplementary information (available from stacks.iop.org/NJP/13/013027/mmedia). For the ellipsoids in our simulation (semi-major axis length a 1 = 5 µm and semi-minor axis length a 2 = a 3 = 100 nm) that are immersed in water of viscosity (8.9 × 10 −4 ) Pa s at 300 K, the translational diffusion coefficient along the major axis is 0.306 µm 2 s −1 and along the minor axis is 0.199 µm 2 s −1 . The rotational diffusion coefficient aboutẑ is 0.0197 rad 2 s −1 . The simulations in section 5 account for ellipsoid dynamics due to both fluid flow and Brownian motion.

Feedback control of the object's position and orientation
In the presence of Brownian motion of the ellipsoid, one can control the motion of the ellipsoid by using a feedback control algorithm. This control algorithm computes the voltages that need to be applied at the electrodes so that the resultant electric field creates a flow in the device, which translates and rotates the object from the currently measured to the desired position and orientation. For any position and orientation of the object, there exists a linear map between the object's velocity and the voltages applied at the electrodes (as explained next). The necessary control electrode voltages can be computed by inverting this map, using least squares. At successive time steps, the object moves to a new position and orientation, there is a new linear map, and we solve another least squares problem to get the next set of electrode voltages. This computation can be done in real time, even for complex situations, as demonstrated in our previous experiments in which we controlled the position of multiple particles at once [5].
The desired trajectory of the ellipsoid is a series of discrete, closely spaced points in the control region with a prescribed desired orientation at each point. At each control update, the difference between the current measured position (orientation) of the object and the desired position (orientation) is multiplied by a proportionality constant called the control gain, yielding the desired (translational and angular) velocity of the object until the next control update. The controller inverts the map between the voltages and the desired translational and rotational velocities to determine the electrode voltages.
We now discuss how to compute this needed linear map, which is a composition of three individual maps related to the three physical processes that control the particle motion. In our previous experiments [5], we did not observe any significant spatial or temporal variations in the surface properties of the PDMS or glass substrates, or the viscosity of the fluid; hence we treat the zeta potential, permittivity and viscosity as spatially and temporally invariant. The first map relates the applied voltages to the resulting electric field (including the gradient of the electric field) in the control region of the device. The second relates the electric field to the fluid flow field in the device. The final map relates the flow field to the object's translational and rotational velocities. We will then show how the composition of these three maps is inverted using least squares.
For the first map, Laplace's equation for the electric potential ∇ · ( ∇( ) = 0 is evaluated, where is the permittivity of water and is the electric potential in the domain shown in figure 5 (which contains the control region as shown in right panel of figure 2).
Eight electric fields are pre-computed, one for each channel. During operation, the electric field (or its gradient) at any point in the control region is a linear superposition of these eight electrical fields. Mathematically, this allows us to write the equations of the first map where η i are the applied voltages and the matrix A(x,ŷ) (of size 4 × 8) is known for each point (x,ŷ) in the control region using the pre-computed fields. From equation (1), the fluid velocity at any point in the control region is directly proportional to the electric field. The flow velocity is assumed to reach steady state instantaneously after the potential difference is applied 8 at the electrodes. The fluid shear (spatial gradient of the velocity) is directly proportional to the gradient of the electric field at that point, which gives the second map, between the electric field in the device and the resulting flow and shear field,û If the ellipsoid's center of mass is at the point (x,ŷ), then the four relations in equation (6) for the fluid's velocity and the shear at (x,ŷ) are the only ones needed for the third map, which relates the fluid's velocity and shear to the object's translational and rotational velocities by equations (2) and (4). After combining the three maps, the object's translational velocities Ux , Uŷ and its rotational velocity ωẑ are given by the final composite linear map where η i are the applied voltages and F c (θ) and F s (θ) are given by F c (θ) = 1−e 2 1+e 2 cos(2θ) and F s (θ) = 1−e 2 1+e 2 sin(2θ ). We now show how best to select the eight electrode voltages to achieve 8 The time T ss required to reach steady state is given by T ss = ( log(100) π 2 ) · ρh 2 µ , where h and ρ are the height of the channel and the density of the fluid, respectively, as derived in [36]. T ss equals 0.32 ms in the device, which is much less than the control update interval of 10 ms that is typically used in our flow control device. No discernible time lag in fluid actuation was observed in previous, position control experiments [5] in our group. Hence, we assume that the flow velocity reaches steady state as soon as the voltages are applied. the desired object velocities Ux , Uŷ and ωẑ that will translate and rotate the object from where it was to where it should be (a similar argument will hold for the general case of N electrodes in an N -channel device).
At every control update, we request the desired translational and rotational velocities of the object, which are uncoupled in the case of orthotropic particles like the ellipsoid, to be, respectively, proportional to the positional and orientational deviations from the desired trajectory. These deviations, or errors, from the desired path are εx =x des −x, εŷ =ŷ des −ŷ and ε θ = θ des − θ , so (εx , εŷ) is the difference between the desired and current object positions and ε θ is the difference between the desired and current object orientations. The desired object velocities Ux , Uŷ and ωẑ are set to be proportional to the errors εx , εŷ and ε θ by the proportionality gain matrix K prop The gains K r and K θ are penalties on the translational and orientational errors, respectively. A higher value of K r forces the controller to select voltages that will translate the object to the desired position more quickly. Similarly, a higher value of K θ forces the controller to select voltages that will rotate the object to the desired orientation quicker. The relative values of K r and K θ decide whether the controller spends more of its control authority on the object's translation or on its rotation. Combining equations (7) and (8), there are more unknowns (the actuator voltages η = (η 1 η 2 · · · η 8 ) T ) than there are known quantities (the desired velocities Ux , Uŷ and ωẑ). A least square solution, which chooses the minimal size control that achieves the desired velocities, is used to find η. Denote the linear map of equation (7) by the matrix P(x,ŷ) The least square fit computes the voltages η i , which minimize the 2-norm [37] of the electrode voltages η 2 . The optimal voltages are given by  where P + (x,ŷ) = (P T (x,ŷ)P(x,ŷ)) −1 P T (x,ŷ) is the pseudo-inverse [37] of the matrix P(x,ŷ) and P T (x,ŷ) is the transpose of the matrix P(x,ŷ). This is the control law-it states how to compute the electrode voltages given the difference between the actual and desired ellipsoid positions and orientations.
As for the control design in our previous experimental work [5], in order to avoid electrolysis (the formation of bubbles at the electrodes that can disrupt the intended flow) we limit the voltages to a maximum value (termed the saturation voltage) η sat = 0.15 V. Hence, the voltages determined in equation (10) are linearly scaled so that this constraint is not violated. If the maximum absolute value of the eight voltages η i , computed by the controller in equation (10), is η max , then the eight scaled voltages η scaled i that are eventually applied at the 12 electrodes are given by Since the ellipsoid velocity is linearly proportional to the applied voltages, this scaling limits the magnitude of the maximum achievable translational and rotational velocities of the object. At every time step this scaling might reduce the magnitude of the object's velocity but not the direction (even after voltage scaling, the particle is steered to correct for the deviation from the desired path, but possibly at a lower speed).

Non-dimensionalized equations for the feedback loop
A non-dimensionalized version of the governing equations using the flow Peclet numbers is presented here. The relevant physical parameters governing the dynamics are the length scales given by the radius of the control region of the device r dev (in figure 5, r dev is the distance from the center of the control region to the midpoint of any straight edge at the boundary), the saturation voltage η sat , the zeta potential ζ at the fluid-PDMS interface and the fluid's permittivity and viscosity µ. In what follows, the non-dimensional parameters are superscripted with an asterisk.
Since the electric field in the device scales as η sat r dev , the translational velocity V EOF due to EOF will scale as V EOF = ( ζ µ )( η sat r dev ). The time needed for the particle to traverse the control region will then scale as t dev = r dev V EOF . The electro-osmotic shear σ EOF generated in the device will scale proportionally to the gradient of the electric field. Since the gradient of the electric field in the control region scales as η sat r 2 dev , the shear will scale as σ EOF = ( ζ µ )( η sat r 2 dev ), so σ EOF = 1 t dev . This scaling of V EOF and σ EOF ignores the contribution of the shape of the control region, which is fixed once the number channels and the ratio c dev r dev are fixed (c dev , the channel width, is the width of the straight edge at the boundary of the control region in figure 5). For the rest of this paper, as shown in figure 5, the number of channels is fixed at 8 and c dev r dev at 0.5. The non-dimensional displacement and time parameters are chosen asx * =x/r dev ,ŷ * = y/r dev , θ * = θ and t * = t/t dev . The translational [û(x,ŷ),v(x,ŷ)] and shear components [ ∂û(x,ŷ) ∂ŷ , ∂v(x,ŷ) ∂ŷ ] of the flow field and their non-dimensional counterparts [û(x,ŷ) * ,v(x,ŷ) * ] and [( ∂û(x,ŷ) ∂ŷ ) * , ( ∂v(x,ŷ) ∂ŷ ) * ] are related byû(x,ŷ) = V EOF ·û(x,ŷ) * ,v(x,ŷ) = V EOF ·v(x,ŷ) * , ∂û(x,ŷ) ∂ŷ = σ EOF · ( ∂û(x,ŷ) ∂ŷ ) * and ∂v(x,ŷ) ∂ŷ = σ EOF · ( ∂v(x,ŷ) ∂ŷ ) * . With the non-dimensional map A(x,ŷ) * chosen as the non-dimensional flow components can then be stated in terms of A(x,ŷ) * as   (10) and (11).
The motion of the ellipsoid is governed by e * 1 (= a 1 r dev ), e (= a 2 a 1 ) and the rotational Peclet number Pe θ = σ EOF /D θ where D θ , is the rotational diffusion coefficient aboutẑ. Since the expected time for the ellipsoid to rotate by 1 radian aboutẑ due to diffusion and the applied actuation are 1 2D θ and 1 σ EOF , respectively, the quantity Pe θ compares the actuation's ability to compensate for the rotational diffusive motion [38]. A larger value of Pe θ signifies a higher ability of the actuation to compensate for the diffusive motion aboutẑ. Denoting the translational diffusion coefficients along the major and minor axes of the ellipsoid by D x and D y , respectively, we define two functions, T θ x (e) = D x a 2 1 D θ and T θ y (e) = D y a 2 2 D θ , both of which depend solely on the parameter e (see section B.1 of the supplementary information, available from stacks.iop.org/NJP/13/013027/mmedia for details). The function T θ x (e) is the ratio of the expected time taken for an ellipsoid to rotate by 1 radian due to rotational diffusion aboutẑ, to the expected time taken by the ellipsoid to diffuse by a distance a 1 along the major axis of the ellipsoid. Similarly, T θ y (e) is the ratio of the expected rotational diffusion time, to the expected time taken to diffuse by a distance a 2 along the minor axis of the ellipsoid.
With the control law given by equations (8) and (10), the controlled dynamics of an ellipsoid, including the effect of Brownian motion, is as follows. The geometric center translates by an amount dx * alongx * , dŷ * alongŷ * , while the ellipsoid rotates by an amount dθ * aboutẑ in time dt * according to the stochastic dynamics update given by dθ * = where dB * x and dB * y are the translational Brownian displacements alongx andŷ, respectively, and dB * θ is the rotational Brownian displacement aboutẑ, which are given by where e * 2 = a 2 /r dev = e * 1 · e. The factors N x (0, 1), N y (0, 1) and N θ (0, 1) denote independent Gaussian random variables with mean 0 and variance 1 which reflect the random character of Brownian motion (see [39] for an introduction to stochastic update formulae). Equations (10)-(15) describe the controlled motion of the ellipsoid in non-dimensional terms. We simulate this motion for different ellipsoid manipulation tasks and Peclet numbers in the next section. Major axis length of ellipsoid 5 µm a 2 (= a 3 ) Minor axis length of ellipsoid 100 nm

Numerical simulations of an ellipsoid's manipulation
In this section, the control law described previously is used to simulate the steering of an ellipsoid along a desired trajectory, while the ellipsoid is being perturbed by Brownian motion (according to the controlled ellipsoid dynamics described in equations (10)- (15)). We show how the ellipsoid can be translated and simultaneously rotated as well as rotated in place while trapped at a given location. We will show the dependence of the root mean square (rms) orientation error (while the ellipsoid is trapped) on Pe θ and the shape parameter e. We explain the source of this error in detail with a simple Fokker-Planck model which can be used to theoretically predict the trapping error in orientation given the size and aspect ratio of the rod and the flow parameters.
For showing control over different trajectories, the simulations were performed with parameters e * 1 = .05, e = .02 and Pe θ = 30. Table 1 states a sample set of values of the physical parameters that reflect the above non-dimensional numbers.
Unless otherwise stated, the non-dimensional simulation results are stated in degrees for rotation, the translational displacements in units of r dev , time in units of t dev and applied voltages in units of η sat . The applied voltages are updated every dt * = 5.9 × 10 −3 , which corresponds to the 10 ms time lag that is expected in experiments 9 . This time lag occurs due to the finite frame rate of the camera and the computational time required by the control and object vision detection algorithm. The ellipsoid's position and orientation are assumed to be perfectly known in these simulations. The choice of gain coefficients K r and K θ (see equation (8)) determines the extent to which the control authority is spent on controlling position versus orientation, respectively. We experimented with different values of the gain coefficients, for achieving the optimal trade-off between positional and orientational errors in simulations and settled onK r = K r * t dev = 1.7 × 10 6 andK θ = K θ * t dev = 8 × 10 4 for all the simulations shown here. 9 A computational time of 3.5 ms, using MATLAB (www.mathworks.com), see footnote 6, was required for previous experiments that implemented position control of spherical particles. This includes the time needed to estimate the particle position as well as the time needed to compute the control voltages. The frame rate of the Pike camera from Allied Vision Technologies (www.alliedvisiontec.com), see footnote 6, is 200 Hz at full resolution. Thus the total feedback loop time of 8.5 ms is within the 10 ms time assumed here. We will discuss the relation between η sat , the rotational diffusion coefficient, and the gain coefficients at the end of this section. The first simulation is of the ellipsoid tracing a square path, shown in figures 6 and 7.
The square path has a side length equal to 0.8. The ellipsoid initially starts off at the bottom left corner of the square and is returned to the same point at the end of the simulation. In the initial part of the simulation (from the bottom left to the top left corner of the square path), the ellipsoid is controlled to move along a straight line, at a constant orientation (θ * = 90 • ). Along the next three sides of the square, the ellipsoid is controlled to rotate by 90 • by the time it reaches the end of that side (see movie M1 in supplementary information, available from stacks.iop.org/NJP/13/013027/mmedia). The positional and orientation errors are shown in figure 8. In figure 9, we demonstrate the ellipsoid being controlled along a more complex 'hour glass' path shape that spans the entire control region. The ellipsoid starts at the bottom left corner of the trajectory at an initial orientation of θ * = 90 • . It is then controlled to move to the top right corner, then to the top left corner, down to the bottom right corner and finally back to the bottom left corner. While translating, the ellipsoid is controlled to simultaneously rotate so that its major axis is aligned with each of the four segments of the trajectory by the time it reaches the end of that segment (see movie M2 in supplementary information, available from stacks.iop.org/NJP/13/013027/mmedia).
The next simulation shows the ellipsoid both trapped in place and being rotated from an initial orientation of 90 • to a final orientation of 0 • by EOF control. The change in orientation of the ellipsoid over time and the associated flow field and electrode voltages are shown in figure 10 (see movie M3 in supplementary information).
In the initial part of the simulation (until t * = 1.17), the controller spends most of its control authority on rotating the ellipsoid from θ = 90 • to θ = 0 • compared to correcting for positional error due to translational Brownian motion. This results in a nonzero mean translational error of 1.9 × 10 −4 in position, with an rms of 1.3 × 10 −4 . In the tail end of the above simulation, when the ellipsoid is nominally trapped at the desired location and orientation, as shown in figure 12, we observe an rms error in the orientation angle. The rms error in orientation, rms sim , is computed from the simulations in the following manner: where θ * (t * i ) is the orientation of the ellipsoid at the ith time step. There are a total of (m + n) time steps in the simulation; the ellipsoid first reaches θ * = 0 • at the mth time step (termed the first passage time) and then θ * = 1 n+1 ( i=m+n i=m θ(t * i )). This rms error is a result of the competition  Denote the coordinates of the ellipsoid at time t * i by (x * (t * i ),ŷ * (t * i ), θ * (t * i )) and the coordinates that were desired at the previous time step t * . For the first side of the square (0 < t * < 1.2), apart from countering rotational Brownian motion, there was no need to rotate the ellipsoid. Hence, the available control authority could be committed more fully to correcting positional errors. At subsequent stages, the need to turn the ellipsoid leads to a decreased capability to correct positional errors.
between Brownian motion and the controller. The rotational Brownian motion of the ellipsoid will tend to rotate the ellipsoid away from θ * = 0 • , while the controller tries to bring it back to θ * = 0 • . Since the rotational dynamics depends on the rotational Peclet number Pe θ and the shape parameter e, we plot the dependence of the rms trapping error as a function of these variables in figure 11. For better visualization, the error is plotted against the aspect ratio 1/e. As one would expect, the orientation error decreases with increasing Pe θ because of comparatively larger actuation compensating for the particle's diffusive motion. For a fixed Pe θ , the plot shows a sharp increase in error as the aspect ratio approaches 1, i.e. as the shape of the ellipsoid approaches that of the sphere (unit aspect ratio). In addition, there is a slow increase in error, for a fixed Pe θ , as the aspect ratio increases (compare errors between aspect ratios 10 and 50). We explain this next.

The Fokker-Planck equation describing the orientation error
It is possible to theoretically estimate the variation of the rms error in orientation, termed rms theor , with respect to the size and the aspect ratio of the ellipsoid as where rms theor is in radians, σ max is the maximum rotational velocity with which one can turn the ellipsoid and D θ ≡ kT µ g(e) a 3 1 (see [35] and expressions for the rotational diffusion coefficient in section B of the supplementary information, available from stacks.iop.org/NJP/13/013027/mmedia) to see that D θ is inversely proportional to a 3 1 ). The value of σ max equals the maximum allowable fluid shear that can be actuated by the electrodes under the constraint that the applied voltages do not exceed η sat . In terms of Pe θ , equation (17) can be rewritten as where σ * max is the non-dimensional maximum shear that only depends on the shape of the device, i.e. on the number of channels and the device geometry parameter c dev r dev . A derivation of the expression in equation (17) makes use of two observations noted in the simulations. Firstly, the controller spends most of its control authority in maintaining the ellipsoid's orientation at θ = 0 • and relatively less authority on position control. Secondly, as seen in the bottom panel of figure 12 (which corresponds to e * 1 = 0.05, e = 0.02 and Pe θ = 30), the controller exerts this authority by maintaining the value of the shear component ( ∂û ∂ŷ ) * at Figure 11. Effect of aspect ratio and rotational Peclet number on rms orientation error: the parameter e * 1 , which corresponds to the semi-major axis length of the ellipsoid, was fixed at the value 0.05. For each (Pe θ , e), the plotted values are the average of 10 runs of the simulation with the ellipsoid being nominally trapped at 0 • . The error bar for each average is plotted on the graph, while a surface is fitted to the data to guide the eye. Increasing values of Pe θ means that more shear is available to compensate for the perturbation due to the rotational Brownian motion of the ellipsoid, thus decreasing the error. For a fixed Pe θ , there is a sharp increase in error as the shape of the ellipsoid approaches that of the sphere (aspect ratio ≈ 1.1). This effect is explained in the text. the maximum of ±1.6 for most time steps after t * ≈ 1.1 (when it is first oriented at θ * = 0 • ). It abruptly switches between ±1.6 as needed in order to counteract the rotational Brownian motion of the ellipsoid. Thus, the controller essentially executes a simple bang-bang-type control law [41] for maintaining the ellipsoid's orientation at θ * = 0 • : it checks whether the ellipsoid has positive orientation θ * > 0 • , or negative orientation, θ * < 0 • , and attempts to apply the maximum allowable shear (σ * max = 1.6), which can rotate the ellipsoid back to θ * = 0 • . The associated Fokker-Planck equation for the probability distribution function of θ that arises from the stochastic differential equation that describes the above simplified control law yields equation (17), as explained in detail in section C of the supplementary information (available from stacks.iop.org/NJP/13/013027/mmedia).
In figure 13, we compare rms theor and rms sim for a range of ellipsoid sizes a 1 and aspect ratios a 1 a 2 (note: e = a 2 a 1 ). The simulations for computing rms sim were performed with the flow parameters stated in table 1, with the ellipsoid being trapped by EOF control at the center of the control region. For plotting rms theor , we set the maximum shear, which is a property of the device/controller and not the ellipsoid, as σ max = 0.94 rad s −1 (corresponding to σ * max = 1.6 shown in figure 12). Given that we disregard any loss of control authority for position control in the simple model of the control law that was used to derive rms theor , the good match between rms sim and rms theor for the range of particle sizes and aspect ratios under consideration verifies that most of the controller's authority is indeed spent in controlling the orientation rather than the position of the ellipsoid.
For a fixed aspect ratio, the rotational diffusion coefficient varies inversely as a 3 1 (as opposed to the translational diffusion coefficient, which is inversely proportional to a 1 ). Hence as a 1 decreases, the rms error should increase as seen in figure 13. The increase in aspect ratio ( a 1 a 2 ), for fixed a 1 , causes a comparatively weaker increase in the diffusion coefficient and consequently in the rms error. For a fixed a 1 , the aspect ratio-dependent term 1+e 2 1−e 2 (in equation (17)) increases the rms error as e → 1. This is observed in figure 13, more noticeably for smaller particles, at a 1 a 2 = 1.1. This increase in rms error reflects the loss in the ability of a curl-free flow to rotate a near-spherical-shaped particle.  (17)) and the simulation-based estimate rms sim (equation (16)) for flow parameters stated in table 1. For each ellipsoid size and aspect ratio, the plotted value of rms sim is the average of 10 runs of the simulation with each run lasting 50 s while the ellipsoid was trapped at the center of the control region. The error bars for the simulation-based estimate rms sim are shown in the inset for a 1 = 2 µm, while the error bars for other values of a 1 are too small to depict on the plot. The plot shows the rms as the ellipsoid size (a 1 ) decreases. For a fixed a 1 , the plot shows an increase in rms error at the highest and lowest values of the aspect ratios ( a 1 a 2 ) with a dip in rms error at a 1 a 2 ≈ 2. The diffusion coefficient increases weakly with aspect ratio. Consequently, for a fixed size a 1 , there is a comparatively slow increase in rms error seen for values a 1 a 2 > 5 in the plot. As the aspect ratio decreases ( a 1 a 2 → 1), there is a sharper increase in rms error due to the term 1+e 2 1−e 2 in equation (17). This effect is noticeable for smaller ellipsoids (a 1 3 µm). The good match between rms theor and rms sim shows that orientation control is harder to achieve and so the controller spends most of its authority in controlling orientation rather than position. The ellipsoid size a 1 = 3 µm and the aspect ratio a 1 a 2 = 1.2 are seen to be lower bounds for maintaining an rms error of 5 • .
In another work by our group [6], we have been able to experimentally show position control of nanoscopic particles by EOF-based position control in a highly viscous fluid. This was possible because the actuation was not saturated while trying to control the particle position but the increased viscosity reduced the translation diffusion coefficient. However, this avenue is not feasible for orientation control because the actuation does saturate while trying to compensate for rotational diffusion. Thus the effectiveness of rotational control is set by actuator saturation. The rotational Peclet number Pe θ is independent of fluid viscosity. Seen in dimensional terms, the rotational diffusion coefficient D θ as well as the maximum shear σ max (which was set as the maximum value of the fluid shear ∂ ∂ŷ (û) observed in figure 12) are both inversely proportional to the fluid viscosity (see equation (11) in section B of the supplementary information, available from stacks.iop.org/NJP/13/013027/mmedia, for an expression of the diffusion coefficient and equation (6) in the main text for the shear field). Hence, increasing the fluid viscosity does not decrease the value of the rms error in orientation (since rms theor ≈ D θ /σ max ). Thus, in the control model described in this paper, for the flow parameters stated in table 1, even though slender particles (where a 2 = 100 nm for example) can be controlled, a 1 = 3 µm and a 1 a 2 = 1.2 are seen to be lower bounds on the particle size and aspect ratio for maintaining rms errors in orientation of around 5 • , irrespective of the fluid viscosity. Since the control voltages are not allowed to exceed η sat penalizing the orientation error in equation (10) by indefinitely increasing the control gain K θ does not reduce the orientation error below what is allowable by the maximum shear σ max . Instead one could reduce electrolysis at the electrodes by protecting them with a film of nafion [42]. This would allow for an increased value of η sat , and hence an increased value of σ max , thus reducing the rms error. This will be investigated in future experiments.

Additional considerations towards experiments
In this section, we consider additional features that will be important for future experiments and applications. We first consider the effect of the size and shape of more general objects for flow control. For other types of orthotropic particles, for example a cylinder, there will be a different constant, instead of 1−e 2 1+e 2 , that determines the rotational velocity of the particle in equation (4). This constant can be computed by simulating the flow around the cylinder with the Stokes flow equations (described in section A of the supplementary information) and integrating the resulting shear and pressure distributions on the surface of the cylinder. The cylinder's motion can be controlled using a similar control law modified with this new constant. Also in experiments, the exact size of the cylinder can only be known to a certain precision. Since the control voltages depends on the object size (through the P(x,ŷ) matrix in equation (9)) as well as on the gain matrix K prop in equation (10), one could compensate for the imprecise knowledge of the object size by experimenting with the value of K prop until a suitable performance is achieved.
For non-orthotropic particles, since translation and rotation are coupled, the structure of the P(x,ŷ) matrix will change. In particular, some of the elements of the matrix that premultiplies A(x,ŷ) in equation (9) will have a nonzero value. However, one can still compute the pseudoinverse of P(x,ŷ) and obtain the control voltages, after one makes an appropriate choice of K prop . In this case, K prop should be chosen to be non-diagonal to ensure that orientation and position errors are penalized in a manner that accounts for the coupling between translation and rotation.
We have shown simulations in which the object is being controlled in two dimensions; however, the translational diffusion of the object in the z-direction and the rotational 'pitching' diffusion about the body-fixed y-axis can bring the object in contact with (and thus might make it stick to) the glass and PDMS surfaces that form the floor and ceiling of the device, respectively. However, the object may have a natural electrostatic repulsion to the wall and, if not, coatings can be applied to the device that will prevent sticking of the object to the device [43]. With such coatings and by considering devices that are thin with respect to the object length, one can mitigate the tendency of the object to pitch out of the plane in which it is being controlled. Even when the object pitches out of plane, one does not completely lose the ability to control the rod in the plane of the device. As long as the pitching angle is small, one should expect to control the rotation aboutẑ (normal to the device plane) with the same control law as proposed in this paper, with a modified rotational gain K θ (see equation (8)).
Our control algorithm assumes that the object velocity exhibits a linear dependence on the electric field (see equations (2) and (4)). However, the object itself might have a charged Debye layer at its interface with the surrounding fluid. In the presence of an applied electric field, the ions in the Debye layer will move the local fluid surrounding the object, which can in turn impart an unintended translational and rotational velocity to the object. Such an electrophoretic motion of the particle, due to its own (uniform) zeta potential, is linearly dependent on the electric field [44] and can be readily accommodated within our control law as explained in section D of the supplementary information (available from stacks.iop.org/NJP/13/013027/mmedia).
However, depending on whether the object is strongly polarizable, one could observe a quadratic dependence of the velocity on the electric field, as opposed to the strictly linear dependence described in the previous section. This quadratic dependence is induced charge electrophoresis (ICEP) [45]- [48]. This effect is explained in section D of the supplementary information where we show that the magnitude of the rotational velocity of the ellipsoid due to ICEP using the approximations used by Saintillan et al [48] is negligible (less than 6% of that due to EOF) due to the low electric field strengths (less than 7 V cm −1 ) at the ellipsoid's location in the control region. The translational velocity of the ellipsoid due to ICEP is identically zero [48]. This matches well with experimental results [49], where the rotational velocity of a comparably sized object (6 µm × 300 nm) in water due to ICEP was noticeable only at electric field strengths that were greater than 30 V cm −1 . Hence this nonlinear effect can be safely neglected in our simulations.

Conclusion
We have described the physics that shows how a planar electro-osmotic flow translates and rotates an object that is immersed in it. The map describing this physics can be inverted to obtain a control law that allows one to manipulate the position and orientation of the object by electrically actuating the flow. Our simulations and theoretical model show how the performance of the algorithm scales with the size and aspect ratio of an ellipsoidal object and we have explained how this approach can be extended to objects of more general shapes. In order to realize what we have discussed experimentally, we need to extend our vision algorithm to also estimate the orientation of objects. The vision algorithm has to be robust to the uncertainty in measurement due to pixelation in the camera sensor, dark noise from the camera, low photon counts if the objects are dimly fluorescent or measurement noise due to imperfect estimation of a small metallic object's orientation from its scattered light. We are currently developing the needed vision algorithms as well as implementing our control methods in hardware, to demonstrate position and orientation control of one and multiple objects simultaneously in experiments. These results will be reported in future publications.