A general perturbative approach for bead-based microswimmers reveals rich self-propulsion phenomena

TL;DR

This paper introduces a versatile perturbative method for analyzing bead-based microswimmers, revealing complex self-propulsion behaviors and extending understanding beyond existing models through validation and broad applicability.

Contribution

A general perturbative calculation scheme for bead-based microswimmers driven by arbitrary forces, incorporating hydrodynamic interactions, and validated against numerical and analytical models.

Findings

Identified qualitative differences in far-field flow patterns.

Revealed richer swimmer behaviors than previous models predicted.

Validated approach with 3-bead assemblies and existing geometries.

Abstract

Studies of model microswimmers have significantly contributed to the understanding of the principles of self-propulsion we have today. However, only a small number of microswimmer types have been amenable to analytic modeling, and further development of such approaches is necessary to identify the key features of these active systems. Here we present a general perturbative calculation scheme for swimmers composed of beads interacting by harmonic potentials, driven by an arbitrary force protocol. The hydrodynamic interactions are treated using the Oseen and Rotne-Pragner approximations. We validate our approach by using 3 bead assemblies and comparing the results with the numerically obtained full-solutions of the governing equations of motion, as well as with the existing analytic models for a linear and a triangular swimmer geometries. While recovering the relation between the force and swimming velocity, our detailed analysis and the controlled level of approximation allow us to find qualitative differences already in the far field flow of the devices. Consequently, we are able to identify a behavior of the swimmer that is richer than predicted in previous models. Given its generality, the framework can be applied to any swimmer geometry, driving protocol and beads interactions, as well as in many swimmers problems.

Equations

$${\varphi }_{ij}({\vec{R}}_i,{\vec{R}}_j)\equiv {\varphi }_{ij}(|{\vec{R}}_i-{\vec{R}}_j|)=\frac{k_{ij}}{2}{\left(|{\vec{R}}_i-{\vec{R}}_j|-L_{ij}\right)}^2.$$

$$\varphi (\mathbf{R})=\sum _{conn.\mathrm{pairs}\mathrm{ij}}{\varphi }_{ij}({\vec{R}}_i,{\vec{R}}_j),$$

$$\nabla p(\vec{r},t)-\eta {\nabla }^2\vec{u}(\vec{r},t)=\vec{F}(\vec{r},t).$$

$$\nabla \vec{u}(\vec{r},t)=0.$$

$$\begin{array}{rl}& \vec{u}(\vec{r})=\int d^3\vec{s}{\hat{T}}^{\mathrm{O}}(\vec{r}-\vec{s})\vec{F}(\vec{s})={\hat{T}}^{\mathrm{O}}(\vec{r}-{\vec{r}}_0){\vec{F}}_0=\\ & \frac{1}{8\pi \eta |\vec{r}-{\vec{r}}_0|}\left(\hat{1}+\frac{(\vec{r}-{\vec{r}}_0)\otimes (\vec{r}-{\vec{r}}_0)}{(\vec{r}-{\vec{r}}_0{)}^2}\right){\vec{F}}_0,\end{array}$$

$$\begin{array}{r}{\hat{T}}^{\mathrm{RP}}(\vec{r}-{\vec{r}}_0)=\frac{1}{8\pi \eta |\vec{r}-{\vec{r}}_0|}\left(\hat{1}+\frac{(\vec{r}-{\vec{r}}_0)\otimes (\vec{r}-{\vec{r}}_0)}{(\vec{r}-{\vec{r}}_0{)}^2}\right)+\\ \frac{a_i^2+a_j^2}{24\pi \eta |\vec{r}-{\vec{r}}_0{|}^3}\left(\hat{1}-3\frac{(\vec{r}-{\vec{r}}_0)\otimes (\vec{r}-{\vec{r}}_0)}{(\vec{r}-{\vec{r}}_0{)}^2}\right).\end{array}$$

$$\begin{split}\underline{\mu}(\bm{R})\equiv\underline{\mu}((\vec{R}_{1},...,\vec{R}_{n}))=\begin{cases}\begin{array}[]{lr}\hat{1}/(6\pi\eta a_{i}),&\text{for }i=j\\ \hat{T}(\vec{R}_{i}-\vec{R}_{j}),&\text{for }i\neq j\end{array}\end{cases},\end{split}$$

$${\vec{E}}_i(t,\mathbf{R})=A_i{\vec{E}}_i(\mathbf{R})\sin (\omega t+{\alpha }_i),$$

$$\frac{d}{dt}\mathbf{R}=\underline{\mu }(\mathbf{R})[\mathbf{E}(t,\mathbf{R})+\mathbf{G}(\mathbf{R})],$$

$$\mathbf{E}(t,\mathbf{R})=({\vec{E}}_1(t,\mathbf{R}),...,{\vec{E}}_n(t,\mathbf{R})),\mathbf{G}(\mathbf{R})=-{\nabla }_{\mathbf{R}}\varphi (\mathbf{R}).$$

$$a_i^{\prime }:=a_i/a,k_{ij}^{\prime }:=k_{ij}/k,L_{ij}^{\prime }:=L_{ij}/L,A_i^{\prime }:=A_i/A,$$

$$\frac{d}{d\tau }{\mathbf{R}}^{\prime }={\underline{\mu }}^{\prime }({\mathbf{R}}^{\prime })[ϵ{\mathbf{E}}^{\prime }({\mathbf{R}}^{\prime },\tau )+{\mathbf{G}}^{\prime }({\mathbf{R}}^{\prime })],$$

$${\mathbf{R}}^{\prime }:=\frac{1}{a}({\vec{R}}_1,...,{\vec{R}}_n);{\mathbf{E}}^{\prime }(\tau ,{\mathbf{R}}^{\prime }):=\frac{1}{A}({\vec{E}}_1(\tau t_V,a{\mathbf{R}}^{\prime }),...,{\vec{E}}_n(\tau t_V,a{\mathbf{R}}^{\prime })),$$

$$\begin{array}{r}{\underline{\mu }}^{\prime }({\mathbf{R}}^{\prime })=6\pi \eta a\underline{\mu }(a{\mathbf{R}}^{\prime }),\end{array}$$

$${\mathbf{G}}^{\prime }({\mathbf{R}}^{\prime })=\frac{1}{ka}\mathbf{G}(a{\mathbf{R}}^{\prime }).$$

