1. Introduction
Some notion of control is widely sought in many contexts. In fluid dynamics, particularly on the microscale, experimental techniques such as optical tweezers, flow modulation and applied magnetic fields are capable of precisely influencing the microfluidic environment, with the potential to realise refined control. The subjects of such control are varied, spanning both synthetic and biological agents, including the canonical spermatozoon (Zaferani, Cheong & Abbaspourrad Reference Zaferani, Cheong and Abbaspourrad2018). Correspondingly, the intended outcomes of control are broad, from the delivery of therapeutics in medicine (Kei Cheang et al. Reference Kei Cheang, Lee, Julius and Kim2014; Yasa et al. Reference Yasa, Erkoc, Alapan and Sitti2018; Tsang et al. Reference Tsang, Demir, Ding and Pak2020) to the trapping and sorting of cells in microdevices (Walker et al. Reference Walker, Ishimoto, Wheeler and Gaffney2018; Zaferani et al. Reference Zaferani, Cheong and Abbaspourrad2018).
Recently, many studies have sought to control microswimmers. Many examples concern the directed motion of magnetotactic helical microswimmers (Mahoney et al. Reference Mahoney, Sarrazin, Bamberg and Abbott2011; Tottori et al. Reference Tottori, Zhang, Qiu, Krawczyk, Franco-Obregón and Nelson2012; Liu et al. Reference Liu, Xu, Guan, Yan, Ye and Wu2017; Wang et al. Reference Wang, Hu, Schurz, De Marco, Chen, Pané and Nelson2018) part of the broader class of magnetic micromachines (Tierno et al. Reference Tierno, Golestanian, Pagonabarraga and Sagués2008; Grosjean et al. Reference Grosjean, Lagubeau, Darras, Hubert, Lumay and Vandewalle2015; Khalil et al. Reference Khalil, Klingner, Hamed, Hassan and Misra2020). Abstracted away from the modality of control, a noteworthy and increasingly popular approach to the design of control schemes is the application of machine learning techniques (Colabrese et al. Reference Colabrese, Gustavsson, Celani and Biferale2017; Schneider & Stark Reference Schneider and Stark2019), with a particularly remarkable example being the reinforcement-learning scheme of Mirzakhanloo, Esmaeilzadeh & Alam (Reference Mirzakhanloo, Esmaeilzadeh and Alam2020) that is used to enact theoretical hydrodynamic cloaking. Many of these modern methodologies are evaluated via a number of simulated or experimental examples, with efficacy justified by successful realisation of control in these cases. Although such a trial-based validation is reassuring, there is often an absence of a complimentary theoretical basis that provides rigorous assurances that control is indeed possible from a given configuration, a desirable if not necessary property when seeking to elicit control and guidance in practice. The topic of such an uncommon analysis is that of controllability, a broad field that, in this context, seeks to determine the conditions in which a given system can be controlled, querying the theoretical existence of a trajectory in the state space that connects given initial and target configurations. This topic is rich, varied and often framed abstractly; thus, as a practical introduction, we aim to concretely summarise elementary but powerful aspects of control theory, with a view to application in the Stokes regime.
In the context of microswimming, there has been extensive exploration of the specific cases of general deformable or propulsive bodies (Martín, Takahashi & Tucsnak Reference Martín, Takahashi and Tucsnak2007; Lohéac, Scheid & Tucsnak Reference Lohéac, Scheid and Tucsnak2013; Lohéac & Munnier Reference Lohéac and Munnier2014; Dal Maso, DeSimone & Morandotti Reference Dal Maso, DeSimone and Morandotti2015; Lohéac & Takahashi Reference Lohéac and Takahashi2020) and model microrobots, formed of connected spheres (Desimone et al. Reference Desimone, Heltai, Alouges and Lefebvre-Lepot2012; Alouges Reference Alouges2013; Alouges et al. Reference Alouges, DeSimone, Heltai, Lefebvre-Lepot and Merlet2013b; Gérard-Varet & Giraldi Reference Gérard-Varet and Giraldi2015; DeSimone Reference DeSimone2020) or links (Alouges et al. Reference Alouges, DeSimone, Giraldi and Zoppello2013a). These studies utilise the gauge field formulation of force- and torque-free microswimming first posed by Shapere & Wilczek (Reference Shapere and Wilczek1987, Reference Shapere and Wilczek1989b), in which the swimming velocity of a deformable body is described in a form linearly related to the surface deformation velocity. As such, the control functions considered in these works are invariably given by the shape. This geometrical phase theory has been used extensively for computing efficient locomotion (Shapere & Wilczek Reference Shapere and Wilczek1989a; Ramasamy & Hatton Reference Ramasamy and Hatton2017; Bittner, Hatton & Revzen Reference Bittner, Hatton and Revzen2018; Ramasamy & Hatton Reference Ramasamy and Hatton2019), and was recently extended to a more-general swimmer including the case with external forces (Koens & Lauga Reference Koens and Lauga2021). Recent studies have also addressed the controllability of magneto-elastic microrobots (Giraldi & Pomet Reference Giraldi and Pomet2017; Moreau Reference Moreau2019), with the control function being an external magnetic field, requiring characteristically different but nevertheless case-specific analysis than the standard geometrical phase theory.
These studies have each assessed the controllability of a single object immersed in fluid, though the applicability of their conclusions to multi-swimmer systems remains unclear. Indeed, the long-range hydrodynamic interactions present in the Stokes limit can lead to surprising and complex behaviours, such as flagellar synchronisation (Brumley et al. Reference Brumley, Wan, Polin and Goldstein2014; Bruot & Cicuta Reference Bruot and Cicuta2016) and the bound swimming states of algae (Drescher et al. Reference Drescher, Leptos, Tuval, Ishikawa, Pedley and Goldstein2009). With this richness of emergent behaviours, the consideration of multi-particle controllability thereby warrants the treatment of non-local interactions, of particular relevance given recent developments in experimental control via optical tweezers (Zou et al. Reference Zou, Zheng, Wu and Lei2020). One such exploration, though notably in the case of inviscid flows, considered the control of a passive particle in a two-dimensional fluid, with control effected by the movement of a cylinder and the induced flow (Or et al. Reference Or, Vankerschaver, Kelly, Murray and Marsden2009), though similar such evaluations in Stokes flow are currently lacking, to the best of the authors’ knowledge.
Thus, as a further significant aim of this work, we investigate the theoretical control and controllability of particles via hydrodynamic interactions, considering both accurate and simplified descriptions of canonical Stokes problems, with such justification currently lacking for even the most simple multi-object settings in the Stokes limit. Further, we seek to highlight additional benefits of conducting such an analysis, with particular regard to the notion of mechanical efficiency.
Hence, and as the principal aim of this work, we seek to demonstrate the broad utility of controllability analysis in the context of particle motion and hydrodynamic interactions in the Stokes limit. In doing so, and as an additional aim of this study, we evaluate and explore the controllability of experimentally motivated and seemingly simple systems of particles in a Stokesian fluid, which, despite their idealised nature, give rise to non-trivial control and hydrodynamic problems. To begin, in § 2, we recount key principles and definitions of control theory, cast in a readily applicable form that may be easily translated to other Stokes problems. Equipped with such a framework, we then investigate the controllability of multi-sphere systems in detail, focusing in particular on two-particle problems. In the first instance, motivated by exemplar modalities of actuation on the microscale, in § 3, we evaluate the degree to which two differently sized spheres can be controlled by the application of forces or torques to one of the particles, resulting in a sufficient controllability condition for a pair of particles in three dimensions in the Stokes regime and an integrable system of torque control. Moving to assess the efficiency of aspects of force-driven control, we make use of both high-accuracy and far-field hydrodynamic approximations to the interactions of these spheres, noting differences in predicted controllability and future extensions to many-sphere systems. Next, in §§ 4 and 5, we consider the same problem setting in the limit where the passive spheres are of vanishing size, moving to establish the controllability of passive tracer particles that are advected by the flow due to a force-driven finite-sized sphere, later considering a further simplification of this system. Finally, in § 6, by highlighting and exemplifying the constructive process of Lafferriere & Sussmann (Reference Lafferriere and Sussmann1992), we demonstrate how these analyses can be used to explicitly formulate controls that effect a desired state change, in this case the targeted motion of spheres along prescribed trajectories by both applied forces and the resulting flows.
2. Controllability
In this study, we consider the motion of 
$N$ spheres of potentially unequal radius in a Newtonian fluid in the regime of vanishing Reynolds number, with the fluid velocity 
$\boldsymbol {u}$ and the pressure 
$p$ being described by the Stokes equations,
\begin{equation} \mu\nabla^2\boldsymbol{u}=\boldsymbol{\nabla} p. \end{equation}
The flow also satisfies the incompressibility condition 
$\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u}=0$ and the fluid viscosity 
$\mu$ is assumed to be constant. We assume that all quantities have been appropriately non-dimensionalised, considering dimensionless mobility coefficients throughout. Let 
$\boldsymbol {x}_i$ (
$i=1,2,\ldots,N$) be the position of the 
$i$th sphere, with 
$\boldsymbol {U}_i$ and 
$\boldsymbol {\varOmega }_i$ being its translational and rotational velocities, respectively. Here and in § 3, we focus on the effects of external forces and torques, denoted by 
$\boldsymbol {f}_i$ and 
$\boldsymbol {m}_i$, respectively, on the 
$N$-sphere system, in contrast to the numerous aforementioned studies of shape-deforming swimmers and motivated by the manipulation of microparticles by optical tweezers and the actuated spinning of magnetic rotors. By the linearity of the Stokes equations, the dynamics of the spheres may be written as (Kim & Karrila Reference Kim and Karrila2005)
\begin{equation} [\begin{matrix}
\boldsymbol{U}_1;\cdots;\boldsymbol{U}_N;
\boldsymbol{\varOmega}_1;\cdots;\boldsymbol{\varOmega}_N
\end{matrix}]=\boldsymbol{\mathsf{A}} [\begin{matrix}\,\boldsymbol{f}_1;\cdots;\boldsymbol{f}_N;
\boldsymbol{m}_1;\cdots;\boldsymbol{m}_N \end{matrix}].
\end{equation}
Here, a semicolon denotes vertical concatenation, with the resulting composite force and velocity vectors being of scalar dimension 
$6N$. The 
$6N\times 6N$ tensor 
$\boldsymbol{\mathsf{A}}$, often termed the grand mobility tensor, is symmetric and positive-definite, noting that any stresslet contributions vanish in the absence of background flow. Further, due to the boundary-value nature of Stokes flow, this mobility tensor is completely determined by the shape of the fluid boundary. Thus, for the systems of spheres considered in this work, this tensor depends only on the positions and relative sizes of the spheres, with individual symmetry naturally removing any rotational dependence. For convenience, let 
$n=3N$ and 
$m=6N$, which represent the composite dimension of the vectors of sphere positions and velocities, respectively. We introduce the 
$n$-dimensional position vector
\begin{equation} \boldsymbol{X}=[\boldsymbol{x}_{1};\boldsymbol{x}_{2};\cdots;\boldsymbol{x}_{N}] \end{equation}
and the 
$m$-dimensional external control vector
\begin{equation} \boldsymbol{F}=[\boldsymbol{f}_{1};\boldsymbol{f}_{2};\cdots; \boldsymbol{f}_{N};\boldsymbol{m}_{1};\boldsymbol{m}_{2};\cdots;\boldsymbol{m}_{N}], \end{equation}so that, by linearity, we may write the translational evolution of the sphere system in the form
\begin{equation} \frac{\mathrm{d}\boldsymbol{X}}{\mathrm{d}t} = \boldsymbol{\mathsf{G}}\boldsymbol{F}. \end{equation}
identifying 
$\boldsymbol{\mathsf{G}}$ with the upper half of 
$\boldsymbol{\mathsf{A}}$ from (2.2) and with sphere orientation now implicit. Representing the 
$n\times m$ tensor 
$\boldsymbol{\mathsf{G}}=\boldsymbol{\mathsf{G}}(\boldsymbol {X})$ by its column vectors 
$\boldsymbol {g}_1,\ldots,\boldsymbol {g}_m$, with 
$\boldsymbol{\mathsf{G}}=[\boldsymbol {g}_1,\boldsymbol {g}_2,\ldots \boldsymbol {g}_{m}]$, leads to the system
\begin{equation} \frac{\mathrm{d}\boldsymbol{X}}{\mathrm{d}t} = \sum_{i=1}^m F_i(t)\boldsymbol{g}_i (\boldsymbol{X}), \end{equation}
which is the standard form of a driftless control-affine system. Here, the 
$n$-dimensional vector 
$\boldsymbol {X}$ is the state of the system, evolving in the state space 
$P \subset \mathbb {R}^n$, and the 
$F_i(t)$ are our control functions, which correspond here to components of applied forces and torques. Of note, the hydrodynamic interactions between the spheres are included in the 
$m$ column vectors 
$\boldsymbol {g}_i$, whose time dependence is only via the state 
$\boldsymbol {X}$ due to the Stokes limit, a property that will generally hold for problems of Stokes flow. Following Coron (Reference Coron2007), to which we refer the interested reader for a detailed summary of the theory of control, we define the term controllable (in any time) for the sphere system of (2.6).
Definition 2.1 The system is controllable within a state set 
$Q \subset P$ if, for any given start state 
$\boldsymbol {X}_0 \in Q$, end state 
$\boldsymbol {X}_1 \in Q$ and duration of time 
$T>0$, there exists a control function 
$\boldsymbol {F}(t)$ that transports the set of spheres from 
$\boldsymbol {X}_0$ at 
$t=0$ to 
$\boldsymbol {X}_1$ at 
$t=T$.
As the field 
$\boldsymbol {g}_i$ associated with the control 
$F_i(t)$ is tangent to the trajectory of the dynamical system, it might seem that the number of scalar controls needs to be as large as the dimension of the system in order to yield a controllable system. However, as two fields do not commute in general, a sequence of controls may generate additional reachable directions and, thus, give rise to controllability in a system of relatively high dimension. This non-commutativity is mathematically represented by a Lie bracket, defined for two fields 
$\boldsymbol {g}_i$ and 
$\boldsymbol {g}_j$ as
\begin{equation} [{\boldsymbol{g}_i},{\boldsymbol{g}_j}](\boldsymbol{X})=(\boldsymbol{\nabla}_{\boldsymbol{X}}\boldsymbol{g}_j)\boldsymbol{g}_i-(\boldsymbol{\nabla}_{\boldsymbol{X}}\boldsymbol{g}_i)\boldsymbol{g}_j, \end{equation}
where 
$\boldsymbol {\nabla }_{\boldsymbol {X}}\boldsymbol {g}_j=[{\partial \boldsymbol {g}_j}/{\partial X _1} \cdots {\partial \boldsymbol {g}_j}/{\partial X _n}]$ denotes the 
$n\times n$ Jacobian matrix. The Lie bracket satisfies 
$[{\boldsymbol {g}_i},{\boldsymbol {g}_j}]=-[{\boldsymbol {g}_j},{\boldsymbol {g}_i}]$, so that 
$[{\boldsymbol {g}_i},{\boldsymbol {g}_i}]=\boldsymbol {0}$, as well as the Jacobi identity.
As an explicit example, for 
$\tau \ll 1$ consider the following sequence of controls on the short time interval 
$t\in [0, 4\tau ]$:
\begin{equation} \boldsymbol{F}(t)=\left\{\begin{array}{@{}ll} F_i= + 1, & \text{for } t\in[0,\tau), \\ F_j= + 1, & \text{for } t\in[\tau,2\tau), \\ F_i={-}1, & \text{for } t\in[2\tau,3\tau), \\ F_j={-}1, & \text{for } t\in[3\tau,4\tau], \end{array}\right. \end{equation}
where all unspecified components of 
$\boldsymbol {F}$ are assumed to be zero. Under this control, the state 
$\boldsymbol {X}(0)$ at 
$t=0$ evolves to
\begin{equation} \boldsymbol{X}(4\tau)=\boldsymbol{X}(0)+\tau^2[\boldsymbol{g}_i, \boldsymbol{g}_j](\boldsymbol{X}(0))+O(\tau^3), \end{equation}
as illustrated in figure 1(a), with a potential lack of commutativity leading to a non-zero final displacement. This reasoning extends to higher-order Lie brackets, with these terms together spanning all the reachable directions around the state 
$\boldsymbol {X}$ via 
$\boldsymbol {F}$. This set of all possible Lie brackets, which includes the original fields 
$\boldsymbol {g}_i$, is denoted by 
$\mathrm {Lie}(\boldsymbol {g}_1,\ldots,\boldsymbol {g}_m)$, and may be written explicitly as
\begin{gather} \mathrm{Lie}(\boldsymbol{g}_1,\ldots,\boldsymbol{g}_m)=\{ \boldsymbol{g}_1, \dots,\boldsymbol{g}_m, [\boldsymbol{g}_1, \boldsymbol{g}_2], \ldots,[\boldsymbol{g}_1, [\boldsymbol{g}_1, \boldsymbol{g}_2]],\ldots,[\boldsymbol{g}_1, [\boldsymbol{g}_1, [\boldsymbol{g}_1, \boldsymbol{g}_2]]]\ldots\}. \end{gather}
In what follows, we refer to the original fields 
$\boldsymbol {g}_i$ as the first-order brackets, with 
$[{\boldsymbol {g}_i},{\boldsymbol {g}_j}]$ being termed a second-order bracket and this convention naturally extending to higher-order terms. The elements of 
$\mathrm {Lie}(\boldsymbol {g}_1,\ldots,\boldsymbol {g}_m)$, each evaluated at a state 
$\boldsymbol {X}\in P$, naturally span a vector space 
$B(\boldsymbol {X})$ when seen as elements of 
$\mathbb {R}^n$. We can then define the controllable state set 
$S\subseteq P$, as illustrated in figure 1(c), to be the set of states from which one can move in any direction in the state space or, perhaps more simply,
\begin{equation} S\colon=\left\{ \boldsymbol{X}\in P: \dim{B(\boldsymbol{X})}=n\right\}. \end{equation}
As intuition suggests, 
$S$ is inherently related to the controllability of the system (Coron Reference Coron2007):
Theorem 2.2 (Rashevski–Chow theorem)
The driftless control-affine system of (2.6) is controllable within any connected subset of the controllable state set 
$S$.
Figure 1. Schematic of (a) the Lie bracket, (b) the control system in 
$\mathbb {R}^n$ and (c) the controllable state set. (a) The Lie bracket, 
$[{\boldsymbol {g}_i},{\boldsymbol {g}_j}]$, can be computed by considering sequences of controls with infinitesimal duration of time. (b) The controllability of the system guarantees that one can effect evolution from an initial state 
$\boldsymbol {X}_0$ to a final state 
$\boldsymbol {X}_1$ at a given time 
$T$. In a driftless control-affine system, the field 
$\boldsymbol {g}_i(\boldsymbol {X})$ associated with a control 
$F_i(t)$ is a tangent of the trajectory but, if the system is controllable, one can nevertheless control the state around a prescribed path (dashed line). The realised path in the state space is shown as a solid curve, approximately coincident with the prescribed path. (c) The controllable state set 
$S$ is a subspace of 
$P (\subset \mathbb {R}^n)$, though the system is only controllable within connected subsets of 
$S$.
Therefore, in what follows, we seek to establish the controllability of our 
$N$-sphere system by determining the dimension of the vector space 
$B(\boldsymbol {X})$ associated with each state 
$\boldsymbol {X}\in P$, thereby constructing the controllable state set 
$S$.
Further, a driftless control-affine system that satisfies this full-rank condition also enjoys a property known as small-time local controllability (STLC) at every 
$\boldsymbol {X} \in S$, which, in the context of this work, guarantees that a prescribed trajectory in the controllable state space can be followed with arbitrary accuracy, as illustrated in figure 1(b). We revisit this property briefly in § 6 and refer the interested reader to the work of Coron (Reference Coron2007) for a thorough definition of STLC.
3. Finite size spheres
3.1. Geometry and hydrodynamics
As a particular case of the 
$N$-sphere system introduced above, we consider the motion of two spheres with the further restriction that only a single sphere may be externally driven, as illustrated in figure 2. With subscripts of 1 and 2 corresponding to the driven control sphere and the other passive sphere, respectively, these assumptions amount to taking 
$\boldsymbol {f}_2=\boldsymbol {m}_2=\boldsymbol {0}$ and we write 
$\boldsymbol {f}_1=\boldsymbol {f}$ and 
$\boldsymbol {m}_1=\boldsymbol {m}$ for convenience. The fixed dimensionless radii of the two spheres are denoted by 
$a_1$ and 
$a_2$, respectively, from which we define their ratio 
$\lambda = a_2/a_1$ and henceforth assume that the radius of the control sphere is unity, without loss of generality. In addition, we proceed by assuming that the spheres do not overlap or make contact with one another, so that the state space is given by 
$P = \{\boldsymbol {X}\in \mathbb {R}^6: \left \lVert \boldsymbol {x}_2-\boldsymbol {x}_1\right \rVert > 1 + \lambda \}$.
Figure 2. Schematic of the control of two spheres by moving one sphere. The force-control problem queries the existence of a forcing function 
$\boldsymbol {f}(t)$ that transports the two spheres between given initial and end positions in a given time 
$T$. Moreover, the full-rank condition in the sphere system allows us to control the spheres to approximately follow prescribed trajectories (dashed lines) in the controllable space.
With these definitions, the general equation of motion (2.5) may be written as
\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t}\begin{bmatrix} \boldsymbol{x}_1 \\ \boldsymbol{x}_2 \end{bmatrix}=\begin{bmatrix} \boldsymbol{\mathsf{M}}\\ \boldsymbol{\mathsf{N}} \end{bmatrix}\boldsymbol{f}- \begin{bmatrix} \boldsymbol{\mathsf{M}}_T\\ \boldsymbol{\mathsf{N}}_T \end{bmatrix}\boldsymbol{m}, \end{equation}
exploiting the symmetry of the two-sphere problem to simplify the form of the hydrodynamic force–velocity relations. Writing 
$\boldsymbol {r}=\boldsymbol {x}_2-\boldsymbol {x}_1$ for the position of the passive sphere relative to that of the control sphere, the resulting 
$3\times 3$ tensors 
$\boldsymbol{\mathsf{M}}$, 
$\boldsymbol{\mathsf{N}}$, 
$\boldsymbol{\mathsf{M}}_T$, and 
$\boldsymbol{\mathsf{N}}_T$ are given exactly by
\begin{gather} \boldsymbol{\mathsf{M}}=M_\parallel{}\frac{\boldsymbol{r}\boldsymbol{r}^{\rm T}}{\left\lVert \boldsymbol{r}\right\rVert^2}+M_\perp{} \left(\boldsymbol{\mathsf{I}}-\frac{\boldsymbol{r}\boldsymbol{r}^{\rm T}}{\left\lVert \boldsymbol{r}\right\rVert^2}\right)\quad \text{and}\quad \boldsymbol{\mathsf{N}}=N_\parallel{}\frac{\boldsymbol{r}\boldsymbol{r}^{\rm T}}{r^2}+N_\perp{} \left(\boldsymbol{\mathsf{I}}-\frac{\boldsymbol{r}\boldsymbol{r}^{\rm T}}{\left\lVert \boldsymbol{r}\right\rVert^2}\right), \end{gather}
\begin{gather}\boldsymbol{\mathsf{M}}_T= M_T\frac{\boldsymbol{\epsilon}\boldsymbol{\cdot}\boldsymbol{r}}{\left\lVert \boldsymbol{r}\right\rVert} \quad\text{and}\quad \boldsymbol{\mathsf{N}}_T= N_T\frac{\boldsymbol{\epsilon}\boldsymbol{\cdot}\boldsymbol{r}}{\left\lVert \boldsymbol{r}\right\rVert}, \end{gather}
respectively, by simple symmetry arguments. Here, 
$\boldsymbol{\mathsf{I}}$ is the 
$3\times 3$ identity tensor and 
$\boldsymbol {\epsilon }$ is the Levi-Civita tensor. The scalar mobility coefficients 
$M_\parallel$, 
$M_\perp$, 
$N_\parallel$, 
$N_\perp$, 
$M_T$ and 
$N_T$ are all functions of only 
$r\colon= \left \lVert \boldsymbol {r}\right \rVert$, owing to the symmetry of the problem. Further, physical interpretation of each of the coefficients suggests that all but 
$N_T$ are strictly positive, and that 
$N_T<0$. However, although positivity may be rigorously deduced for 
$M_\parallel$ and 
$M_\perp$ via calculations of energy dissipation, which we later present in § 3.3, we are not aware of similar reasoning that applies to the remaining coefficients. Therefore, to support this intuition, we numerically evaluate these coefficients in Appendix A, shown in figure 9, with these findings being consistent with expected signs of these mobility coefficients. Of note, the convergence of the calculation of 
$M_T$ described in Appendix A is not sufficiently fast to conclude that 
$M_T$ is strictly positive for 
$\lambda \ll 1$, though this does not affect our later explorations.
If we additionally introduce 
$\boldsymbol{\mathsf{m}} = \boldsymbol{\mathsf{M}} - \boldsymbol{\mathsf{N}}$ and 
$\boldsymbol{\mathsf{m}}_T = \boldsymbol{\mathsf{M}}_T - \boldsymbol{\mathsf{N}}_T$, each functions of 
$\boldsymbol {r}$, we may further simplify the above into the control-affine form
\begin{equation} \frac{\mathrm{d}\tilde{\boldsymbol{X}}}{\mathrm{d}t}=\sum_{j} f_j(t)\boldsymbol{g}_j(\tilde{\boldsymbol{X}})+\sum_{j} m_j(t)\boldsymbol{h}_j(\tilde{\boldsymbol{X}}), \end{equation}
where 
$j\in \{x,y,z\}$ and we henceforth adopt these natural Latin subscripts to refer to the components of forces and torques that act in the directions of the right-handed orthonormal triad 
$\{\boldsymbol {e}_x,\boldsymbol {e}_y,\boldsymbol {e}_z\}$, which forms the fixed basis of the laboratory frame. We similarly index 
$\boldsymbol {g}_j$ and 
$\boldsymbol {h}_j$, with 
$\boldsymbol {g}_j$ and 
$\boldsymbol {h}_j$ being the fields associated with 
$f_j$ and 
$m_j$, respectively. Here, the modified state and the fields are defined by
\begin{equation}
\tilde{\boldsymbol{X}}=\begin{bmatrix} \boldsymbol{x}_1 \\
\boldsymbol{r} \end{bmatrix},\quad [\boldsymbol{g}_x,\boldsymbol{g}_y,
\boldsymbol{g}_z]=\begin{bmatrix}
\boldsymbol{\mathsf{M}}(\boldsymbol{r}) \\
-\boldsymbol{\mathsf{m}}(\boldsymbol{r})
\end{bmatrix},\quad \text{and}\quad [\boldsymbol{h}_x,\boldsymbol{h}_y,\boldsymbol{h}_z]=\begin{bmatrix}
-\boldsymbol{\mathsf{M}}_T(\boldsymbol{r}) \\
\boldsymbol{\mathsf{m}}_T(\boldsymbol{r})
\end{bmatrix}. \end{equation}
More verbosely, defining 
$m_\parallel = M_\parallel - N_\parallel$, 
$m_\perp = M_\perp - N_\perp$ and 
$m_T = M_T - N_T \neq M_T$, each functions of 
$r$, the control-affine system of (3.4) may be written as
\begin{equation} \frac{\mathrm{d}\tilde{\boldsymbol{X}}}{\mathrm{d}t} =\frac{\mathrm{d}}{\mathrm{d}t} \begin{bmatrix} \boldsymbol{x}_1 \\ \boldsymbol{r} \end{bmatrix}= \begin{bmatrix} M_\parallel{}(r)\dfrac{\boldsymbol{r}\boldsymbol{r}^{\rm T}}{\left\lVert \boldsymbol{r}\right\rVert^2}+M_\perp{}(r) \left(\boldsymbol{\mathsf{I}}-\dfrac{\boldsymbol{r}\boldsymbol{r}^{\rm T}}{\left\lVert \boldsymbol{r}\right\rVert^2}\right) \\ -m_\parallel{}(r)\dfrac{\boldsymbol{r}\boldsymbol{r}^{\rm T}}{\left\lVert \boldsymbol{r}\right\rVert^2}-m_\perp{}(r) \left(\boldsymbol{\mathsf{I}}-\dfrac{\boldsymbol{r}\boldsymbol{r}^{\rm T}}{\left\lVert \boldsymbol{r}\right\rVert^2}\right) \end{bmatrix}\boldsymbol{f} + \begin{bmatrix} -M_T(r)\dfrac{\boldsymbol{\epsilon}\boldsymbol{\cdot}\boldsymbol{r}}{r}\\ m_T(r)\dfrac{\boldsymbol{\epsilon}\boldsymbol{\cdot}\boldsymbol{r}}{r} \end{bmatrix}\boldsymbol{m}. \end{equation}
In the remainder of this section, we consider only this relative two-sphere control system and, therefore, drop the tilde on the relative composite state vector for notational convenience, hereafter representing the relative state by 
$\boldsymbol {X}$.
3.2. Establishing controllability
3.2.1. Control by external force
We first consider the controllability of the two-sphere system when 
$\boldsymbol {m}=\boldsymbol {0}$, so that only an external force acts on the control sphere through the fields 
$\boldsymbol {g}_x$, 
$\boldsymbol {g}_y$ and 
$\boldsymbol {g}_z$. As we are therefore seeking to control a six-dimensional system with three scalar controls, computation of the Lie brackets of the 
$\boldsymbol {g}_i$ is warranted. Writing 
$\boldsymbol {r} = (r_x,r_y,r_z)$ with respect to 
$\{\boldsymbol {e}_x,\boldsymbol {e}_y,\boldsymbol {e}_z\}$, we compute
\begin{equation}[{\boldsymbol{g}_i},{\boldsymbol{g}_j}]\left(\tilde{\boldsymbol{X}}\right)=\begin{bmatrix} \left(m_\parallel{} M'_\perp{-}m_\perp{} \dfrac{M_\parallel{}-M_\perp{}}{r} \right) \dfrac{r_j\boldsymbol{e}_i-r_i\boldsymbol{e}_j}{r}\\ -\left( m_\parallel{} m'_\perp{-}m_\perp{} \dfrac{m_\parallel{}-m_\perp{}}{r} \right) \dfrac{r_j\boldsymbol{e}_i-r_i\boldsymbol{e}_j}{r} \end{bmatrix}, \end{equation}
where a prime indicates a derivative with respect to 
$r$.
The evaluation of the fields 
$\boldsymbol {g}_i$ and the second-order brackets 
$[{\boldsymbol {g}_i},{\boldsymbol {g}_j}]$ at a general state 
$\boldsymbol {X}\in P\subset \mathbb {R}^6$ would yield notationally cumbersome expressions for the Lie bracket, which would be a barrier to further analysis. However, owing to the symmetry of the two-sphere problem, we may, without loss of generality, evaluate the system at a particular state 
$\boldsymbol {X}^{\ast }{}$, in which the control sphere is located at the origin and the target sphere lies on the 
$\boldsymbol {e}_x{}$ axis at a distance 
$r$, i.e. 
$\boldsymbol {x}_1 = \boldsymbol {0}$ and 
$\boldsymbol {r} = r\boldsymbol {e}_x{}$. This reduced configuration is clearly equivalent to a general state up to rotation and translation, with controllability invariant under such transformations in this context. In this parameterised state, the fields 
$\boldsymbol {g}_i$ reduce to
\begin{gather} \boldsymbol{g}_x(\boldsymbol{X}^{{\ast}}{})= [ M_\parallel{}, 0, 0, -m_\parallel{}, 0, 0]^{\rm T}, \end{gather}
\begin{gather}\boldsymbol{g}_y(\boldsymbol{X}^{{\ast}}{})= \left[ 0, M_\perp{}, 0, 0, -m_\perp{}, 0\right]^{\rm T}, \end{gather}
\begin{gather}\boldsymbol{g}_z(\boldsymbol{X}^{{\ast}}{})= \left[ 0, 0, M_\perp{}, 0, 0, -m_\perp{}\right]^{\rm T}. \end{gather}
Evaluating the second-order brackets at 
$\boldsymbol {X}^{\ast }{}$, we have
\begin{align}
[{\boldsymbol{g}_x},{\boldsymbol{g}_y}](\boldsymbol{X}^{{\ast}}{})&=\begin{bmatrix}
0\\ -m_\parallel{} M'_\perp + m_\perp{}
\dfrac{M_\parallel{}-M_\perp{}}{r} \\ 0 \\ 0 \\
m_\parallel{} m'_\perp - m_\perp{}
\dfrac{m_\parallel{}-m_\perp{}}{r} \\ 0 \end{bmatrix},\\
[{\boldsymbol{g}_x},{\boldsymbol{g}_z}](\boldsymbol{X}^{{\ast}}{})&=\begin{bmatrix}
0 \\ 0\\ -m_\parallel{} M'_\perp + m_\perp{}
\dfrac{M_\parallel{}-M_\perp{}}{r} \\ 0 \\ 0 \\
m_\parallel{} m'_\perp - m_\perp{}
\dfrac{m_\parallel{}-m_\perp{}}{r} \end{bmatrix}, \end{align}
whereas 
$[{\boldsymbol {g}_y},{\boldsymbol {g}_z}](\boldsymbol {X}^{\ast }{})=\boldsymbol {0}$. Clearly, additional elements are required in order to span a space of dimension six. Hence, we consider the third-order Lie brackets given by
\begin{equation} [{\boldsymbol{g}_i},{[{\boldsymbol{g}_j},{\boldsymbol{g}_k}]}]=\begin{pmatrix} \boldsymbol{L}_{ijk}\\ \boldsymbol{\ell}_{ijk} \end{pmatrix}, \end{equation}
where the three-dimensional field 
$\boldsymbol {L}_{ijk}$ is explicitly given by
\begin{align} \boldsymbol{L}_{ijk} &=
\frac{M_{{\parallel}}-M_\perp{}}{r^3} \left(m_\parallel{}
m_\perp{}'-m_\perp{}\frac{m_\parallel{}-m_\perp{}}{r}\right)
(r_jr_k\boldsymbol{e}_i+r_j\delta_{ik}\boldsymbol{r}-r_ir_k\boldsymbol{e}_j-r_i\delta_{jk}\boldsymbol{r})
\nonumber\\ &\quad -\frac{m_\perp{}}{r} \left(m_\parallel{}
M_\perp{}'-m_\perp{}\frac{M_\parallel{}-M_\perp{}}{r}\right)
\left[
(\delta_{kj}\boldsymbol{e}_i-\delta_{ik}\boldsymbol{e}_j) +
\left(
\frac{r_ir_k}{r^2}\boldsymbol{e}_j-\frac{r_jr_k}{r^2}\boldsymbol{e}_i
\right)\right] \nonumber\\ &\quad
-\frac{m_\parallel{}}{r^2}\left[\left(m_\parallel{}
M_\perp{}'-m_\perp{}\frac{M_\parallel{}-M_\perp{}}{r}\right)'
\right] (r_jr_k\boldsymbol{e}_i-r_ir_k\boldsymbol{e}_j)
\end{align}
and 
$\boldsymbol {\ell }_{ijk}$ is obtained by replacing 
$M_\parallel {}$ and 
$M_\perp {}$ with 
$-m_\parallel {}$ and 
$-m_\perp {}$, respectively, in the above expression. Here, the Kronecker delta is such that 
$\delta _{ij}=\boldsymbol {e}_i\boldsymbol {\cdot }\boldsymbol {e}_j$ and primes again denote differentiation with respect to 
$r$, with 
$i,j,k\in \{x,y,z\}$. As before, this cumbersome expression simplifies significantly when evaluating at 
$\boldsymbol {X}^{\ast }{}$. In particular, 
$[{\boldsymbol {g}_y},{[{\boldsymbol {g}_x},{\boldsymbol {g}_y}]}]$ conditionally generates the additional dimension required to yield controllability and is explicitly given by
\begin{align} &
[{\boldsymbol{g}_y},{[{\boldsymbol{g}_x},{\boldsymbol{g}_y}]}](\boldsymbol{X}^{{\ast}}{})
\nonumber\\ &\quad =\begin{bmatrix}
-(M_\parallel{}-M_\perp{})\left(\dfrac{m_\parallel{}
m'_\perp}{r}-\dfrac{m_\perp{}(m_\parallel{}-m_\perp{})}{r^2}\right)
-m_\perp{}\left(\dfrac{m_\parallel{}
M'_\perp}{r}-\dfrac{m_\perp{}(M_\parallel{}-M_\perp{})}{r^2}\right)\\
0 \\ 0 \\ m_\parallel{}\left(\dfrac{m_\parallel{}
m'_\perp}{r}-\dfrac{m_\perp{}(m_\parallel{}-m_\perp{})}{r^2}\right)\\
0 \\ 0 \end{bmatrix}. \end{align}
Additional third-order brackets are simply linear combinations of the lower-order brackets, with the exception of 
$[{\boldsymbol {g}_z},{[{\boldsymbol {g}_x},{\boldsymbol {g}_z}]}]$, which is parallel to 
$[{\boldsymbol {g}_y},{[{\boldsymbol {g}_x},{\boldsymbol {g}_y}]}]$ by the symmetry of the two-sphere problem. Together, the six terms of (3.10), (3.11a,b) and (3.14) in general span a space of dimension six, with this property holding unless these fields are linearly dependent. This condition may be succinctly investigated by considering the determinant of the 
$6\times 6$ controllability matrix
\begin{equation}
\boldsymbol{\mathsf{C}}= [\boldsymbol{g}_x,
\boldsymbol{g}_y, \boldsymbol{g}_z,
[{\boldsymbol{g}_x},{\boldsymbol{g}_y}],
[{\boldsymbol{g}_x},{\boldsymbol{g}_z}],
[{\boldsymbol{g}_y},{[{\boldsymbol{g}_x},{\boldsymbol{g}_y}]}]], \end{equation}
with the full-rank controllability condition being satisfied if 
$\det {\boldsymbol{\mathsf{C}}}\neq 0$. Of particular note, 
$\det {\boldsymbol{\mathsf{C}}}=0$ only implies that the controllability condition is not satisfied by this particular controllability matrix, with the dimension of the Lie algebra being bounded below by the rank of any controllability matrix. Here, this determinant is explicitly given by
\begin{equation} \det{\boldsymbol{\mathsf{C}}} = \frac{m_\parallel{}^4 m_\perp{}^6}{r} \left[ \frac{1}{r}\left( \frac{M_\parallel{}}{m_\parallel{}} -\frac{M_\perp{}}{m_\perp{}}\right) - \left(\frac{M_\perp{}}{m_\perp{}}\right)' \right]^3, \end{equation}
where we recall that all mobility coefficients are functions of the separation 
$r$ for a given sphere size ratio 
$\lambda$.
In general, closed-form expressions for the mobility coefficients are unavailable, though the coefficients may be obtained as infinite series (Jeffrey & Onishi Reference Jeffrey and Onishi1984). Therefore, as detailed in Appendix A, we numerically evaluate these coefficients and compute the value of 
$\det {\boldsymbol{\mathsf{C}}}$, which we show in figure 3 as solid black lines for a range of separations 
$d=r - (a_1+a_2) = r - (1+\lambda )$ and relative sphere radii 
$\lambda$. For 
$d>0$, we observe that the determinant does not vanish. Thus, we conclude that the controllability matrix 
$\boldsymbol{\mathsf{C}}$ is of rank six throughout the entire admissible state space 
$P= \{\boldsymbol {X}\in \mathbb {R}^6: d>0\}$, so that 
$\dim {\boldsymbol{\mathsf{B}}(\boldsymbol {X})}=6$ on all of 
$P$ and, hence, 
$S = P$. Therefore, by the Rashevski–Chow theorem, the force-driven two-sphere problem is controllable everywhere, excluding the trivial case when the spheres intersect with one another.
Figure 3. Values of the determinant 
$\det {\boldsymbol{\mathsf{C}}}$ as a function of separation 
$d = r - (1+\lambda )$ for different relative sphere radii 
$\lambda$. The determinants corresponding to the far-field approximation are shown as dotted red curves, whereas those corresponding to the full coefficients are shown in black, each evaluated numerically to a precision beyond the resolution of these plots. The strict positivity of each of these curves as 
$d\rightarrow \infty$ is confirmed by the leading-order expression of (3.18). The dotted line in (c) corresponds to 
$\det {\boldsymbol{\mathsf{C}}}=0$.
Though we have used accurate expressions for the mobility coefficients in the above analysis, the widespread use of leading-order far-field approximations in the study of motion in Stokes flow motivates repeating the above calculation with the far-field analogues of the mobility coefficients. In this case, having non-dimensionalised forces relative to the mobility coefficient 
$M_\perp {}$, the far-field approximations to the mobility coefficients are simply
\begin{equation} M_\parallel{} = M_\perp{} = 1,\quad m_\parallel{} = 1- \frac{3}{2r},\quad m_\perp{} = 1 - \frac{3}{4r},\end{equation}
independent of 
$\lambda$ in the far-field, with the far-field approximation to the determinant being
\begin{equation} \det{\boldsymbol{\mathsf{C}}}=\frac{27}{512}\frac{(2r-3)(8r-9)^3}{(4r-3)r^9}, \end{equation}
shown as red dotted lines in figure 3. Should this determinant vanish, this would entail that there is a linear dependence in the constituent Lie brackets and a non-trivial nullspace, corresponding to the existence of directions that are not controllable via these Lie brackets. Hence, as can be seen from the expression of (3.18), and as is evident in the figure, this far-field approximation does not yield an everywhere-controllable system for 
$\lambda >1/2$ and this choice of 
$\boldsymbol{\mathsf{C}}$, with the determinant vanishing at 
$r=3/2$ and 
$r=9/8$. In these cases, the image of 
$\boldsymbol{\mathsf{C}}$ corresponds to the reachable directions. When 
$r=3/2$, the direction 
$[0,0,0,1,0,0]^{\rm T}$ is unreachable, corresponding to increasing sphere separation, whereas, when 
$r=9/2$, there are two unreachable directions, 
$[0,1,0,0,3,0]^{\rm T}$ and 
$[0,0,1,0,0,3]^{\rm T}$, which each correspond to particular weighted motions of both the control sphere and passive sphere perpendicular to the displacement vector 
$\boldsymbol {r}$. The zeros of the determinant can be seen in figure 3(c) and correspond to 
$d = 1/2 - \lambda$ and 
$d = 1/8 - \lambda$, respectively. Despite this, as should be expected, good agreement between the two determinant calculations is evident for 
$d\gtrsim 2$, though the far-field theory fails to predict the local maxima in the determinant seen in figures 3(a) and 3(b). Of note, when the determinant is non-zero but near vanishing, this can be associated with some notion of poor controllability in certain directions, which we discuss and explore further in § 3.3.
3.2.2. Control via an external torque
Complimentary to the above enquiry, we now seek to evaluate the controllability of a torque-driven system, taking 
$\boldsymbol {f}=\boldsymbol {0}$ instead of 
$\boldsymbol {m}=\boldsymbol {0}$. As in the analysis of the force-driven system, we exploit the symmetry of the two-sphere problem to evaluate 
$B(\boldsymbol {X})$ on a reduced set of states, which we again denote by 
$\boldsymbol {X}^{\ast }{}$, without loss of generality. First, however, we compute the second-order Lie brackets, which may be written as
\begin{equation}
[{\boldsymbol{h}_i},{\boldsymbol{h}_j}]=
\begin{bmatrix}
\dfrac{-2m_TM_T}{r^2}(r_j\boldsymbol{e}_i-r_i\boldsymbol{e}_j)\\
\dfrac{2m_T^2}{r^2}(r_j\boldsymbol{e}_i-r_i\boldsymbol{e}_j)
\end{bmatrix}, \end{equation}
where 
$i,j\in \{x,y,z\}$. Evaluated at a state in the reduced space, these become
\begin{gather}
{[}{\boldsymbol{h}_x},{\boldsymbol{h}_y}](\boldsymbol{X}^{{\ast}}{})=
\left[0, \frac{2m_TM_T}{r}, 0, 0, -\frac{2m_T^2}{r}, 0\right]^{\rm T} \end{gather}
\begin{gather}
{[}{\boldsymbol{h}_x},{\boldsymbol{h}_z}](\boldsymbol{X}^{{\ast}}{})=
\left[0, 0, \frac{2m_TM_T}{r}, 0, 0, -\frac{2m_T^2}{r}\right]^{\rm T}, \end{gather}
and 
$[{\boldsymbol {h}_y},{\boldsymbol {h}_z}](\boldsymbol {X}^{\ast }{})=\boldsymbol {0}$, whereas the fields corresponding to the control are simply
\begin{gather} \boldsymbol{h}_x(\boldsymbol{X}^{{\ast}}{})= \left[ 0, 0, 0, 0, 0, 0\right]^{\rm T} \end{gather}
\begin{gather}\boldsymbol{h}_y(\boldsymbol{X}^{{\ast}}{})= \left[ 0, 0, M_T, 0, 0, -m_T\right]^{\rm T}, \end{gather}
\begin{gather}\boldsymbol{h}_z(\boldsymbol{X}^{{\ast}}{})= \left[ 0, -M_T, 0, 0, m_T, 0 \right]^{\rm T}. \end{gather}
Immediately, and perhaps surprisingly, we see that 
$[{\boldsymbol {h}_x},{\boldsymbol {h}_y}]$ is parallel to 
$\boldsymbol {h}_z$, whereas 
$[{\boldsymbol {h}_x},{\boldsymbol {h}_z}]$ is parallel to 
$\boldsymbol {h}_y$, so that these Lie brackets do not generate additional directions in this example. Thus, so far, we have demonstrated only that 
$\dim {B(\boldsymbol {X}^{\ast }{})}\geq 2$, spanned by 
$\boldsymbol {h}_y$ and 
$\boldsymbol {h}_z$, noting that 
$m_T\neq M_T$. In seeking higher-dimensional control, one might be tempted to consider higher-order Lie brackets, such as 
$[{\boldsymbol {h}_x},{[{\boldsymbol {h}_x},{\boldsymbol {h}_y}]}]$.
However, this lower bound on 
$\dim {B(\boldsymbol {X}^{\ast }{})}$ is in fact sharp, though drawing such a conclusion from direct consideration of 
$\text {Lie}(\boldsymbol {h}_x,\boldsymbol {h}_y,\boldsymbol {h}_z)$ can be cumbersome in general, in this case requiring inductive arguments on the iterated Lie brackets; hence, we reason directly from the dynamical system of (3.6). First, it is clear from symmetry that the torque-driven system is not controllable along the 
$\boldsymbol {r}$ direction, with direct calculation of 
$\mathrm {d}\left \lVert \boldsymbol {r}\right \rVert ^2/\mathrm {d}t$ via (3.6) highlighting that the distance between the spheres is unchanged throughout the motion. Explicitly, recalling (3.6), we have
\begin{equation} \frac{1}{2}\frac{\mathrm{d}\left\lVert \boldsymbol{r}\right\rVert^2}{\mathrm{d}t} = \boldsymbol{r}^{\rm T}\frac{\mathrm{d}\boldsymbol{r}}{\mathrm{d}t}= \frac{m_T(r)}{r}\boldsymbol{r}^{\rm T}(\boldsymbol{\epsilon}\boldsymbol{\cdot}\boldsymbol{r})\boldsymbol{m}=0 \end{equation}
by the antisymmetry of the Levi-Civita tensor. Thus, the dimension of 
$B(\boldsymbol {X})$ is at most five. Further, introducing the weighted point 
$\boldsymbol {x}_h=(m_T\boldsymbol {x}_1+M_T\boldsymbol {r})$, simple calculation yields that 
$\mathrm {d}\boldsymbol {x}_h/\mathrm {d}t=\boldsymbol {0}$, irrespective of the applied torque, with the mobility coefficients being explicit functions of only geometry. In more detail, noting from (3.6) that 
$\mathrm {d}\boldsymbol {r}/\mathrm {d}t = - (m_T/M_t)\,\mathrm {d}\boldsymbol {x}/\mathrm {d}t$ when 
$\boldsymbol {f}=\boldsymbol {0}$, we compute
\begin{equation} \frac{\mathrm{d}\boldsymbol{x}_h}{\mathrm{d}t} = m_T \frac{\mathrm{d}\boldsymbol{x}}{\mathrm{d}t} + M_T \frac{\mathrm{d}\boldsymbol{r}}{\mathrm{d}t} = m_T \frac{\mathrm{d}\boldsymbol{x}}{\mathrm{d}t} + M_T \left(-\frac{m_T}{M_T}\frac{\mathrm{d}\boldsymbol{x}}{\mathrm{d}t}\right)= \boldsymbol{0}, \end{equation}
recalling that we have shown that 
$r$ is constant in time in this scenario, so that the mobility coefficients are also constant in time. Hence, with this last constraint imposing three scalar conditions, we can view the system as evolving in a smaller state space 
$P'$ of dimension two, namely
\begin{equation} P'=\{\boldsymbol{X}\in\mathbb{R}^6:\ r=\textrm{const.},\ \boldsymbol{x}_h=\textrm{const.}\}, \end{equation}
and, further, note that the controllable state set within this new state space is found to be 
$S=P'$. Therefore, by the Rashevski–Chow theorem applied to this subspace, the torque-driven system is controllable in 
$S$, indicating that one may control the two spheres around the point 
$\boldsymbol {x}_h$ with a fixed distance between them.
3.3. Efficient control
As demonstrated in the previous section, control in the entire state space cannot be effected by a simple torque control. Hence, with controllability naturally being a desirable property, we consider only the controllable force-driven problem in further detail. In particular, the presence of local maxima in figure 3 is suggestive of some notion of optimally efficient control, with the determinant providing an estimate for the volume of the state space explored by the columns of the controllability matrix 
$\boldsymbol{\mathsf{C}}$. This is a measure that we can conceptually relate, in the framework of geometric control theory, to the size of the sub-Riemannian ball estimated by the ball-box theorem (Bellaïche Reference Bellaïche1997), though notably in a different context. Motivated by this, we now look to examine the mechanical efficiency of the force-driven control problem of § 3.2.1.
Let us consider the energy dissipation 
$E$ that occurs whilst attempting to control the spheres in a particular direction for a small time 
$\tau$ from a state 
$\boldsymbol {X}^{\ast }{}$. We focus on the six fields that make up the columns of 
$\boldsymbol{\mathsf{C}}$ in (3.15), as they span the space around 
$\boldsymbol {X}^{\ast }{}$, whereas higher-order brackets are generally associated with lower efficiencies, an example of which we show later. For a given motion, the energy dissipation rate 
$\dot {E}$ is given by
\begin{equation} \dot{E}=\boldsymbol{U}_1\boldsymbol{\cdot}\boldsymbol{f}_1,\end{equation}which can be written as the quadratic form
\begin{equation} \dot{E}(t)= \frac{\mathrm{d}\boldsymbol{X}_1}{\mathrm{d}t}\boldsymbol{\cdot}\boldsymbol{f}=\boldsymbol{f}^{\rm T}\left[ M_\parallel{}\frac{\boldsymbol{r}\boldsymbol{r}^{\rm T}}{\left\lVert \boldsymbol{r}\right\rVert^2}+ M_\perp{}\left(\boldsymbol{\mathsf{I}}-\frac{\boldsymbol{r}\boldsymbol{r}^{\rm T}}{\left\lVert \boldsymbol{r}\right\rVert^2}\right)\right]\boldsymbol{f}\geq0. \end{equation}
The total energy dissipation 
$E$ can then be calculated by the time integral of 
$\dot {E}$.
Let 
$\gamma$ denote the magnitude of the external force and, in the first instance, set 
$\boldsymbol {f} = [\gamma,0,0]^{\rm T}$, independent of time. The six-dimensional displacement after time 
$\tau$ under this control is simply given by 
$\boldsymbol {\varDelta }_x = \varepsilon \boldsymbol {g}_x(\boldsymbol {X}^{\ast }{})$, where we have neglected higher-order terms in 
$\varepsilon \colon= \gamma \tau \ll 1$. Defining the projection map 
$\left \langle \cdot \right \rangle _{\boldsymbol {r}}$ such that 
$\left \langle [\boldsymbol {a};\boldsymbol {b}]\right \rangle _{\boldsymbol {r}} = \boldsymbol {b}\in \mathbb {R}^3$ for 
$[\boldsymbol {a};\boldsymbol {b}]\in \mathbb {R}^6$, the change in relative position during this motion may be written simply as 
$\left \langle \boldsymbol {\varDelta }_x\right \rangle _{\boldsymbol {r}}$. The energy 
$E_x$ corresponding to this process is given by the integral of (3.29), explicitly
\begin{equation} E_x=\int_0^\tau \gamma^2 M_\parallel{}\,\mathrm{d}t=\gamma^2\tau M_\parallel{}=\varepsilon \gamma M_\parallel{} \end{equation}
to leading order in 
$\varepsilon$. We now define the mechanical efficiency for this process as the relative change in displacement 
$\boldsymbol {r}$ per unit energy consumption, the reciprocal of the widely used ‘cost of transport’, written
\begin{equation} \eta_{x}=
\frac{\langle \boldsymbol{\varDelta}_x\rangle_{\boldsymbol{r}}\boldsymbol{\cdot}\langle\boldsymbol{g}_x\rangle_{\boldsymbol{r}}}{E_x
\lVert \langle\boldsymbol{g}_x\rangle_{\boldsymbol{r}}\rVert},
\end{equation}
where displacement is being measured in the 
$\boldsymbol {g}_x$ direction and we are suppressing the argument 
$\boldsymbol {X}^{\ast }{}$ of 
$\boldsymbol {g}_x$ for brevity. Making use of the expression of (3.10), the unit vector 
$\langle \boldsymbol {g}_x\rangle _{\boldsymbol {r}}/\lVert \langle \boldsymbol {g}_x\rangle _{\boldsymbol {r}}\rVert$ is simply equal to 
$-\boldsymbol {e}_x{}$, though we note that this may not always be the case when considering the far-field approximation of the mobility coefficients.
Similarly, we may calculate the displacement and energy consumption for the other control fields 
$\boldsymbol {g}_y$ and 
$\boldsymbol {g}_z$, giving
\begin{equation} \boldsymbol{\varDelta}_y = \varepsilon \boldsymbol{g}_y(\boldsymbol{X}^{{\ast}}{}),\quad \boldsymbol{\varDelta}_z =\varepsilon \boldsymbol{g}_z(\boldsymbol{X}^{{\ast}}{}),\quad E_y=E_z=\varepsilon \gamma M_\perp{}. \end{equation}
We may then define 
$\eta _{y}$ and 
$\eta _{z}$ as the relative displacement in the 
$\boldsymbol {g}_y(\boldsymbol {X}^{\ast }{})$ and 
$\boldsymbol {g}_z(\boldsymbol {X}^{\ast }{})$ directions, respectively, per unit mechanical energy consumption by the control, analogously to (3.31). These definitions can be easily extended to the other elements of 
$B(\boldsymbol {X})$, recalling that Lie brackets represent piecewise sequences of applied external forces. The displacements and energies associated with 
$[{\boldsymbol {g}_x},{\boldsymbol {g}_y}]$ and 
$[{\boldsymbol {g}_x},{\boldsymbol {g}_z}]$ are, to leading order in 
$\varepsilon$,
\begin{gather}
\boldsymbol{\varDelta}_{xy}=\varepsilon^2
[{\boldsymbol{g}_x},{\boldsymbol{g}_y}](\boldsymbol{X}^{{\ast}}{}),\quad
\boldsymbol{\varDelta}_{xz} = \varepsilon^2
[{\boldsymbol{g}_x},{\boldsymbol{g}_z}](\boldsymbol{X}^{{\ast}}{}),\quad
E_{xy}=E_{xz}=2\varepsilon \gamma (M_\parallel{} +
M_\perp{}), \end{gather}
respectively. The iterated Lie bracket 
$[{\boldsymbol {g}_y},{[{\boldsymbol {g}_x},{\boldsymbol {g}_y}]}]$ gives rise to a displacement on the order of 
$\varepsilon ^3$,
\begin{equation} \boldsymbol{\varDelta}_{yxy}=\varepsilon^3[{\boldsymbol{g}_y},{[{\boldsymbol{g}_x},{\boldsymbol{g}_y}]}](\boldsymbol{X}^{{\ast}}{}), \end{equation}with brackets of higher order giving rise to displacements of even higher order. The energy consumption of this term is given by
\begin{equation} E_{yxy}=2\varepsilon \gamma (2M_\parallel{}+3M_\perp{}), \end{equation}
which we note is the same order of 
$\varepsilon$ as the other Lie brackets.
We now similarly define the mechanical efficiencies of each of the considered Lie brackets, denoting them by 
$\eta _{xy}$, 
$\eta _{xz}$ and 
$\eta _{yxy}$ for each of 
$[{\boldsymbol {g}_x},{\boldsymbol {g}_y}]$, 
$[{\boldsymbol {g}_x},{\boldsymbol {g}_z}]$ and 
$[{\boldsymbol {g}_y},{[{\boldsymbol {g}_x},{\boldsymbol {g}_y}]}]$, respectively, the latter of which is given explicitly by
\begin{equation}
\eta_{yxy}=\varepsilon^2\frac{\langle\boldsymbol{\varDelta}_{yxy}\rangle_{\boldsymbol{r}}\boldsymbol{\cdot}
\langle\boldsymbol{g}_{yxy}\rangle_{\boldsymbol{r}}}{\gamma
(4M_\parallel{}+6M_\perp{})\lVert \langle\boldsymbol{g}_{yxy}\rangle_{\boldsymbol{r}}\rVert}, \end{equation}
making use of the compact notation 
$\boldsymbol {g}_{yxy}\equiv [{\boldsymbol {g}_y},{[{\boldsymbol {g}_x},{\boldsymbol {g}_y}]}]$. In figure 4, we compute 
$\eta _{yxy}$ for various relative radii 
$\lambda$, each normalised by 
$\varepsilon ^2/\gamma$. For completeness, we also repeat the above calculations using the far-field approximation for the hydrodynamics, with the resulting efficiencies plotted as dotted red curves in figure 4. Of note, the direction of 
$[{\boldsymbol {g}_y},{[{\boldsymbol {g}_x},{\boldsymbol {g}_y}]}](\boldsymbol {X}^{\ast }{})$ reverses around 
$d=0.5$ in figure 4(c) for the far-field case, giving rise to a non-smooth efficiency. Returning to the full mobility coefficients, we also note the presence of a local maximum in each of these plots, a property shared with 
$\eta _{xy}$ (not shown), which corresponds to the energy-optimal distance at which this control should be applied.
Figure 4. Values of the mechanical efficiency 
$\eta _{yxy}$ normalised by 
$\varepsilon ^2/\gamma$ for different sphere radius ratios 
$\lambda$. The results of the full calculations are shown in black, whereas the analogous results for the far-field hydrodynamic approximation are shown as dotted red curves.
We compute this optimal separation 
$d^{\ast {}}$ as a function of the relative radius 
$\lambda$, shown in figure 5(a), from which we observe that the optimal separation increases to 
$d^{\ast }\approx 1$ as the size of the passive particle decreases (
$\lambda \rightarrow 0$) and attains a local minimum around 
$\lambda =1$, when the spheres are of equal size. Further, as 
$\eta _{yxy}=O(\varepsilon ^2)$, whereas the other efficiencies are 
$O(1)$ or 
$O(\varepsilon )$, the optimisation of 
$\eta _{yxy}$ corresponds to improving the minimally efficient component of a control scheme utilising the considered elements of 
$\text {Lie}(\boldsymbol {g}_x,\boldsymbol {g}_y,\boldsymbol {g}_z)$, which may be desirable in the context of seeking out efficient schemes for sphere control. One might also attempt to optimise the other efficiencies, for example 
$\eta _{x}$ and 
$\eta _{y}$, though, in fact, these do not possess the local maxima of 
$\eta _{xy}$ and 
$\eta _{yxy}$ (not shown). Instead, these efficiencies are monotonically increasing functions of the sphere separation 
$d$, highlighting that simple controls, such as 
$\boldsymbol {f}=[\gamma,0,0]^{\rm T}$, are relatively inefficient methods for altering the separation 
$\boldsymbol {r}$ of the spheres when they are close to one another, as might be intuitively expected, which warrants the use of more-complex Lie bracket controls.
Figure 5. Efficiency, determinant and optimisation. (a) Optimal distance in terms of the mechanical efficiency 
$\eta _{yxy}$ as a function of relative sphere radius 
$\lambda$, with the 
$\lambda$ axis scaled logarithmically. (b) The relationship between the determinant of the controllability matrix and the mechanical efficiencies, plotted for 
$\lambda =1$. All quantities in (b) are normalised by their respective maxima.
Finally, in order to investigate a potential link between mechanical efficiency and the determinant of figure 3, we plot the various efficiencies against the corresponding values of 
$\det {\boldsymbol{\mathsf{C}}}$ in figure 5(b), fixing 
$\lambda =1$ and normalising each quantity with respect to its maximum. We observe a surprising correlation between the efficiency 
$\eta _{yxy}$ and 
$\det {\boldsymbol{\mathsf{C}}}$, with their maxima being approximately coincident, though this agreement is not present to the same extent in the other computed efficiencies. Thus, this partially supports the hypothesised link between efficiency and the determinant of the controllability matrix, though it remains pertinent to consider alternative definitions of efficiency and thoroughly explore the design of efficient controls.
4. Tracer limit
4.1. Far-field hydrodynamics
A natural limit of the force-controlled finite-size sphere scenario is for 
$\lambda \ll 1$, in which the passive spheres are of negligible size relative to the control sphere. These dynamics are of tracer particles being advected by the flow induced by a translating sphere, which, to a first approximation, is simply that of a Stokeslet. Explicitly, with the position of the control sphere as 
$\boldsymbol {x}$ and the relative positions of the passive tracers as 
$\boldsymbol {r}_1,\ldots,\boldsymbol {r}_{N-1}$, mirroring the previous notation, we have
\begin{equation} \frac{\mathrm{d}\boldsymbol{r}_k}{\mathrm{d}t} = \frac{3}{4}\left(\frac{\boldsymbol{\mathsf{I}}}{\left\lVert \boldsymbol{r}_k\right\rVert} + \frac{\boldsymbol{r}_k\boldsymbol{r}_k^{\rm T}}{\left\lVert \boldsymbol{r}_k\right\rVert^3}\right) \boldsymbol{f} - \boldsymbol{f}\end{equation}
for 
$k=1,\ldots,N-1$, where the final term arises from subtracting 
$\dot {\boldsymbol {x}}=\boldsymbol {f}$ to give the relative velocity of the tracers, noting that the dimensionless mobility coefficient is simply unity in the tracer limit. Written together, we have the 
$3N$-dimensional control system
\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t}\begin{bmatrix} \boldsymbol{r}_1 \\ \vdots \\ \boldsymbol{r}_{N-1} \\ \boldsymbol{x} \end{bmatrix}= \begin{bmatrix} \dfrac{3}{4}\left(\dfrac{\boldsymbol{\mathsf{I}}}{\left\lVert \boldsymbol{r}_1\right\rVert} + \dfrac{\boldsymbol{r}_1\boldsymbol{r}_1^{\rm T}}{\left\lVert \boldsymbol{r}_1\right\rVert^3}\right) - \boldsymbol{\mathsf{I}} \\ \vdots \\ \dfrac{3}{4}\left(\dfrac{\boldsymbol{\mathsf{I}}}{\left\lVert \boldsymbol{r}_{N-1}\right\rVert} + \dfrac{\boldsymbol{r}_{N-1}\boldsymbol{r}_{N-1}^{\rm T}}{\left\lVert \boldsymbol{r}_{N-1}\right\rVert^3}\right) - \boldsymbol{\mathsf{I}} \\ \boldsymbol{\mathsf{I}} \end{bmatrix}\boldsymbol{f} = \sum_{i\in\{x,y,z\}} f_i(t) \boldsymbol{g}_i(\boldsymbol{X}), \end{equation}
having concatenated the equations for tracer evolution from (4.1) along with the simple equation of motion 
$\mathrm {d}\boldsymbol {x}/\mathrm {d}t=\boldsymbol {f}$ for the control sphere, whose position is denoted by 
$\boldsymbol {x}$. We particularly note that the position of the control sphere can only be controlled directly, i.e. not with any Lie brackets that are second order or higher, as the 
$\boldsymbol {x}$ components of the fields 
$\boldsymbol {g}_i$ commute with one another, corresponding to the lower block of (4.2). Therefore, the problem of controllability reduces to asking if we can control the relative position of the tracers using only second-order or higher Lie brackets. For brevity, we define the restricted three-dimensional fields 
$\tilde {\boldsymbol {g}}_{i,k}$ as the 
$\boldsymbol {r}_k = (r_{x,k},r_{y,k},r_{z,k})$ components of 
$\boldsymbol {g}_i$. Fixing 
$k\in \{1,\ldots,N-1\}$, straightforward calculation shows that an iterated Lie bracket 
$\boldsymbol {l}$ of order 
$p$ of the fields 
$\tilde {\boldsymbol {g}}_{i,k}$ (in particular, 
$p=1$ corresponds to the original fields) has components of the form 
$l_j = L_j r_k^{-5(p-1)-3}$. Here, 
$L_j$ is a homogeneous polynomial of degree 
$3(p-1)+2$ in 
$r_{x,k}$, 
$r_{y,k}$, and 
$r_{z,k}$, whose coefficients depend on 
$(\frac {3}{4}-r_k)$, where 
$r_k=\left \lVert \boldsymbol {r}_k\right \rVert$. Moreover, given the polynomials 
$L_j$ of 
$\boldsymbol {l}$, the components of the bracket 
$[\tilde {\boldsymbol {g}}_{i,k},\boldsymbol {l}]$ can be obtained recursively via the relation
\begin{align} r_k^{5p+3} [\tilde{\boldsymbol{g}}_{i,k},\boldsymbol{l}]_j &= r_k^2 \left( \frac{3}{4}-r_k \right) \left (\frac{\partial L_j}{\partial r_{i,k}} r_k^2 - p r_{i,k} L_j \right ) \nonumber\\ &\quad + \frac{3}{4} \sum_{q=1}^{3} L_q (3 r_{i,k} r_{j,k} r_{q,k} - r_k^2 (\delta_{iq} r_{j,k} + \delta _{jq} r_{i,k} - \delta_{ij} r_{q,k})). \end{align}After some simplifications, this gives the order two brackets as
\begin{equation}[{\tilde{\boldsymbol{g}}_{i,k}},{\tilde{\boldsymbol{g}}_{j,k}}] = \frac{3 (9-8r_k)}{16 r_k^4} ( r_{j,k} \boldsymbol{e}_i - r_{i,k} \boldsymbol{e}_j). \end{equation}
The expression of (4.4) allows us to immediately conclude that the second-order Lie brackets span a space of dimension two, except when 
$r_k=9/8$, which we recall also prohibited controllability in the far-field case of § 3.2.1. Here, however, the inclusion of third-order brackets, in particular 
$[{\tilde {\boldsymbol {g}}_x},{[{\tilde {\boldsymbol {g}}_x},{\tilde {\boldsymbol {g}}_y}]}]$ and 
$[{\tilde {\boldsymbol {g}}_x},{[{\tilde {\boldsymbol {g}}_x},{\tilde {\boldsymbol {g}}_z}]}]$, is sufficient to ensure that the Lie brackets indeed span a space of dimension at least two.
In order to make further progress without lengthy calculation, we now restrict to the case with a single tracer, taking 
$N=2$. Hence, we suppress the second subscript of 
$k$ in 
$\tilde {\boldsymbol {g}}_{i,k}$, with 
$k$ simply being unity. With this simplification, we may again make use of the state-space reduction of § 3, exploiting the invariance of the controllability of this system to translation and rotation, evaluating the Lie brackets at a reduced configuration in which 
$\boldsymbol {x}=\boldsymbol {0}$ and 
$\boldsymbol {r}_1\equiv \boldsymbol {r} = r\boldsymbol {e}_x$, without loss of generality. In this reduced state, the space spanned by 
$[{\tilde {\boldsymbol {g}}_x},{[{\tilde {\boldsymbol {g}}_x},{\tilde {\boldsymbol {g}}_y}]}]$ and 
$[{\tilde {\boldsymbol {g}}_x},{[{\tilde {\boldsymbol {g}}_x},{\tilde {\boldsymbol {g}}_z}]}]$ is spanned by 
$\boldsymbol {e}_y$ and 
$\boldsymbol {e}_z$, with the 
$\boldsymbol {e}_x$ direction seemingly uncontrollable.
In an attempt to yield controllability in this final direction, we consider the 
$\boldsymbol {e}_x$ components of further Lie brackets, with these components having been identically zero for those brackets considered thus far. The lowest-order brackets with a non-zero 
$\boldsymbol {e}_x$ component are 
$[{\tilde {\boldsymbol {g}}_y},{[{\tilde {\boldsymbol {g}}_x},{\tilde {\boldsymbol {g}}_y}]}]$ and 
$[{\tilde {\boldsymbol {g}}_z},{[{\tilde {\boldsymbol {g}}_x},{\tilde {\boldsymbol {g}}_z}]}]$, with this shared component being explicitly given as
\begin{equation} \frac{3}{32}\frac{(2r-3)(8r-9)}{r^5}. \end{equation}
Vanishing when 
$r=3/2$ and 
$r=9/8$, this closely mirrors the structure of (3.18), with these values seemingly being an emergent feature of the far-field approximation. From these third-order Lie brackets, we can therefore conclude that the two-particle tracer-limit system is controllable by the Rashevski–Chow theorem, except when 
$r=3/2$ or 
$r=9/8$. Although these both correspond to the spheres being in close proximity, with the far-field approximation almost certainly being inappropriate here, it is instructive to further consider the former of these two cases, with controllability able to be readily established at 
$r=9/8$ via either of the fourth-order Lie brackets 
$[{\tilde {\boldsymbol {g}}_y},{[{\tilde {\boldsymbol {g}}_x},{[ {\tilde {\boldsymbol {g}}_x},{\tilde {\boldsymbol {g}}_y}}]}]]$ and 
$[{\tilde {\boldsymbol {g}}_x},{[{\tilde {\boldsymbol {g}}_y},{[ {\tilde {\boldsymbol {g}}_x},{\tilde {\boldsymbol {g}}_y}}]}]]$.
In order to best examine the scenario where 
$r=3/2$, we consider the evolution of the squared separation distance 
$\left \lVert \boldsymbol {r}\right \rVert ^2$, as in § 3.2.2. This is simply given by
\begin{equation} \frac{\mathrm{d}\left\lVert \boldsymbol{r}\right\rVert^2}{\mathrm{d}t} = 2\boldsymbol{r}^{\rm T}\frac{\mathrm{d}\boldsymbol{r}}{\mathrm{d}t} = \frac{3-2\left\lVert \boldsymbol{r}\right\rVert}{\left\lVert \boldsymbol{r}\right\rVert}\boldsymbol{r}^{\rm T}\boldsymbol{f}, \end{equation}
vanishing and independent of the control 
$\boldsymbol {f}$ when 
$\left \lVert \boldsymbol {r}\right \rVert =3/2$. Hence, if the tracer is located at this distance from the control sphere, the particle will remain trapped at this separation, irrespective of the applied control, thus is not controllable from this configuration. It should be noted, however, that confounding factors in a physical system, such as the presence of noise or other objects in the flow, may act to move the particle away from such a zero-measure configuration, with the particle rendered controllable.
Therefore, in summary, we have found that the force-driven sphere-tracer system, with the flow due to the sphere being approximated with simply a Stokeslet and only one tracer being considered, is controllable unless the object centres are separated by a distance of 
$3/2$, with the controllable state set formally written as
\begin{equation} S = \{ \boldsymbol{X}\in P : \left\lVert \boldsymbol{r}\right\rVert\neq 3/2 \}. \end{equation}4.2. Including the potential dipole
It is well known that the exact solution for Stokes flow around a translating sphere includes a potential-dipole correction to the Stokeslet term of (4.2) (Kim & Karrila Reference Kim and Karrila2005), which captures finite-volume effects and satisfies the no-slip condition on the surface of the sphere. Hence, noting the difference in controllability found in § 3 when using hydrodynamic representations of differing accuracy, we modify the system of (4.2), again taking 
$N=2$ for simplicity, to include the weighted potential dipole contribution, which we write as
\begin{equation} \frac{\mathrm{d}\boldsymbol{r}}{\mathrm{d}t} = \frac{3}{4}\left(\frac{\boldsymbol{\mathsf{I}}}{\left\lVert \boldsymbol{r}\right\rVert} + \frac{\boldsymbol{r}\boldsymbol{r}^{\rm T}}{\left\lVert \boldsymbol{r}\right\rVert^3}\right) \boldsymbol{f}+ \frac{1}{4}\left(\frac{\boldsymbol{\mathsf{I}}}{\left\lVert \boldsymbol{r}\right\rVert^3} - \frac{3\boldsymbol{r}\boldsymbol{r}^{\rm T}}{\left\lVert \boldsymbol{r}\right\rVert^5}\right)\boldsymbol{f} -\boldsymbol{f}.\end{equation}
This can readily be seen to satisfy 
$\mathrm {d}\boldsymbol {r}/\mathrm {d}t=\boldsymbol {0}$ when 
$\left \lVert \boldsymbol {r}\right \rVert =1$. For completeness, the full control system is now given by
\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t}\begin{bmatrix} \boldsymbol{r} \\ \boldsymbol{x} \end{bmatrix}= \begin{bmatrix} \dfrac{3}{4}\left(\dfrac{\boldsymbol{\mathsf{I}}}{\left\lVert \boldsymbol{r}\right\rVert} + \dfrac{\boldsymbol{r}\boldsymbol{r}^{\rm T}}{\left\lVert \boldsymbol{r}\right\rVert^3}\right) + \dfrac{1}{4}\left(\dfrac{\boldsymbol{\mathsf{I}}}{\left\lVert \boldsymbol{r}\right\rVert^3} - \dfrac{3\boldsymbol{r}\boldsymbol{r}^{\rm T}}{\left\lVert \boldsymbol{r}\right\rVert^5}\right) - \boldsymbol{\mathsf{I}} \\ \boldsymbol{\mathsf{I}} \end{bmatrix}\boldsymbol{f} = \sum_{i\in\{x,y,z\}} f_i \boldsymbol{g}_i. \end{equation}
We now seek to repeat the analysis of § 4.1 with this improved representation of the hydrodynamics and, by doing so, assess the result of finite-volume effects on tracer controllability. By the same reasoning as above, we seek to establish the controllability of the tracer using only Lie bracket terms, and write 
$\tilde {\boldsymbol {g}}_i$ for the 
$\boldsymbol {r}$ components of 
$\boldsymbol {g}_i$. The second-order brackets corresponding to this augmented control system are given by
\begin{equation} [{\tilde{\boldsymbol{g}}_i},{\tilde{\boldsymbol{g}}_j}] = \frac{-3(r-1)^2 (8r^3 + 7r^2+6r+3)}{16r^8}(r_j \boldsymbol{e}_i - r_i \boldsymbol{e}_j). \end{equation}
As with the far-field approximation, these span a space of dimension two, but now without the condition 
$r=9/8$. Indeed, the polynomial factor in (4.10) is non-zero for 
$r>1$, unsurprisingly vanishing when the tracer touches the control sphere at 
$r=1$. When evaluated in the reduced configuration, where 
$\boldsymbol {x}=\boldsymbol {0}$ and 
$\boldsymbol {r}=r\boldsymbol {e}_x$, without loss of generality, these Lie brackets unconditionally span the 
$\boldsymbol {e}_y{}\boldsymbol {e}_z{}$ plane from any point in the admissible state space 
$P\colon= \{\boldsymbol {X}\in \mathbb {R}^6: \left \lVert \boldsymbol {r}\right \rVert >1\}$. Analogously to § 4.1, the final direction, 
$\boldsymbol {e}_x$, is generated by the additional brackets 
$[{\tilde {\boldsymbol {g}}_y},{[{\tilde {\boldsymbol {g}}_x},{\tilde {\boldsymbol {g}}_y}]}]$ and 
$[{\tilde {\boldsymbol {g}}_z},{[{\tilde {\boldsymbol {g}}_x},{\tilde {\boldsymbol {g}}_z}]}]$, with their shared 
$\boldsymbol {e}_x$ entry being
\begin{equation} \frac{3(r - 1)^4 (2 r + 1) (8 r^3 + 7 r^2 + 6 r + 3)}{32 r^{11}}. \end{equation}
However, in contrast to the previous analysis, this component is non-zero for all admissible 
$r$, vanishing at 
$r=1$ and two points inside the control sphere, with the direction generated unconditionally. Hence, having considered a tracer that is advected by the exact velocity field around a translating sphere, though still neglecting any disturbance by the tracer particle, we can now conclude that the sphere-tracer system is, in fact, controllable everywhere, in stark contrast to the partial controllability result obtained when considering only the Stokeslet flow induced by the control sphere. This highlights that some degree of care should be taken in accurately specifying the hydrodynamics that govern the evolution of the control system. Nevertheless, we do note some level of robustness, at least in this case, with the difference in theoretical controllability between the Stokeslet and exact flow representations only being for a single separation, though questions of practical controllability, such as mechanical efficiency, may be differently affected and warrant future exploration.
5. Tracers in Stokeslet flow
5.1. Geometry and hydrodynamics
In order to facilitate progress in the tracer limit above, we simplified the sphere–tracer system to only a single tracer particle. Seeking to consider the controllability of multiple tracer particles in a shared domain, we now make an alternative simplification that allows for ready evaluation of multi-tracer controllability. Further, we employ and exemplify a systematic approach for establishing controllability using multiple controllability matrices, which can be applied in more-general contexts.
In the far-field limit, with the separation of the spheres large compared with their sizes, all interaction terms between the particles vanish, as can be seen explicitly for the two-sphere case in (3.17a–c), where 
$M_\parallel {}$, 
$M_\perp {}$, 
$m_\parallel {}$ and 
$m_\perp {}$ all approach unity in this limit. With the grand mobility tensor of (2.2) thereby rendered diagonal, it is clear that the control of one particle by another is no longer plausible. Taking this large-separation limit and removing any driving forces or torques on the particles, a scenario that corresponds to a dilute suspension of passive tracers, we consider a different modality of control, one in which tracer motion is effected by an externally imposed flow. Motivated by the far-field approximation used in § 4.1 and recent experimental trappings of self-propelling microswimmers (Zou et al. Reference Zou, Zheng, Wu and Lei2020), we investigate the motion of tracers in the flow field produced by a dimensionless Stokeslet of strength 
$\boldsymbol {f}(t)$, though other flows may be similarly considered. In contrast to the exploration of the previous section, here we control the strength of a fixed-position Stokeslet directly, leading to a characteristically different control system to that considered above, which used a moving Stokeslet as the far-field approximation to the flow induced by a moving sphere.
Concretely, subject to the above assumptions and indexing each tracer by 
$k\in \{1,\ldots,N\}$, the evolution of the 
$k$th tracer in the laboratory frame is simply given by
\begin{equation} \frac{\mathrm{d}\boldsymbol{x}_k}{\mathrm{d}t} = \left(\frac{\boldsymbol{\mathsf{I}}}{\left\lVert \boldsymbol{x}_k\right\rVert} + \frac{\boldsymbol{x}_k\boldsymbol{x}_k^{\rm T}}{\left\lVert \boldsymbol{x}_k\right\rVert^3}\right) \boldsymbol{f},\end{equation}
where 
$\boldsymbol {x}_k$ is the swimmer position and the fixed Stokeslet is located at the origin, without loss of generality. Concatenating the 
$\boldsymbol {x}_k$ into 
$\boldsymbol {X}\in \mathbb {R}^{3N}$ as in § 2 and forming composite 
$\boldsymbol {g}_x$, 
$\boldsymbol {g}_y$ and 
$\boldsymbol {g}_z$ in the same manner, these equations of motion combine to yield the 
$N$-tracer driftless control-affine system
\begin{equation} \frac{\mathrm{d}\boldsymbol{X}}{\mathrm{d}t} = \sum_{i\in\{x,y,z\}} f_i(t) \boldsymbol{g}_i(\boldsymbol{X}). \end{equation} We now restrict to the case where 
$N=2$, corresponding to the motion of two inert tracers in Stokeslet flow, and again make use of compact notation for the left-iterated Lie bracket, writing 
$\boldsymbol {g}_{xyz}= [{\boldsymbol {g}_x},{[{\boldsymbol {g}_y},{\boldsymbol {g}_z}]}]$ and similarly for other iterated brackets, with 
$\boldsymbol {g}_{yxxy} = [{\boldsymbol {g}_y},{[{\boldsymbol {g}_x},{[ {\boldsymbol {g}_x},{\boldsymbol {g}_y}}]}]]$, for example. With 
$i,j\in \{x,y,z\}$, the second-order brackets are explicitly given by
\begin{equation} \tilde{\boldsymbol{g}}_{ij,k} = \frac{3}{\left\lVert \boldsymbol{x}_k\right\rVert^2}[(\boldsymbol{x}_k\boldsymbol{\cdot}\boldsymbol{e}_j)\boldsymbol{e}_i - (\boldsymbol{x}_k\boldsymbol{\cdot}\boldsymbol{e}_i)\boldsymbol{e}_j], \end{equation}
where, for 
$k\in \{1,2\}$, 
$\tilde {\boldsymbol {g}}_{\ast {},k}$ denotes the 
$\boldsymbol {x}_k$ components of 
$\boldsymbol {g}_{\ast {}}$, analogous to the notation of § 4.
As in § 3, we compute the determinants of controllability matrices, seeking to identify the points of the six-dimensional state space 
$P$ at which these matrices are of full rank, with 
$P=\{[\boldsymbol {x}_1;\boldsymbol {x}_2]\in \mathbb {R}^6 : \boldsymbol {x}_1,\boldsymbol {x}_2\neq \boldsymbol {0}\}$. However, the evaluation of these expressions would prove to be cumbersome should we retain the six-dimensional representation of the control system. Previously, noting the invariance of controllability to rotations and translations in these systems, we successfully reduced the state space to facilitate simple analysis, which remains valid in this case. Indeed, we may now simplify even further, with this tracer system being self-similar with respect to rescalings in space, compensated for by rescalings in time. Thus, we rotate and rescale the spatial description of the two-tracer problem, taking 
$\boldsymbol {x}_1=\boldsymbol {e}_x$ and 
$\boldsymbol {x}_2=x\boldsymbol {e}_x + y\boldsymbol {e}_y$ and omitting the subscripts on the components of 
$\boldsymbol {x}_2$ for notational convenience. We also remark that we may take 
$y\geq 0$, without loss of generality. This significantly reduced configuration is illustrated in figure 6.
Figure 6. A reduced two-tracer state space for evaluating controllability. The original six-dimensional state space (a), with 
$\boldsymbol {x}_1$ and 
$\boldsymbol {x}_2$ being general points in 
$\mathbb {R}^3$ and the Stokeslet situated at the origin, labelled 
$\boldsymbol {f}$, is instantaneously equivalent to the reduced planar configuration (b) up to rotation and rescaling, noting the invariance of controllability properties to such transformations in the context of the Stokeslet-driven control system of (5.2). Of note, the reduced system still admits unrestricted motion of each tracer in 
$\mathbb {R}^3$, with the instantaneous controllability properties merely being evaluated in this notationally reduced but nevertheless general configuration.
Of particular note, although we later evaluate controllability matrices at points in this reduced two-dimensional space, with elements denoted by 
$\boldsymbol {X}^{\ast }{}$, motion out of the plane of this instantaneous configuration is still permitted. Indeed, 
$\boldsymbol {f}$ need not lie within this plane. Hence, we duly seek to identify controllability matrices of rank six, despite evaluating their entries in this notationally and computationally convenient two-dimensional setting. We refer to this reduced space as 
$\tilde {P}=\{(x,y)\in \mathbb {R}^2: y\geq 0,\ (x,y)\neq (0,0)\}$, noting that controllability on 
$\tilde {P}$ is equivalent to controllability on 
$P$.
5.2. Establishing controllability
Seeking a controllability matrix 
$\boldsymbol{\mathsf{C}}$ of rank six, the natural and perhaps minimal choice of 
$\boldsymbol{\mathsf{C}}=[\boldsymbol {g}_x,\boldsymbol {g}_y,\boldsymbol {g}_z,\boldsymbol {g}_{xy},\boldsymbol {g}_{xz},\boldsymbol {g}_{yz}]$ leads to 
$\det {\boldsymbol{\mathsf{C}}}=0$ identically, with higher-order Lie brackets therefore warranted in order to yield controllability. In what follows, we consider the three controllability matrices
\begin{gather} \boldsymbol{\mathsf{C}}_1 = [\boldsymbol{g}_x,\boldsymbol{g}_y,\boldsymbol{g}_z,\boldsymbol{g}_{xy},\boldsymbol{g}_{yz},\boldsymbol{g}_{xxy}], \end{gather}
\begin{gather}\boldsymbol{\mathsf{C}}_2 = [\boldsymbol{g}_x,\boldsymbol{g}_y,\boldsymbol{g}_z,\boldsymbol{g}_{xy},\boldsymbol{g}_{xz},\boldsymbol{g}_{yxy}], \end{gather}
\begin{gather}\boldsymbol{\mathsf{C}}_3 = [\boldsymbol{g}_x,\boldsymbol{g}_y,\boldsymbol{g}_z,\boldsymbol{g}_{xy},\boldsymbol{g}_{xz},\boldsymbol{g}_{yxxy}], \end{gather}
which differ only in the final two columns. For completeness, the full expressions for 
$\boldsymbol{\mathsf{C}}_1$, 
$\boldsymbol{\mathsf{C}}_2$, and 
$\boldsymbol{\mathsf{C}}_3$ are given in Appendix B along with their determinants. As also further explained in Appendix B, we now identify the controllable sets 
$S_1,S_2,S_3\subseteq \tilde {P}$ where each of 
$\boldsymbol{\mathsf{C}}_1$, 
$\boldsymbol{\mathsf{C}}_2$, and 
$\boldsymbol{\mathsf{C}}_3$ have non-vanishing determinant, respectively, noting that the system is controllable at a point 
$\boldsymbol {p}\in \tilde {P}$ in the reduced configuration space if 
$\boldsymbol {p}\in S_1\cup S_2\cup S_3$ by the Rashevski–Chow theorem.
The set 
$Z_1\colon= \tilde {P}\setminus S_1$, on which 
$\det {\boldsymbol{\mathsf{C}}_1}$ vanishes, is found to be a subset of the upper right quadrant and the x axis. The analogously defined 
$Z_2$ intersects with 
$Z_1$ at precisely two points 
$(x,y)\in \{(-1,0),(1,0)\}$, with the origin excluded from 
$\tilde {P}$. These sets and their points of intersection are shown in figure 7, with controllability established everywhere except at these intersections. Notably, the process of identifying potentially uncontrollable states may be streamlined once 
$Z_1$ is found, with only the intersection 
$Z_1\cap Z_2$, rather than the whole of 
$Z_2$, being needed to further interrogate controllability, as demonstrated explicitly in Appendix B.
Figure 7. Sets 
$Z_1$ and 
$Z_2$ on which 
$\det {\boldsymbol{\mathsf{C}}_1}$ and 
$\det {\boldsymbol{\mathsf{C}}_2}$ in turn vanish, shown as black and grey curves, respectively. Their intersections are marked as black points, corresponding to three cases: (i) 
$\boldsymbol {x}_1=-\boldsymbol {x}_2$, the tracers are mirror images in the Stokeslet; (ii) 
$\boldsymbol {x}_2=\boldsymbol {0}$, the tracer is at the singular point of the flow, excluded from 
$P$ (shown open); (iii) 
$\boldsymbol {x}_1=\boldsymbol {x}_2$, the tracers share location and the system is degenerate.
The two points of intersection each correspond to a simple configuration. For 
$(x,y)=(1,0)$ the tracers overlap and share position for all time, so that 
$\boldsymbol {x}_1=\boldsymbol {x}_2$, with independent control of both 
$\boldsymbol {x}_1$ and 
$\boldsymbol {x}_2$ clearly impossible. The second case, however, entails only that 
$\boldsymbol {x}_1=-\boldsymbol {x}_2$, so that the two tracers are mirror images with respect to the location of the Stokeslet, a configuration that does not obviously prohibit controllability. Indeed, consideration of the matrix 
$\boldsymbol{\mathsf{C}}_3$ is sufficient to provide controllability in this circumstance, which may be easily established by noting that 
$\det {\boldsymbol{\mathsf{C}}_3}$ is non-zero at 
$(x,y)=(-1,0)$. Thus, 
$(-1,0)\notin Z_3$, so that 
$Z_1\cap Z_2 \cap Z_3 = \{(1,0)\}$. Therefore, 
$S_1\cup S_2\cup S_3 = \tilde {P}\setminus \{(1,0)\}$ and we have demonstrated controllability on all of the reduced space, excluding the degenerate point where the tracers overlap. Though controllable, the requirement of a fourth-order bracket when 
$\boldsymbol {x}_1=-\boldsymbol {x}_2$ entails that control is somewhat more challenging in this configuration, with higher-order brackets needing more-intricate controls to realise in practice.
Recalling that the three considered controllability matrices differed only in their final two columns, a succinct summary of the above calculations is that the composite controllability matrix
\begin{equation} \boldsymbol{\mathsf{C}} =[\boldsymbol{g}_x,\boldsymbol{g}_y,\boldsymbol{g}_z,\boldsymbol{g}_{xy},\boldsymbol{g}_{yz},\boldsymbol{g}_{xz},\boldsymbol{g}_{xxy},\boldsymbol{g}_{yxy},\boldsymbol{g}_{yxxy}] \end{equation}
is of rank six everywhere in 
$\tilde {P}$ apart from when the tracers are coincident. Therefore, the two-tracer system is controllable from any such configuration.
6. Constructing control schemes
Controllability analysis theoretically predicts that driving a system from any state 
$\boldsymbol {X}_0$ to any target state 
$\boldsymbol {X}_1$ is possible within the controllable state space, though, as yet, it may be unclear how such a transition can be achieved. In fact, exploiting the directions given by the Lie brackets found to span in the previous sections, it is possible to explicitly design control policies that realise this goal, as described for the first time in the celebrated work of Lafferriere & Sussmann (Reference Lafferriere and Sussmann1992), which contributed a general motion planning algorithm for non-holonomic driftless control-affine systems. Since then, additional schemes for motion planning have been devised, such as the so-called sinusoidal methodology of Murray & Sastry (Reference Murray and Sastry1993) and the approach of Melli, Rowley & Rufat (Reference Melli, Rowley and Rufat2006), the latter being presented in the context of articulated bodies. In this section, we briefly describe and numerically exemplify the algorithm of Lafferriere and Sussmann as applied to the two-sphere system of § 3, inheriting its motivation and context. For the sake of legibility, we reduce the dimension of the system and focus only on planar motion; carrying out the algorithm in higher dimensions is strictly analogous.
Let 
$X_0$ be the initial state and 
$X_1$ a target state for the two-sphere system. We have proved in § 3.2.1 that force control of one sphere allows controllability of the system, owing to the controllability matrix composed of four fields, 
$\boldsymbol{\mathsf{C}} = ( \boldsymbol {g}_x, \boldsymbol {g}_y, [{\boldsymbol {g}_x},{\boldsymbol {g}_y}],[{\boldsymbol {g}_y},{[{\boldsymbol {g}_x},{\boldsymbol {g}_y}]}])$, having full rank everywhere in the state space. The principal idea of the motion planning scheme is to utilise these fields and consider the extended control system
\begin{equation} \dot{\boldsymbol{X}} = v_1 \boldsymbol{g}_x + v_2 \boldsymbol{g}_y + v_3 [{\boldsymbol{g}_x},{\boldsymbol{g}_y}] + v_4 [{\boldsymbol{g}_y},{[{\boldsymbol{g}_x},{\boldsymbol{g}_y}]}] = \boldsymbol{\mathsf{C}} \boldsymbol{v}, \end{equation}
with a new set of four controls 
$\boldsymbol {v} = (v_1,v_2,v_3,v_4)$. Now having as many controls as equations, we can build a control 
$\boldsymbol {v}$ that drives the system along a set trajectory 
$\boldsymbol {\bar {X}}(t)$ that links 
$\boldsymbol {X}_0$ to 
$\boldsymbol {X}_1$ on a set time interval 
$[0,T]$ by simply inverting system (6.1):
\begin{equation} \boldsymbol{v}(t) = \boldsymbol{\mathsf{C}}^{{-}1} \dot{\bar{\boldsymbol{X}}}(t). \end{equation}
We typically choose 
$\boldsymbol {\bar {X}}$ as a straight line between 
$\boldsymbol {X}_0$ and 
$\boldsymbol {X}_1$, in which case 
$\boldsymbol {\bar {X}}(t) = t \boldsymbol {X}_1 + (T-t) \boldsymbol {X}_0$ and 
$\boldsymbol {v}(t) = \boldsymbol {v} = \boldsymbol{\mathsf{C}}^{-1} (\boldsymbol {X}_1 - \boldsymbol {X}_0)$ is constant. Note that, for other control systems, such as that of § 5, the overall controllability matrix 
$\boldsymbol{\mathsf{C}}$ may not be square, in which case 
$\boldsymbol {v}$ may be defined with a pseudoinverse of 
$\boldsymbol{\mathsf{C}}$ instead of 
$\boldsymbol{\mathsf{C}}^{-1}$.
By construction, with 
$\boldsymbol {v}$ defined as in (6.2), we have 
$\boldsymbol {X}(T) = \boldsymbol {X}_1$. Next, we wish to use this control 
$\boldsymbol {v}$ to build controls 
$(f_x,f_y)$ that will drive the original two-sphere system to the same target. This can be done by following an elegant, albeit technical, process of calculating so-called Philip Hall coordinates of the solution 
$\boldsymbol {X}$, which requires formal calculation with non-commutative indeterminates 
$X$ and 
$Y$ that represent the fields 
$\boldsymbol {g}_x$ and 
$\boldsymbol {g}_y$. For more detail on this formalism, we refer the reader to Lafferriere & Sussmann (Reference Lafferriere and Sussmann1992) or the more comprehensive introduction to geometric control theory of Kawski & Sussmann (Reference Kawski and Sussmann1997). A Philip Hall basis of the polynomials in 
$(X,Y)$ up to order three is defined as 
$(X,Y,[X,Y],[Y,[X,Y]],[X,[X,Y]])$, with the first four elements corresponding to the brackets used to define the controllability matrix 
$\boldsymbol{\mathsf{C}}$. Then, one can associate the Philip Hall coordinates 
$(h_1,h_2,h_3,h_4,h_5)$ to the formal trajectory 
$\tilde {\boldsymbol {X}}$ in such a way that
\begin{equation} \tilde{\boldsymbol{X}}(t) = {\rm e}^{h_1(t) X}\,{\rm e}^{h_2(t) Y}\,{\rm e}^{h_3(t) [X,Y]} \,{\rm e}^{h_4(t) [Y,[X,Y]]}\,{\rm e}^{h_5(t) [X,[X,Y]]}, \end{equation}
decomposing the trajectory with respect to the directions given by each of the (formal) brackets in the Philip Hall basis. Moreover, via the well-known Baker–Campbell–Hausdorff formula, these coordinates 
$(h_1,h_2,h_3,h_4,h_5)$ may be easily computed as algebraic expressions of the solution to a simple ordinary differential equation (ODE) system, which is, in turn, dependent on 
$\boldsymbol {v}$:
\begin{equation} \left\{\begin{array}{@{}l} h_1 = g_1, \\ h_2 = g_2, \\ h_3 = g_3 - g_1 g_2, \\ h_4 = g_4 - g_2 g_3 + \frac{1}{2}g_1 g_2^2 - \frac{1}{6} g_1^2 g_2, \\ h_5 = g_4 - g_1 g_3 + \frac{1}{2}g_1^2 g_2 + \frac{1}{6} g_1 g_2^2, \end{array}\right.\quad \text{where} \left\{\begin{array}{@{}l} \dot{g}_1 = v_1, \\ \dot{g}_2 = v_2, \\ \dot{g}_3 = v_3 + g_1 v_2, \\ \dot{g}_4 = v_4 + g_1 g_2 v_2, \\ \dot{g}_5 = \frac{1}{2} g_1^2 v_2 + g_1 v_3,\\ g_i(0) = 0,\quad i = 1\ldots 5. \end{array}\right. \end{equation}
Finally, the controls 
$f_x$ and 
$f_y$ can be used to intuitively generate each term in (6.3) at 
$t=T$. For example, the first term 
${\rm e}^{h_1(T) X}$ can be generated by setting 
$f_x (t) = h_1(T)$ for 
$t\in [0,1]$. Terms associated with higher-order brackets are obtained through concatenations of 
$f_x$ and 
$f_y$, as demonstrated in (2.8) of § 2. As an explicit example in this context, the second-order bracket direction 
${\rm e}^{h_3(T) [X,Y]}$ can be achieved on 
$[0,4]$ with the control policy
\begin{equation} \boldsymbol{f}(t)=\left\{\begin{array}{@{}ll} f_x=\zeta, & \text{for } t\in[0,1), \\ f_y=\zeta, & \text{for } t\in[1,2), \\ f_x={-}\zeta, & \text{for } t\in[2,3), \\ f_y={-}\zeta, & \text{for } t\in[3,4], \end{array}\right. \end{equation}
with 
$\zeta = h_3 (T)^{1/2}$ if 
$h_3(T)\geq 0$; if 
$h_3(T) < 0$, we simply interchange the roles of 
$f_x$ and 
$f_y$ in this description. In this example, the controls are specified as piecewise constant functions, but one may instead utilise smoother controls, as further remarked upon in Appendix C. Repeating this process for the brackets of third order and performing some minor simplifications, we can define a control policy made of 23 steps that drives 
$\boldsymbol {X}_0$ to 
$\boldsymbol {X}_1$, which we give explicitly in Appendix C. However, it should be noted that this control policy only reaches the approximate location of the target, as we have neglected the effects of brackets of order four and higher in the above calculations. To reduce the effects of these higher-order terms, one can subdivide the trajectory into smaller steps by defining intermediate targets, a process with guaranteed convergence (Lafferriere & Sussmann Reference Lafferriere and Sussmann1992). In more detail, given a step size of 
$\varDelta$ and a final error tolerance of 
$\delta$, the algorithm can be summarised as
We numerically apply this algorithm to control our two-sphere system using the MATLAB ODE solver ode45 (Shampine & Reichelt Reference Shampine and Reichelt1997), with the results displayed in figure 8 and Movies 1 and 2 of the Supplementary Material available at https://doi.org/10.1017/jfm.2022.253. Figures 8(a)–8(c) and Movie 1 of the Supplementary Material show the character of the trajectory produced by the algorithm on a simple example, with step size 
$\varDelta = 0.5$. The trajectory of each sphere after one iteration is displayed in figure 8(b), from which we note that generating motion in the target direction requires moving the sphere around a ‘loop’ whose length is significantly greater than the final displacement. This serves to validate the study of efficiency in § 3.3, which suggested that the third-order brackets require high-energy input due the factor 
$\varepsilon ^2$ appearing in (3.35). With this value of 
$\varDelta$, the algorithm converges after seven iterations, with the results shown in figure 8(c). The size of the ‘loop’ can be reduced by decreasing the value of 
$\varDelta$, which naturally entails that the target will be reached after more iterations. This is shown in figure 8(e) and Movie 2 of the Supplementary Material, in which the same target is reached with 
$\varDelta = 0.1$ after 
$27$ iterations. Finally, this algorithm may be applied repeatedly, allowing control of the spheres over more-complex trajectories, as illustrated in figure 8(f).
Figure 8. Application of the motion planning algorithm to the finite-size two-sphere control problem. The initial configuration is shown in (a), with the control sphere in red centred at 
$(0,1)$ and the passive sphere in green centred at 
$(3,3)$, both having a radius of unity. The target positions of the spheres, 
$(0,0)$ and 
$(1,3)$, are represented on the graph by black and blue crosses, respectively. (b) The position of the spheres and their trajectory after one iteration of the motion planning algorithm, whereas (c) shows the full trajectory that results in the spheres reaching their targets after seven iterations. (d) The controls 
$f_x$ and 
$f_y$ to be applied to generate the trajectories on (b). (e) The trajectory after one iteration for a smaller value of 
$\varDelta$ (
$\varDelta = 0.1$), reducing the size of the loops made by the spheres. The black and blue crosses indicate the positions of the spheres after each iteration of the algorithm until they reach the target positions. (f) A more elaborate trajectory where the control and passive spheres are assigned to reach successive targets, along a circle and a square, respectively, each shown as a red line. The full trajectories of the spheres are represented as light grey lines. The respective initial positions of the control and passive sphere are 
$(-0.5,1.5)$ and 
$(-1,-1)$ and they traverse their target trajectories anticlockwise and clockwise, respectively. In addition to this figure, Movies 1 and 2 of the Supplementary Material display the trajectories presented on (c) and (f), respectively.
Thus, via this concrete example, we have seen how a controllability analysis may be leveraged practically to design and explicitly construct motion planning policies to realise a control task, with the presented approach also being immediately generalisable. This basic description may be readily augmented and refined, with smoother, more-direct trajectories being obtained via a combination of adaptive step sizing and the considered design of control combinations when seeking to generate the Lie bracket directions, of likely pertinence to potential applications and experimental set-ups. This exploration has also served to validate the controllability analysis of this two-sphere control system, with control able to be effected using the controllability matrix of § 3.2.1.
7. Discussion
In this work, we have presented an introduction to control and controllability in the context of multiple particles in Stokes flow, illustrated by a number of explicit examples of quite general control problems involving hydrodynamic particle–particle interactions, with this framework contrasting with other explorations of controllability in this regime (Shapere & Wilczek Reference Shapere and Wilczek1989b; Martín et al. Reference Martín, Takahashi and Tucsnak2007; Desimone et al. Reference Desimone, Heltai, Alouges and Lefebvre-Lepot2012; Alouges et al. Reference Alouges, DeSimone, Giraldi and Zoppello2013a; Lohéac & Munnier Reference Lohéac and Munnier2014; Hatton & Choset Reference Hatton and Choset2015; Lohéac & Takahashi Reference Lohéac and Takahashi2020). In doing so, we have seen how rigorous controllability can be established with elementary, if cumbersome, calculation, highlighting the potential for the systematic justification of controllability in other studies, which is often lacking in contemporary works. This broad applicability draws from the time stationarity of many Stokes flow problems, with explicit reference to time being absent from the governing equations of motion. As we have seen in the cases presented in this study, this can give rise to systems in control-affine form, with simple yet powerful results of geometric control theory then being readily wielded to interrogate the controllability. Absent from the problems considered here, the presence of drift terms, those with no dependence on any control, such as propulsive velocities or prescribed background flows, somewhat complicate the analysis of controllability, though similar tools are able to query the property of STLC, leading to results pertaining to local rather than global control. Such systems are of particular relevance to self-propelled microswimmers, with an analysis of such a swimmer in a controlled background flow being presented by Moreau et al. (Reference Moreau, Ishimoto, Gaffney and Walker2021).
Though illustrative of a general approach to controllability in the Stokes regime, the particular examples explored in this study may be interpreted in physical or biophysical contexts. For instance, the force-driven controls that we have studied can be interpreted as being effected in practice by optical tweezers, a methodology that may also be used to generate the Stokeslet flow considered in § 5 (Zou et al. Reference Zou, Zheng, Wu and Lei2020). Facilitating many of the control modalities explored in this work, hydrodynamic interactions between particles have also been studied in various biophysical contexts, such as the enhancement of hydrodynamic capturing of swimming sperm near an egg (Ishimoto, Cosson & Gaffney Reference Ishimoto, Cosson and Gaffney2016), the colonisation of bacteria near a nutrient source (Desai, Shaik & Ardekani Reference Desai, Shaik and Ardekani2019) and bacterial predator–prey dynamics (Ishimoto, Gaffney & Walker Reference Ishimoto, Gaffney and Walker2020). Further, flow-mediated particle control has been discussed in the context of nutrition uptake in microorganisms and ciliated cells (Riisgård & Larsen Reference Riisgård and Larsen2010; Kiørboe et al. Reference Kiørboe, Jiang, Gonçalves, Nielsen and Wadhwa2014). Hydrodynamic interactions between particles are also utilised as a mechanism of cargo transport in an active colloid system (Demirörs et al. Reference Demirörs, Akan, Poloni and Studart2018; Wang & Simmchen Reference Wang and Simmchen2019; Yang & Bevan Reference Yang and Bevan2020; Zou et al. Reference Zou, Zheng, Wu and Lei2020), with the controllability results established in the present study rigorously guaranteeing that this particle transportation is possible. However, this analysis has also indicated the relative inefficiency of particle repulsion and attraction via hydrodynamic interactions when objects are in very close proximity.
As a consequence of the explicit manner in which we have established control throughout this study, identifying spanning elements of the problem-specific Lie algebra, we have been able to constructively and systematically design schemes for control towards a particular goal, employing the approach of Lafferriere & Sussmann (Reference Lafferriere and Sussmann1992). Indeed, by utilising the spanning sets found in our analysis, this methodology constructs an appropriate control, with the resulting control policy demonstrating the key role that non-commutativity can play in realising the manipulation of high-dimensional systems with low-dimensional control modalities. This inherently constructive approach is starkly different to many modern methods for identifying suitable controls, with popular and powerful techniques such as machine learning and high-dimensional nonlinear optimisation being in common use. In contrast to these ‘black-box’ methods, the highlighted explicit approach of Lafferriere & Sussmann (Reference Lafferriere and Sussmann1992) has the potential to readily develop our understanding of the physics of the control systems considered, with the aforementioned importance of non-commutative controls being a significant example of this.
Bordering on the related field of optimal control, we have also explored the mechanical efficiency of aspects of control, reporting on the efficiency-optimal distances at which particular controls should be applied. The dependence of these optimal distances on the sphere-size ratio 
$\lambda$ suggests further investigation in the noted biophysical interpretations of the two-sphere problem, such as the efficient capture of prey and the realisation of sperm–egg encounters (Jabbarzadeh & Fu Reference Jabbarzadeh and Fu2018). There also remains significant scope for the consideration of optimality in further generality, with possible extensions of the explicit method of Lafferriere & Sussmann (Reference Lafferriere and Sussmann1992) being a likely topic of further investigation. In particular, informed by the analysis of the force-driven sphere system in § 3.3, we expect that optimal controls for maximising efficiency would make minimal use of terms that correspond to higher-order Lie brackets, with these typically requiring significant additional control input in order to realise comparable changes in state to lower-order brackets. In a practical context, alternative definitions of efficiency, or indeed optimality, may be appropriate, with the investigation of these settings being an additional direction for exploration.
A further theoretical limitation present in this study lies in the geometric nature of the steady Stokes equations. If the oscillatory frequency of the control force were to be large, the commonplace assumption of negligible inertial effects may no longer be valid. For such a rapidly oscillating force actuation, one might seek to apply the unsteady Stokes equations and the associated Stokeslet, called the oscillet, which leads to a correction in the far-field (Wei et al. Reference Wei, Dehnavi, Aubin-Tam and Tam2021). In this work, we have exploited symmetry to improve the computational tractability of the systems that we have considered, a property not present in many-sphere settings or in more-complex geometries. Indeed, with confinement having been noted to enhance the controllability of a microrobot (Alouges & Giraldi Reference Alouges and Giraldi2013), assessing the control properties of torque-driven spheres near a boundary may be of great interest. With the boundary modifying the mobility coefficients of the spheres, one can imagine that they may be sufficiently modified to enhance controllability, which we have seen is otherwise limited in an unbounded domain. The extension to many-sphere systems, be they finite-size or tracer particles, represents a likely direction for future work; indeed, the principle of the analysis presented in this work is readily applicable to higher-dimensional systems, with the caveat of increased algebraic and notational complexity. However, if such multi-sphere systems are driven with relatively few degrees of freedom, as we have done in this study, their control will almost certainly require the use of high-order Lie brackets, with practical control perhaps limited by this likely necessity.
In summary, we have explored the controllability of multiple physically inspired problems in the Stokes regime, highlighting how explicit assessments of controllability can both facilitate the constructive design of controls and deepen our understanding of underlying fluid mechanics. In examining these example cases, which has also led to an evaluation of mechanical efficiency and its optimisation, we have summarised and demonstrated the application of key principles of geometric control theory such that they may be readily translated to other research contexts, overall serving to illustrate how these simple yet powerful tools can be utilised to rigorously evaluate control in the Stokes limit.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2022.253.
Funding
B.J.W. is supported by the UK Engineering and Physical Sciences Research Council (EPSRC), grant EP/N509711/1. K.I. is supported by JSPS-KAKENHI for Young Researchers (18K13456) and JST, PRESTO grant number JPMJPR1921. C.M. is a JSPS International Research Fellow (PE20021).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Mobility coefficients
The mobility coefficients, 
$M_\parallel {}$, 
$N_\parallel {}$, 
$m_\parallel {}$, 
$M_\perp {}$, 
$N_\perp {}$, 
$m_\perp {}$, 
$M_T$, 
$N_T$ and 
$m_T$ of § 3 are obtained by the inversion of a corresponding resistance matrix. This matrix is constructed from the widely used expressions of Jeffrey & Onishi (Reference Jeffrey and Onishi1984), in which they connect the far-field expansion of the hydrodynamics, up to 
$O(1/r^{11})$, with results from lubrication theory. These intricate expressions are evaluated numerically, with select cases plotted in figure 9. Of note, for 
$\lambda \ll 1$, the series expression for 
$M_T$ does not converge sufficiently fast to enable us to conclude that it is non-negative, even if the series is summed up to 
$O(1/r^{300})$.
Figure 9. Mobility coefficients for various relative sphere sizes, as computed from the expressions of Jeffrey & Onishi (Reference Jeffrey and Onishi1984). For 
$\lambda =0.1$, 
$M_T$ is not distinguishable from zero at the resolution of this plot.
Appendix B. Explicit analysis of two tracers in flow
With reference to the set-up of § 5.2, the controllability matrices 
$\boldsymbol{\mathsf{C}}_1$, 
$\boldsymbol{\mathsf{C}}_2$ and 
$\boldsymbol{\mathsf{C}}_3$ are given explicitly by
\begin{gather} \boldsymbol{\mathsf{C}}_1 = \begin{bmatrix} 2 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & -3 & 0 & 21 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ \dfrac{x^2}{r^{3}}+\dfrac{1}{r} & \dfrac{x y}{r^{3}} & 0 & \dfrac{3 y}{r^{4}} & 0 & -\dfrac{27 x y}{r^{7}} \\ \dfrac{x y}{r^{3}} & \dfrac{y^2}{r^{3}}+\dfrac{1}{r} & 0 & -\dfrac{3 x}{r^{4}} & 0 & \dfrac{3 \left(7 x^2-2 y^2\right)}{r^{7}} \\ 0 & 0 & \dfrac{1}{r} & 0 & -\dfrac{3 y}{r^{4}} & 0 \end{bmatrix}, \end{gather}
\begin{gather}\boldsymbol{\mathsf{C}}_2 = \begin{bmatrix} 2 & 0 & 0 & 0 & 0 & 6 \\ 0 & 1 & 0 & -3 & 0 & 0 \\ 0 & 0 & 1 & 0 & -3 & 0 \\ \dfrac{x^2}{r^{3}}+\dfrac{1}{r} & \dfrac{x y}{r^{3}} & 0 & \dfrac{3 y}{r^{4}} & 0 & \dfrac{3 \left(2 x^2-7 y^2\right)}{r^{7}} \\ \dfrac{x y}{r^{3}} & \dfrac{y^2}{r^{3}}+\dfrac{1}{r} & 0 & -\dfrac{3 x}{r^{4}} & 0 & \dfrac{27 x y}{r^{7}} \\ 0 & 0 & \dfrac{1}{r} & 0 & -\dfrac{3 x}{r^{4}} & 0 \end{bmatrix}, \end{gather}
\begin{gather} \boldsymbol{\mathsf{C}}_3 = \begin{bmatrix} 2 & 0 & 0 & 0 & 0 & -21 \\ 0 & 1 & 0 & -3 & 0 & 0 \\ 0 & 0 & 1 & 0 & -3 & 0 \\ \dfrac{x^2}{r^{3}}+\dfrac{1}{r} & \dfrac{x y}{r^{3}} & 0 & \dfrac{3 y}{r^{4}} & 0 & \dfrac{27 x y^2 \left(y^2-2 x^2\right)+6 x y^2 \left(x^2-2 y^2\right)-21 x^3 \left(x^2-2 y^2\right)}{r^{12}} \\ \dfrac{x y}{r^{3}} & \dfrac{y^2}{r^{3}}+\dfrac{1}{r} & 0 & -\dfrac{3 x}{r^{4}} & 0 & -\dfrac{3 \left(9 x^2 y \left(x^2+4 y^2\right)+7 x^2 y \left(x^2-2 y^2\right)+2 y^3 \left(2 y^2-x^2\right)\right)}{r^{12}} \\ 0 & 0 & \dfrac{1}{r} & 0 & -\dfrac{3 x}{r^{4}} & 0 \end{bmatrix}, \end{gather}
where 
$r=\sqrt {x^2+y^2}$ and these matrices are evaluated at 
$\boldsymbol {X}^{\ast }{}$ in the reduced two-tracer space 
$\tilde {P}$. The determinants may be computed as
\begin{gather} \det{\boldsymbol{\mathsf{C}}_1} ={-}\frac{108 y^2 \left(7r^4 -8 xr + 1\right)}{r^{13}}, \end{gather}
\begin{gather}\det{\boldsymbol{\mathsf{C}}_2} ={-}\frac{108 \left(r^{10} -2 x r^{7} +7 y^2 r^4 + 2 x r^{3}-8 x y^2 r-x^2\right)}{r^{13}}\end{gather}and
\begin{align} \det{\boldsymbol{\mathsf{C}}_3} &= 54r^{{-}16} (7 r^{13} -14 x r^{10} +7 r^{9} -7 y^2 r^{7} -7 x r^{6} +14 r^{5} \nonumber\\ &\quad +21 x y^2 r^{4} -31 y^2 r^{3} -7 x r^{2}+21 y^4 r+3 x y^2). \end{align} We next identify 
$Z_1\colon= \tilde {P}\setminus S_1$, the set on which 
$\det {\boldsymbol{\mathsf{C}}_1}$ vanishes. This trivially includes all points with 
$y=0$, with the remaining elements 
$\boldsymbol {X}^{\ast }{}\in \tilde {P}$ necessarily satisfying
\begin{equation} 8xr = 7r^4+1. \end{equation}
This latter relation admits solutions only when 
$x>0$ and, thus, defines a curve in the upper right quadrant of the 
$xy$ plane. With controllability thereby already broadly established, though still lacking on 
$Z_1$, we now consider the zeros of 
$\det {\boldsymbol{\mathsf{C}}_2}$. Notably, for the purpose of establishing controllability, we may restrict our rootfinding to 
$Z_1$, computing the intersection 
$Z_1\cap Z_2$, which significantly simplifies the necessary calculation. First seeking solutions with 
$y=0$, the condition 
$\det {\boldsymbol{\mathsf{C}}_2}=0$ reduces to
\begin{equation} x^{10} - 2x\left\lvert{x}\right\rvert^7 + 2x\left\lvert{x}\right\rvert^3 - x^2 = 0, \end{equation}
which holds precisely when 
$x\in \{-1,0,1\}$. Hence, 
$\{(-1,0),(1,0)\}\subseteq Z_1\cap Z_2$, recalling that 
$\{(0,0)\}\not \in \tilde {P}$. The remaining elements of 
$Z_1\cap Z_2$ are solutions of the simultaneous equations
\begin{gather} 0 = 7r^4+1 - 8xr, \end{gather}
\begin{gather}0 = r^{10} -2 x r^{7} +7 y^2 r^4 + 2 x r^{3}-8 x y^2 r-x^2. \end{gather}
Recalling that the first equation admits solutions only for 
$x>0$, this system can be solved analytically via elementary manipulations, leading to the single solution 
$(x,y) = (1,0)$. Therefore, we in fact have that 
$Z_1\cap Z_2 = \{(-1,0),(1,0)\}$ precisely, with controllability of the two-tracer system thereby established except at these three points of the reduced space. Both 
$Z_1$ and 
$Z_2$ are illustrated in figure 7, computed numerically and with their three points of intersection shown as black dots. The controllability analysis is now readily completed by evaluating 
$\det {\boldsymbol{\mathsf{C}}_3}$ on 
$Z_1\cap Z_2$, which vanishes only at 
$(1,0)$. Hence, 
$S_1\cup S_2\cup S_3 = \tilde {P}\setminus \{(1,0)\}$.
Appendix C. Equations for constructing control schemes
The piecewise-constant control policies built from the Philip Hall coordinates 
$(h_1,h_2,h_3,h_4,h_5)$ are defined on 
$t=[0,24]$ as follows. We define
\begin{equation} \left.\begin{gathered} \alpha = h_1(T),\quad \beta = h_2(T),\quad \gamma = \sqrt{| h_3(T) |},\\ \delta = \left ( h_5(T) - \tfrac{1}{2} \gamma^3 \right )^{{1}/{3}},\quad \epsilon = \left ( h_4(T) + \tfrac{1}{2} \gamma^3 \right)^{{1}/{3}}. \end{gathered}\right\} \end{equation}
The terms 
$\frac {1}{2} \gamma ^3$ in the definition of 
$\delta$ and 
$\epsilon$ are used to remove the third-order residuals generated by the control policy generating the direction 
$[\boldsymbol {g}_x,\boldsymbol {g}_y]$. Then, if 
$h_3(T) \geqslant 0$, we define
\begin{gather} f_x(t) = \left\{\begin{array}{@{}ll} \alpha & \text{on}\ [0,1), \\ \gamma & \text{on}\ [2,3), \\ -\gamma & \text{on}\ [4,5), \\ \delta & \text{on}\ [6,8) \cup [13,14), \\ -\delta & \text{on}\ [9,10) \cup [11,12) \cup [15,16), \\ \epsilon & \text{on}\ [17,18) \cup [21,22), \\ -\epsilon & \text{on}\ [19,20) \cup [23,24), \end{array}\right. f_y(t) = \left\{\begin{array}{@{}ll } \beta & \text{on}\ [1,2), \\ \gamma & \text{on}\ [3,4), \\ -\gamma & \text{on}\ [5,6), \\ \delta & \text{on}\ [8,9) \cup [12,13), \\ -\delta & \text{on}\ [10,11) \cup [14,15),\\ \epsilon & \text{on}\ [16,17) \cup [18,19), \\ -\epsilon & \text{on}\ [20,21) \cup [22,23), \end{array}\right. \end{gather}
and, if 
$h_3(T) < 0$, we instead define
\begin{gather} f_x(t) = \left\{\begin{array}{@{}ll} \alpha & \text{on}\ [0,1), \\ \gamma & \text{on}\ [3,4), \\ -\gamma & \text{on}\ [5,6), \\ \delta & \text{on}\ [8,9) \cup [12,13), \\ -\delta & \text{on}\ [10,11) \cup [14,15),\\ \epsilon & \text{on}\ [16,17) \cup [18,19), \\ -\epsilon & \text{on}\ [20,21) \cup [22,23), \end{array}\right. f_y(t) = \left\{\begin{array}{@{}ll} \beta & \text{on}\ [1,2), \\ \gamma & \text{on}\ [2,3), \\ -\gamma & \text{on}\ [4,5), \\ \delta & \text{on}\ [6,8) \cup [13,14), \\ -\delta & \text{on}\ [9,10) \cup [11,12) \cup [15,16), \\ \epsilon & \text{on}\ [17,18) \cup [21,22), \\ -\epsilon & \text{on}\ [19,20) \cup [23,24), \end{array}\right. \end{gather}
Smoother controls can be obtained by multiplying the constant value set on an interval 
$[T_0,T_0+1]$ by any function 
$\varphi$ supported on this interval whose integral is unity. For the control policies displayed in figure 8(d), we have made use of 
$\varphi (t) = 2 \sin ^2 ({\rm \pi} (t-T_0))$.
Open access
