# Lattice Boltzmann simulations of two linear microswimmers using the immersed
boundary method

D. Geyer  1  S. Ziegler  2  A. Sukhov  1  M. Hubert  2  A.-S. Smith  2  O.
Aouane  1  P. Malgaretti and J. Harting  1  1,3  11affiliationmark:  Helmholtz
Institute Erlangen-Nürnberg for Renewable Energy (IEK-11), Forschungszentrum
Jülich, Cauerstraße 1, 91058 Erlangen, Germany
22affiliationmark:  PULS Group, Department of Physics, Interdisciplinary
Center for Nanostructured Films, Friedrich-Alexander-Universität Erlangen-
Nürnberg, Cauerstraße 3, 91058 Erlangen, Germany
33affiliationmark:  Department of Chemical and Biological Engineering and
Department of Physics, Friedrich-Alexander-Universität Erlangen-Nürnberg,
Fürther Straße 248, 90429 Nürnberg, Germany<EMAIL_ADDRESS>(J.
Harting)

###### Abstract

The performance of a single or the collection of microswimmers strongly
depends on the hydrodynamic coupling among their constituents and themselves.
We present a numerical study for a single and a pair of microswimmers based on
lattice Boltzmann method (LBM) simulations. Our numerical algorithm consists
of two separable parts. Lagrange polynomials provide a discretization of the
microswimmers and the lattice Boltzmann method captures the dynamics of the
surrounding fluid. The two components couple via an immersed boundary method.
We present data for a single swimmer system and our data also show the onset
of collective effects and, in particular, an overall velocity increment of
clusters of swimmers.

###### keywords:

Immersed boundary method, lattice Boltzmann method, finite element method,
microswimmer, collective motion

##  1 Introduction

In his seminal work, Purcell pointed out that the properties of microscopic
objects placed in fluids are significantly different from their macroscopic
counterparts [1]. In particular, Purcell showed that in order to swim in the
low Reynolds number regime, a micrometric swimmer has to move its parts in
such a manner as to break the time inversion symmetry. This fact led to the
well-known “scallop theorem”, which states that in order to attain self-
propulsion in the low Reynolds number regime, at least two degrees of freedom
are needed. Since then, numerous attempts have been done to elucidate the
dynamics of microswimmers by means of theoretical models [2, 3, 4, 5, 6, 7, 8,
9, 10, 11, 12, 13, 14], experimental setups [15, 16, 17, 18, 19, 20, 21, 22,
23, 24, 25] and numerical simulations [26, 27, 28, 29, 30, 31, 32]. In
particular, in Ref. [5] Najafi and Golestanian proposed a theoretical model
that precisely fulfills the requirement of Purcell. Indeed, they offered a
very simple swimmer composed of three aligned solid spherical particles
suspended in a viscous fluid and actuated internally by changing the distances
between neighboring particles. In a Newtonian fluid (see Refs. [33, 34] for a
theoretical and experimental extension of the problem in the case of the
underdamped regime and non-Newtonian fluids), the dynamics of such a swimmer
is fully determined by the two degrees of freedom of the swimmer, namely, the
two distances among subsequent beads [6]. Later, a variation of this swimmer
has been proposed, where harmonic springs connect neighboring particles and
external forces drive the swimmer [3, 9]. Once actuated with a proper
protocol, the motion of the three beads at the swimmer results to be non-
reciprocal and hence leads to a net displacement.

Recently, the focus of theoretical research has shifted towards swimming in
complex environments like channels [35, 36] or near fluid interfaces [31, 37],
and the question of collective swimmer dynamics has become of major interest.
The present work aims at understanding the collective dynamics of several
microswimmers [38, 39]. For this purpose, we consider a relatively simple
situation of two linear microswimmers in different configurations.
Additionally, we exploit the suitability of the employed simulation method for
this task. The structure of the manuscript is as follows: In Sec. 2, we
present the theoretical model and the numerical implementation. In Sec. 3 we
present our numerical results for a single and two microswimmers in diverse
arrangements, which are then discussed in Sec. 4. Finally, in Sec. 5 we
provide some concluding remarks.

##  2 Model

We focus on investigating the behavior of a single and a pair of bead-spring
microswimmers in a resting Newtonian fluid. Each microswimmer consists of
three aligned equal beads connected with springs which is a modification [9]
of the model proposed by Najafi and Golestanian [5]. The distances between the
beads represent the two degrees of freedom necessary to attain self-
propulsion, as per the scallop theorem [1] and depicted in Fig. 1.

Figure 1: Scheme of the linear, three bead swimmer connected with two springs.
The motion direction $v_{\text{swim}}$ depends on the interplay of the
amplitudes and phases of external forces as detailed in [17]. Distances
between beads are measured from the respective centers of mass.

The microswimmer particles consist of rigid spherical shells filled with a
Newtonian fluid with the same viscosity and density as the external fluid. The
fluid-particle interactions are incorporated into the model through the
immersed boundary method (IBM) [40]. The flow field in the entire
computational domain is computed using the lattice Boltzmann method (LBM).

###  2.1 Fluid dynamics

In this study we use the single relaxation time LBM with a Bhatnagar, Gross,
Krook [41] collision operator to solve for the flow field on an Eulerian frame
in the weakly compressible limit. The macroscopic fluid dynamics are recovered
from the mesoscopic Boltzmann equation as detailed in Ref. [42, 43]. The LBM
consists of two steps, namely collision and advection. We adopt a three-
dimensional 19 velocities lattice (D3Q19) model and represent lengths in units
of the lattice spacing $\Delta x$ and times in units of the integration time
step $\Delta t$. The time evolution of the distribution functions
$f_{i}(\textbf{x},t)$ is obtained by solving the discrete lattice Boltzmann
equation in velocity space such that

$f_{i}({\mathbf{x}}+\mathbf{c}_{i}\Delta t,t+\Delta
t)-f_{i}({\mathbf{x}},t)=-\frac{\Delta
t}{\tau}[f_{i}({\mathbf{x}},t)-f_{i}^{eq}({\mathbf{x}},t)]$ (1)