$${\mathbf{\xi }}^{\prime }:={\mathbf{R}}^{\prime }-{\mathbf{R}}^{\prime \mathrm{eq}}=ϵ{\mathbf{\xi }}^{\prime (1)}+{ϵ}^2{\mathbf{\xi }}^{\prime (2)}+{ϵ}^3{\mathbf{\xi }}^{\prime (3)}+...,$$

$$\begin{array}{rl}& \frac{d}{d\tau }({\mathbf{R}}^{\prime \mathrm{eq}}+{\mathbf{\xi }}^{\prime })=\\ & \left[{\underline{\mu }}^{\prime }({\mathbf{R}}^{\prime \mathrm{eq}})+{\partial }_{\alpha }{\underline{\mu }}^{\prime }({\mathbf{R}}^{\prime \mathrm{eq}}){\mathbf{\xi }}^{\prime \alpha }+\frac{1}{2}{\partial }_{\alpha }{\partial }_{\beta }{\underline{\mu }}^{\prime }({\mathbf{R}}^{\prime \mathrm{eq}}){\mathbf{\xi }}^{\prime \alpha }{\mathbf{\xi }}^{\prime \beta }+...\right]\times \\ & [\left(ϵ{\mathbf{E}}^{\prime }({\mathbf{R}}^{\prime \mathrm{eq}},\tau )+ϵ{\partial }_{\alpha }{\mathbf{E}}^{\prime }({\mathbf{R}}^{\prime \mathrm{eq}},\tau ){\mathbf{\xi }}^{\prime \alpha }+\frac{ϵ}{2}{\partial }_{\alpha }{\partial }_{\beta }{\mathbf{E}}^{\prime }({\mathbf{R}}^{\prime \mathrm{eq}},\tau ){\mathbf{\xi }}^{\prime \alpha }{\mathbf{\xi }}^{\prime \beta }+...\right)\\ & +\left({\mathbf{G}}^{\prime }({\mathbf{R}}^{\prime \mathrm{eq}})+{\partial }_{\alpha }{\mathbf{G}}^{\prime }({\mathbf{R}}^{\prime \mathrm{eq}}){\mathbf{\xi }}^{\prime \alpha }+\frac{1}{2}{\partial }_{\alpha }{\partial }_{\beta }{\mathbf{G}}^{\prime }({\mathbf{R}}^{\prime \mathrm{eq}}){\mathbf{\xi }}^{\prime \alpha }{\mathbf{\xi }}^{\prime \beta }+...\right)].\end{array}$$

$$\frac{d}{d\tau }{\mathbf{\xi }}^{\prime (p)}={\underline{K}}^{\prime }{\mathbf{\xi }}^{\prime (p)}+{\mathbf{S}}^{\prime (p)}(\tau ),$$

$${\underline{K}}^{\prime }={\underline{\mu }}^{\prime }({\mathbf{R}}^{\prime \mathrm{eq}})\circ ({\nabla }_{{\mathbf{R}}^{\prime }}{\mathbf{G}}^{\prime })({\mathbf{R}}^{\prime \mathrm{eq}}),$$

$${\mathbf{S}}^{\prime (1)}(\tau )={\underline{\mu }}^{\prime }({\mathbf{R}}^{\prime \mathrm{eq}}){\mathbf{E}}^{\prime }({\mathbf{R}}^{\prime \mathrm{eq}},\tau ),$$

$$\begin{gathered} {\mathbf{S}}^{\prime (2)}= \\ {\underline{\mu }}^{\prime }({\mathbf{R}}^{\prime \mathrm{eq}})\left(\frac{1}{2}{\partial }_{\alpha }{\partial }_{\beta }{\mathbf{G}}^{\prime }({\mathbf{R}}^{\prime \mathrm{eq}}){\mathbf{\xi }}^{\prime (1)\alpha }{\mathbf{\xi }}^{\prime (1)\beta }+{\partial }_{\alpha }{\mathbf{E}}^{\prime }({\mathbf{R}}^{\prime \mathrm{eq}},\tau ){\mathbf{\xi }}^{\prime (1)\alpha }\right)+ \\ {\partial }_{\alpha }{\underline{\mu }}^{\prime }({\mathbf{R}}^{\prime \mathrm{eq}}){\mathbf{\xi }}^{\prime (1)\alpha }\left({\mathbf{E}}^{\prime }({\mathbf{R}}^{\prime \mathrm{eq}},\tau )+{\partial }_{\beta }{\mathbf{G}}^{\prime }({\mathbf{R}}^{\prime \mathrm{eq}}){\mathbf{\xi }}^{\prime (1)\beta }\right). \end{gathered}$$

$${\mathbf{\xi }}_{hom}^{\prime }=\sum _{\alpha }{\mathbf{X}}_{\alpha }\exp ({\lambda }_{(\alpha )}\tau ),$$

$${\vec{u}}_{\mathrm{fluid}}^{\prime }({\vec{r}}^{\prime })=\sum _{i=1}^n{\hat{T}}^{\prime }({\vec{r}}^{\prime }-{\vec{R}}_i^{\prime }){\left(ϵ{\mathbf{E}}^{\prime }({\mathbf{R}}^{\prime },\tau )+{\mathbf{G}}^{\prime }({\mathbf{R}}^{\prime })\right)}^i,$$

$$\begin{array}{r}{\vec{u}}_{\mathrm{fluid}}^{\prime }({\vec{r}}^{\prime })=ϵ\sum _{i=1}^n{\hat{T}}^{\prime }({\vec{r}}^{\prime }-{\vec{R}}_i^{\prime \mathrm{eq}}){\left({\mathbf{E}}^{\prime }({\mathbf{R}}^{\prime \mathrm{eq}},\tau )+{\partial }_{\alpha }{\mathbf{G}}^{\prime }({\mathbf{R}}^{\prime \mathrm{eq}}){\mathbf{\xi }}^{\prime (1)\alpha }\right)}^i+\\ {ϵ}^2[\sum _{i=1}^n{\mathbf{\xi }}_i^{\prime (1)}\cdot {\nabla }_{{\vec{R}}_i^{\prime }}{\hat{T}}^{\prime }({\vec{r}}^{\prime }-{\vec{R}}_i^{\prime }){|}_{{\vec{R}}_i^{\prime }={\vec{R}}_i^{\prime \mathrm{eq}}}{\left({\mathbf{E}}^{\prime }({\mathbf{R}}^{\prime \mathrm{eq}},\tau )+{\partial }_{\alpha }{\mathbf{G}}^{\prime }({\mathbf{R}}^{\prime \mathrm{eq}}){\mathbf{\xi }}^{\prime (1)\alpha }\right)}^i+\\ \sum _{i=1}^n{\hat{T}}^{\prime }({\vec{r}}^{\prime }-{\vec{R}}_i^{\prime \mathrm{eq}})({\partial }_{\alpha }{\mathbf{E}}^{\prime }({\mathbf{R}}^{\prime \mathrm{eq}},\tau ){\mathbf{\xi }}^{\prime (1)\alpha }+{\partial }_{\alpha }{\mathbf{G}}^{\prime }({\mathbf{R}}^{\prime \mathrm{eq}}){\mathbf{\xi }}^{\prime (2)\alpha }+\\ {\frac{1}{2}{\partial }_{\alpha }{\partial }_{\beta }{\mathbf{G}}^{\prime }({\mathbf{R}}^{\prime \mathrm{eq}}){\mathbf{\xi }}^{\prime (1)\alpha }{\mathbf{\xi }}^{\prime (1)\beta })}^i]+{ϵ}^3...,\end{array}$$

