In this section, we describe our methods for steering particles using both electrokinetic and diffusiophoretic controllers. We begin by describing our BD simulations, and then, we describe the control approaches for each method. Finally, we briefly discuss the physical parameters used in the simulations.

A. Brownian dynamics simulations

The simulations are designed assuming that a single colloidal particle will be controlled within a 2D region, as shown in Fig. 2. Bordering the region are either electrodes or reactive chemical “probes.” For the electric field simulations, four electrodes are used, one on each of the north, south, east, and west sides of the control region [shown in Fig. 2(a)]. Similarly, for the chemical gradient simulations, four chemical reaction “probes” are placed on each side of the control region [shown in Fig. 2(b)]. The electrodes/probes are separated by a distance L, and the control region (the region in which the particles are allowed to move) is the area between the electrodes/probes. The simulation domain extends beyond the control region to a length of 3L. This makes the boundary conditions of the simulation distant enough that it does not disturb the motion of the particle.We use Brownian dynamics (BD) to simulate the motion of a colloidal particle moving in an externally applied field. We simulate a single particle (implying that there are no interparticle interactions) that feels either electrokinetic or diffusiophoretic forces, with other forces assumed to be zero. With the assumption that the motion of the particle is overdamped, the velocity of the particle is described by39–4139. K. D. Dorfman, D. Gupta, A. Jain, A. Muralidhar, and D. R. Tree, Eur. Phys. J. Spec. Top. 223, 3179 (2014). https://doi.org/10.1140/epjst/e2014-02326-440. H. C. Öttinger, Stochastic Processes in Polymeric Fluids (Springer, Berlin, 1996).41. D. T. Gillespie, Am. J. Phys. 64, 225 (1996). https://doi.org/10.1119/1.18210 X˙=1γ(Fdiff+Felec)+2Dcξ,(1)where X=[xparticle,yparticle]T is the 2D position of the colloidal particle; Fdiff and Felec are the diffusiophoretic and electrokinetic forces, respectively; and ξ is a Gaussian white noise term defined as ξ≡limdt→0N(0,dt−1),(2)where N(0,dt−1) is a normal distribution with a mean of 0 and a variance of dt−1 and dt is the time differential. Additionally, γ=6πηRc is the Stokes friction coefficient, and Dc=kbT/γ is the diffusion coefficient of the colloidal particle, which introduces Boltzmann’s constant kb, the system temperature T, the viscosity of the surrounding fluid η, and the radius of the colloidal particle Rc.The electrokinetic force Felec is a combination of the electrophoretic force Fep on the particle and the force imparted by electroosmotic flows Feo. In general, a particle can be electrokinetically steered using electroosmotic flows if it is neutral or a combination of electrophoresis and electroosmosis if it is charged.88. M. D. Armani, S. V. Chaudhary, R. Probst, and B. Shapiro, J. Microelectromech. Syst. 15, 945 (2006). https://doi.org/10.1109/JMEMS.2006.878863 However, electroosmosis complicates the comparison between electrokinetic steering and diffusiophoretic steering. Accordingly, since our primary purpose is to examine particle steering via diffusiophoresis, we will limit our analysis to the case where electroosmotic flows are suppressed (Feo=0), and a charged particle is steered solely using electrophoresis (Felec=Fep).The electrophoretic force is given by4242. R. F. Probstein, Physicochemical Hydrodynamics: An Introduction (John Wiley & Sons, Inc., Hoboken, NJ, 1994), Vol. 1, pp. 211–236. Felec(X,t)=γμeE(X,t)=−γμe∇Φ(X,t),(3)where μe is the electrophoretic mobility, E is the electric field, and Φ is the electric potential. As emphasized in Eq. (3), the electric potential is a function of position and time and is determined by the voltages ϕi(t) of the i electrodes. We calculate the electric potential by solving Laplace’s equationwith boundary conditions Φ(Xprobe,i,t)=ϕi(t) at the electrode locations Xprobe,i and Φ = 0 at the edge of the computational domain as indicated in Fig. 2. We discretize the domain using a 2D finite difference method and solve Laplace’s equation using numerical relaxation via the Gauss–Seidel algorithm.The diffusiophoretic force for a dilute, non-ionic solute is given by3838. J. Anderson, Annu. Rev. Fluid Mech. 21, 61 (1989). https://doi.org/10.1146/annurev.fluid.21.1.61 Fdiff(X,t)=γμd∇C(X,t),(5)where C is the concentration of the solute and μd is the diffusiophoretic mobility. The sign and the magnitude of the mobility depend on interactions between the particle, the solute, and the solvent,43–4543. J. T. Ault, P. B. Warren, S. Shin, and H. A. Stone, Soft Matter 13, 9015 (2017). https://doi.org/10.1039/c7sm01588g44. S. Shin, E. Um, B. Sabass, J. T. Ault, M. Rahimi, P. B. Warren, and H. A. Stone, Proc. Natl. Acad. Sci. U.S.A. 113, 257 (2016). https://doi.org/10.1073/pnas.151148411245. J. S. Paustian, C. D. Angulo, R. Nery-Azevedo, N. Shi, A. I. Abdel-Fattah, and T. M. Squires, Langmuir 31, 4402 (2015). https://doi.org/10.1021/acs.langmuir.5b00300 and the units for the mobility are m2/(M s), where M means moles per liter.Analogous to the need to know the electric potential gradient in the case of electrophoresis, one must determine the concentration gradient near the colloid to compute the diffusiophoretic force. Assuming the diffusion coefficient of the solute is constant (which is appropriate for a dilute solution) and that there is no coupling between hydrodynamics and diffusion (i.e., that there are no Marangoni effects or density-induced flows and that the colloid velocity is slow relative to diffusion), we determine the concentration field using a reaction–diffusion equation ∂C(x,y,t)∂t=Ds∇2C(x,y,t)+G(x,y,t),(6)where Ds is the diffusion coefficient of the solute and G(x,y,t) is the reaction rate of a local time-varying chemical reaction caused by reactive “probes” located at various points in the simulation domain as shown in Fig. 2. We used Dirichlet conditions C=0 at the boundary of the computational domain. As previously discussed, the zero concentration boundary acts as a concentration sink that keeps the average concentration from monotonically increasing with time. We calculate the concentration field using a finite difference method by discretizing Eq. (6) in two dimensions using a forward time centered space (FTCS) scheme on a regular grid. We chose a time step of Δtsim≤0.5(Δx)2/Ds to satisfy the Neumann stability criterion for the FTCS method,4646. A. Greenbaum and T. P. Chartier, Numerical Methods: Design, Analysis, and Computer Implementation of Algorithms (Princeton University Press, Princeton, NJ, 2012). where Δx is the distance between grid nodes.