where $f_{i}$ describes the discrete probability of finding a fluid particle
at position $\mathbf{x}$ and time $t$ moving with velocity $\mathbf{c}_{i}$
for $i=1,\ldots,19$. $\tau$ is the relaxation time or rate at which the system
relaxes toward a local equilibrium distribution function $f_{i}^{eq}$
corresponding to the truncated expansion of the Maxwell-Boltzmann distribution
for the velocities in an ideal gas. $f_{i}^{eq}$ is expressed as

$f_{i}^{eq}=\omega_{i}\rho\bigg{[}1+\frac{\mathbf{c}_{i}\cdot\mathbf{u}^{eq}}{c_{s}^{2}}-\frac{\left(\mathbf{u}^{eq}\cdot\mathbf{u}^{eq}\right)}{2c_{s}^{2}}+\frac{\left(\mathbf{c}_{i}\cdot\mathbf{u}^{eq}\right)^{2}}{2c_{s}^{4}}\bigg{]},$
(2)

where $\rho$ and $\mathbf{u}^{eq}$ are the macroscopic number density and the
equilibrium velocity, and $c_{s}=1/\sqrt{3}\Delta x/\Delta t$ is the lattice
speed of sound. $\omega_{i}$ are the lattice weights which read as $1/3$,
$1/18$ and $1/36$ for $i=1$, $i=2\dots 7$, and $i=8\dots 19$, respectively.
$\rho$ and $\mathbf{u}^{eq}$ are obtained from the moments of the particle
distribution function $f_{i}$ such that

$\rho=\sum_{i=1}^{19}f_{i}({\mathbf{x}},t),\quad\text{and}\hskip
7.11317pt\rho\mathbf{u}^{eq}=\sum_{i=1}^{19}f_{i}(\mathbf{x},t)\mathbf{c}_{i}+\tau\mathbf{F}$
(3)

where $\mathbf{F}$ is the total force applied by the particles on the fluid.
Note that since we are using the Shan and Chen scheme [44] to include external
forces on the fluid, the equilibrium velocity $(\mathbf{u}^{eq})$ used in Eq.
(2) is different from the physical velocity of the fluid nodes $(\mathbf{u})$
which is defined as

$\rho\mathbf{u}=\sum_{i=1}^{19}f_{i}({\mathbf{x}},t)\mathbf{c}_{i}+\frac{\Delta
t}{2}\mathbf{F}.$ (4)

Finally, the dynamic viscosity $\mu=\nu/\rho_{0}$ with $\nu$ being the
kinematic viscosity, reads as

$\mu=\rho_{0}c_{s}^{2}\left(\tau-\frac{\Delta t}{2}\right),$ (5)

where $\rho_{0}=m\rho$ is the mass density with $m$ here being the molecular
mass of the fluid particles. For convenience, $\Delta x$, $\Delta t$, and
$\rho_{0}$ are set to unity in this study.

###  2.2 Microswimmer model

We consider a force-driven microswimmer consisting of three beads of equal
size with two harmonic springs as depicted in Fig. 1. This is one of the
simplest possible microswimmer models [3]. A driving force $\textbf{f}^{d}$
exerted on each bead of the microswimmer is prescribed as [9]

$\displaystyle\textbf{f}_{\text{1}}^{d}(t)$ $\displaystyle=A\ \sin(\omega t)\
\hat{\mathbf{x}},$ $\displaystyle\textbf{f}_{\text{2}}^{d}(t)$
$\displaystyle=-\mathbf{f}_{\text{1}}^{d}(t)-\mathbf{f}_{\text{3}}^{d}(t)\text{
and }$ (6) $\displaystyle\textbf{f}_{\text{3}}^{d}(t)$ $\displaystyle=B\
\sin(\omega t+\alpha)\ \hat{\mathbf{x}}\text{ with }\alpha\in[-\pi,\pi].$

Here $A$ and $B$ are non-negative amplitudes of the time-dependent driving
forces $\mathbf{f}_{1}^{d}(t)$ and $\mathbf{f}_{3}^{d}(t)$ applied along the
$x$-axis to the outer beads at the frequency $\omega=2\pi/T$. A crucial
condition for self-propelled objects is that the sum over all external driving
forces equals zero at all times. For simplicity, we assume that the two
harmonic springs connecting the beads are identical with a stiffness $k$ and
an equilibrium length $L_{0}$ which implies that the total length of a single
swimmer is $\sim 2L_{0}$. The spring forces on the beads are given by

$\displaystyle\textbf{f}_{\text{1}}^{s}(t)$
$\displaystyle=-k(|\mathbf{R}_{1}(t)-\mathbf{R}_{2}(t)|-L_{0})\frac{\mathbf{R}_{1}(t)-\mathbf{R}_{2}(t)}{|\mathbf{R}_{1}(t)-\mathbf{R}_{2}(t)|},$
$\displaystyle\textbf{f}_{\text{2}}^{s}(t)$
$\displaystyle=-\mathbf{f}_{\text{1}}^{s}(t)-\mathbf{f}_{\text{3}}^{s}(t),$
(7) $\displaystyle\textbf{f}_{\text{3}}^{s}(t)$
$\displaystyle=-k(|\mathbf{R}_{3}(t)-\mathbf{R}_{2}(t)|-L_{0})\frac{\mathbf{R}_{3}(t)-\mathbf{R}_{2}(t)}{|\mathbf{R}_{3}(t)-\mathbf{R}_{3}(t)|},$

with $\mathbf{R}_{i}$ the position vector of bead $i$.

Our model is limited to small forces such that the assumptions for the Stokes
regime remain valid. For this purpose, it is critical to differentiate between
two Reynolds numbers

$\displaystyle
Re_{\text{swim}}=\frac{{\rho_{0}}\cdot|\textbf{v}_{\text{swim}}|\cdot 2\
L_{0}}{\mu}\quad\text{ and }\quad
Re_{\text{bead}}=\frac{{\rho_{0}}\cdot|\textbf{v}_{\text{max fastest
bead}}|\cdot 2\ r_{\text{bead}}}{\mu},$ (8)

and ensure that both $Re_{\text{swim}}$ and $Re_{\text{bead}}$ are small
enough.

Each spherical bead in our microswimmer is generated from an icosahedron that
is refined recursively until obtaining a sufficiently smooth surface described
by a triangular mesh [45]. Our fluid-filled particles are modeled using a
strain-hardening constitutive law known as the Skalak strain energy [46, 47],
which is written as