$${\vec{E}}_1(t)=A\sin (\omega t),{\vec{E}}_2(t)=A\sin (\omega t+\alpha ),{\vec{E}}_3(t)=-{\vec{E}}_1(t)-{\vec{E}}_2(t),$$

$${\mathbf{S}}^{\prime (1)}(\tau )={\mathbf{S}}_{s1}^{\prime (1)}\sin (\Gamma \tau )+{\mathbf{S}}_{c1}^{\prime (1)}\cos (\Gamma \tau ),$$

$${\mathbf{S}}^{\prime }(\tau )={\mathbf{S}}_{sf}^{\prime }\sin (f\Gamma \tau )+{\mathbf{S}}_{cf}^{\prime }\cos (f\Gamma \tau ),$$

$$\begin{array}{r}{{\mathbf{\xi }}^{\prime }}_{sf}={\left(f^2{\Gamma }^2\hat{\underline{1}}+{\underline{K}}^{\prime 2}\right)}^{-1}(f\Gamma {\mathbf{S}}_{cf}^{\prime }-{\underline{K}}^{\prime }{\mathbf{S}}_{sf}^{\prime }),\\ {{\mathbf{\xi }}^{\prime }}_{cf}=-{\left(f^2{\Gamma }^2\hat{\underline{1}}+{\underline{K}}^{\prime 2}\right)}^{-1}(f\Gamma {\mathbf{S}}_{sf}^{\prime }+{\underline{K}}^{\prime }{\mathbf{S}}_{cf}^{\prime }).\end{array}$$

$${\mathbf{S}}^{\prime (2)}(\tau )={\mathbf{S}}_{s2}^{\prime (2)}\sin (2\Gamma \tau )+{\mathbf{S}}_{c2}^{\prime (2)}\cos (2\Gamma \tau )+{\mathbf{S}}_{const}^{\prime (2)}.$$

$$\frac{d}{d\tau }{\overline{\mathbf{\xi }}}^{\prime (2)\alpha }(\tau )={\lambda }_{(\alpha )}{\overline{\mathbf{\xi }}}^{\prime (2)\alpha }(\tau )+{\overline{\mathbf{S}}}_{const}^{\prime (2)\alpha },$$

Full text

††: New J. Phys.

Keywords: microswimmers, perturbation theory, active matter, low Reynolds number dynamics, self-propulsion

1 Introduction

The locomotion of swimmers at small scales has been an active area of research in recent years [1], with a variety of microswimmer models being proposed, both experimental [2, 3, 4, 5, 6, 7, 8, 9, 10, 11] and theoretical [12, 13, 14, 15, 16, 17, 18, 19, 20, 21]. A number of these models aims at understanding the propulsion mechanisms of small organisms such as bacteria or algae cells, or at designing artificial microswimmers. Due to the time-independence of the Stokes equation, modelling microswimmers has turned out to be a tradeoff between as little degrees of freedom as possible and enough degrees to break the time-reversal symmetry [22].

A milestone in the field has been the minimalistic swimmer consisting of three spherical beads arranged in a linear fashion, introduced by Najafi and Golestanian [12]. In their model, each neighbouring pair of beads is connected by an extendible arm of a length that is prescribed as a function of time. Calculating the swimming velocity to leading order in the extension of the arms gives rise to a simple intuition for the swimmer’s speed: The displacement of the swimmer corresponding to one swimming stroke is proportional to the area enclosed by the swimmer’s trajectory in the conformation space [13]. This model has been used to investigate the hydrodynamic properties of microswimmers, including the flow fields they produce and their mutual interaction [23, 24, 25], as well as the interaction of a swimmer with a wall [26, 27].

Despite its immense usefulness, this model suffers from the constriction of all internal degrees of freedom by the swimming stroke, namely the internal dynamical behavior of the swimmer cannot react to its surrounding. This was overcome by replacing the arms with springs and prescribing the forces acting on the beads rather than the stroke itself [28]. In this so-called bead-spring swimmer model, the existence of a viscosity maximising the swimming velocity [29] and synchronization effects of the stroke have been discovered [20]. Recently, an altered version of the bead-spring model has been proposed, where the swimmer was driven by periodic changes in the equilibrium lengths of the springs [30, 31].

The boundedness of the linear swimmer to one dimension is broken in a triangular swimmer geometry, allowing for translational as well as rotational motion [32, 14, 16]. This geometry has also been used to model Chlamydomonas reinhardtii and investigate in particular the synchronization between the beating of its two flagella [18, 33, 19]. Experimentally, a triangular swimmer with intrinsic elasticity has been realized by placing ferromagnetic beads subject to an oscillating magnetic field at a water air interface [10, 9] and a similar system has been investigated by means of lattice Boltzmann simulations [34]. The full controllability of a triangular swimmer in the 2D space has been shown recently analytically [21], expanding on the use of the bead spring model [20]. However, the perturbative calculations in [20, 29, 21] hold only in the limit of very large bead separations where the swimming velocity becomes extremely small. Especially since this limit is inaccessible in experiments (see [5, 9]) as external disturbances become exuberant, an investigation of swimmers with bead separation of a few radii is required and still lacking.

In this paper, we fill this gap by presenting a general perturbative framework to calculate the full behavior of arbitrarily shaped bead-spring swimmers, i.e. not only their stroke and swimming velocity, but also the average deformation and the produced flow field with tunable precision of the results in terms of the bead separation. Our calculation is essentially different from previous calculations as we split the equation of motion by orders of the driving force and systematically solve the long-term limit of each order [28, 29, 21]. Doing so significantly increases the precision of the result as the assumption of very large bead separation is not made in our calculation. The so-obtained correction is found to be small for the symmetric linear swimmer, but grows strongly when the swimmer becomes asymmetric. We investigate the dynamics of the triangular swimmer, and find a transient phase of rotation towards a stable steady state in which the swimmer propagates without rotating. In the state of purely translational motion, the swimmer produces a dipolar flow field, which has not been reported for the triangular bead-spring swimmer yet.