B. Electrokinetic micro-actuation control algorithm

Our electric field micro-actuation control loop includes the following four steps: (1) A vision system measures the position of a particle relative to a set of electrodes arranged regularly around the edge of the device. (2) Using the position information, a dynamic model is used to relate the motion of the particle to the voltages of the electrodes. (3) A user-defined reference velocity is chosen and then modified with a feedback term to correct for Brownian motion and other inaccuracies in sensing, actuation, and modeling. (4) Electrode voltages are chosen by minimizing the difference between the dynamic model velocity and the reference velocity. Since these same steps will be used for the chemical gradient controller, with slight variations, we will explain each step in detail in the following paragraphs.

In step (1), we assume that a vision system is able to measure the position of a colloidal particle in real time. The proposed vision system would include a camera attached to a microscope, along with any image processing algorithms necessary to extract the x and y position of the particle relative to the electrodes.

In step (2), we model the dynamics of a particle that moves with the drift velocity due to the electrophoretic force only, neglecting any other forces or Brownian motion. As seen from Eqs. (1) and (3), this produces a velocity that is only a function of the gradient of the electric potential Φ at the location of the particle. What remains is to relate the electric potential Φ to the voltages of the electrodes. As mentioned previously, this can be done using a numerical solution to the Laplace equation. However, we found it more convenient to instead use the charge simulation method47–4947. N. H. Malik, IEEE Trans. Electr. Insul. 24, 3 (1989). https://doi.org/10.1109/14.1986148. M. M. A. Salama and R. Hackam, IEEE Power Eng. Rev. PER-4, 25 (1984). https://doi.org/10.1109/MPER.1984.552621149. H. Singer, H. Steinbigler, and P. Weiss, IEEE Trans. Power Appar. Syst. PAS-93, 1660 (1974). https://doi.org/10.1109/TPAS.1974.293898 to relate the electrode voltage ϕi to an effective charge qeff,i at the electrode location. There are two primary advantages to this approach. (1) We can make use of an analytical relation between the effective charge and the potential that simplifies the calculation of the electric field. (2) There is a useful analogy that we wish to highlight between the effective charge qeff,i of an electrode in electrokinetic micro-actuation and the reaction rate at a point in space in diffusiophoretic micro-actuation. As detailed in Sec. III of the supplementary material, the electrophoretic velocity using the charge simulation method is given by v^elec(t)=−μe2πϵ∑i=0Nprobeqeff,i(t)ri2ri,(7)where ϵ is the electric permittivity of the medium; ri=X−Xprobe, where Xprobe is the location of the effective charge qeff,i; and ri is the magnitude of ri.In step (3), we define a reference trajectory and modify it using a feedback term to correct for thermal noise and sensing and modeling inaccuracies. First, we pre-define a reference trajectory xref(t)=[xref(t),yref(t)]T for the colloidal particle to follow. We also take the derivative of the pre-defined trajectory, x˙ref(t), which is the open-loop velocity that our controller must produce. However, to track the trajectory in the presence of noise, measurement errors, and process errors, we add a feedback term kelec(X−xref), or the difference between the measured position and the reference position, to ensure that the particle actually follows the desired trajectory. Thus, we define vdes≡x˙ref−kelec(X−xref)(8)as the desired velocity to produce using the optimization algorithm.In step (4), we minimize the difference between the predicted velocity v^elec from Eq. (7) and the desired velocity vdes from Eq. (8). To solve this optimization problem, we used a nonlinear optimization algorithm to calculate the electrode voltages by minimizing the cost function, Jelec(uelec)=welec,1‖v^elec−vdes‖2v02+welec,2‖uelec‖2qscale2,(9)where v0 is a velocity scale that will be presented in Sec. , qscale is a scale for the electric charge derived in Sec. I of the supplementary material, and uelec is the vector of control inputs defined as uelec=[qeff,1,qeff,2,…,qeff,Nprobe]T.(10)The weights welec,1 and welec,2 are tuning parameters of the method. The optimization was implemented in Python using the least_squares function from the scipy.optimize5050. P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, I. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, A. Vijaykumar, A. P. Bardelli, A. Rothberg, A. Hilboll, A. Kloeckner, A. Scopatz, A. Lee, A. Rokem, C. N. Woods, C. Fulton, C. Masson, C. Häggström,C. Fitzgerald, D. A. Nicholson, D. R. Hagen, D. V. Pasechnik, E. Olivetti, E. Martin, E. Wieser, F. Silva, F. Lenders, F. Wilhelm, G. Young, G. A. Price, G. L. Ingold, G. E. Allen, G. R. Lee, H. Audren, I. Probst, J. P. Dietrich, J. Silterra, J. T. Webber, J. Slavi.c, J. Nothman, J. Buchner, J. Kulick, J. L. Schönberger, J. V. de Miranda Cardoso, J. Reimer, J. Harrington, J. L. C. Rodríguez, J. Nunez-Iglesias, J. Kuczynski, K. Tritz, M. Thoma, M. Newville, M. K. Nummerer, M. Bolingbroke, M. Tartre, M. Pak, N. J. Smith, N. Nowaczyk, N. Shebanov, O. Pavlyk, P. A. Brodtkorb, P. Lee, R. T. McGibbon, R. Feldbauer, S. Lewis, S. Tygier, S. Sievert, S. Vigna, S. Peterson, S. More, T. Pudlik, T. Oshima, T. J. Pingel, T. P. Robitaille, T. Spura, T. R. Jones, T. Cera, T. Leslie, T. Zito, T. Krauss, U. Upadhyay, Y. O. Halchenko, and Y. V. Lazquez-Baeza, Nat. Methods 17, 261 (2020). https://doi.org/10.1038/s41592-019-0686-2 library, which uses the “Trust Region Reflective” algorithm. The voltages were constrained to be positive to make a better comparison with our chemical controller results. The simulations were performed with a variety of initial positions, domain sizes, and paths. A subset of these simulations will be presented in Sec. .There are two major differences between our methods and the methods for electric field steering used previously by Chaudhary and Shapiro.55. S. Chaudhary and B. Shapiro, IEEE Trans. Control Syst. Technol. 14, 669 (2006). https://doi.org/10.1109/TCST.2006.876636 First, in step (2), they used a finite element solution to the Laplace equation to create a dynamic model instead of using the charge simulation method. Our choice to use the charge simulation method added a small amount of modeling error, but it also made the dynamic model less expensive to compute and will provide a clearer analogy with our chemical gradient controller. Second, in step (4), instead of using nonlinear optimization, they used a minimum norm least-squares algorithm to calculate the voltages of the electrodes in a linearized system of equations, and they dealt with ill-conditioning using a singular value decomposition. While their minimum norm solution is simpler and faster to compute, we used nonlinear optimization to be consistent with our chemical controller, which cannot be computed using a minimum norm solution.