$\mathcal{W}_{sk}=\frac{\kappa_{s}}{4}\oint[I_{1}^{2}+2I_{1}-2I_{2}+CI_{2}^{2}]dA,$
(9)

where $I_{1}=\lambda_{1}^{2}+\lambda_{2}^{2}-2$ and
$I_{2}=\lambda_{1}^{2}\lambda_{2}^{2}-1$ are the deformation invariants, and
$C$ is a constant parameter controlling the extensibility of the membrane. The
area dilatation modulus $\kappa_{a}$ is defined such as
$\kappa_{a}/\kappa_{s}=1+2C$, with $\kappa_{s}$ being the shear elastic
modulus. The deformations are evaluated using a linear finite element method
[48]. In addition to resistances to shear elasticity and area dilatation, our
particle membrane can withstand out of plane deformations (i.e. bending). The
curvature energy is accounted for via the Helfrich free energy

$\mathcal{W}_{b}=\frac{\kappa_{b}}{2}\oint[2H]^{2}dA+{\kappa_{g}}\oint K\,dA,$
(10)

where $H=\frac{1}{2}\sum_{k=1}^{2}\bar{C}_{k}$, and
$K=\prod_{k=1}^{2}\bar{C}_{i}$ are the mean and Gaussian curvatures.
$\bar{C}_{1}$ and $\bar{C}_{2}$ are the two principal curvatures. $\kappa_{b}$
and $\kappa_{g}$ are the bending and Gaussian curvatures modulii. The
discretization of the bending energy follows the approach of Kantor and Nelson
[49] for flat triangulated meshes. The volume conservation of the capsule is
enforced using a penalty function reading as

$\mathcal{W}_{v}=\frac{\kappa_{v}}{2}\frac{[V-V_{0}]^{2}}{V_{0}},$ (11)

where $V_{0}$ is the reference volume of the stress-free particle, and
$\kappa_{v}$ is a constant parameter. The $i$-th particle nodal forces
$\mathbf{f}_{i}^{part}$ resulting from bending, shear elasticity and the
constraint on the volume are evaluated using the principle of virtual work
such that
$\mathbf{f}_{i}^{part}(\mathbf{r}_{i,j})=-\frac{\partial\mathcal{W}_{i}}{\partial\mathbf{r}_{i,j}}$
with $\mathcal{W}=\mathcal{W}_{sk}+\mathcal{W}_{b}+\mathcal{W}_{v}$. The
driving and spring forces are distributed over the mesh nodes of the
corresponding particle such that
$\hat{\mathbf{f}}_{i}^{\alpha}(\mathbf{r}_{i,j})\equiv\mathbf{f}_{i}^{\alpha}/N_{v}$,
with $N_{v}$ being the total number of mesh nodes on particle $i$, and the
index $i$ runs over the number of particles, here from $1$ to $3$. The
superscript $\alpha$ stands for either $d$ or $s$. $\mathbf{r}_{i,j}$ is the
position of the $j$-th mesh node belonging to the $i$-th particle. The total
nodal force ($\mathbf{f}_{i}^{tot}$) on particle $i$, which includes the
contribution of the external forces; is defined as

$\mathbf{f}_{i}^{tot}(\mathbf{r}_{i,j})=\mathbf{f}_{i}^{part}(\mathbf{r}_{i,j})+\hat{\mathbf{f}}_{i}^{d}(\mathbf{r}_{i,j})+\hat{\mathbf{f}}_{i}^{s}(\mathbf{r}_{i,j}).$
(12)

By choosing the appropriate values of $\kappa_{s}$, $C$, $\kappa_{b}$ and
$\kappa_{v}$ it is possible to work in the small deformation regime where the
particles are quasi-rigid.

###  2.3 Fluid-Particle interaction

The immersed boundary method (IBM) is a fluid-structure coupling method that
was first introduced by Peskin in the early seventies of the last century to
model the flow patterns around heart valves [50]. The IBM involves both
Eulerian and Lagrangian quantities. The Eulerian variables exist on a
Cartesian grid representing the fluid region while the Lagrangian variables
are based on a moving curvilinear mesh representing the interface. A smoothed
approximation of the Dirac delta function is used to transfer data from one
mesh to the other. The distribution of the particle nodal forces to the
neighboring fluid nodes reads as

$\mathbf{F}(\mathbf{x},t)=\sum_{j}\sum_{i}\mathbf{f}_{i,j}^{\textnormal{tot}}\Delta(\mathbf{x}-\mathbf{r}_{i,j}),$
(13)

where $\sum_{j}$ is a sum over all the membrane nodes $\mathbf{r}_{i,j}$
located within an interpolation range from the fluid node $\mathbf{x}$. At
this point, we reason in term of interactions between Lagrangian mesh nodes
(not in term of individual particles) and Eulerian fluid nodes. $\Delta$ is
the smoothed discrete Dirac delta function, and $\mathbf{F}(\mathbf{x},t)$ is
the force density acting on the fluid at the Eulerian node
$\mathbf{x}(x_{1},x_{2},x_{3})$ due to the contributions of the membrane total
force ($\mathbf{f}^{\textnormal{part}}$), the external driving force
($\hat{\mathbf{f}}^{\textnormal{d}}$) and the spring force
($\hat{\mathbf{f}}^{\textnormal{s}}$). Similarly, the interpolation of the
velocity of the neighboring fluid nodes onto a membrane Lagrangian node
$\mathbf{r}(r_{1},r_{2},r_{3})$ is performed as

$\dot{\mathbf{r}}_{i,j}=\sum_{x}\mathbf{u}(\mathbf{x})\Delta(\mathbf{x}-\mathbf{r}_{i,j}),$
(14)

where $\sum_{x}$ is a sum over all the fluid nodes within an interpolation
range from the membrane node $\mathbf{r}_{i,j}$.
$\Delta(\mathbf{x}-\mathbf{r}_{i,j})$ is then replaced by a two-point linear
interpolation function as detailed in [51, 52]. Eqs. (13) and (14) describe
the spreading of the interfacial forces to the surrounding fluid nodes and the
interpolation of the fluid velocity to the deformable interface. The particle
forces are included directly in the fluid node velocity using the method
proposed by Shan and Chen [44] as described in Eq. (3).