The remainder of the article is structured as follows: in the next section we introduce the general model of bead-spring microswimmers. We proceed by analyzing the equation of motion by means of perturbation theory and present the calculation of the swimmer’s velocity, flow field and the beads’s trajectories step by step in section 3. Subsequently, in section 4, we apply this framework to the linear swimmer as a benchmark. Finally, in section 5, we use it to investigate the triangular swimmer in both external and internal driving. Section 6 contains the discussion and conclusion.

2 Model

We consider a microswimmer composed of $n$ spherical beads of radius $a_{i}$ in the $d$-dimensional space. Some pairs of beads are connected by linear springs which we assume here to be harmonic corresponding to the interaction potential

[TABLE]

Here, $\vec{R}_{i}$ and $\vec{R}_{j}$ are the positions of the beads, $k_{ij}$ the stiffness and $L_{ij}$ the length of the spring connecting bead $i$ and $j$ in the swimmer’s mechanical equilibrium. Our approach can also easily be adopted to more complex interaction potentials, e.g. magneto-capillary potentials (see [35]). Note that the springs are not tied to a certain direction but can freely rotate as the beads move in the fluid. The total spring potential of the device is given by

[TABLE]

where we sum over all pairs of beads that are connected and define $\bm{R}=(\vec{R}_{1},...,\vec{R}_{n})$ as a vector with $n\cdot d$ components.

We assume that the Reynolds number of the beads in the fluid is small and that the relaxation of the fluid takes place sufficiently fast [36] such that the fluid dynamics can be described by the Stokes equation

[TABLE]

Here, $p(\vec{r},t)$ denotes the pressure in the fluid, $\vec{u}(\vec{r},t)$ the velocity of the fluid, $\vec{F}(\vec{r},t)$ the force density applied to the fluid and $\eta$ its dynamic viscosity. $\vec{r}$ and $t$ denote position vector and time, respectively. The fluid is assumed to be incompressible:

[TABLE]

A point force $\vec{F}(\vec{r})=\vec{F}_{0}\delta(\vec{r}-\vec{r}_{0})$, acting on the fluid at position $\vec{r}_{0}$, induces a fluid flow at position $\vec{r}$ given by the Oseen tensor $\hat{T}^{\mathrm{O}}(\vec{r}-\vec{r}_{0})$ as

[TABLE]

with $\otimes$ the tensor product and $\hat{1}$ the $d$-dimensional unity matrix. For a large separation between two beads in terms of their radii, the Oseen tensor is a sufficient description for their interaction. Being interested in swimmers with small bead separations, we make use of the Rotne-Prager approximation [37], which is given by

[TABLE]

The mobility matrix of an ensemble of $n$ spheres is a $(n\cdot d)\times(n\cdot d)$ matrix defined in terms of $d\times d$ blocks as

[TABLE]

with $\hat{T}$ either the Oseen or the Rotne-Prager tensor and $i,j=1,...,n$. In (7), we account for both, the Stokes drag as well as the interaction of the beads due to the fluid. The springs are assumed to be not interacting with the fluid.

An oscillating force $\vec{E}_{i}(t,\bm{R})$ with fixed frequency $\omega$ is acting on each bead $i$ as

[TABLE]

with $A_{i}$ encoding the amplitude of the driving for each bead $i$ and $\alpha_{i}$ the phase shift associated to each bead. The vector $\vec{E}_{i}(\bm{R})$ is dimensionless and subject to a suitable normalization that will be chosen specifically for each geometry and driving protocol. As indicated above, we allow for a dependence of the driving forces on the current configuration of the swimmer. This is required if e.g. the driving forces result from spatially dependent magnetic or electric fields, or if certain demands to the driving protocol shall be fulfilled, like a constantly zero total torque for an effectively rotating swimmer.

The temporal evolution of the system is governed by the constitutive equation of the mobility matrix on the $(n\cdot d)$-dimensional configuration space of the beads as

[TABLE]

with

[TABLE]

$\nabla_{\bm{R}}$ denotes the gradient of a function with respect to all $n\cdot d$ components of $\bm{R}$. In view of rescaling of the equations of motion, we introduce for each family of parameters $a_{i},k_{ij},L_{ij}$ and $A_{i}$ a characteristic value ($a,k,L$ and $A$, respectively) and define dimensionless parameters by

[TABLE]

which become 1 in the case of equal parameters of one type.

3 Analysis

To calculate how the swimmer behaves around a stable mechanical equilibrium, we develop a perturbative approach that allows to split the equation of motion by orders in the driving force. We firstly rescale the equation of motion using the characteristic length $a$ and the characteristic time $t_{V}=(6\pi\eta a)/k$ [29], denoted as viscous time. We find the effective parameters to be $\epsilon=A/(ka)$ for the rescaled driving force amplitude, $\nu=a/L$ the aspect ratio between bead radius and separation, and $\Gamma=\omega t_{V}$ the rescaled driving frequency. Rescaled variables are marked with an additional dash and the rescaled time is $\tau:=t/t_{V}$. The equation of motion can then be re-expressed as

[TABLE]

with

[TABLE]

[TABLE]

and

[TABLE]

The equation of motion (12) can be solved with a perturbative approach in the vicinity of $\epsilon=0$. Therefore, we employ a suitable power series ansatz in $\epsilon$ for the displacement out of the equilibrium

[TABLE]

where $\bm{R}^{\prime\mathrm{eq}}$ is the rescaled equilibrium configuration of the swimmer. A Taylor expansion of all $\bm{R}^{\prime}$-dependent parts of (12) around $\bm{R}^{\prime\mathrm{eq}}$ yields

[TABLE]

$\alpha,\beta$ are summed over when appearing repeatedly and go from 1 to $n\cdot d$, $\bm{\xi}^{\prime\alpha}$ denotes the $\alpha$-th component of $\bm{\xi}^{\prime}$, and $\partial_{\alpha}$ the derivative with respect to the $\alpha$-th component of $\bm{R}^{\prime}$. The $\tau$-derivative of $\bm{R}^{\prime\mathrm{eq}}$ is zero as well as the spring forces evaluated in the equilibrium $\bm{G}^{\prime}(\bm{R}^{\prime\mathrm{eq}})$.

