1. Introduction
Manipulation of microparticles suspended in fluids has relevance to several applications, including targetted drug delivery (Nelson, Kaliakatsos & Abbott Reference Nelson, Kaliakatsos and Abbott2010; Li et al. Reference Li, Esteban-Fernández de Ávila, Gao, Zhang and Wang2017; Zhou et al. Reference Zhou, Mayorga-Martinez, Pané, Zhang and Pumera2021; Ezike et al. Reference Ezike2023), environmental remediation (Wang, Khezri & Pumera Reference Wang, Khezri and Pumera2016), cell sorting (Bhagat et al. Reference Bhagat, Bow, Hou, Tan, Han and Lim2010; Wang, Jalikop & Hilgenfeldt Reference Wang, Jalikop and Hilgenfeldt2011), assisted fertilization (Fishel et al. Reference Fishel, Timson, Green, Hall, Dowell and Klentzeris1993) and microassembly (Ghadiri et al. Reference Ghadiri, Weigel, Esen and Ostendorf2012; Agnus et al. Reference Agnus2013).
Some common mechanisms for transporting large collections of particles in microfluidic devices are using pressure-driven fluid flow along channels, electrokinetic effects and acoustic streaming (Chakraborty & Chakraborty Reference Chakraborty and Chakraborty2010; Wu et al. Reference Wu, Ozcelik, Rufo, Wang, Fang and Jun Huang2019). It has also been shown that spatially and temporally patterned fluid flows can be generated in microfluidic chambers through buoyancy and electrokinetic effects associated with chemical reactions (Sengupta et al. Reference Sengupta, Patra, Ortiz-Rivera, Agrawal, Shklyaev, Dey, Córdova-Figueroa, Mallouk and Sen2014; Ortiz-Rivera et al. Reference Ortiz-Rivera, Shum, Agrawal, Sen and Balazs2016; Niu et al. Reference Niu, Kreissl, Brown, Rempfer, Botin, Holm, Palberg and de Graaf2017; Shum & Balazs Reference Shum and Balazs2018) or by harnessing bacterial or artificial cilia carpets (Darnton et al. Reference Darnton, Turner, Breuer and Berg2004; Kim et al. Reference Kim, Cheang, Kim, Ali and Kim2015). The motion of individual particles subject to these effects could be dependent on the size or other properties of the particle, so there is some control over at least the direction and speed of transport, but these mechanisms are not typically used for fine manipulation of individual particles along specific, arbitrary paths.
Precise manipulation of individual particles suspended in fluids is challenging, but can be achieved mechanically using micropipettes (Zhang et al. Reference Zhang, Yuan, Yue, Liu, Liu and Sun2024). Rotating externally applied magnetic fields are commonly used as a non-invasive method to produce rotational motion of magnetic particles, possibly leading to translational propulsion depending on the shape of the particle or proximity to a wall (Ghosh & Fischer Reference Ghosh and Fischer2009; Khalil et al. Reference Khalil, Keuning, Abelmann and Misra2012; Martinez-Pedrero et al. Reference Martinez-Pedrero, Navarro-Argemí, Ortiz-Ambriz, Pagonabarraga and Tierno2018). Optical tweezers use focused laser beams to exert optical forces on particles, allowing cells and other microparticles to be held in a trap, which can be moved to manipulate particles along arbitrary paths (Ghadiri et al. Reference Ghadiri, Weigel, Esen and Ostendorf2012; Bradac Reference Bradac2018; Polimeno et al. Reference Polimeno2018; Bunea & Glückstad Reference Bunea and Glückstad2019; Jamil, Pokharel & Park Reference Jamil, Pokharel and Park2022).
The problem of swimming by body deformations, without the application of external forces or torques, has been studied under the framework of geometric control theory (see, e.g. Agrachev & Sachkov Reference Agrachev and Sachkov2004) in low- and high-Reynolds-number flows (Shapere & Wilczek Reference Shapere and Wilczek1989; San Martín, Takahashi & Tucsnak Reference San Martín, Takahashi and Tucsnak2007; Chambrion & Munnier Reference Chambrion and Munnier2012; Hatton & Choset Reference Hatton and Choset2013; Koens & Lauga Reference Koens and Lauga2021). Similar approaches have been used to prove controllability of a single non-spherical swimmer by an externally applied torque in an unbounded fluid (Buzhardt & Tallapragada Reference Buzhardt and Tallapragada2020). For spherical particles in unbounded fluid, an external torque produces only rotational motion so translational motion is not attainable. For a pair of such microrotors, each can generate translational motion of the other, but Buzhardt, Fedonyuk & Tallapragada (Reference Buzhardt, Fedonyuk and Tallapragada2018b) showed that their positions are still not controllable in an unbounded two-dimensional fluid domain, while they are in a confined domain (see also Buzhardt et al. Reference Buzhardt, Fedonyuk, Sudarsanam and Tallapragada2018a).
In addition to controlling the motion of magnetic particles through magnetic fields, it has been shown theoretically (Buzhardt & Tallapragada Reference Buzhardt and Tallapragada2021) that such particles can be used as robots or active agents to arbitrarily manipulate another, passive particle through hydrodynamic interactions. Or et al. (Reference Or, Vankerschaver, Kelly, Murray and Marsden2009) showed that arbitrarily small bounded motions of a cylinder could be used to control the motion of a passive tracer particle in a two-dimensional domain, noting the same geometric structure for inviscid fluids as for Stokes flow. Walker et al. (Reference Walker, Ishimoto, Gaffney and Moreau2022) studied a similar system where an active spherical particle controlled by either external forces or external torques is used to direct the motion of a passive spherical particle in three-dimensional Stokes flow. In particular, they showed that the Lie brackets generated by the vector fields controlling the velocity of the active particle spanned the full six-dimensional configuration space for one active and one passive particle; controllability follows owing to the Rashewsky–Chow theorem (Agrachev & Sachkov Reference Agrachev and Sachkov2004, Theorem 5.9).
In the current work, we build on the previous theoretical studies, exploring a strategy for manipulating passive particles suspended in a viscous fluid through hydrodynamic interactions with one active particle directly controlled by other means (for example, optical tweezers or magnetic forces).
The general set-up for our model is that a collection of neutrally buoyant particles is suspended in an incompressible Newtonian fluid. Since we are primarily interested in microfluidic systems, it is natural to adopt the equations of incompressible Stokes flow to characterize the behaviour of the fluid in low-Reynolds-number applications. For particles separated by a distance 
$r=100\ \mathrm {\mu }{\rm m}$ and moving with characteristic speeds 
$U=100\ \mathrm {\mu }{\rm m}\ {\rm s}^{-1}$ in a fluid with the kinematic viscosity 
$\nu = 1\ {\rm mm}^2\ {\rm s}^{-1}$ (comparable to that of water at 
$20\,^\circ {\rm C}$), the Reynolds number is 
$Re =Ur/\nu =0.01$, for example.
We consider the case where one of these particles is directly controlled by external forces, so that its velocity and position are prescribed functions of time, and the remaining particles move passively in the fluid flow generated by the actively controlled one. We use the far-field expressions for the Stokes flow field produced by a translating rigid, spherical particle to determine the velocity and trajectory of passive particles in the fluid.
Using explicit constructions inspired by the work of Dal Maso, DeSimone & Morandotti (Reference Dal Maso, DeSimone and Morandotti2015), we prove total controllability of systems consisting of one active and either one or two passive particles. That is, such particles can be moved from arbitrary initial positions to arbitrary final positions in unbounded three-dimensional space, provided that the particles are far apart from each other in these configurations so that the far-field approximation is valid.
We provide a strategy that breaks down the desired displacements into a sequence of steps that can be achieved by iteratively applying elementary moves, in which a passive particle moves in the radial direction or in the polar direction around the active particle, for example. We extend the controllability result to two passive particles, proposing a separate strategy for this case. Using numerical solutions, we test the accuracy of asymptotic expressions for displacements predicted for our elementary moves and also assess the errors associated with adopting the far-field hydrodynamic approximation.
Our results are a step towards the more general problem of independently manipulating an arbitrary number of passive particles using a small number of control variables. Although the strategies we discuss in the current work do not readily extend to larger numbers of particles, the framework can be generalized to describe such systems in a straightforward manner.
To motivate the study of controlling passive particles through one or more actively manipulated particles, we note that optical traps can be used for multiple particles simultaneously, but this requires more complicated experimental procedures, making it inconvenient as a method for controlling many particles. Moreover, optical tweezers face challenges and limitations from heating effects, the dependence on the size and dielectric properties of the particle being trapped, and degraded laser focus when particles are located deeper in the fluid, far from the bounding glass surface (Melzer & McLeod Reference Melzer and McLeod2018; Jamil et al. Reference Jamil, Pokharel and Park2022; Malinowska et al. Reference Malinowska, van Mameren, Peterman, Wuite and Heller2024).
In contrast, exploiting the hydrodynamic interaction has several advantages: forces propagate throughout the whole fluid and influence particles that are far away from one that is moving (or being moved), overcoming both the spatial limitation of optical tweezers and the need to have one tweezer for each particle to be displaced. Another advantage is that, in the far-field approximation, the size of the passive particles does not enter the equations of motion, so, in principle, one can control arbitrarily large particles (provided other modelling assumptions remain valid).
The paper is organized as follows: in § 2, we present the mathematical formulation of the problem; in § 3, we construct the elementary and compound moves that will be used in § 4 to prove controllability of our system of one active particle and one or two passive ones. In § 5, we discuss the errors due to finite amplitudes and separations in comparison to the theoretically predicted trajectories. Finally, in § 6, we offer an overview of the results we obtained and an outlook for future research.
2. Mathematical formulation
We introduce the general set-up for a system of spherical particles immersed in a viscous fluid. Of these, one is an active particle, whose velocity is directly prescribed, and the remaining 
$M$ are passive particles, whose motions are determined by the interaction with the active one.
We let 
$a>0$ be the radius of the active particle and we denote by 
$t\mapsto \boldsymbol {x}(t)$ its position in space at time 
$t$ and by 
$t\mapsto \dot {\boldsymbol {x}}(t)$ its velocity at time 
$t$. Analogously, we denote by 
$a_j>0$, 
$t\mapsto \boldsymbol {y}_j(t)$ and 
$t\mapsto \dot {\boldsymbol {y}}_j(t)$, for every 
$j=1,\ldots,M$, the radius, position and velocity, respectively, of the 
$j$th passive particle.
Assuming that the Reynolds number is small enough that inertial effects may be neglected, the fluid flow is governed by the equations of incompressible Stokes flow,
\begin{equation} \boldsymbol{\nabla} p -\mu \nabla^2\boldsymbol{u}=\boldsymbol{0}, \quad \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} = 0, \end{equation}
where 
$\boldsymbol {u}$ is the velocity field, 
$p$ is the pressure field and 
$\mu$ is the dynamic viscosity of the fluid. We assume that the velocity field vanishes at infinity and satisfies the no-slip boundary conditions on the surfaces of the particles, namely,
\begin{equation} \boldsymbol{u}(\boldsymbol{z}) = \begin{cases} \boldsymbol{U}_0 + \boldsymbol{\varOmega}_0\times (\boldsymbol{z} - \boldsymbol{x}) & \text{for }\ \lVert\boldsymbol{z}-\boldsymbol{x}\rVert = a,\\ \boldsymbol{U}_j + \boldsymbol{\varOmega}_j\times (\boldsymbol{z} - \boldsymbol{y}_j) & \text{for }\ \lVert\boldsymbol{z}-\boldsymbol{y}_j\rVert = a_j,\ j= 1,\ldots,M, \end{cases} \end{equation}
where 
$\boldsymbol {U}_j$ and 
$\boldsymbol {\varOmega }_j$ (
$\kern 1.5pt j=0,\ldots,M$) are the translational and rotational velocities, respectively, of the active and passive particles. We denote by 
$t\mapsto \boldsymbol {F}(t)$ the force that the active particle exerts, at time 
$t$, on the surrounding fluid and we assume that all passive particles are force-free. All particles, whether active or passive, are torque-free.
By linearity of the equations of Stokes flow, the relationship between active forces and the velocities of the particles, in the absence of background flows, are generically described by (Kim & Karrila Reference Kim and Karrila1991)
\begin{equation} \begin{pmatrix}
\boldsymbol{U}_j \\ \boldsymbol{\varOmega}_j \end{pmatrix}
= \begin{pmatrix} \boldsymbol{\mathsf{M}}_{j} \\ \boldsymbol{\mathsf{N}}_{j}
\end{pmatrix} \boldsymbol{F}\quad\text{for }\
j=0,\ldots,M. \end{equation}
The quantities 
$\boldsymbol{\mathsf{M}}_{j}$ (
$\kern 1.5pt j=0,\ldots,M$) are the mobility tensors for the translational velocities of the active and passive particles, respectively, due to the force on the active particle, and 
$\boldsymbol{\mathsf{N}}_{j}$ (
$\kern 1.5pt j=0,\ldots,M$) are the mobility tensors for the rotational velocities of the active and passive particles, respectively, due to the force on the active particle. In general, the mobility tensors depend on the relative positions of all particles and it is not possible to obtain a closed-form expression for them. By symmetry of the spherical particles, the mobility tensors are independent of the orientations of the particles. We focus on the problem of controlling the particle positions without regard to their orientations, hence, the rotational velocities need not be considered.
Let 
$\boldsymbol {d}_{j}\unicode{x2254} \boldsymbol {y}_j-\boldsymbol {x}$ (
$\kern 1.5pt j=1,\ldots,M$) be the displacement vector of the 
$j$th passive particle from the active one. Assuming that all of the pairwise distances 
$r_{j}\unicode{x2254} \lVert \boldsymbol {d}_{j}\rVert$ are large compared with all of the particle radii 
$a$ and 
$a_j$, as well as the mutual distances between the passive particles, we use the far-field approximation for the translational mobility tensors (Zuk et al. Reference Zuk, Wajnryb, Mizerski and Szymczak2014; Graham Reference Graham2018) given by
\begin{equation} \boldsymbol{\mathsf{M}}_{j}=\begin{cases} \displaystyle \frac{1}{6{\rm \pi}\mu a}\operatorname{\mathbb{1}} & \text{for }\ j=0,\\[10pt] \displaystyle \left(1+\frac{a^2+a_j^2}{6}\nabla^2\right)\mathcal{G}(\boldsymbol{d}_j) & \text{for }\ j=1,\ldots,M, \end{cases} \end{equation}
where the function 
$\mathcal {G}\colon \mathbb {R}^3\setminus \{\boldsymbol {0}\}\to \mathbb {R}^{3\times 3}_{sym}$ defined by
\begin{equation} \mathcal{G}(\boldsymbol{d}) \unicode{x2254} \frac{1}{8{\rm \pi} \mu}\left(\frac{\operatorname{\mathbb{1}}}{r}+\frac{\boldsymbol{d}\otimes \boldsymbol{d}}{r^3}\right) \end{equation}
is the Stokeslet Green's function for the Stokes equation (2.1a,b) such that 
$\mathcal {G}(\boldsymbol {d})\boldsymbol {F}$ is the velocity at 
$\boldsymbol {d}$ due to a singular force 
$\boldsymbol {F}$ applied at the origin. Equation (2.4) is accurate up to 
$O(r_j^{-3})$ and can be extended to higher orders of accuracy by the method of reflections (Kim & Karrila Reference Kim and Karrila1991). Note that, to this order of accuracy, the passive particles do not affect the velocities of other particles (active or passive).
The matrix 
$\boldsymbol{\mathsf{M}}_0$ is evidently invertible and its inverse 
$\boldsymbol{\mathsf{R}}_0=\boldsymbol{\mathsf{M}}_0^{-1}=6{\rm \pi} \mu a\operatorname {\mathbb{1}}$ is the resistance matrix describing the linear relationship between forces applied to the fluid and the velocities of the particles, 
$\boldsymbol {F} = \boldsymbol{\mathsf{R}}_0\boldsymbol {U}_0$. The equations of motion for our system of active and passive particles are then
\begin{equation} \begin{cases}
\dot{\boldsymbol{x}}(t)=\boldsymbol{u}(t), & \\
\displaystyle \dot{\boldsymbol{y}}_j(t)=
\boldsymbol{\mathsf{M}}_{j}(t)\boldsymbol{F}(t) = \bar{\boldsymbol{\mathsf{M}}}_{j}(t)\boldsymbol{u}(t) &
\text{for }\ j=1,\ldots,M, \end{cases}
\end{equation}
where 
$\bar {\boldsymbol{\mathsf{M}}}_{j}(t) = \boldsymbol{\mathsf{M}}_{j}(t) \boldsymbol{\mathsf{R}}_0(t)=6{\rm \pi} \mu a\boldsymbol{\mathsf{M}}_j(t)$. We further simplify the equations by retaining only the leading order terms, namely, those of order 
$O(r_j^{-1})$. Notice that, with this approximation, the radii of the passive particles do not enter the system of equations, which becomes
\begin{equation} \begin{cases} \dot{\boldsymbol{x}}(t)=\boldsymbol{u}(t), \\[4pt] \displaystyle \dot{\boldsymbol{y}}_j(t)=\frac{3a}{4}\left(\frac{\operatorname{\mathbb{1}}}{r_{j}(t)}+\frac{\boldsymbol{d}_{j}(t)\otimes \boldsymbol{d}_{j}(t)}{r_{j}(t)^3}\right)\boldsymbol{u}(t) = \boldsymbol{\mathsf{G}}(\boldsymbol{d}_j(t))\boldsymbol{u}(t) & \text{for }\ j=1,\ldots,M, \end{cases} \end{equation}
where 
$\boldsymbol{\mathsf{G}}(\boldsymbol {d}) \unicode{x2254} ({3a}/{4})({\operatorname {\mathbb{1}}}/{r}+{\boldsymbol {d}\otimes \boldsymbol {d}}/{r^3})$ is a 
$3\times 3$ matrix. We remark that the passive particles move as tracers or point particles in the flow field induced by the moving active particles. Equation (2.7) can be written in matrix form as
\begin{equation} \begin{pmatrix}
\dot{\boldsymbol{x}} \\ \dot{\boldsymbol{y}}_1 \\ \vdots \\
\dot{\boldsymbol{y}}_M
\end{pmatrix}=\boldsymbol{\mathsf{H}}\boldsymbol{u},\quad
\text{where}\ \boldsymbol{\mathsf{H}}=\begin{pmatrix}
\operatorname{\mathbb{1}}_{3\times3} \\ \boldsymbol{\mathsf{G}}_{1}
\\ \vdots \\ \boldsymbol{\mathsf{G}}_{M}\end{pmatrix}, \quad
\boldsymbol{\mathsf{G}}_{j} = \boldsymbol{\mathsf{G}}(\boldsymbol{d}_j)\quad \text{for
}\ j=1,\ldots,M. \end{equation}
3. Elementary and compounds moves
In this section, we construct the elementary and compound moves, which are the building blocks for proving the controllability results in § 4. Below, we use 
$\boldsymbol {e}_k, k=1,2,3$ to denote the standard basis vectors in 
$\mathbb {R}^3$ and we use 
$\boldsymbol {0}$ to denote the zero vector in 
$\mathbb {R}^3$.
3.1. Elementary moves and Lie brackets of vector fields
We describe three elementary classes of control functions from which we will construct strategies to move the active and passive particles from arbitrary initial positions to arbitrary target positions. Since the equations for the passive particles are decoupled, i.e. the velocity of the 
$j$th passive particle does not depend on the 
$i$th passive particle with 
$i\neq j$, the action of the active particle is the same on all the passive particles. For this reason, in what follows, the three elementary classes of control functions will be described for the case 
$M=1$ in (2.7) and (2.8). Let 
$\boldsymbol {h}_k\in \mathbb {R}^{6}$ be the 
$k$th column of 
$\boldsymbol{\mathsf{H}}$ in (2.8) for 
$M=1$ and 
$k=1,2,3$. The first three components of 
$\boldsymbol {h}_k$ represent the velocity of the active particle and the last three components represent the velocity of the passive particle when the control 
$\boldsymbol {u}=\boldsymbol {e}_k$ is applied.
3.1.1. Zeroth-order control
Given an amplitude 
$\varepsilon >0$ and time duration 
$\varDelta t>0$, consider the zeroth-order (constant) control
\begin{equation} \boldsymbol{u}(t) = \boldsymbol{u}_{\Delta t}^{\varepsilon,\,\boldsymbol{h}_k}(t) \unicode{x2254}\begin{cases} \displaystyle \frac{\varepsilon}{\Delta t} \boldsymbol{e}_k, & 0\leqslant t < \Delta t,\\[6pt] \boldsymbol{0}, & \text{otherwise.} \end{cases} \end{equation}
The net displacements of the active and passive particles over the time interval 
$[0,\Delta t]$ with initial positions 
$\boldsymbol {x}(0)=\boldsymbol {x}^\circ$ and 
$\boldsymbol {y}(0)=\boldsymbol {x}^\circ +\boldsymbol {d}$ (with 
$\boldsymbol {d}\in \mathbb {R}^3\setminus \{\boldsymbol {0}\}$), respectively, are
\begin{align} \varDelta^{\varepsilon,\boldsymbol{h}_k}(\boldsymbol{d}) &= \begin{pmatrix} \varDelta_{\boldsymbol{x}}^{\varepsilon,\boldsymbol{h}_k}(\boldsymbol{d})\\ \varDelta_{\boldsymbol{y}}^{\varepsilon,\boldsymbol{h}_k}(\boldsymbol{d}) \end{pmatrix}\unicode{x2254} \begin{pmatrix} \boldsymbol{x}(\Delta t)-\boldsymbol{x}(0)\\ \boldsymbol{y}(\Delta t)-\boldsymbol{y}(0) \end{pmatrix} = \begin{pmatrix} \varepsilon\,\boldsymbol{e}_k \\[4pt] \dfrac{3a\varepsilon}{4r}\left(\boldsymbol{e}_k+\dfrac{d_k \boldsymbol{d}}{r^2}\right) + \boldsymbol{O}(\varepsilon^2) \end{pmatrix} \nonumber\\ &= \varepsilon \boldsymbol{h}_k(\boldsymbol{d}) + \begin{pmatrix} \boldsymbol{0}\\ \boldsymbol{O}(\varepsilon^2) \end{pmatrix}, \end{align}
where 
$\boldsymbol {O}(\boldsymbol {\varepsilon ^2})$ denotes a vector in 
$\mathbb {R}^3$ with components that are 
$O(\varepsilon ^2)$ as 
$\varepsilon \to 0$. This zeroth-order control is primarily used to control the position of the active particle since its velocity is directly prescribed by the control. The passive particle will also move, due to the flow field generated by the translating active particle and we can compute the trajectory of the passive particle according to (3.2), see figure 1(a).
Figure 1. Streamline plots of the displacement vector field of a passive particle, given by components 4–6 of the Lie brackets, for (a) the 
$\boldsymbol {h}_1$ zeroth-order control, (b) the 
$[\boldsymbol {h}_1,\boldsymbol {h}_2]$ first-order control, (c) the 
$x$–
$y$ plane of the 
$[\boldsymbol {h}_2,[\boldsymbol {h}_1,\boldsymbol {h}_2]]$ second-order control and (d) the 
$x$–
$z$ plane of the 
$[\boldsymbol {h}_2,[\boldsymbol {h}_1,\boldsymbol {h}_2]]$ second-order control. In all cases, we use an active particle of radius 
$a=1$ and the plots are in the same plane as the active particle, which is located at the origin. By symmetry, there are no out-of-plane components of displacements for the passive particle. The colour scale indicates the base 10 logarithm of the magnitude.
To displace the passive particle without any net displacement of the active particle, we construct higher-order controls.
3.1.2. First-order control
Consider a control that moves the active particle around a square loop with sides of length 
$\varepsilon$ in the 
$\boldsymbol {e}_k$ and 
$\boldsymbol {e}_l$ directions, namely,
\begin{equation} \boldsymbol{u}_{\Delta t}^{\varepsilon, [\boldsymbol{h}_k,\boldsymbol{h}_l]}(t) \unicode{x2254} \boldsymbol{u}_{\Delta t/4}^{\varepsilon, \boldsymbol{h}_k}(t) + \boldsymbol{u}_{\Delta t/4}^{\varepsilon, \boldsymbol{h}_l}\left(t-\frac{\Delta t}4\right) - \boldsymbol{u}_{\Delta t/4}^{\varepsilon, \boldsymbol{h}_k}\left(\frac{3\Delta t}4-t\right) - \boldsymbol{u}_{\Delta t/4}^{\varepsilon, \boldsymbol{h}_l}(\Delta t - t) \end{equation}
for 
$k,l=1,2,3$. The inversion of the time variable in the last two terms represents performing the reverse control of the first two terms in the sum. The controls here are constant over each subinterval but in later examples, it will be necessary to make this distinction. Explicitly, the function in (3.3) can be expressed as
\begin{equation} \boldsymbol{u}_{\Delta t}^{\varepsilon, [\boldsymbol{h}_k,\boldsymbol{h}_l]}(t) = \begin{cases} \displaystyle\frac{4\varepsilon}{\Delta t}\boldsymbol{e}_k, & \displaystyle0\leqslant t < \frac{\Delta t}4,\\[10pt] \displaystyle\frac{4\varepsilon}{\Delta t}\boldsymbol{e}_l, & \displaystyle\frac{\Delta t}4 \leqslant t < \frac{\Delta t}2,\\[10pt] \displaystyle-\frac{4\varepsilon}{\Delta t}\boldsymbol{e}_k, & \displaystyle\frac{\Delta t}2 \leqslant t < \frac{3\Delta t}4,\\[10pt] \displaystyle-\frac{4\varepsilon}{\Delta t}\boldsymbol{e}_l, & \displaystyle\frac{3\Delta t}4 \leqslant t < {\Delta t},\\ 0, & \text{otherwise.} \end{cases} \end{equation}
The net displacements, for small 
$\varepsilon$, are
\begin{equation} \varDelta^{\varepsilon,[\boldsymbol{h}_k,\boldsymbol{h}_l]}(\boldsymbol{d}) = \begin{pmatrix} \varDelta_{\boldsymbol{x}}^{\varepsilon,[\boldsymbol{h}_k,\boldsymbol{h}_l]}(\boldsymbol{d})\\ \varDelta_{\boldsymbol{y}}^{\varepsilon,[\boldsymbol{h}_k,\boldsymbol{h}_l]}(\boldsymbol{d}) \end{pmatrix} = \varepsilon^2 [\boldsymbol{h}_k, \boldsymbol{h}_l](\boldsymbol{d}) + \begin{pmatrix} \boldsymbol{0}\\ \boldsymbol{O}(\varepsilon^3) \end{pmatrix}, \end{equation}
where 
$[\boldsymbol {h}_k,\boldsymbol {h}_l] \unicode{x2254} \boldsymbol {h}_k\boldsymbol {\cdot }\boldsymbol {\nabla }\boldsymbol {h}_l - \boldsymbol {h}_l\boldsymbol {\cdot }\boldsymbol {\nabla }\boldsymbol {h}_k$ is the first-order Lie bracket of the vector fields 
$\boldsymbol {h}_k$ and 
$\boldsymbol {h}_l$. The Lie bracket evaluates to
\begin{equation} [\boldsymbol{h}_k,\boldsymbol{h}_l](\boldsymbol{d}) = \begin{pmatrix} \boldsymbol{0}\\ \dfrac{3a}{2r^3}\left(1-\dfrac{9a}{8r}\right)(d_k\boldsymbol{e}_l - d_l\boldsymbol{e}_k) \end{pmatrix} = \begin{pmatrix} \boldsymbol{0}\\ \pmb{\omega}\times\boldsymbol{d} \end{pmatrix}, \end{equation}
where 
$\pmb {\omega }(r) = ({3a}/{2r^3})(1-{9a}/{8r})\boldsymbol {e}_k\times \boldsymbol {e}_l$ and 
$\boldsymbol {0}$ is the three-dimensional zero vector. Hence, for 
$k\neq l$ and small 
$\varepsilon$, the control 
$\boldsymbol {u}_{t_0, \Delta t}^{\varepsilon,[\boldsymbol {h}_k,\boldsymbol {h}_l]}$ results approximately in a rotation of the passive particle by an angle
\begin{equation} \Delta\theta^\varepsilon(r) = \frac{3a\varepsilon^2}{2r^3}\left(1-\frac{9a}{8r}\right) \end{equation}
about the axis passing through the active particle and perpendicular to 
$\boldsymbol {e}_k$ and 
$\boldsymbol {e}_l$. By construction, the active particle returns to its initial position, 
$\boldsymbol {x}(\Delta t)=\boldsymbol {x}(0)$. If 
$k=l$, then the net displacements are exactly zero as this is a time-reciprocal motion.
The interpretation of this result is that a particle forced to move around in a closed, square loop produces a net displacement field that is, to leading order, equivalent to a rotlet, see figure 1(b). Indeed, the time-averaged distribution of forces applied to the fluid over the interval 
$[0,\Delta t]$ corresponds to the sum of a Stokeslet dipole with force in the 
$\boldsymbol {e}_l$ direction and displacement in the 
$\boldsymbol {e}_k$ direction and a Stokeslet dipole with force in the 
$-\boldsymbol {e}_k$ direction and displacement in the 
$\boldsymbol {e}_l$ direction.
We refer to the control 
$\boldsymbol {u}_{\Delta t}^{\varepsilon, [\boldsymbol {h}_k,\boldsymbol {h}_l]}$ as the first-order control (and assume that 
$k\neq l$) since it corresponds to a first-order Lie bracket.
Since rotations preserve the distance 
$r$, we require another class of controls: one that generates net displacements of the passive particle in the radial direction, with respect to the active particle, without a net displacement of the active particle.
3.1.3. Second-order control
Consider second-order control functions of the form 
$\boldsymbol {u}(t) = \boldsymbol {u}_{\Delta t}^{\varepsilon, [\boldsymbol {h}_k,[\boldsymbol {h}_l, \boldsymbol {h}_m]]}(t)$ defined as in (3.3), replacing 
$\boldsymbol {h}_l$ with 
$[\boldsymbol {h}_l,\boldsymbol {h}_m]$. Since 
$[\boldsymbol {h}_l, \boldsymbol {h}_m]$ corresponds to a rotlet-like flow field if 
$l\neq m$, the Lie bracket 
$[\boldsymbol {h}_k,[\boldsymbol {h}_l, \boldsymbol {h}_m]]$ has the approximate form of a rotlet dipole, with axis 
$(\boldsymbol {e}_l\times \boldsymbol {e}_m)$ and displacement in the 
$\boldsymbol {e}_k$ direction, acting on the passive particle,
\begin{align}
&[\boldsymbol{h}_k,[\boldsymbol{h}_l,\boldsymbol{h}_m]](\boldsymbol{d})
\nonumber\\ &\quad= \begin{pmatrix}
\boldsymbol{0}\\
\dfrac{3a}{2r^3}\left(1-\dfrac{9a}{8r}\right)(\delta_{mk}\boldsymbol{e}_l
- \delta_{lk}\boldsymbol{e}_m) +
\dfrac{9a}{2r^3}\left(1-\dfrac{3a}{2r}\right)^2\dfrac{d_k}{r}\left(\dfrac{d_l\boldsymbol{e}_m
- d_m\boldsymbol{e}_l}{r}\right) \\ \quad -
\dfrac{9a^2}{8r^4}\left(1-\dfrac{9a}{8r}\right)\delta_{mk}\left(\boldsymbol{e}_l+\dfrac{d_l\boldsymbol{d}}{r^2}\right)
\end{pmatrix} \end{align}
\begin{align} &\quad= \begin{pmatrix} \boldsymbol{0}\\ \dfrac{3a}{2r^3}\left\{(\boldsymbol{e}_k\times(\boldsymbol{e}_l\times\boldsymbol{e_m}) +\dfrac{3d_k(\boldsymbol{e}_l\times\boldsymbol{e}_m)\times\boldsymbol{d}}{r^2}\right\} + \boldsymbol{O}(1/r^4) \end{pmatrix}, \end{align}see figure 1(c).
The net displacements of the active and passive particles are given, for small 
$\varepsilon$, by
\begin{equation} \varDelta^{\varepsilon,[\boldsymbol{h}_k,[\boldsymbol{h}_l,\boldsymbol{h}_m]]}(\boldsymbol{d}) = \begin{pmatrix} \varDelta_{\boldsymbol{x}}^{\varepsilon,[\boldsymbol{h}_k,[\boldsymbol{h}_l,\boldsymbol{h}_m]]}(\boldsymbol{d})\\ \varDelta_{\boldsymbol{y}}^{\varepsilon,[\boldsymbol{h}_k,[\boldsymbol{h}_l,\boldsymbol{h}_m]]}(\boldsymbol{d}) \end{pmatrix} = \varepsilon^3 [\boldsymbol{h}_k,[\boldsymbol{h}_l,\boldsymbol{h}_m]](\boldsymbol{d}) + \begin{pmatrix} \boldsymbol{0}\\ \boldsymbol{O}(\varepsilon^4) \end{pmatrix}. \end{equation}
Notice that if, for a given relative position of the passive particle with respect to the active particle, we choose a right-handed reference frame in which 
$\boldsymbol {d} = d_1\boldsymbol {e}_1$, then the control corresponding to 
$[\boldsymbol {h}_2,[\boldsymbol {h}_1,\boldsymbol {h}_2]]$ results in a passive particle displacement
\begin{equation} \varDelta_{\boldsymbol{y}}^{\varepsilon,[\boldsymbol{h}_2,[\boldsymbol{h}_1,\boldsymbol{h}_2]]}(d_1\boldsymbol{e}_1) = \frac{3a\varepsilon^3}{2r^3}\left(1-\frac{3a}{2r}\right)\left(1-\frac{9a}{8r}\right)\boldsymbol{e}_1 + \boldsymbol{O}(\varepsilon^4), \end{equation}
where 
$r=|d_1|$. To leading order in 
$\varepsilon$, this produces a displacement in the 
$\boldsymbol {e}_1$ direction.
In the more general configuration assuming only that 
$d_2 = 0$, we have the result that
\begin{equation} \varDelta_{\boldsymbol{y}}^{\varepsilon,[\boldsymbol{h}_2,[\boldsymbol{h}_1,\boldsymbol{h}_2]]}(d_1\boldsymbol{e}_1+d_3\boldsymbol{e}_3) = \frac{3a\varepsilon^3}{2r^3}\boldsymbol{e}_1 + \boldsymbol{O}(1/r^4) + \boldsymbol{O}(\varepsilon^4), \end{equation}
which implies that, to leading order in 
$1/r$ and 
$\varepsilon$, the displacement of passive particles in the plane 
$d_2=0$ is purely in the 
$\boldsymbol {e}_1$ direction and the magnitude of the displacement depends on the magnitude but not the direction of the vector 
$\boldsymbol {d}$. This is illustrated in figure 1(d).
3.2. Compound moves
In this section, we establish a set of manipulations that can be performed on a system of one active and two passive particles (i.e. 
$M=2$), assuming that their motion is governed by system (2.7), using the elementary moves discussed in § 3.1. Since this system is based on the far-field approximation, we will ensure that the particles always remain well separated, according to the following definition.
Definition 3.1 Let 
$R>0$ be the minimum separation we wish to maintain between any two particles. Given a number 
$M\in \mathbb {N}$ of passive particles, we define the state space of the system as the set of configurations that respect this minimum distance,
\begin{align} \mathcal{S}_R^{1+M}=\{(\boldsymbol{x},\boldsymbol{y}_1,\ldots,\boldsymbol{y}_M)\in\mathbb{R}^{3(1+M)}: \lVert {\boldsymbol{y}_i-\boldsymbol{x}} \rVert>R, \lVert {\boldsymbol{y}_j-\boldsymbol{y}_{i}} \rVert>R\ \text{for all}\; i,j=1,\ldots,M\}. \end{align}
We say that a configuration 
$(\boldsymbol {x}(t),\boldsymbol {y}_1(t),\ldots,\boldsymbol {y}_M(t))$ at time 
$t$ is well separated if it is in 
$\mathcal {S}_R^{1+M}$. We say that a solution, or trajectory, 
$(\boldsymbol {x},\boldsymbol {y}_1,\ldots,\boldsymbol {y}_M)\in (AC([0,T];\mathbb {R}^3))^{1+M}$ of system (2.7) is well separated if the configuration 
$(\boldsymbol {x}(t),\boldsymbol {y}_1(t),\ldots,\boldsymbol {y}_M(t)) \in \mathcal {S}_R^{1+M}$ for all times 
$t\in [0,T]$.
By 
$AC([0,T];\mathbb {R}^{3})$, we denote the class of 
$\mathbb {R}^3$-valued absolutely continuous functions defined on 
$[0,T]$.
For brevity, we will not reiterate the well-separated conditions in Propositions 3.2–3.7 that follow, but these conditions will be implied in all cases, namely, we will always move the particles from an initial configuration 
$(\boldsymbol {x}^\circ,\boldsymbol {y}_1^\circ,\boldsymbol {y}_2^\circ )\in \mathcal {S}_R^3$ to a final configuration 
$(\boldsymbol {x}^f,\boldsymbol {y}_1^f,\boldsymbol {y}_2^f)\in \mathcal {S}_R^3$ ensuring that the particles stay well separated at all times.
Proposition 3.2 (Equidistant to non-equidistant configurations)
The particles can be moved from any non-collinear initial configuration 
$(\boldsymbol {x}^\circ,\boldsymbol {y}_1^\circ,\boldsymbol {y}_2^\circ )$ with the two passive particles equidistant from the active particle, i.e. 
$\lVert {\boldsymbol {y}_1^\circ - \boldsymbol {x}^\circ } \rVert = \lVert {\boldsymbol {y}_2^\circ - \boldsymbol {x}^\circ } \rVert$, to a final configuration 
$(\boldsymbol {x}^f,\boldsymbol {y}_1^f,\boldsymbol {y}_2^f)$ satisfying 
$\lVert {\boldsymbol {y}_1^f - \boldsymbol {x}^f} \rVert \neq\, \lVert {\boldsymbol {y}_2^f - \boldsymbol {x}^f} \rVert$.
Proof. The goal can be achieved using the second-order control. We choose a right-handed orthonormal reference frame 
$\{\boldsymbol {e}_1,\boldsymbol {e}_2,\boldsymbol {e}_3\}$ in which 
$\boldsymbol {e}_3 = \boldsymbol {d}_1^\circ /\lVert {\boldsymbol {d}_1^\circ } \rVert$, 
$\boldsymbol {d}_1$ and 
$\boldsymbol {d}_2$ belong to the span of 
$\boldsymbol {e}_1$ and 
$\boldsymbol {e}_3$ so that 
${d}_{i,2}=0$ for 
$i=1,2$, and 
$\boldsymbol {e}_1\boldsymbol {\cdot } \boldsymbol {d}_2 > 0$ (see figure 2a). Applying the second-order control 
$\boldsymbol {u}_{\Delta t}^{\varepsilon, [\boldsymbol {h}_2,[\boldsymbol {h}_1,\boldsymbol {h}_2]]}$ causes no net displacement of the active particle and displaces each passive particle by the same distance in the 
$\boldsymbol {e}_1$ direction (to leading order in both 
$1/r$ and 
$\varepsilon$), according to (3.12). Moving both passive particles in the 
$\boldsymbol {e}_1$ direction breaks the symmetry and results in a final configuration with 
$\lVert {\boldsymbol {d}_2^f} \rVert > \lVert {\boldsymbol {d}_1^f} \rVert$. To leading order, the distance between passive particles is unchanged and the distance between each passive particle and the active particle is increased so the configuration remains well separated.
Figure 2. Schematic illustration of the compound moves for two passive particles. The objectives of these moves are to (a) produce a non-equidistant configuration, (b) produce a collinear configuration, (c) make collinear particles equidistant, (d) rotate equidistant collinear particles to a new orientation, (e) translate equidistant collinear particles together and (f) adjust the distance between equidistant collinear particles.
Proposition 3.3 (Arbitrary to collinear configurations)
The particles can be moved from any initial configuration 
$(\boldsymbol {x}^\circ,\boldsymbol {y}_1^\circ,\boldsymbol {y}_2^\circ )$ to a collinear final configuration with the active particle between the two passive particles, i.e. a configuration 
$(\boldsymbol {x}^f,\boldsymbol {y}_1^f,\boldsymbol {y}_2^f)$ satisfying 
$\boldsymbol {d}_1^f = -\alpha \boldsymbol {d}_2^f$ for some 
$\alpha > 0$.
Proof. Suppose that the initial configuration does not satisfy the desired final properties. We may also suppose without loss of generality that 
$\lVert {\boldsymbol {d}_1^\circ } \rVert > \lVert {\boldsymbol {d}_2^\circ } \rVert$. If this is not the case (i.e. if the two passive particles are exactly the same distance from the active particle), then we apply Proposition 3.2 to produce a new configuration in which the first passive particle is farther from the active particle.
We choose a right-handed orthonormal reference frame 
$\{\boldsymbol {e}_1,\boldsymbol {e}_2,\boldsymbol {e}_3\}$ in which 
$\boldsymbol {e}_1 = \boldsymbol {d}_1/\lVert {\boldsymbol {d}_1} \rVert$, 
$\boldsymbol {d}_1$ and 
$\boldsymbol {d}_2$ belong to the span of 
$\boldsymbol {e}_1$ and 
$\boldsymbol {e}_2$, and 
$\boldsymbol {e}_2\boldsymbol {\cdot }\boldsymbol {d}_2 \geqslant 0$ (see figure 2b).
Ignoring the 
$O(\varepsilon ^3)$ terms, applying the first-order control 
$\boldsymbol {u}_{\Delta t}^{\varepsilon, [\boldsymbol {h}_1,\boldsymbol {h}_2]}$ with small 
$\varepsilon$ causes no net displacement of the active particle and rotates each of the passive particles around the active particle about the 
$\boldsymbol {e}_3$-axis in the anticlockwise direction by small angles 
$\Delta \theta _1$ and 
$\Delta \theta _2$ given by (3.7). In the far field, 
$\Delta \theta$ decreases with 
$r$, so 
$\Delta \theta _1 < \Delta \theta _2$ so the angle 
$\theta _{12}=\arccos ({\boldsymbol {d}_1\boldsymbol {\cdot }\boldsymbol {d}_2}/({\lVert {\boldsymbol {d}_1} \rVert \times \lVert {\boldsymbol {d}_2} \rVert }))$ between 
$\boldsymbol {d}_1$ and 
$\boldsymbol {d}_2$ increases by
\begin{equation} \Delta\theta_{12}=\Delta\theta_2-\Delta\theta_1 = \frac{3a\varepsilon^2}{2}\left[\frac{1}{r_2^3}\left(1-\frac{9a}{8r_2}\right)-\frac{1}{r_1^3}\left(1-\frac{9a}{8r_1}\right)\right] \unicode{x2255} \widetilde{\Delta\theta}_{12}\varepsilon^2 \end{equation}upon each application of the first-order control.
Repeating the application of this control with a suitable choice of 
$\varepsilon$, we can rotate the particles to the final desired collinear configuration with the two passive particles on opposite sides of the active particle, corresponding to 
$\theta _{12}^f={\rm \pi}$. Specifically, if we choose to apply the control 
$N$ times, we will choose
\begin{equation} \varepsilon=\left(\frac{{\rm \pi}-\theta_{12}^\circ}{N\widetilde{\Delta\theta}_{12}}\right)^{1/2}. \end{equation}The distances from the active particle to the passive ones are unchanged by repeated applications of the first-order control while the distance between the two passive particles increases so the system remains well separated.
Proposition 3.4 (Collinear to equidistant collinear configurations)
The particles can be moved from any collinear initial configuration 
$(\boldsymbol {x}^\circ,\boldsymbol {y}_1^\circ,\boldsymbol {y}_2^\circ )$ with 
$\boldsymbol {d}_1^\circ = -\alpha \boldsymbol {d}_2^\circ$ for some 
$\alpha > 0$ to a collinear final configuration with the two passive particles equidistant from the active particle, i.e. a configuration 
$(\boldsymbol {x}^f,\boldsymbol {y}_1^f,\boldsymbol {y}_2^f)$ satisfying 
$\boldsymbol {d}_1^f = -\boldsymbol {d}_2^f$.
Proof. Suppose without loss of generality that 
$\lVert {\boldsymbol {d}_1^\circ } \rVert > \lVert {\boldsymbol {d}_2^\circ } \rVert$ and that 
$\boldsymbol {e}_1 = \boldsymbol {d}_1^\circ /\lVert {\boldsymbol {d}_1^\circ } \rVert$ (see figure 2c). Applying the zeroth-order control 
$\boldsymbol {u}_{\Delta t}^{\varepsilon, \boldsymbol {h}_1}$ moves all three particles in the 
$\boldsymbol {e}_1$ direction according to (3.2). The active particle translates by the largest magnitude, 
$\varepsilon$, and the second passive particle translates by a larger distance than the first particle because the displacement decreases with 
$r$. Hence, 
$\lVert {\boldsymbol {d}_1} \rVert$ decreases and 
$\lVert {\boldsymbol {d}_2} \rVert$ increases with each application of the zeroth-order control. Note that, within the far-field assumption, the active particle can approach the first passive particle arbitrarily closely by repeated applications of the control. Hence, by continuity, a point can be reached at which 
$\lVert {\boldsymbol {d}_1^f} \rVert = \lVert {\boldsymbol {d}_2^f} \rVert$. The configurations are always well separated if the initial configuration is.
Proposition 3.5 (Reorienting equidistant collinear configurations)
The particles can be moved from any equidistant collinear initial configuration 
$(\boldsymbol {x}^\circ,\boldsymbol {y}_1^\circ,\boldsymbol {y}_2^\circ )$ with 
$\boldsymbol {d}_1^\circ = -\boldsymbol {d}_2^\circ$ to any other equidistant collinear final configuration 
$(\boldsymbol {x}^f=\boldsymbol {x}^\circ,\boldsymbol {y}_1^f,\boldsymbol {y}_2^f)$ with the same position of the active particle and the same distance between the active and passive particles, i.e. a configuration with 
$\boldsymbol {d}_1^f = -\boldsymbol {d}_2^f$ and 
$\lVert {\boldsymbol {d}_1^f} \rVert = \lVert {\boldsymbol {d}_1^\circ } \rVert$.
Proof. We choose a right-handed orthonormal reference frame 
$\{\boldsymbol {e}_1,\boldsymbol {e}_2,\boldsymbol {e}_3\}$ in which 
$\boldsymbol {d}_1^\circ$ and 
$\boldsymbol {d}_1^f$ lie in the span of 
$\boldsymbol {e}_1$ and 
$\boldsymbol {e}_2$. Then, the transformation from the initial to the desired final configuration can be described as a rotation about the 
$\boldsymbol {e}_3$ axis through the active particle by an angle 
$\theta = \arccos (\boldsymbol {d}_1^\circ \boldsymbol {\cdot } \boldsymbol {d}_1^f/\lVert {\boldsymbol {d}_1^\circ } \rVert ^2)$.
As in the proof of Proposition 3.3, we repeatedly apply the first-order control 
$\boldsymbol {u}_{\Delta t}^{\varepsilon, [\boldsymbol {h}_1,\boldsymbol {h}_2]}$, but now 
$r_1 = r_2$ so the particles rotate by equal angles and remain collinear.
Proposition 3.6 (Translating a group of equidistant collinear particles)
Let the particles be initially equidistant and collinear, i.e. the initial configuration 
$(\boldsymbol {x}^\circ,\boldsymbol {y}_1^\circ,\boldsymbol {y}_2^\circ )$ satisfies 
$\boldsymbol {d}_1^\circ = -\boldsymbol {d}_2^\circ$. Given any scalar 
$\ell >0$ and unit vector 
$\boldsymbol {e}_{\perp } \perp \boldsymbol {d}_1^\circ$, the active and two passive particles can be translated by the same vector 
$\boldsymbol {\varDelta } = \ell \boldsymbol {e}_{\perp }$ to the final configuration 
$(\boldsymbol {x}^f,\boldsymbol {y}_1^f,\boldsymbol {y}_2^f) =(\boldsymbol {x}^\circ +\boldsymbol {\varDelta },\boldsymbol {y}_1^\circ +\boldsymbol {\varDelta },\boldsymbol {y}_2^\circ +\boldsymbol {\varDelta })$.
Proof. We consider a right-handed orthonormal reference frame 
$\{\boldsymbol {e}_1,\boldsymbol {e}_2,\boldsymbol {e}_3\}$ in which 
$\boldsymbol {e}_1 = \boldsymbol {e}_\perp$ and 
$\boldsymbol {e_2} = \boldsymbol {d}_2^\circ /\lVert {\boldsymbol {d}_2^\circ } \rVert$. In this reference frame, 
$\boldsymbol {d}_1^\circ = (0, -d^\circ, 0)^\top$ and 
$\boldsymbol {d}_2^\circ = (0, d^\circ, 0)^\top$, where 
$d^\circ = \lVert {\boldsymbol {d}_1^\circ } \rVert$.
Our strategy is to use the zeroth-order control to move the particles along the 
$\boldsymbol {e}_1$ direction. Since the passive particles move less than the active particle, they gradually fall behind. We use the first-order control to bring the passive particles in front of the active particle as needed and use symmetry to arrive with the same relative configuration as the initial state (see figure 2e).
More specifically, the first stage of our strategy involves applying the zeroth-order control 
$\boldsymbol {u}_{\Delta t}^{\ell /2, \boldsymbol {h}_1}$ to move the active particle to 
$\bar {\boldsymbol {x}} = \boldsymbol {x}^\circ + ({\ell }/{2})\boldsymbol {e}_1$. For small 
$\ell$, the displacements of the passive particles are approximated by (3.2). To leading order in 
$\ell$, the two passive particles undergo the same displacement, which is in the 
$\boldsymbol {e}_1$ direction. For finite (possibly large) 
$\ell$, the exact displacements of the two passive particles are constrained by symmetry to have the form 
$\boldsymbol {\varDelta }_{\boldsymbol {y}_1} = ({\varDelta }_1, \varDelta _2, 0)^\top$ and 
$\boldsymbol {\varDelta }_{\boldsymbol {y}_2} = (\varDelta _1, -\varDelta _2, 0)^\top$. Hence, the relative position vectors of the passive particles in this configuration are of the form 
$\bar {\boldsymbol {d}}_1=\bar {\boldsymbol {y}}_1-\bar {\boldsymbol {x}}=(-\bar {d}_{1},-\bar {d}_{2}, 0)^\top$ and 
$\bar {\boldsymbol {d}}_2=\bar {\boldsymbol {y}}_2-\bar {\boldsymbol {x}}=(-\bar {d}_{1}, \bar {d}_{2}, 0)^\top$. Since the component of the velocity in the 
$\boldsymbol {e}_1$ direction is always larger for the active particle than for the passive particles, we expect 
$\bar {d}_1 > 0$. We note, however, that this observation is not necessary for our proof.
In the second stage, we use Proposition 3.5 to rotate the passive particles by the angle 
${\rm \pi}$ about the 
$\boldsymbol {e}_2$ axis through 
$\bar {\boldsymbol {x}}$ to achieve the configuration 
$(\bar {\boldsymbol {x}}, \tilde {\boldsymbol {y}}_1,\tilde {\boldsymbol {y}}_2)$ with relative position vectors 
$\tilde {\boldsymbol {d}}_1 = (\bar {d}_1, -\bar {d}_{2}, 0)^\top$ and 
$\tilde {\boldsymbol {d}}_2 = (\bar {d}_1, \bar {d}_2, 0)^\top$.
In the third stage, we apply the zeroth-order control 
$\boldsymbol {u}_{\Delta t}^{\ell /2, \boldsymbol {h}_1}$, bringing the active particle to 
$\boldsymbol {x}^f = \bar {x} + ({\ell }/{2})\boldsymbol {e}_1 = \boldsymbol {x}^\circ + \ell \boldsymbol {e}_\perp$. This is equivalent to a time-reversal of the control applied at the beginning of our strategy in a coordinate frame that has been rotated by 
${\rm \pi}$ about the 
$\boldsymbol {e}_2$ axis through 
$\bar {\boldsymbol {x}}$. Hence, the displacements of the passive particles are the negative of the displacements 
$\boldsymbol {\varDelta }_{\boldsymbol {y}_1}$ and 
$\boldsymbol {\varDelta }_{\boldsymbol {y}_2}$ described earlier, rotated about the 
$\boldsymbol {e}_2$ direction. The final positions of the passive particles are, therefore, 
$\boldsymbol {y}_1^f = \boldsymbol {x}^f + \boldsymbol {d}_1^\circ = \boldsymbol {y}_1^\circ + \ell \boldsymbol {e}_\perp$ and 
$\boldsymbol {y}_2^f = \boldsymbol {x}^f + \boldsymbol {d}_2^\circ = \boldsymbol {y}_2^\circ + \ell \boldsymbol {e}_\perp$.
The procedure described above achieves the desired outcome provided that the configuration remains well separated at all times. Note that the second stage of the strategy does not alter distances between particles and the third stage is a rotated reversal of the first stage. Hence, we need only consider the changes in distances during the first stage of our strategy. In this stage, the distances 
$r_1=\lVert {\boldsymbol {y}_1 - \boldsymbol {x}} \rVert$ and 
$r_2=\lVert {\boldsymbol {y}_2-\boldsymbol {x}} \rVert$ increase from their initial values 
$r_1^\circ$ and 
$r_2^\circ$, respectively, because the active particle moves faster than the passive particles in a direction away from them. In contrast, the distance between the passive particles, 
$\lVert {\boldsymbol {y}_2-\boldsymbol {y}_1} \rVert =2\bar {d}_2$, decreases during this stage because the passive particles move towards each other in the 
$\boldsymbol {e}_2$-direction and always have the same position in the 
$\boldsymbol {e}_1$-direction. Hence, it is possible for the distance between the passive particles to decrease to the well-separated limit 
$R$.
Suppose that this would occur at the point where the active particle has travelled a distance 
$\ell ^*$. To avoid reaching this point, we break up the motion into a number 
$N_\ell$ of shorter segments of length 
$\bar {\ell }=\ell /N_\ell < \ell ^*$, so that we may accomplish translations of the particles by a displacement vector 
$\bar {\ell }\boldsymbol {e}_1$ without violating the well-separated condition. The repetition of this shortened motion achieves the desired final outcome and maintains the desired separation.
Proposition 3.7 (Adjusting distances between particles in an equidistant collinear configuration)
The distance between the active and passive particles in an equidistant collinear configuration can be changed arbitrarily. That is, given an initial configuration 
$(\boldsymbol {x}^\circ,\boldsymbol {y}_1^\circ = \boldsymbol {x}^\circ -r^\circ \boldsymbol {e}_2, \boldsymbol {y}_2^\circ = \boldsymbol {x}^\circ + r^\circ \boldsymbol {e}_2)$ with 
$r^\circ > R$, there is a control that achieves the final configuration 
$(\boldsymbol {x}^f=\boldsymbol {x}^\circ,\boldsymbol {y}_1^f = \boldsymbol {x}^\circ -r^f\boldsymbol {e}_2, \boldsymbol {y}_2^f = \boldsymbol {x}^\circ + r^f \boldsymbol {e}_2)$ with 
$r^f > R$.
Proof. We first describe how to achieve a final distance 
$r^f < r^\circ$. Repeated applications of the second-order control 
$\boldsymbol {u}_{\Delta t}^{\varepsilon, [\boldsymbol {h}_2,[\boldsymbol {h}_1,\boldsymbol {h}_2]]}$ leave the active particle at the initial position and move the passive particles along curves in the 
$\boldsymbol {e}_1$–
$\boldsymbol {e}_2$ plane that lead to the active particle, shown as streamlines in figure 1(c). By symmetry, the passive particles maintain equal and opposite displacements in the 
$\boldsymbol {e}_2$-direction. We stop applying this control when we reach relative positions 
$\boldsymbol {e}_2\boldsymbol {\cdot }\boldsymbol {d}_2 = -\boldsymbol {e}_2\boldsymbol {\cdot } \boldsymbol {d}_1 = r^f$.
We then apply the second-order control 
$\boldsymbol {u}_{\Delta t}^{\varepsilon, [\boldsymbol {h}_3,[\boldsymbol {h}_1,\boldsymbol {h}_3]]}$. By (3.12), replacing the index 2 with 3, this control moves the two passive particles in the 
$\boldsymbol {e}_1$-direction. We repeat this control until we achieve 
$\boldsymbol {e}_1\boldsymbol {\cdot } \boldsymbol {d}_1 = \boldsymbol {e}_1\boldsymbol {\cdot }\boldsymbol {d}_2 = 0$, which then satisfies the desired final configuration. The complete process is illustrated in figure 2(f).
By the assumption that 
$r^f > R$, we have that 
$\lVert {\boldsymbol {d}_j} \rVert \geqslant \boldsymbol {e}_2\boldsymbol {\cdot }\boldsymbol {d}_j \geqslant r^f >R$ for 
$j=1,2$ and 
$\lVert {\boldsymbol {y}_2-\boldsymbol {y}_1} \rVert = 2\lVert {\boldsymbol {d}_1} \rVert > R$ throughout this strategy, so the configuration remains well separated.
In the case where 
$r^f > r^\circ$, we apply the control strategy above in reverse.
4. Controllability for one or two passive particles
In this section, we prove controllability results for systems of one active particle and one or two passive particles, based on system (2.7). We first prove the case for two passive particles. Note that controllability with one passive particle follows from the controllability with two passive particles, since the two passive particles act as tracers and do not affect the dynamics of each other or of the active particle. The control strategy, however, can be simplified for a single passive particle, so we will present a separate proof.
Theorem 4.1 (Controllability with 
$M=2$ passive particles)
An active particle and two passive particles can be moved from any well-separated initial configuration 
$(\boldsymbol {x}^\circ, \boldsymbol {y}_1^\circ,\boldsymbol {y}_2^\circ )$ to any well-separated final configuration 
$(\boldsymbol {x}_1^f, \boldsymbol {y}_1^f, \boldsymbol {y}_2^f)$ along a well-separated trajectory. That is, given 
$(\boldsymbol {x}^\circ, \boldsymbol {y}_1^\circ,\boldsymbol {y}_2^\circ ), (\boldsymbol {x}_1^f, \boldsymbol {y}_1^f, \boldsymbol {y}_2^f) \in \mathcal {S}_R^3$, there exist 
$T\in (0,+\infty )$ and a control map 
$\boldsymbol {u}\in L^\infty ([0,T];\mathbb {R}^3)$ such that system (2.7) with 
$M=2$, namely,
\begin{equation} \begin{cases} \dot{\boldsymbol{x}}(t)=\boldsymbol{u}(t), \\[4pt] \displaystyle \dot{\boldsymbol{y}}_1(t)= \frac{3a}{4}\left(\frac{1}{\lVert {\boldsymbol{d}_1(t)} \rVert}\operatorname{\mathbb{1}}+\frac{1}{\lVert {\boldsymbol{d}_1(t)} \rVert^3}\boldsymbol{d}_1(t)\otimes \boldsymbol{d}_1(t)\right){\boldsymbol{u}}(t),\\[10pt] \displaystyle \dot{\boldsymbol{y}}_2(t)= \frac{3a}{4}\left(\frac{1}{\lVert {\boldsymbol{d}_2(t)} \rVert}\operatorname{\mathbb{1}}+\frac{1}{\lVert {\boldsymbol{d}_2(t)} \rVert^3}\boldsymbol{d}_2(t)\otimes \boldsymbol{d}_2(t)\right){\boldsymbol{u}}(t), \end{cases} \end{equation}
admits a solution 
$(\boldsymbol {x},\boldsymbol {y}_1,\boldsymbol {y}_2)\in AC([0,T];\mathcal {S}_R^{3})$ (depending on 
$\boldsymbol {u}$), such that 
$(\boldsymbol {x}(0),\boldsymbol {y}_1(0), \boldsymbol {y}_2(0))=(\boldsymbol {x}^\circ,\boldsymbol {y}_1^\circ, \boldsymbol {y}_2^\circ )$ and 
$(\boldsymbol {x}(T),\boldsymbol {y}_1(T),\boldsymbol {y}_2(T))=(\boldsymbol {x}^f,\boldsymbol {y}_1^f,\boldsymbol {y}_2^f)$.
Proof. We first describe three parts of the control strategy and then explain how they are combined.
Part 1. Using Propositions 3.3 and 3.4, we can bring the particles from the arbitrary initial configuration 
$(\boldsymbol {x}^\circ, \boldsymbol {y}_1^\circ,\boldsymbol {y}_2^\circ )$ to an intermediate configuration 
$(\tilde {\boldsymbol {x}}^\circ, \tilde {\boldsymbol {y}}_1^\circ,\tilde {\boldsymbol {y}}_2^\circ )$ that is equidistant and collinear.
Part 2. Likewise, we can bring the particles from the final configuration 
$(\boldsymbol {x}_1^f, \boldsymbol {y}_1^f, \boldsymbol {y}_2^f)$ to a configuration 
$(\tilde {\boldsymbol {x}}^f, \tilde {\boldsymbol {y}}_1^f,\tilde {\boldsymbol {y}}_2^f)$ that is equidistant and collinear.
Part 3. Using Propositions 3.5, 3.6 and 3.7, the particles can be moved from 
$(\tilde {\boldsymbol {x}}^\circ, \tilde {\boldsymbol {y}}_1^\circ,\tilde {\boldsymbol {y}}_2^\circ )$ to 
$(\tilde {\boldsymbol {x}}^f, \tilde {\boldsymbol {y}}_1^f,\tilde {\boldsymbol {y}}_2^f)$.
The sequence of steps in our complete strategy is, therefore: (i) apply Part 1; (ii) apply Part 3; (iii) apply Part 2 in reverse. Indeed, by the time-reversibility property of the Stokes equation, we can reverse the control for Part 2 to bring particles from 
$(\tilde {\boldsymbol {x}}^f, \tilde {\boldsymbol {y}}_1^f,\tilde {\boldsymbol {y}}_2^f)$ to the given final configuration. This strategy takes inspiration from Dal Maso et al. (Reference Dal Maso, DeSimone and Morandotti2015).
Each of the Parts 1–3 above can be achieved in finite time. Defining 
$T$ as the sum of these times, a control map 
$\boldsymbol {u}\in L^\infty ([0,T];\mathbb {R}^3)$ can be constructed by concatenating the controls associated with the parts above. This control map steers the system from the initial conditions to the final conditions along a well-separated trajectory. The theorem is proved.
Theorem 4.2 (Controllability with 
$M=1$ passive particle)
An active particle and a single passive particle can be moved from any well-separated initial configuration 
$(\boldsymbol {x}^\circ, \boldsymbol {y}^\circ )$ to any well-separated final configuration 
$(\boldsymbol {x}^f, \boldsymbol {y}^f)$ along a well-separated trajectory. That is, given 
$(\boldsymbol {x}^\circ, \boldsymbol {y}^\circ ), (\boldsymbol {x}^f, \boldsymbol {y}^f) \in \mathcal {S}_R^2$, there exist 
$T\in (0,+\infty )$ and a control map 
$\boldsymbol {u}\in L^\infty ([0,T];\mathbb {R}^3)$ such that system (2.7) with 
$M=1$, namely,
\begin{equation} \begin{cases} \dot{\boldsymbol{x}}(t)=\boldsymbol{u}(t), \\[3pt] \displaystyle \dot{\boldsymbol{y}}(t)= \frac{3a}{4}\left(\frac{1}{\lVert {\boldsymbol{d}(t)} \rVert}\operatorname{\mathbb{1}}+\frac{1}{\lVert {\boldsymbol{d}(t)} \rVert^3}\boldsymbol{d}(t)\otimes \boldsymbol{d}(t)\right){\boldsymbol{u}}(t), \end{cases} \end{equation}
admits a solution 
$(\boldsymbol {x},\boldsymbol {y})\in AC([0,T];\mathcal {S}_R^{2})$ (depending on 
$\boldsymbol {u}$), such that 
$(\boldsymbol {x}(0),\boldsymbol {y}(0))=(\boldsymbol {x}^\circ,\boldsymbol {y}^\circ )$ and 
$(\boldsymbol {x}(T),\boldsymbol {y}(T))=(\boldsymbol {x}^f,\boldsymbol {y}^f)$.
Proof. We provide a constructive proof of controllability, which is achieved by an appropriate composition of the zeroth-, first- and second-order controls. We apply the propositions from § 3.2, which concerned systems with two passive particles; by simply neglecting the second passive particle, those propositions describe possible moves for an active particle and a single passive particle.
Let us denote by 
$\varPi ^\circ$ the plane containing 
$\boldsymbol {x}^\circ$, 
$\boldsymbol {y}^\circ$ and 
$\boldsymbol {x}^f$, and let us notice that it is not restrictive to assume that 
$\boldsymbol {x}^\circ$ is the origin. We may also choose the reference frame such that 
$\boldsymbol {e}_3$ is perpendicular to 
$\varPi ^\circ$ and 
$\boldsymbol {x}^f=x^f\boldsymbol {e}_1$ lies on the positive 
$x$-axis (see figure 3).
Figure 3. A schematic illustration of the four steps for moving one active and one passive particle from arbitrary initial positions 
$\boldsymbol {x}^\circ$, 
$\boldsymbol {y}^\circ$ to arbitrary final positions 
$\boldsymbol {x}^f$, 
$\boldsymbol {y}^f$. Steps 1 and 2 are shown in panel (a), while Steps 3 and 4 are shown in panel (b) with a change of reference frame.
Step 1: rotation about 
$\boldsymbol {x}^\circ$. By Proposition 3.5, the passive particle can be rotated about 
$\boldsymbol {x}^\circ$ to lie on the negative 
$\boldsymbol {e}_1$ axis (see figure 3a); its position at the end of this step will be 
$\hat {\boldsymbol {y}}=-d^\circ \boldsymbol {e}_1$, where 
$d^\circ =\lVert {\boldsymbol {y}^\circ } \rVert$.
Step 2: translation of active particle. Using the direct control 
$\boldsymbol {u}(t)=\boldsymbol {u}^{x^f, \boldsymbol {h}_1}_1(t) = \boldsymbol {x}^f-\boldsymbol {x}^\circ$, we move the active particle from 
$\boldsymbol {x}^\circ$ to 
$\boldsymbol {x}^f$. The passive particle moves in the 
$\boldsymbol {e}_1$ direction to 
$\bar {\boldsymbol {y}}$. Since the velocity of the active particle is greater than that of the passive particle, the particles remain well separated during this step.
Step 3: rotation about 
$\boldsymbol {x}^f$. We consider the plane 
$\varPi ^f$ containing 
$\bar {\boldsymbol {y}}$, 
$\boldsymbol {x}^f$ and 
$\boldsymbol {y}^f$, and change the reference frame, using orthonormal vectors 
$\boldsymbol {e}_1^f = \boldsymbol {d}^f/\lVert {\boldsymbol {d}^f} \rVert$ and 
$\boldsymbol {e}_2^f$ in this plane. By Proposition 3.5, we can rotate the passive particle around the active one until it reaches a position 
$\bar {\boldsymbol {y}}^f=\bar {d}\boldsymbol {e}_1^f$ on the positive 
$x$-axis analogously to Step 1.
Step 4: translation of passive particle (distance adjustment). Note that 
$\boldsymbol {y}^f$ and 
$\bar {\boldsymbol {y}}^f$ both lie on the positive 
$\boldsymbol {e}_1$-axis relative to the active particle 
$\boldsymbol {x}^f$. By Proposition 3.7, we can adjust the distance between the active and passive particle to achieve the desired final configuration 
$(\boldsymbol {x}^f, \boldsymbol {y}^f)$.
Conclusion of the proof. Each of the Steps 1–4 above can be achieved in finite time. Defining 
$T$ as the sum of these times, a control map 
$\boldsymbol {u}\in L^\infty ([0,T];\mathbb {R}^3)$ can be constructed by concatenating the controls associated with the steps above. This control map steers the system from the initial conditions to the final conditions along a well-separated trajectory.
Remark 4.3 (An alternative proof strategy for Theorems 4.1 and 4.2)
Controllability can also be proved by appealing to the Rashewsky–Chow theorem (Agrachev & Sachkov Reference Agrachev and Sachkov2004, Theorem 5.9). For the case of 
$M=1$, it is sufficient that the six-dimensional vector fields 
$\boldsymbol {h}_1,\boldsymbol {h}_2,\boldsymbol {h}_3$ defined in § 3.1 generate a Lie algebra 
$\mathfrak {Lie}\{\boldsymbol {h}_1,\boldsymbol {h}_2,\boldsymbol {h}_3\}$ that satisfies the Lie algebra rank condition
\begin{equation} \dim\left(\mathfrak{Lie}\{\boldsymbol{h}_1(\boldsymbol{x},\boldsymbol{y}),\boldsymbol{h}_2(\boldsymbol{x},\boldsymbol{y}),\boldsymbol{h}_3(\boldsymbol{x},\boldsymbol{y})\}\right)=6\quad\text{for all } (\boldsymbol{x},\boldsymbol{y})\in\mathcal{S}_R^2. \end{equation}
Since the system is translationally and rotationally invariant, we may consider a reference frame in which the active particle is at the origin and the passive particle lies along the positive 
$x$-axis. The Lie algebra rank condition can then be verified by direct computation of the determinant of the 
$6\times 6$ matrix whose columns are (in this reference frame) 
$\boldsymbol {h}_1,\boldsymbol {h}_2,\boldsymbol {h}_3,[\boldsymbol {h}_1,\boldsymbol {h}_2],[\boldsymbol {h}_3,\boldsymbol {h}_1],[\boldsymbol {h}_2,[\boldsymbol {h}_1,\boldsymbol {h}_2]]$, evaluated at the point 
$(\boldsymbol {x},\boldsymbol {y})=(0,0,0,y_1,0,0)^\top$. To see why we expect the span of these vectors to be six-dimensional, note that 
$\boldsymbol {h}_1,\boldsymbol {h}_2,\boldsymbol {h}_3$ move the active particle precisely in the 
$\boldsymbol {e}_1,\boldsymbol {e}_2,\boldsymbol {e}_3$ directions, respectively. All first- and higher-order Lie brackets result in zero velocity of the active particle. The components of the first-order Lie bracket 
$[\boldsymbol {h}_1,\boldsymbol {h}_2]$ corresponding to the passive particle are 
$({3a}/{2y_1^2})(1-{9a}/{8y_1})(0,1,0)^\top$ for the configuration 
$\boldsymbol {x}=(0,0,0)^\top$, 
$\boldsymbol {y}=(y_1,0,0)^\top$ according to (3.6), implying that the velocity of the passive particle is in the 
$\boldsymbol {e}_2$-direction. Similarly, 
$[\boldsymbol {h}_1,\boldsymbol {h}_3]$ produces a passive particle velocity in the 
$\boldsymbol {e}_3$-direction. For the second-order Lie bracket 
$[\boldsymbol {h}_2,[\boldsymbol {h}_1,\boldsymbol {h}_2]]$, the velocity of the passive particle is in the 
$\boldsymbol {e}_1$-direction, according to (3.11). It follows that we may generate motion in any direction independently for the active and the passive particle. Similar arguments can be made for the case 
$M=2$.
5. Errors due to finite amplitudes and separations
In the control strategies for compound moves and the general controllability theorems of § 4, we used the far-field hydrodynamic flow field associated with a moving particle and we used Lie brackets to generate the necessary directions of motion of the passive particles in the asymptotic limit 
$\varepsilon \to 0$. Since the displacement per cycle decreases as 
$\varepsilon$ decreases, it may be preferable in practice to use a relatively large value of 
$\varepsilon$. In this section, we present numerical results for solutions of system (2.7), applying the first- and second-order controls to achieve a fixed target angular or linear displacement of a single passive particle with various values of 
$\varepsilon$. Numerical results for the zeroth-order control are omitted because the strategies we propose do not rely on the leading-order form of displacements due to this control. We characterize the error between the intended exact (target) displacement and the numerically computed displacement with finite 
$\varepsilon$. Numerical trajectories were obtained using the solve_ivp function with the RK45 ODE solver from the SciPy Python library. Additionally, we characterize the differences in displacements using the far-field approximation (2.7) compared with applying the same controls to system (2.6), which includes the 
$O(1/r^3)$ potential dipole in the velocity field, represented by the 
$\nabla ^2\mathcal {G}$ term in (2.4). We intentionally consider an initial separation that is only a few times the particle diameter and in the following two subsections, we show that the leading-order expressions based on far-field hydrodynamics from § 3.1 give good estimates for the particle displacements; we expect that errors would be reduced if particles are further apart.
5.1. Angular displacements
To characterize the errors associated with finite 
$\varepsilon$ when rotating a passive particle about the active particle, we consider a target angular displacement of 
$\theta ={\rm \pi} /6$ about the 
$z$-axis for a passive particle initially at position 
$\boldsymbol {y}^\circ = (5, 0, 0)^\top$ and an active particle of radius 
$a=1$ at the origin. The target position for the passive particle is, accordingly, 
$\boldsymbol {y}_{target} = (5\sqrt {3}/2, 5/2, 0)^\top$.
For a range of choices of integers 
$N$, we use (3.7) to define corresponding values of 
$\varepsilon$ such that each of 
$N$ applications of the first-order control 
$\boldsymbol {u}_{\Delta t}^{\varepsilon, [\boldsymbol {h}_1, \boldsymbol {h}_2]}$ is expected to produce a rotation by the angle 
$\Delta \theta ^\varepsilon = \theta /N$ for the relevant values 
$a=1, r=5$. We then numerically solve system (2.7) for 
$N$ applications of the control, obtaining the final position 
$\boldsymbol {y}^\varepsilon _{num}$.
Equation (3.7) neglects terms of order 
$\varepsilon ^3$ from (3.5). Hence, we can expect that the displacement generated with a finite value of 
$\varepsilon$ deviates from the desired (leading-order) motion. By symmetry, however, the velocity of the passive particle remains in the 
$x$–
$y$ plane at all times, while the active particle moves in this plane. Hence, the exact displacement (without truncation to finite order in 
$\varepsilon$) of the passive particle has zero component in the 
$z$-direction. The two components of interest are the error in the angular (polar) displacement, which in our case is
\begin{equation} \eta_\theta^\varepsilon = \frac{1}{\theta} \arccos\left(\displaystyle\frac{\boldsymbol{y}^\varepsilon_{num}}{\lVert {\boldsymbol{y}^\varepsilon_{num}} \rVert} \boldsymbol{\cdot}\frac{\boldsymbol{y}_{target}}{\lVert {\boldsymbol{y}_{target}} \rVert}\right), \end{equation}and the radial error, which we define as
\begin{equation} \eta_r^\varepsilon = \frac{\left| \lVert {\boldsymbol{y}^\varepsilon_{num}} \rVert-r^\circ \right|}{r^\circ}. \end{equation} The convergence of these two components of error is shown in figure 4 (labelled ‘one-cycle polar’ and ‘one-cycle radial’ in the legend). Both components of error appear to converge to zero with first-order rate of convergence with 
$\varepsilon$, which we could anticipate since (3.7) neglects terms of order 
$\varepsilon ^3$ and the number of required iterations 
$N$ scales as 
$1/\varepsilon ^2$.
Figure 4. Convergence of the passive particle displacement with 
$N$ repeated applications of the first-order control corresponding to 
$[\boldsymbol {h}_1,\boldsymbol {h}_2]$ with control amplitude 
$\varepsilon$ for a fixed target rotation by the angle 
${\rm \pi} /6$ about the origin.
We can, however, improve the rate of convergence by modifying the first-order control. The general first-order control (3.3) applied with any of the following 
$(k,l)$ combinations of indices – 
$(1,2), (-1, -2), (2,-1), (-2,1)$ – all result in the same leading-order term in the displacement (3.5), but the higher-order terms differ. Here, a negative index indicates that we apply the control in the negative direction. As shown in figure 4, we find that a strategy of alternating between the first two choices, which we refer to as a two-cycle control, produces errors that decay quadratically with 
$\varepsilon$. A four-cycle strategy, in which we cycle through all four of the listed pairs of 
$(k,l)$, results in radial errors that decay cubically and polar errors that decay quadratically with 
$\varepsilon$. Note that values of 
$\varepsilon$ as large as 1 can be used for rotations with errors of 
$10^{-3}$ or less. The paths of the active and passive particles over one application of the four-cycle control are illustrated in figure 5.
Figure 5. A multi-cycle control corresponding to the first-order Lie bracket 
$[\boldsymbol {h}_1,\boldsymbol {h}_2]$. Initially, the active particle is at 
$\boldsymbol {x}^\circ = (0,0,0)^\top$ and the passive particle is at 
$\boldsymbol {y}^\circ = (5,0,0)^\top$. The panel on the right is a magnification of the path of the passive particle. The two-cycle consists of portions 1 and 2, whereas the four-cycle consists of portions 1 to 4, in sequence. The amplitude of the control is 
$\varepsilon =2$.
When we include the potential dipole term in the flow field, we find that the radial error is indistinguishable from the case without the potential dipole. In contrast, the polar error does not converge to zero with the four-cycle strategy but it remains below 
$10^{-2}$ for 
$\varepsilon < 1$.
We assert that it is more important to reduce the radial error than it is to guarantee a small polar error because we can easily adjust 
$\varepsilon$ or the number of iterations to compensate for errors in the angular component, whereas we require a different type of control to correct for changes in radial distance. In particular, the second-order control could generate a corrective radial displacement for one passive particle, but this may not be able to correct the trajectories of two passive particles simultaneously.
Note that we considered particles that are initially close together, with 
$\lVert {\boldsymbol {y}^\circ - \boldsymbol {x}^\circ } \rVert = 5a$. This was chosen to illustrate that even when the far-field regime is not strictly observed, the first-order control is an effective strategy for achieving circular motion of a passive particle around an active one. Radial deviations are small and angular displacements are well approximated by the leading-order expression given by (3.7). The analytical formula could be modified to include the effect of the potential dipole if greater accuracy were required.
5.2. Linear displacements
Errors for the second-order control are analysed similarly, using a target displacement from 
$\boldsymbol {y}^\circ = (5,0,0)^\top$ to 
$\boldsymbol {y}_{target}=(5.1,0,0)^\top$. For a given number 
$N$ of applications of the second-order control, we numerically determine the value of 
$\varepsilon$ that would result in the target displacement according to (3.11), noting that 
$r$ in this equation changes with each application of the control. The polar component of error is defined as
\begin{equation} \eta_\theta^\varepsilon = \arccos\left(\frac{y^\varepsilon_{{num}, 1}}{\lVert {\boldsymbol{y}^\varepsilon_{num}} \rVert}\right), \end{equation}and the radial component of error for this motion is
\begin{equation} \eta_r^\varepsilon = \frac{\left| y^\varepsilon_{{num},1}-y_{{target},1} \right|}{\left|y_{{target},1}-y^\circ_1\right|} = 10\left| y^\varepsilon_{{num},1}-y_{{target},1} \right|. \end{equation} Following the description in § 3.1.3, the one-cycle second-order control is obtained by applying the general pattern (3.3) with 
$\boldsymbol {h}_k$ corresponding to the first-order control 
$[\boldsymbol {h}_1,\boldsymbol {h}_2]$ and 
$\boldsymbol {h}_l = -\boldsymbol {h}_2$. The two-cycle control alternates between this and the control with 
$\boldsymbol {h}_k$ corresponding to the first-order control 
$[\boldsymbol {h}_1,-\boldsymbol {h}_2]$ and 
$\boldsymbol {h}_l = \boldsymbol {h}_2$. The motion of the two particles generated by the two-cycle control is illustrated in figure 6.
Figure 6. A multi-cycle control corresponding to the second-order Lie bracket 
$[[\boldsymbol {h}_1,\boldsymbol {h}_2],-\boldsymbol {h}_2]=[\boldsymbol {h}_2,[\boldsymbol {h}_1,\boldsymbol {h}_2]]$. Initially, the active particle is at 
$\boldsymbol {x}^\circ = (0,0,0)^\top$ and the passive particle is at 
$\boldsymbol {y}^\circ = (5,0,0)^\top$. The panel on the right is a magnification of the path of the passive particle. The direction of travel along the curves is indicated by arrows and the portions of the cycles correspond to applying the controls: (1a) 
$\boldsymbol {u}^{\varepsilon, [\boldsymbol {h}_1,\boldsymbol {h}_2]}$, (1b) 
$\boldsymbol {u}^{\varepsilon, -\boldsymbol {h}_2}$, (2a) 
$\boldsymbol {u}^{\varepsilon, [\boldsymbol {h}_1,-\boldsymbol {h}_2]}$ and (2b) 
$\boldsymbol {u}^{\varepsilon, \boldsymbol {h}_2}$ all with 
$\varepsilon =2$.
As shown in figure 7, the errors in the radial (linear) direction are similar for the one- and two-cycle controls over the range of 
$\varepsilon$ considered, decaying approximately linearly with 
$\varepsilon$. Polar errors decay quadratically with 
$\varepsilon$ with the one-cycle strategy and quintically for the two-cycle strategy. Both polar and radial errors are essentially unchanged when the potential dipole terms are included, as shown in figure 7.
Figure 7. Convergence of the passive particle displacement with 
$N$ repeated applications of the second-order control corresponding to 
$[\boldsymbol {h}_2,[\boldsymbol {h}_1,\boldsymbol {h}_2]]$ with control amplitude 
$\varepsilon$ for a fixed target displacement from 
$\boldsymbol {y}^\circ =(5,0,0)^\top$ to 
$\boldsymbol {y}_{target} = (5.1,0,0)^\top$.
Since polar errors decay rapidly as 
$\varepsilon$ decreases, we can readily achieve linear motion of a passive particle using the second-order control. The errors in the radial direction, which are larger in magnitude than those in the polar direction, can be corrected either by considering higher-order terms in the analytic expression for the displacement or by reducing 
$\varepsilon$, perhaps incrementally as the target position is approached.
6. Conclusions
In this paper, we have presented the motion planning problem for a system of one active and 
$M$ passive spherical particles immersed in a viscous fluid, in the far-field approximation. We explicitly constructed elementary moves that, suitably concatenated, resulted in strategies to achieve total controllability in the specific cases 
$M=1$ and 
$M=2$. Moreover, the strategies we proposed ensure that the particles can maintain an arbitrary minimum separation compatible with their initial and final configurations. The elementary and compound moves were expressed in terms of zeroth-, first- and second-order controls characterized by an amplitude parameter 
$\varepsilon$, with asymptotic expressions valid in the limit 
$\varepsilon \to 0$. We showed that in this limit, the passive particle displacements resulting from the zeroth-, first- and second-order controls correspond to the Stokeslet, rotlet and rotlet dipole singularity solutions of Stokes flow, respectively. Higher-order singularities can be generated by extension of the controls to higher orders.
Through numerical solutions of the equations of motion, we demonstrated that the two key components of our motion planning strategy, namely, moving passive particles in a circular orbit around an active one and translating a passive particle without a net displacement of the active particle, could be achieved to a high accuracy even with 
$\varepsilon \approx 1$ and with the particles as close as a few diameters apart.
This research contributes to the growing literature on ensembles of microparticles subject to hydrodynamic interactions in low-Reynolds-number flows, including the possibly chaotic behaviour of sedimenting particles (Hocking Reference Hocking1964; Jánosi et al. Reference Jánosi, Tél, Wolf and Gallas1997), mixing and transport in suspensions of microswimmers (Katija & Dabiri Reference Katija and Dabiri2009; Lauga & Powers Reference Lauga and Powers2009; Pushkin, Shum & Yeomans Reference Pushkin, Shum and Yeomans2013), idealized models of swimmers such as Purcell's scallop or three-link swimmers (Purcell Reference Purcell1977) and three linked spheres (Najafi & Golestanian Reference Najafi and Golestanian2004). Mathematical treatments of control of model swimmers started with the seminal paper by Shapere & Wilczek (Reference Shapere and Wilczek1989) and have since been applied in many contexts (Alouges, DeSimone & Lefebvre Reference Alouges, DeSimone and Lefebvre2008; Chambrion & Munnier Reference Chambrion and Munnier2012; Dal Maso et al. Reference Dal Maso, DeSimone and Morandotti2015; Chambrion, Giraldi & Munnier Reference Chambrion, Giraldi and Munnier2019; Lohéac & Takahashi Reference Lohéac and Takahashi2020; Zoppello, Morandotti & Bloomfield-Gadêlha Reference Zoppello, Morandotti and Bloomfield-Gadêlha2022; Attanasi, Zoppello & Napoli Reference Attanasi, Zoppello and Napoli2024).
The present contribution sets the basis for further investigations from multiple viewpoints. Six areas of future research that could be of interest to the mathematical, physical and engineering communities are as follows.
(1) Using periodic controls for the active particles to produce flow fields that act as hydrodynamic traps (Lutz, Chen & Schwartz Reference Lutz, Chen and Schwartz2006; Chamolly, Lauga & Tottori Reference Chamolly, Lauga and Tottori2020). Rather than moving a passive particle from one specific location to another, we may want to attract all nearby particles to a target and hold them there, possibly against other effects such as a background flow or gravity.
(2) Considering cases in which we have
$N$ active particles and $M$ passive ones, with both $N$ and $M$ large. Generalizing the formulation (2.8) to arbitrary numbers of active and passive particles is relatively straightforward but the task of effectively controlling $M$ passive particles with a minimal number of active particles is challenging. It could also be of interest to investigate whether partial controllability results can be proved for an even lower number of controllers. We stress that, even in the case $N=1$ and $M=3$, the strategies proposed in our proofs would have to be substantially modified, since the presence of a third passive particle disrupts the symmetry that has been exploited in some of the moves. For example, the strategy used in Proposition 3.6 (translating a group of equidistant collinear particles) does not work even if the particles are not all required to be collinear, because we cannot guarantee that the symmetry is preserved for the third passive particle.(3) Manipulating particles in bounded domains. As proved by, e.g. Alouges & Giraldi (Reference Alouges and Giraldi2013), Buzhardt et al. (Reference Buzhardt, Fedonyuk, Sudarsanam and Tallapragada2018a), the presence of boundaries can enhance controllability. Confining walls are especially relevant to microfluidic systems so studying such domains is a natural extension.
(4) Accounting for near-field interactions. One limitation of our setting is the far-field approximation that we use to simplify the equations. In real-world systems, especially spatially confined ones with many particles, it may be necessary to accurately account for the mutual influence of particles that are close to one another, in which case the mobility tensor of each particle depends on the positions of all (nearby) particles in the system.
(5) Mixing fluids at low Reynolds number. This is known to be challenging in microfluidic devices (Ward & Fan Reference Ward and Fan2015); one proposed technique involves using magnetic particles in rotating magnetic fields (Munaz et al. Reference Munaz, Kamble, Shiddiky and Nguyen2017), which corresponds to our model of actively actuated particles but with applied torques and rotations of the active particles playing significant roles.
(6) Accounting for stochastic terms (i.e. Brownian motion) in the dynamics of the passive particles (Graham Reference Graham2018). In the current work, we assumed that particles were large enough that Brownian motion could be neglected, but this may not be valid if the particles are small (more precisely, when the Péclet number is small). Including random diffusion effects would be particularly interesting when the number of passive particles is very large (ideally, diverging to infinity), to the point that a description in terms of the distribution of particles would be preferable. In this context, it is customary to study the partial differential equation (PDE) arising for the limiting distribution, which, in this context, is expected to be a Fokker–Planck-type equation featuring a transport term coming from the action of the active particles, with the diffusion term resulting from the Brownian motion. While this is a very promising and fertile research field, it is beyond the scope of the present paper.
Acknowledgements
The authors thank AnhadSingh Bagga for assistance with some of the figures. M.Z. is a member of the Gruppo Nazionale di Fisica Matematica of the Istituto Nazionale di Alta Matematica. M.M. is a member of the Gruppo Nazionale per l'Analisi Matematica, la Probabilità e le loro Applicazioni of the Istituto Nazionale di Alta Matematica. This study was carried out within the GNAMPA2024 project Analisi asintotica di modelli evolutivi di interazione (CUP E53C23001670001).
Funding
H.S. and M.A. acknowledge the support of the Natural Sciences and Engineering Research Council of Canada (NSERC) [funding reference number RGPIN-2018-04418].
Cette recherche a été financée par le Conseil de recherches en sciences naturelles et en génie du Canada (CRSNG) [numéro de référence RGPIN-2018-04418].
Funding from the Mathematics for Industry 4.0 2020F3NCPX PRIN2020 (M.M. and M.Z.), funded by the Italian MUR, and the Geometric-Analytic Methods for PDEs and Applications 2022SLTHCE (M.M.) is gratefully acknowledged. This manuscript reflects only the authors’ views and opinions and the Ministry cannot be considered responsible for them.
Declaration of interests
The authors report no conflict of interest.
Open access