The IBM belongs to the class of the front-tracking methods where the sharp
interface is explicitly known through a set of Lagrangian marker points acting
on the fluid via body forces and moving with the same velocity as the ambient
fluid, thus enforcing the no-slip boundary condition. The deformation of the
interface is governed by the chosen strain-stress constitutive law and not by
the IBM itself, making it a very popular method to simulate biological fluid
dynamics (e.g., blood flow) [53, 54, 55, 56, 57].

###  2.4 Analytical calculation

To compare the results of the lattice Boltzmann method, we make use of an
analytical framework [39] where the hydrodynamic interactions between the
particles are calculated using either the Oseen or Rotne-Prager approximation
[58]. The equation of motion of the bead-spring system reads

$\dot{\mathbf{R}}_{i}(t)=\frac{1}{6\pi\mu
r_{\text{bead}}}(\mathbf{f}^{d}_{i}(t)+\mathbf{f}^{s}_{i}(t))+\sum_{j\neq
i}\hat{T}(\mathbf{R}_{i}(t)-\mathbf{R}_{j}(t))\cdot(\mathbf{f}^{d}_{j}(t)+\mathbf{f}^{s}_{j}(t)),$
(15)

with particle indices $i,j$, the position $\mathbf{R}_{i}$ of particle $i$,
the bead radius $r_{\text{bead}}$ and the spring force $\mathbf{f}^{s}_{i}$ on
bead $i$. The tensor $\hat{T}$ is given by either the Oseen tensor,

$\hat{T}(\mathbf{r})=\frac{1}{8\pi\mu|\mathbf{r}|}\left(\hat{1}+\frac{\mathbf{r}\otimes\mathbf{r}}{\mathbf{r}^{2}}\right),$
(16)

or the Rotne-Prager tensor [30]

$\hat{T}(\mathbf{r})=\frac{1}{8\pi\mu|\mathbf{r}|}\left(\hat{1}+\frac{\mathbf{r}\otimes\mathbf{r}}{\mathbf{r}^{2}}\right)+\frac{r_{\text{bead}}^{2}}{12\pi\mu|\mathbf{r}|^{3}}\left(\hat{1}-3\frac{\mathbf{r}\otimes\mathbf{r}}{\mathbf{r}^{2}}\right).$
(17)

Here, $\hat{1}$ denotes the $3\times 3$ unit matrix and $\otimes$ the tensor
product. The equation of motion (15) is then solved numerically, and the
swimmer velocities extracted from the particle trajectories are readily
compared to the velocities obtained in the LBM simulations.

##  3 Results

We consider a fully periodic simulation domain of size $1024\times 256\times
256$. Each microswimmer is formed by three beads of radius $r_{\text{bead}}=5$
and connected with two springs with an equilibrium length of $L_{0}=36$ and a
spring constant $k=2\cdot 10^{-2}$. The amplitudes of the time-dependent
driving forces $\mathbf{f}_{1}^{d}(t)$ and $\mathbf{f}_{3}^{d}(t)$ are chosen
as $A=10^{-1}$ and $B=5\cdot 10^{-3}$, respectively. The period and the phase
shift are set to $T=2\cdot 10^{4}$ and $\alpha=\pi/2$. Since the size of the
beads is quite smaller than the total length of the microswimmer, then, as
mentioned earlier, we can approximate the total length of the swimmer by the
sum of the lengths of the two springs at equilibrium $\sim 2L_{0}$. We choose
$\kappa_{s}=5\cdot 10^{-2}$, $\kappa_{a}=1$, $\kappa_{b}=1.25\cdot 10^{-3}$
and $\kappa_{v}=1$ such that the particles remain spherical during the swimmer
motion.

We calculate the swimmer velocity by the average position of the center of
mass (cm)

$\displaystyle\textbf{v}_{\text{cm}}=\frac{\textbf{x}_{\text{cm}}(t_{0}+T)-\textbf{x}_{\text{cm}}(t_{0})}{T},\quad\text{
with }\textbf{x}_{\text{cm}}=\frac{1}{3}\sum_{\text{i =
1}}^{3}\textbf{x}_{\text{i}}.$ (18)

###  3.1 Single microswimmer

Figure 2: Left: Velocity $\textbf{v}_{\text{ref}}$ of a single microswimmer
(reference swimmer) for different external driving force periods $T$ (see
legend). Points stand for LBM data, solid (dashed) lines for analytical
predictions based on the Oseen [59] (Rotne-Prager [10]) tensor. Right:
relative differences
$\Delta\textbf{v}_{\text{ref}}=|\textbf{v}_{\text{ref}}-\textbf{v}_{\text{theo}}|/\textbf{v}_{\text{theo}}$
between numerical data and analytical, Oseen (open points) and Rotne-Prager
(solid points) predictions. The red lines report the magnitude of the Reynolds
number for the fastest bead (solid) and the swimmer in total (dashed) in the
case of $T=20000$ [l.u.]. The figure shows a better match between LBM data and
Rotne-Pragner tensor than LBM data and Oseen tensor. The Reynolds number of
the fastest bead plays an essential role in comparing. The higher external
driving force period implies a decreased Reynolds number, so the tensor
approach fits better to our simulation data. Figure 3: Averaged fluid velocity
field for a single microswimmer for dynamic viscosity $\mu=0.1\bar{6}$. The
fluid flow (black arrows) surrounding a single microswimmer (black beads) is
shown whereby it is averaged over one swimming period in the stationary case.
The absolute values of the fluid velocity u (black arrows) are shown by the
different colors.

First, we characterize the performance of a single microswimmer. Long-time
simulations over 50 swimming periods $T$ show that the stationary swimming
velocity is reached after $\sim 12$ swimming periods. Accordingly, the average
velocities are defined via Eq. (18) after such a threshold. In particular, we
focus on the dependence of the swimmer’s velocity on the dynamic viscosity of
the fluid $\mu$.