We proceed by replacing $\bm{\xi}^{\prime}$ in (17) by its power series in $\epsilon$ (16). Ordering and splitting the resulting equation by powers of $\epsilon$ yields a vectorial equation for each order $p=1,2,...$ in $\epsilon$. One finds that each of them is of the generic form

[TABLE]

with

[TABLE]

where $\circ$ denotes the matrix multiplication and $(\nabla_{\bm{R}^{\prime}}\bm{G}^{\prime})(\bm{R}^{\prime\mathrm{eq}})$ is the Jacobian matrix of the spring forces, evaluated at the swimmer’s equilibrium. $\bm{S}^{\prime(p)}$ is a term that only depends on the displacements $\bm{\xi}^{\prime(1)},...,\bm{\xi}^{\prime(p-1)}$ and on the derivatives (of first and higher order) of $\underline{\mu}^{\prime},\bm{E}^{\prime}$ and $\bm{G}^{\prime}$, evaluated at $\bm{R}^{\prime\mathrm{eq}}$. Since $\bm{S}^{\prime(p)}$ does not depend on $\bm{\xi}^{\prime(p)}$, it can be considered as a source term that is known assuming (18) are solved in ascending order in $p$. Consequently, the first two source terms read

[TABLE]

[TABLE]

The matrix $\underline{K}^{\prime}$ maps the displacement vector $\bm{\xi}^{\prime}$ to the velocity vector $(d/dt)\bm{R}^{\prime}=(d/dt)\bm{\xi}^{\prime}=\underline{K}^{\prime}\bm{\xi}^{\prime}$ that emerges from the displacement in the situation of zero external forces. We assume the mobility matrix $\underline{\mu}^{\prime}(\bm{R}^{\prime\mathrm{eq}})$ to be positive definite and symmetric. The gradient of the spring forces, $(\nabla_{\bm{R}^{\prime}}\bm{G}^{\prime})(\bm{R}^{\prime\mathrm{eq}})$, has to be negative semi-definite due to the stability of the swimmer’s equilibrium and symmetric as it is the Hessian of the spring potential. For the latter matrix, we distinguish between eigenvalues being zero, associated to translations and rotations of the whole swimmer, and negative eigenvalues, associated to deformations of the swimmer. Then, it is easy to show that, similarly to $(\nabla_{\bm{R}^{\prime}}\bm{G}^{\prime})(\bm{R}^{\prime\mathrm{eq}})$, the matrix $\underline{K}^{\prime}$, defined by (19), is diagonalizable and has only non-positive eigenvalues, where zero as an eigenvalues is associated to translations/rotations and negative eigenvalues are associated to internal degrees of freedom [38].

The explicit way to solve a set of differential equations like (18), accounting for certain initial conditions, is to split the initial conditions by orders of $\epsilon$, find the full solution for each order, and adjust it to these initial conditions by a suitable choice of the parameters in the homogeneous solution. Knowing that $\underline{K}^{\prime}$ has only non-positive eigenvalues $\lambda_{(\alpha)}$, the solution to the homogeneous counterpart to (18) ($\bm{S}^{\prime(p)}=0$) can be written as

[TABLE]

with $\bm{X}_{\alpha}$ a suitably scaled eigenvector of $\underline{K}^{\prime}$ corresponding to the eigenvalue $\lambda_{(\alpha)}$. For translational and rotational degrees of freedom, one has $\lambda_{(\alpha)}=0$ and hence a constant solution. The displacements corresponding to the internal degrees of freedom are exponentially decaying and hence go to zero for large $\tau$. Therefore, the homogeneous solution describes the relaxation of the swimmer in the absence of driving from arbitrary initial conditions. Consequently, all solutions to the inhomogeneous equations ($\bm{S}^{\prime(p)}\neq 0$) differ only by the constant homogeneous solution for large $\tau$. It suffices to find a single arbitrary solution to the full (18) in order to determine the swimming velocity and the deformation of the swimmer. For a fixed source term, the swimmer’s behavior is hence independent of the initial conditions for large $\tau$. The eigenvalues corresponding to internal degrees of freedom are of the order of 1, such that the displacements corresponding to the internal degrees of freedom decay exponentially with a characteristic time of the order of $t_{V}$. In the following calculations, we are interested in the behavior of the swimmer on times scales larger than $t_{V}$. Therefore, we will neglect all terms decaying exponentially at time scale $t_{V}$ arising from the initial conditions in the displacements $\bm{\xi}^{\prime(1)},...,\bm{\xi}^{\prime(p-1)}$ when calculating the source term of each order $p$.

We point out here that if the system is not invariant under the translational or rotational degree of freedom, e.g. because the driving forces are not invariant under these transformations, the source term $\bm{S}^{\prime(p)}(\tau)$ explicitly depends on those degrees of freedom. Via that pathway, the initial conditions may actually have an impact on the swimmer’s behavior. This will be the case for the triangular swimmer in external driving as it will be discussed in section 5.1.

The rescaled flow field $\vec{u}^{\prime}_{\mathrm{fluid}}(\vec{r}^{\prime})$ generated by the swimmer, expressed in dependence of the rescaled position $\vec{r}^{\prime}$, is given by

[TABLE]

with $\hat{T}^{\prime}(\vec{x}^{\prime}):=6\pi\eta a\hat{T}(a\vec{x}^{\prime})$ the rescaled Oseen/Rotne-Prager tensor and $(\bullet)^{i}$ the $i$-th part in a decomposition of $\bullet$ into $n$ partial vectors of length $d$, i.e. the components associated to the $i$-th bead. Given the solutions of (18) up to order $p$, the flow field can be calculated up to the same order $p$ in $\epsilon$ by expanding (23). We here state the explicit expression up to the second order in $\epsilon$:

[TABLE]

with $\nabla_{\vec{R}^{\prime}_{i}}$ the gradient with respect to $\vec{R}^{\prime}_{i}$. Note that having solved for the displacements in advance, we simply need to insert them into (24) to obtain the flow fields produced by the swimmer.

4 Linear three-bead swimmer

The simplest one-dimensional bead-spring swimmer able to swim at low Reynolds number is the linear three-bead swimmer (figure 1). A swimmer with two beads comes with only one internal degree of freedom which is not sufficient to break the time reversal symmetry [22]. It has been studied in detail in previous analytical works [28, 20, 29], where two similar perturbative calculations were performed. In the latter works [20, 29], the results for the swimming velocity were expanded and truncated in orders of $\nu$. With the precision of the Oseen tensor, this calculation does not allow for a predictive result at higher order than the leading order $\nu^{2}$.