C. Diffusiophoretic micro-actuation control algorithm

Similar to the case of electrokinetic micro-actuation, we use optimal control with diffusiophoretic micro-actuation to minimize the difference between the diffusiophoretic velocity v^diff from a dynamic model developed in the following paragraphs and the desired velocity vdes from Eq. (8). However, unlike the near-instantaneous changes in electric potentials, changes in concentration are due to diffusive transport processes, which take time to propagate through the solution. This leads to a significant delay between the chemical signal from the source and the motion of the particle, necessitating a more sophisticated control strategy.To deal with this delay, we use nonlinear model predictive control (NMPC). In NMPC, one specifies a cost function J that includes information about the current state and model predictions of the future states of the particle.36–5236. E. Camacho, C. Bordons, and C. Alba, Model Predictive Control, Advanced Textbooks in Control and Signal Processing (Springer, London, 2004).51. J. Pannek and K. Worthmann, IFAC Proc. Vol. 44, 7969 (2011). https://doi.org/10.3182/20110828-6-IT-1002.0091652. S. Qin and T. A. Badgwell, Control Eng. Pract. 11, 733 (2003). https://doi.org/10.1016/S0967-0661(02)00186-7 One then uses a nonlinear optimization algorithm to minimize this cost function and derive control commands for the concentration sources. The model prediction is necessarily limited to a finite future time, called the event or prediction horizon, resulting in an optimization over a set number of future steps Npred. The choice of the size of the prediction horizon results in a tradeoff between speed and accuracy: a larger prediction horizon may give better performance, but it also increases the computational cost of the optimization problem.5151. J. Pannek and K. Worthmann, IFAC Proc. Vol. 44, 7969 (2011). https://doi.org/10.3182/20110828-6-IT-1002.00916Figure 3 shows an example of how NMPC can be used to steer a colloid. Panel (a) shows the trajectory X(t) and dynamic model prediction x^diff(t) of a colloid as it is steered along a circular reference trajectory xref(t). Panel (b) shows the optimal inputs (i.e., reaction rates) for four probes that surround the particle, obtained by using an optimization algorithm to minimize the cost function J. The key aspect of NMPC demonstrated here is that the input reaction rates in part (b) are determined such that the dynamic model prediction in part (a) lies along the target trajectory. Also, note that noise from Brownian motion may move the particle off of the target trajectory, but feedback will correct such errors.Given this overview, we now describe the specifics of this method, starting with the cost function. We would like a cost function that, when minimized, drives the particle to the reference trajectory while simultaneously minimizing the amount of source concentration needed to actuate motion. A cost function based on the commonly used concept of a squared error3636. E. Camacho, C. Bordons, and C. Alba, Model Predictive Control, Advanced Textbooks in Control and Signal Processing (Springer, London, 2004). that accomplishes this goal is given by Jdiff()=∑k=1Npred+1wdiff,1‖v^diff,k−vdes,k‖2v02+∑k=1Npred+1wdiff,2‖uk‖2gscale2,(11)where v^diff,k is the velocity predicted by the dynamic model for k time steps into the future, vdes,k is the desired velocity for time step k, given in Eq. (12), v0 is a velocity scale that will be presented in Sec. , gscale is a reaction rate scale derived in Sec. II of the supplementary material, and uk=[u1k,u2k,…,uNprobek]T is a vector of Nprobe reaction rates at time step k. The weights wdiff,1 and wdiff,2 are tuning parameters of the method.3636. E. Camacho, C. Bordons, and C. Alba, Model Predictive Control, Advanced Textbooks in Control and Signal Processing (Springer, London, 2004). Generating these tuning parameters is a key step in developing optimal control of the system and will be discussed below. As desired, the cost function has two terms: the first drives the dynamic model to the desired trajectory and the second minimizes the concentration needed to actuate the system.For the desired velocity vdes,k, we modify Eq. (8) to be a function of the time step k, giving vdes,k=x˙ref(tk)−kdiff(x^diff,k−xref(tk)).(12)The main modification we have made is that we now use x^diff,k, the position predicted from the dynamic model for time step k, instead of the measured position X. Since the initial prediction is set equal to the measured position, or x^diff,1=X, this equation reduces to Eq. (8) when k=1.