As expected, we retrieve a non-monotonous velocity dependence of the swimmer
on $\mu$, as reported in Ref. [59], for diverse values of the period of the
forcing $T$. In particular, the velocity of a single swimmer obtained from the
LB simulations compares well to analytical results (Fig. 2). To assess the
origin of the quantitative mismatch between our simulation and analytical
results using the Oseen tensor, we look at the value of the Reynolds number
associated with the swimmer as a whole and with a single bead at its maximum
speed. Indeed, the Oseen approach of Ref. [59] assumes vanishing Reynolds
number. In contrast, the solid and dotted lines of Fig. 2 show that while the
Reynolds number of the swimmer is relatively small ($Re_{swimm}\sim 10^{-4}$),
the Reynolds number associated with the beads at maximum speed is larger by
two orders of magnitude ($Re_{bead}\sim 10^{-2}$). Indeed, recent works have
shown that when the Reynolds number approaches unity, swimming protocols that
exploit inertia can become operational [33, 60]. However, the relative
mismatch between the Oseen predictions and the numerical results suggests that
fluid inertia may not be the only cause of such a discrepancy. An additional
assumption of the analytical model is that the fluid velocity adjusts
instantaneously to the force, whereas in the LBM simulations, such a time
characterized by $\bar{\tau}=L_{0}^{2}/\nu$ is finite and it is represented by
the relaxation time $\tau$ in the BGK collision operator [41]. On the other
hand, also the analytical predictions are derived under approximations. In
fact, the Oseen approach assumes point particles and it does not properly
account for the spherical shape of the particles. Accordingly, we compare our
data also with a more refined approach based on the Rotne-Prager tensor [10].
Interestingly, Fig. 2 shows a better agreement between the LBM simulations and
the predictions based on the Rotne-Prager as with the Oseen model [59].

Finally, the analytical-numerical mismatch can be due to the periodic boundary
conditions that we implement in the LBM simulations compared to the unbound
fluid considered in the analytical models. The periodic boundary conditions
induce spurious interactions of the swimmer with its image. Accordingly, to
assess this issue, we have checked that the self-interaction across the
boundaries is lower than one percent, ruling it out as a possible reason for
the mismatch with the analytical results. This crucial point is also valid for
the case of two microswimmers considered below.

Fig. 3 shows the velocity field, averaged over a period of the external forces
acting on the swimmer. The velocity profile of this swimmer resembles that of
a force dipole in a far-field expansion.

###  3.2 Two aligned microswimmers

Figure 4: Non-scaled sketch of the in-row configuration for two swimmers. The
swimmers move along the x-axis to the right. The leading swimmer swims in
front of the trailing one. The initial distance between the swimmers is
measured as a difference of the center beads within the swimmers
$(p_{2}-p_{1})$. Figure 5: The velocity increase/decrease
$\textbf{v}_{\Delta}=(\textbf{v}-\textbf{v}_{\text{ref}})/\textbf{v}_{\text{ref}}$
of two aligned microswimmers for different configurations. The velocities are
measured during the 12th swimming cycle (stationary case) with simulation
parameters given in Sec. 2. The collective motion of two swimmers increases
the velocities of both swimmers in certain configurations. Figure 6: Averaged
fluid velocity field for two aligned microswimmers with dynamic viscosity
$\mu=0.1\bar{6}$. The fluid flow (black arrows) surrounding the microswimmers
(black beads) is shown whereby it is averaged over one swimming period in the
stationary case. The absolute values of the fluid velocity u (black arrows)
are depicted by the different colors.

Next, we analyze the case of two identical aligned microswimmers, as sketched
in Fig. 4. In this case the crucial observable is the relative velocity
increase/decrease

$\displaystyle\displaystyle\textbf{v}_{\Delta}=\frac{\textbf{v}-\textbf{v}_{\text{ref}}}{\textbf{v}_{\text{ref}}}\,,$
(19)

in comparison to the velocity $\textbf{v}_{\text{ref}}$ of a single
microswimmer.

Fig. 5 shows $\textbf{v}_{\Delta}$ as a function of the dynamic viscosity
$\mu$ for different distances among the centers of mass of the swimmers. As in
the single swimmer case, we observe a dependency of the velocity on the
fluid’s dynamic viscosity. Interestingly, for a distance of $4L_{0}$, both
swimmers move faster than a single swimmer, with a velocity increase of up to
$36\%$ for the leading swimmer (see Fig. 5). For larger initial distances
between the swimmers, the effect decreases. Since we observe that the leading
swimmer moves faster than the trailing one, such configurations are unstable,
i.e., the trailing one will be left behind until the distance between their
centers of mass reaches $\gtrsim 216$ lattice units. Since the body length of
the swimmers amounts to $L_{swim}=72$ lattice units, it implies that the two
swimmers effectively do not affect each other when they are $216/72\sim 3$
body length of the swimmers away from each other (see Fig. 5). Finally, Fig. 6
shows the fluid velocity field, which is similar to that of a single swimmer
shown in Fig. 3.

###  3.3 Two parallel microswimmers

Figure 7: Non-scaled schematics of two microswimmers next to each other.
Considering the $\hat{\textbf{y}}$ direction one swimmer swims on the top of
the other one. Both swimmers are centered in $\hat{\textbf{y}}$,
$\hat{\textbf{z}}$ and swim in $\hat{\textbf{x}}$ direction. The distance
between both is measured from the center of mass ($p_{1}$) of the center bead
to the center of mass of the middle bead ($p_{2}$) in $\hat{\textbf{z}}$
direction. Figure 8: The velocity increase/decrease
$\textbf{v}_{\Delta}=(\textbf{v}-\textbf{v}_{\text{ref}})/\textbf{v}_{\text{ref}}$
of two microswimmer next to each other for different configurations. The
velocities are measured inside the stationary swimming cases, and collective
motion implies different velocity increase/decrease effects. The transition
from decreasing to increasing velocities for larger distances implies an
optimal distance for maximal velocity increment. Figure 9: Averaged fluid
velocity field for two microswimmers next to each other with dynamic viscosity
$\mu=0.1\bar{6}$. The fluid flow (black arrows) surrounding the microswimmers
(black beads) is shown whereby it is averaged over one swimming period in the
stationary case. The absolute values of the fluid velocity u (black arrows) is
shown by the different colors.

Finally, we discuss the case of two parallel identical swimmers as shown in
Fig. 7. As expected by symmetry, we find that the two swimmers always have the
same velocity for all the distances we investigated among them, i.e., all
investigated configurations are stationary on the time scales accessible to
lattice Boltzmann simulations. However, theoretical work [61, 39] has shown
that parallel swimmers typically also experience sideways interaction as well
as rotation, which will in general break stationarity.