We will show that the result of our approach for the swimming velocity employing the Oseen tensor is correct up to order $\nu^{3}$ and that using the Rotne-Prager tensor hydrodynamics has an impact at orders $\nu^{4}$ and higher. Furthermore, we find that our approach coincides at leading order $\nu^{2}$ to the aforementioned previous results but differs at order $\nu^{3}$, explaining why our results also hold at order $\nu^{3}$.

The swimmer consists of three beads for which we choose radii $a_{1}=a_{2}=a$ and $a_{3}$, respectively [28]. The beads are connected by two identical harmonic springs of stiffness $k$ and equilibrium length $L$. The driving forces $\vec{E}_{i}(t)$, specified as

[TABLE]

act on the beads, ensuring that the total force vanishes. The second derivative of $\bm{G}^{\prime}(\bm{R}^{\prime})$ vanishes, because the swimmer is restricted to one dimension. Furthermore, the driving forces are not spatially dependent and hence all their spatial derivatives vanish too.

We here calculate the displacement $\bm{\xi}^{\prime}$ up to second order in $\epsilon$, pointing out that the displacement to higher order can be obtained analogously. The source term $\bm{S}^{\prime(1)}(\tau)$, (20), composes of purely oscillating contributions with frequency $\Gamma$,

[TABLE]

with the indices $s1$ and $c1$ denoting first the correspondence to $\sin$ or $\cos$ and second indicating the argument of the trigonometric function in multiples of $\Gamma\tau$. To safe efforts later, we calculate here the solution for a more general source term given by

[TABLE]

with $f$ an arbitrary positive integer. Given the linear nature of (18), a suitable ansatz for the displacement is $\bm{\xi}^{\prime}(\tau)=\bm{\xi}^{\prime}_{sf}\sin(f\Gamma\tau)+\bm{\xi}^{\prime}_{cf}\cos(f\Gamma\tau)$. The resulting solutions to (18) read

[TABLE]

These results are in full agreement with [28].

Having calculated the first order displacement in $\epsilon$, we proceed by calculating the second order source term (21). We find that it contains oscillating terms of the frequency $2\Gamma$ and a constant contribution:

[TABLE]

Again, the linearity of (18) allows us to calculate its solution for each summand in (29) separately and to add up the results to obtain a full solution. Firstly, the oscillating parts in the source term yield oscillating contributions to the $\epsilon^{2}$-displacement (see (27),(28)), which contribute to the stroke of the swimmer. Secondly, (18) with the constant source term alone can be easily solved by expressing the source term in the eigenbasis of $\underline{K}^{\prime}$, where this matrix becomes diagonal and the equations separate. We find for each eigenvalue $\lambda_{(\alpha)}$, $\alpha=1,...,n\cdot d$:

[TABLE]

with the overline indicating the expression of a vector in the eigenbasis of $\underline{K}^{\prime}$. For the translational and rotational degrees of freedom with respect to $\underline{K}^{\prime}$ we have $\lambda_{(\alpha)}=0$ and the solution is a linear function in time $\tau$ plus a constant which we neglect here, as it is determined by the choice initial conditions:

[TABLE]

For the internal degrees of freedom, $\lambda_{(\alpha)}<0$, the solution to (30) is constant in time plus an exponentially decaying term that we neglect since we are only interested in the limit of $\tau\gg t_{V}$:

[TABLE]

In effect, this describes a deformation of the swimmer such that the beads do not oscillate around their mechanical equilibrium but around a different, deformed configuration.

Analyzing the eigensystem of $\underline{K}^{\prime}$, we find that the linear three-bead spring swimmer with equal radii and restricted to one dimension has one translational mode $X_{1}$ and two internal eigenmodes $X_{2}$ and $X_{3}$ (see figure 2 a) with

[TABLE]

The mode $X_{3}$ is orthogonal (with respect to the standard scalar product) to $X_{1}$ and $X_{2}$. The modes $X_{1}$ and $X_{2}$ are in general not orthogonal to each other, but become orthogonal for $\nu\to 0$. In our calculation, the swimming velocity can be read off from the component parallel to $X_{1}$ in the decomposition of $\bm{S}^{\prime(2)}_{const}$ in terms of the eigenvectors of $\underline{K}$ (green arrow in figure 2 a).

In a previous work [28], the swimming velocity was effectively calculated as

[TABLE]

with $\bm{S}^{\prime(2)i}_{const}$ being the part of the constant source term associated to bead $i$. In the case of equal radii, this calculation is equivalent to an orthogonal projection of the source term onto the vector $X_{1}$ and reading off the velocity from the projected source term in multiples of $X_{1}$ (blue arrow in figure 2 a). Due to $X_{2}\perp X_{1}$ for $\nu\to 0$, the axis projected onto is in this limit orthogonal to the two internal modes and both the projection and the decomposition of the source term yield the same result for the swimming velocity. This explains why both calculations agree for $\nu\to 0$, but also why they differ for finite values of $\nu$. This can also be seen in the explicit ratio between our result $u_{Z}$, given by

[TABLE]

with $B=(a/t_{V})\,\epsilon^{2}\Gamma\nu^{2}\,(9\nu(21\nu-22)+56)$, and the result $u_{0}$ stated in [28]. This ratio (36) converges to 1 for $\nu\to 0$:

[TABLE]

The differences discussed above stem in part from the fact that the original perturbative approach [28] calculates the oscillations of the beads around the undisturbed swimmer shape. In contrast, in the current scheme, the swimmer’s average shape and the oscillations of the beads around it are obtained simultaneously. Actually, both the mean deformation and the swimming velocity (the translation mode) arise from the constant contribution to the second order source term $\bm{S}^{\prime(2)}_{const}$ (29). The distances $d_{1}$ between bead 1 and 2 and $d_{2}$ between bead 2 and 3 in the swimmer’s mean configuration, around which the beads oscillate harmonically, are given up to order $\epsilon^{2}$ and for $a_{3}=a$ by

[TABLE]

and

[TABLE]

Hence, the ratio between deformation and swimming velocity, both of order $\epsilon^{2}$, is a simple geometrical factor. Also, the amplitude of the deformation obeys a similar frequency dependence as the swimming velocity itself and decays similarly as $1/\nu^{2}$ for large bead separations.