In addition to specifying a cost function, our NMPC method requires a dynamic model for the motion of the particle. Rather than use a computationally expensive numerical solution (e.g., with finite differences), we develop here a method based on Green’s function for the diffusion equation that is analogous to the charge simulation method for the Laplace equation described previously. As we will see, a primary advantage of this method is that the concentration gradient ∇C only needs to be evaluated at the colloid location, rather than finding the solution to the unsteady diffusion problem in all space. However, the efficiency of this approach is attenuated by the history dependence of Green’s function, which causes the method to become increasingly costly with time.

For the dynamic model, we assume that the only forces acting on the particle are diffusiophoretic forces v^diff=Fdiffγ=μd∇C(X,t),(13)neglecting thermal noise. The concentration C(X,t) is again given by the reaction–diffusion equation provided in Eq. (6), where the reaction rate term G(x,y,t) is a function of both time and space.In an infinite domain, Eq. (6) has an analytical Green’s function solution,5353. R. Haberman, Applied Partial Differential Equations: With Fourier Series and Boundary Value Problems, 4th ed. (Prentice Hall, 2004). C(X,t)=∫0t∫1[4πDs(t−τ)]n/2e−‖X−σ‖24Ds(t−τ)G(σ,τ)dnσdτ+∫1(4πDst)n/2e−‖X−σ‖24DstC(X,0)dnσ,(14)where n is the number of spatial dimensions and σ and τ are dummy integration variables for space and time, respectively. Boundary correction terms can be added using the method of images, but we chose to neglect these corrections in favor of a simpler and faster dynamic model and anticipate that feedback control will correct for this model simplification.We are interested in the case of a 2D domain, which allows us to specify that n=2. We also assume a solution initially devoid of solute, i.e., C(