Moreover, similarly to the case of aligned swimmers, we observe a dependence
of the relative velocity change on the dynamic viscosity of the fluid, as
shown in Fig. 8. Indeed, there is a transition in Fig. 8 from decreased
velocity to increased velocity between an initial distance of one arm length
$L_{0}$ and $3\,L_{0}$. This transition implies that swimmers that are close
to each other mutually disturb, while swimmers that are further away benefit
from the presence of a nearby (yet not too close) companion. Finally, Fig. 9
represents the averaged fluid velocity field.

##  4 Discussion

Comparison of lattice Boltzmann simulations with theoretical approaches based
on the Oseen or Rotne-Pragner tensor is generally tricky. The theory holds in
the Stokes limit, in which the fluid velocity profile adjusts instantaneously
to the change of position of the beads. However, in our IBM+LBM+FEM
simulations, momentum propagation over a length $L$ occurs on a finite time
determined by $\bar{\tau}=L^{2}/\nu$ with $\nu=\mu/\rho_{0}$ the kinematic
viscosity. Therefore, the velocity field of one microswimmer needs a finite
amount of time to attain its stationary profile and to affect the other beads’
motion. In our simulations we have $L\simeq 36$ and $\nu\in[0.1,1.2]$ (with
$\rho_{0}=1$) and hence we have $\bar{\tau}\in[1000,12000]$ timesteps.
Clearly, for larger values of the viscosity $\bar{\tau}\simeq 1000$, it is
reasonable to have a time scale separation between $\bar{\tau}$ and the period
of the force $T=20000$. In contrast, this is no longer the case for smaller
viscosity values, and hence larger deviations between numerical and analytical
results are to be expected. We have numerically solved the equation of motion
for the swimmers using the Rotne-Prager tensor where the delay time has been
introduced as a fitting parameter. Fig. 10 shows the dependence of the
velocity enhancement of the leading and trailing swimmer as a function of the
viscosity. In order to discuss the magnitude of the fitted delay time,
$\tau_{fit}$, we have compared it with $\tau=L_{0}^{2}/\nu$.

Figure 10: Left: leading (triangle) and trailing (circle) swimmers relative
velocity, $v_{\Delta}$, for diverse proposed delay times while the black
points are the relative velocities of the leading (smooth) and trailing
(dashed) swimmer as obtained from LBM simulations. Right: delay time as
obtained by fitting the numerical integration of Eq. (20) to the data from the
LBM simulations in units of the diffusion time of the momentum,
$\bar{\tau}=L^{2}/\nu$.

Some preliminary results for the finite relaxation time of the velocity field
show an improvement in the analytical/numerical comparison. Here, the equation
of motion in the analytical calculations has been adapted to

$\dot{\mathbf{R}}_{i}(t)=\frac{1}{6\pi\mu
r_{\text{bead}}}(\mathbf{f}^{d}_{i}(t)+\mathbf{f}^{s}_{i}(t))+\sum_{j\neq
i}\hat{T}(\mathbf{R}_{i}(t)-\mathbf{R}_{j}(t-\tau_{fit}))\cdot(\mathbf{f}^{d}_{j}(t-\tau_{fit})+\mathbf{f}^{s}_{j}(t-\tau_{fit})),$
(20)

with $\tau_{fit}$ being the time delay. The left panel of Fig. 10 shows that a
good agreement between the LBM and the analytical model can be obtained by
properly tuning the fitted delay time $\tau_{fit}$. The right panel of Fig. 10
shows that the fitted delay time, $\tau_{fit}$, grows with the viscosity,
$\nu$ and it is comparable to the time, $\tau$, that the momentum takes to
propagate over half of the length of the swimmer.

##  5 Conclusions

We reported on LBM simulations of one and two bead-spring swimmers. As
expected, for the single swimmer, we found that its velocity is sensitive to
the dynamic viscosity of the fluid, $\mu$, and it shows a non-monotonous
dependence on $\mu$ as shown in Fig. 2. Our numerical results show a good
quantitative agreement with the theoretical predictions based on both the
Oseen (Ref. [9, 59]) or Rotne-Pragner (Ref. [10, 39]) tensor as shown in Fig.
2.

Next, we analyzed the collective motion of two microswimmers in the case in
which they are aligned or alongside. In the case of aligned swimmers, we found
that when the swimmers are close, the speedup due to hydrodynamic interactions
of the leading swimmer is more significant than that of the trailing one hence
leading to unstable clusters. At variance, for the case of alongside swimmers,
we found that all the distances among their centers of mass are stationary so
that swimmers can proceed together. In the case of alongside swimmers, we
found that when swimmers are too close, they slow down, whereas they speed up
at larger (yet finite) distances. Since when the separation distance diverges
we expect no speed up, our results suggest the existence of an optimal length
at which the speed up is maximized.

Finally, we have commented on the mismatch of the velocity profile as obtained
from the LBM simulations and the analytical calculations. Supported by some
preliminary results, we argue that this is due to the finite relaxation time
of the fluid velocity in the LBM as compared to the Stokes regime assumed in
the analytical calculations.

## Acknowledgments

The DFG Priority Programme SPP 1726 „Microswimmers—From Single Particle Motion
to Collective Behaviour“ (HA 4382/5-1) and SFB 1411 (Project-ID 416229255)
financially supported this work. We also acknowledge the Jülich Supercomputing
Centre (JSC) for the allocation of computing time.