The comparison to numerical calculations, done by numerically integrating the equation of motion (12), shows a very neat agreement between our result and the numerics with errors below 0.1 % (figure 2 b). Comparing our result with the one obtained in [28], we find for $a_{3}=a$ (i.e. equal radii) and $\nu\approx 1/10$ small differences in the range of percents, but the difference increases drastically for increasing values of $a_{3}$ (figure 2 b), i.e. when the swimmer becomes asymmetric. We observe that for $a_{3}\neq a$, all pairs of eigenvectors are in general no more orthogonal, even in the limit $\nu\to 0$. Also, (34) describes no more a projection onto $X_{1}$ but onto $(1,1,a_{3})$. This explains why the difference between $u_{0}$ and $u_{Z}$ grows for increasing $a_{3}$ in figure 2 b, yet in the limit $\nu\to 0$ both still agree independently of $a_{3}$.

Despite this quantitative difference to previous results [28], we recover the typical dependencies for bead-spring swimmers which have been reported previously [28, 20, 29, 34]: The swimming velocity scales with the square of the driving forces for small amplitudes, the swimming speed becomes maximal in the vicinity of $t_{V}^{-1}$ and decays as $1/L^{2}$ for large bead separations.

The Rotne-Prager approximation (see Appendix) has only an impact onto the swimming velocity at orders $\nu^{4}$ and higher: Using the Rotne-Prager tensor instead of the Oseen tensor yields an additional term to the mobility matrix scaling as $1/r^{3}$ (see (6)), with $r$ the distance between the beads. This results in additional terms to $\underline{\mu}^{\prime}(\bm{R}^{\prime\mathrm{eq}})$ and $\partial_{\alpha}\underline{\mu}^{\prime}(\bm{R}^{\prime\mathrm{eq}})$, scaling with $\nu^{3}$ and $\nu^{4}$, respectively. A closer investigation of the second order source term (21) shows that the factor multiplied to $\underline{\mu}^{\prime}(\bm{R}^{\prime\mathrm{eq}})$ scales linearly in $\nu$, such that the additional term due to the Rotne-Prager extension has an impact at order $\nu^{4}$ and higher on the source term $\bm{S}^{\prime(2)}$ and likewise on the swimming velocity.

5 Triangular swimmer

5.1 External driving

Triangular bead-spring swimmers have been studied in detail recently [21, 39], where the employed driving protocol prescribes forces on each bead parallel to the adjoining sides of the triangle. By varying the amplitudes of and the phase shifts between the driving forces, the swimmer can be steered on arbitrary trajectories. Both, translational and rotational motion were shown to scale with the square of the driving force. Also here, the perturbative approach used in [21] does only hold in the limit $\nu\to 0$.

The triangular swimmer is composed of three spherical beads (radius $a$) connected by identical springs of equilibrium length $L$ and spring constant $k$. All beads are placed in the $x$-$z$-plane with orientation $\theta$ the angle between the connection of bead 3 to the middle between bead 1 and 2 and the $x$-axis (figure 3). Numerous experimental realizations of microswimmers rely on an external field, commonly an electric or magnetic one [11, 5, 2]. Therefore, we first consider here a toy model swimmer that is subject to an external force field which shall act in one direction only (without restriction of generality the $x$-direction) for the sake of simplicity. For the swimmer to be self-propelled, we demand that all forces acting on the three beads sum up to zero and also have vanishing total torque. We determine the remaining degree of freedom by prescribing that the sum of the squares of all forces is equal to a constant, $(2A)^{2}$, such that the forces explicitly are given by

[TABLE]

with

[TABLE]

Notably, the forces depend explicitly on the configuration of the swimmer, in order to satisfy the force-free and torque-free condition throughout the whole swimming stroke.

In numerical studies, we observe that the swimmer typically undergoes a transient phase during which it both rotates and translates. It finally reaches a steady state in which the motion is purely translational (figure 4 a).

A closer numerical investigation shows that the angular velocity of the swimmer depends in sinusoidal fashion on the instantaneous orientation of the swimmer (figure 4 b). The rotational dynamics of the externally driven triangular swimmer hence is equivalent to the one of an overdamped pendulum. Perturbation theory and numerics consistently show that the angular velocity $\Omega$ at fixed angle scales as

[TABLE]

and attains its maximum close to the inverse viscous time, similarly to the translational velocity. The time scale of the rotational relaxation hence can be estimated as $T_{relax}\sim\Omega^{-1}\sim\epsilon^{-4}\,\nu^{-4}$. In the parameter range for which the perturbation approach applies, this time scale is several orders of magnitude larger than $t_{V}$.

We find stable steady states of the swimmer at $\theta=p\cdot 60^{\circ}$ and unstable steady states at $\theta=p\cdot 60^{\circ}-30^{\circ},\ p=1,...,6$ (see figure 4 c). Obviously, the swimmer in external driving is invariant under a 3-fold rotation. Furthermore, a $180^{\circ}$ rotation inverts the swimming direction but does not affect the stroke in internal coordinates, showing that the stability of states is invariant under a 6-fold rotation. Hence, all stable steady states can be considered equivalent and likewise all unstable states. The swimmer is found to always rotate towards the stable steady state closest to the initial orientation as shown in figure 4 b. The existence of stable and unstable steady states results from prescribing the driving forces with respect to the laboratory frame compared to their prescription in the swimmer’s frame of reference [21]. In the latter protocol, the forces are held constant in the internal coordinates throughout the rotation of the swimmer and hence the swimmer undergoes constant translation and rotation.

In the steady states of the externally driven swimmer, we find for the swimming velocity in the Oseen approximation the following expression

[TABLE]

Performing an expansion and truncation of (42) to leading power of $\nu$ shows that this result is in agreement with the result presented in [21]. We recover the typical $u\sim\epsilon^{2}$ dependence and $u\sim\nu^{2}$ in the limit $\nu\to 0$, which seem characteristic for bead-spring swimmers. We can understand the latter dependence from the fact that swimming emerges from the interplay between the hydrodynamic interaction of parts of the swimmer and variations in their distance, suggesting a swimming speed scaling with the gradient of the hydrodynamic interaction and hence with $\nu^{2}$. In the range of $L>3a$, we find that expression (42) is negative meaning that the swimmer swims towards the base (with respect to the symmetry) of the triangle (see figure 5 a). For sufficiently small amplitudes of the driving, we observe that the beads 1 and 2 move on ellipsoidal trajectories that are tilted towards the middle axis of the swimmer. In the orientation $\theta=\pi/2$, exemplaric for the unstable steady states, the swimmer swims at the same speed as in a stable state (42). We point out that this property is sensible to the way the normalization of the forces is done (see  (39)), i. e. it holds only if the sum of the squares is fixed. In contrast to the stable states, the swimming direction is here pointing towards the tip of the triangle and also the trajectories of the beads 1 and 2 are in this case tilted towards the base of the triangle (figure 5 b).

Analysis of the analytically computed flow fields shows that in contrast to previous works, which reported that triangular swimmers produce in average neutral flow fields at second order in the driving force [21, 39], we here find a non-vanishing average dipolar flow field at order $\epsilon^{2}$ in both states (see figure 5 c, d). Going from stable to the unstable state comes with an exact inversion of the flow field transforming the puller-like swimmer in the stable state into a pusher-like swimmer in the unstable one.

5.2 Internal driving and flow field

In order to provide deeper insights into the swimmer’s pusher or puller character, we repeat the calculation of the flow field for a purely translational driving protocol presented similarly in [21]. We assume forces $F_{ij}$ between each pair of beads $i,j$ that are defined as

[TABLE]

The application of the perturbative approach presented here and the analysis of the resulting flow field shows that the fields at order $\epsilon^{2}$ (in dependence of the rescaled position $\vec{x}^{\prime}$) can be approximated to leading order as a superposition of two force dipoles

[TABLE]

with

[TABLE]

$\vec{e}$ denotes the direction in which the forces act, $\vec{d}$ the separation vector of the two forces of the dipole and $\nabla^{\prime}$ the gradient with respect to rescaled coordinates. In the direction orthogonal to the swimmer plane, the two force dipoles cancel up to order $1/y^{4}$, such that in this direction only the quadrupolar flow field, scaling as $1/y^{3}$, is reminiscent. Investigation of the magnitude of the dipolar flow field shows that a swimmer can be tuned from pusher to puller by changing $\alpha_{2}$ and $\gamma_{2}$ (figure 6, black curve).

Comparison of the strength of the force dipole $f(\alpha_{2},\gamma_{2})$ produced by the swimmer and the swimming velocity $u_{swim}$ shows that their ratio is a function of $\nu$ only:

[TABLE]

As this ratio is positive for physical values of $\nu$, a swimming motion towards the tip of the triangle is always associated to a pusher flow field whereas motion towards the base is associated to a puller character. For $\gamma_{2}=1/2$ and $\alpha_{2}=\pi$, we recover the driving of the stable steady state, showing that the swimmer behaves as a pusher. The unstable steady state is recovered for $\gamma_{2}\to\infty$ making the corresponding swimmer a pusher. Plots of the flow fields in the swimmer’s plane and in the orthogonal plane through the symmetry axis of the swimmer (figure 5 c - f) illustrate the dipolar character of the flow in the plane of the swimmer and its fast decay orthogonal to the swimmer’s plane.

The complexity of the flow field at order $\epsilon^{2}$ increases when the radius of bead 1 is chosen differently, in particular the form (44) does not hold and the fast decay of the magnitude in $y$-direction is lost. Also the curve determining the pusher/puller behavior in dependence of $\alpha_{2}$ and $\gamma_{2}$ is sensitive to changes in the parameter $a_{1}$ (figure 6). Still, the labels of the pusher and puller area are positioned such that they do not only hold for $a_{1}=a$ (solid line) but also for $a_{1}=a/3$ (dot-dashed line) and $a_{1}=3$ (dashed line).

6 Discussion and Conclusions

We presented a general perturbative approach to calculate the trajectories, internal dynamics and the flow fields of swimmers consisting of beads interacting by harmonic potentials organized in an arbitrary geometry and subject to a force protocol of choice. We first applied our method to the linear swimmer as a benchmark. Comparison to previous results shows that the qualitative behavior of the swimmer, as presented in earlier works [20], is recovered. However strong quantitative differences are found in the case of asymmetric swimmers. We showed that for both, the symmetric and asymmetric linear swimmer, our approach yields a swimming velocity that is correct up to order $\nu^{3}$, whereas previous approaches are only correct up to order $\nu^{2}$. This is due to the different ways of extracting the swimming velocity from the constant contribution to the source term. The current formulation maintains consistency and provides corrections that were to the best of our knowledge unaccounted for in previous models. The consequences of these terms are negligible in the case of the symmetric swimmer for which the swimming modes are orthogonal. However, since the deviations from orthogonality increase with enhancing the asymmetry of the linear design, important quantitative differences can be observed in the swimming velocity.

We further investigated the dynamics of an externally driven triangular swimmer with equal radii and discovered the existence of stable and unstable rotational steady states. We showed that the swimmer propagates at the same speed in both states, but in opposite directions with respect to the symmetry axis of the swimmer. In contrast to previous results [21], the average flow field produced by the swimmer was shown to be puller- or pusher-like, depending on the forces driving the swimmer. Actually, the character of the dipolar flow field can be directly associated with the swimming direction, which is, unlike in the stroke-based models, a consequence of the driving. Interestingly, the average flow field of the triangular swimmer with equal radii is strong near the swimmer plane, but shows a quick decay in the direction orthogonal to the swimmer plane, in contrast to the linear swimmer which produces a rotationally symmetric flow field. This may prove to be important when considering assemblies of swimmers.

Apart from investigating planar swimmers composed of more than three beads or three-dimensional bead spring swimmers, the framework presented here is also applicable to an ensemble of interacting microswimmers. With small reformulations, our approach allows for the calculation of relative and absolute translational and angular velocities of a swimmer in the presence of an arbitrary configuration of other ones. The model is, therefore, well suited to provide a comprehensive study of the mechanisms underlying the interaction of microswimmers. While this issue has been addresed in part for dumbbell-shaped swimmers [40, 41], the stroke-based linear swimmer [24, 25] and recently, for a force-based linear swimmer [31], a full discussion of the interaction between bead-based microswimmers is still pending. The tools presented herein can be used as the foundation of such analysis, which is a task that we intend to address in future work.

Acknowledgements The authors acknowledge founding by the European Research Council through the grant MembranesAct ERC Stg 2013-337283 and by the DFG Priority Programme “Microswimmers – From Single Particle Motion to Collective behavior” (SPP 1726). We also thank Alexander Sukhov, Galien Grosjean, Oleg Trosman and Jayant Pande for valuable discussions.

Appendix A Velocity of the linear swimmer in the Rotne-Prager approximation

The result for the swimming velocity of the linear swimmer in the Rotne-Prager approximation is given by

[TABLE]

with

[TABLE]