## References

  * [1] E. M. Purcel, Life at low Reynolds number, American Journal Of Physics, 45 (1977), 3-11.
  * [2] J. R. Blake, A spherical envelope approach to cilary propulsion, Journal Of Fluid Mechanics, 46 (1971), 199-208.
  * [3] B. U. Felderhof, The swimming of animalcules, Physics Of Fluids, 18 (2006), 063101.
  * [4] D. J. Earl and C. M. Pooley and J. F. Ryder and I. Bredberg and J. M. Yeoman, Modeling microscopic swimmers at low Reynolds number, Journal Of Chemical Physics, 126 (2007), 064703.
  * [5] A. Najafi and R. Golestanian, Simple swimmer at low Reynolds number: Three linked spheres, Physical Review E, 69 (2004), 6.
  * [6] R. Golestanian and A. Ajdari, Analytic results for the three-sphere swimmer at low Reynolds number, Physical Review E, 77 (2008), 036308.
  * [7] J. E. Avron and O. Kenneth and D. H. Oakni, Pushmepullyou: An efficient micro-swimmer, New Journal Of Physics, 7 (2005), 234-234.
  * [8] M. T. Downton and H. Stark, Simulation of a model microswimmer, Journal Of Physics: Condensed Matter, 21 (2009), 204101.
  * [9] J. Pande and A.-S. Smith, Forces and shapes as determinants of micro-swimming: effect on synchronisation and the utilisation of drag, Soft Matter, 11 (2015), 2364-2371.
  * [10] S. Ziegler and M. Hubert and N. Vandewalle and J. Harting, and A.-S. Smith, A general perturbative approach for bead-based microswimmers reveals rich self-propulsion phenomena, New Journal Of Physics, (2019,11) 21 (2019), 113017.
  * [11] Q. Wang, Optimal strokes of low Reynolds number linked-sphere swimmers, Applied Sciences (Switzerland), 9 (2019), 4023.
  * [12] A. Daddi-Moussa-Ider and M. Lisicki and C. Hoell and H. Löwen Swimming trajectories of a three-sphere microswimmer near a wall, Journal Of Chemical Physics, 148 (2018), 134904.
  * [13] A. Daddi-Moussa-Ider and C. Kurzthaler and C. Hoell and A. Zött and M. Mirzakhanloo and M. R.Alam and A. M. Menzel and H. Löwen and S. Gekle, Frequency-dependent higher-order Stokes singularities near a planar elastic boundary: Implications for the hydrodynamics of an active microswimmer near an elastic interface, Physical Review E, 100 (2019), 032610.
  * [14] M. S. Rizvi and A. Farutin and C. Misba, Three-bead steering microswimmers, Physical Review E, 97 (2018), 023102.
  * [15] R. Dreyfus and J. Baudry and M. L. Roper and M. Fermigier and H. A. Stone and J. Bibett, Microscopic artificial swimmers, Nature, 437 (2005), 862-865.
  * [16] D. Ahmed and M. Lu and A. Nourhani and P. E. Lammert and Z. Stratton and H. S. Muddana and V. H.Crespi and T. J. Huang, Selectively manipulable acoustic-powered microswimmers, Scientific Reports, 5 (2015).
  * [17] G. Grosjean and G. Lagubeau and A. Darras and M. Hubert and G. Lumay and N. Vandewalle, Remote control of self-assembled microswimmers, Scientific Reports, 5 (2015), 16035.
  * [18] G. Grosjean and M. Hubert and G. Lagubeau and N. Vandewalle, Realization of the Najafi-Golestanian microswimmer, Physical Review E, 94 (2016), 021101.
  * [19] S. Palagi and A. G. Mark and S. Y. Reigh and K. Melde and T. Qiu and H. Zeng and C. Parmeggiani and D.Martella and A. Sanchez-Castillo and N. Kapernaum and F. Giesselmann and D. S. Wiersma and E.Lauga and P. Fische, Nature Materials, 15 (2016), 647–653.
  * [20] G. Grosjean and M. Hubert and Y. Collard and S. Pillitteri and N. Vandewalle, Surface swimmers harnessing the interface to self-propel, European Physical Journal E, 41 (2018), 137.
  * [21] J. Zheng and B. Dai and J. Wang and Z. Xiong and Y. Yang and J. Liu and X. Zhan and Z. Wan and J. Tang, Orthogonal navigation of multiple visible-light-driven artificial microswimmers, Nature Communications, 8 (2017), 1438.
  * [22] J. K. Hamilton and P. G. Petrov and C. P. Winlove and A. D. Gilbert and M. T. Bryan and F. Y.Ogrin, Magnetically controlled ferromagnetic swimmers, Scientific Reports, 7 (2017), 44142.
  * [23] M. T. Bryan and S. R. Shelley and M. J. Parish and P. G. Petrov and C. P. Winlove and A. D. Gilbert and F. Y. Ogrin, Emergent propagation modes of ferromagnetic swimmers in constrained geometries, Journal Of Applied Physics, 121 (2017), 073901.
  * [24] J. K. Hamilton and A. D. Gilbert and P. G. Petrov and F. Y. Ogrin, Torque driven ferromagnetic swimmers, Physics Of Fluids, 30 (2018), 092001.
  * [25] Y. Collard and G. Grosjean and N. Vandewalle, Magnetically powered metachronal waves induce locomotion in self-assemblies, Communications Physics, 3 (2020), 112.
  * [26] Gerhard Gompper and Roland G Winkler and Thomas Speck and Alexandre Solon and Cesare Nardini and Fernando Peruani and Hartmut Löwen and Ramin Golestanian and U Benjamin Kaupp and Luis Alvarez and Thomas Kiørboe and Eric Lauga and Wilson C K Poon and Antonio DeSimone and Santiago Muiños-Landin and Alexander Fischer and Nicola A Söker and Frank Cichos and Raymond Kapral and Pierre Gaspard and Marisol Ripoll and Francesc Sagues and Amin Doostmohammadi and Julia M Yeomans and Igor S Aranson and Clemens Bechinger and Holger Stark and Charlotte K Hemelrijk and François J Nedelec and Trinish Sarkar and Thibault Aryaksama and Mathilde Lacroix and Guillaume Duclos and Victor Yashunsky and Pascal Silberzan and Marino Arroyo and Sohan Kale, Journal Of Physics: Condensed Matter, 32 (2020), 193001.
  * [27] L. E. Becker and S. A. Koehler and H. A. Stone, On self-propulsion of micro-machines at low Reynolds number: Purcell’s three-link swimmer, Journal Of Fluid Mechanics, 490 (2003), 15-35.
  * [28] A. Zöttl and H. Star, Nonlinear dynamics of a microswimmer in Poiseuille flow, Physical Review Letters, 108 (2012), 218104.
  * [29] K. Pickl and J. Götz and K. Iglberger and J. Pande and K. Mecke and A.-S. Smith and U. Rüde, All good things come in threes-Three beads learn to swim with lattice Boltzmann and a rigid body solver, Journal Of Computational Science, 3 (2012), 374-387.
  * [30] K. Pickl and J. Pande and H. Köstler, U. Rüde, and A.-S. Smith, Lattice Boltzmann simulations of the bead-spring microswimmer with a responsive stroke - From an individual to swarms, Journal Of Physics Condensed Matter, 29 (2017), 124001.
  * [31] A. Sukhov and S. Ziegler and Q. Xie and O. Trosman and J. Pande and G. Grosjean and M. Hubert and N.Vandewalle and A.-S. Smith and J. Harting, Optimal motion of triangular magnetocapillary swimmers, The Journal Of Chemical Physics, 151 (2019), 124707.
  * [32] T. Peter and P. Malgaretti and N. Rivas and A. Scagliarini and J. Harting and S. Dietrich, Numerical simulations of self-diffusiophoretic colloids at fluid interfaces, Soft Matter, 16 (2020), 3536-3547.
  * [33] M. Hubert amd O. Trosman and Y. Collard and A. Sukhov and J. Harting and N. Vandewalle and A.-S.Smith, Scallop Theorem and Swimming at the Mesoscale, Physical Review Letters, 126 (2021), 224501.
  * [34] T. Qiu and T.-C. Lee and A. G. Mark and K. I. Morozov and R. Münster and O. Mierka and S. Turek and A. M.Leshansky and P. Fischer, Swimming by reciprocal motion at low Reynolds number, Nature Communications, 5 (2014), 5119.
  * [35] J. Elgeti and R. Winkler and G. Gompper, Physics of microswimmers–single particle motion and collective behavior: A review, Reports On Progress In Physics, 78 (2015), 056601.
  * [36] P. Malgaretti and H. Stark, Model microswimmers in channels with varying cross section, Journal Of Chemical Physics, 146 (2017), 174901.
  * [37] P. Malgaretti and J. Harting, Phoretic colloids close to and trapped at fluid interfaces, ChemNanoMat, 7 (2021), 1073-1081.
  * [38] J. Buzhardt and P. Tallapragad, Dynamics of groups of magnetically driven artificial microswimmers, Physical Review E, 100 (2019), 033106.
  * [39] S. Ziegler and T. Scheel and M. Hubert and J. Harting and A.-S. Smith, Theoretical framework for two-microswimmer hydrodynamic interactions, New Journal Of Physics, 23 (2021,7), 073041.
  * [40] C. Peskin, The immersed boundary method, Acta Numerica, 11 (2002), 479-517.
  * [41] P. L. Bhatnagar and E. P. Gross and M. Krook, A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems, Physical Review, 94 (1954), 511-525.
  * [42] R. Benzi and S. Succi and M. Vergasso, The lattice Boltzmann equation: theory and applications, Physics Reports, 222 (1992), 145-197.
  * [43] S. Succi, The lattice Boltzmann equation for fluid dynamics and beyond, Clarendon Press Oxford University Press, 2001.
  * [44] X. Shan and H. Chen, Lattice Boltzmann model for simulating flows with multiple phases and components, Physical Review E, 47 (1993), 1815-1819.
  * [45] O. Aouane and Q. Xie and A. Scagliarini and J. Harting, Mesoscale simulations of Janus particles and deformable capsules in flow, High Performance Computing In Science And Engineering’17 (2018), 369-385.
  * [46] R. Skalak and A. Tozeren and R. P. Zarda and S. Chien, Strain Energy Function of Red Blood Cell Membranes, Biophysical Journal, 13 (1973), 245-264.
  * [47] R. Skalak, Modelling the mechanical behavior of red blood cells, Biorheology, 10 (1973), 229-238.
  * [48] T. Krüger and F. Varnik and D. Raabe, Efficient and accurate simulations of deformable particles immersed in a fluid using a combined immersed boundary lattice Boltzmann finite element method, Computers & Mathematics With Applications, 61 (2011), 3485-3505.
  * [49] Y. Kantor and D. R. Nelson, Crumpling transition in polymerized membranes, Physical Review Letters, 58 (1987), 2774-2777.
  * [50] C. Peskin, Flow patterns around heart valves: a digital computer method for solving the equations of motion, IEEE Transactions On Biomedical Engineering, (1973), 316-317.
  * [51] T. Krüger, Computer simulation study of collective phenomena in dense suspensions of red blood cells under shear, Springer Science & Business Media, 2012.
  * [52] O. Aouane and A. Scagliarini and J. Harting, Structure and rheology of suspensions of spherical strain-hardening capsules, Journal Of Fluid Mechanics, 911 (2021), A11.
  * [53] J. Zhang and P. C. Johnson and A. S. Popel, An immersed boundary lattice Boltzmann approach to simulate deformable liquid capsules and its application to microscopic blood flows, Physical Biology, 4 (2007), 285-295.
  * [54] L. M. Crowl and A. L. Fogelson, Computational model of whole blood exhibiting lateral platelet motion induced by red blood cells, International Journal For Numerical Methods In Biomedical Engineering, 26 (2010), 471-487.
  * [55] B. Kaoui and J. Harting and C. Misbah, Two-dimensional vesicle dynamics under shear flow: Effect of confinement, Physical Review E, 83 (2011), 066319.
  * [56] M. Thiébaud and Z. Shen and J. Harting and C. Misbah, Prediction of anomalous blood viscosity in confined shear flow, Physical Review Letters, 112 (2014), 238304.
  * [57] T. Krüger and D. Holmes and P. V. Coveney, Deformability-based red blood cell separation in deterministic lateral displacement devices—A simulation study, Biomicrofluidics, 8 (2014), 054114.
  * [58] J. Dhont, An Introduction to Dynamics of Colloids, Elsevier, (1996), 1-67.
  * [59] J. Pande and L. Merchant and T. Krüger and J. Harting and A.-S. Smith, Setting the pace of microswimmers: when increasing viscosity speeds up self-propulsion, New Journal Of Physics, 19 (2017), 053024.
  * [60] A. Choudhary and S. Paul and F. Rühle and H. Star, How inertial lift affects the dynamics of a microswimmer in Poiseuille flow, Communications Physics, 5 (2022), 14.
  * [61] C. M. Pooley and G. P. Alexander and J. M. Yeoman, Hydrodynamic interaction between two swimmers at low Reynolds number, Physical Review Letters, 99 (2007), 228103.