February 2021
EPL, 133 (2021) 34001
doi: 10.1209/0295-5075/133/34001

www.epljournal.org

Autonomous elastic microswimmer
Katsutomo Era1 , Yuki Koyano2 , Yuto Hosaka1 , Kento Yasuda1 , Hiroyuki Kitahata3
and Shigeyuki Komura1
1

Department of Chemistry, Graduate School of Science, Tokyo Metropolitan University - Tokyo 192-0397, Japan
Department of Physics, Graduate School of Science, Tohoku University - Sendai 980-8578, Japan3
Department of Physics, Graduate School of Science, Chiba University - Chiba 263-8522, Japan2

received 6 October 2020; accepted in final form 15 December 2020
published online 30 March 2021
47.63.Gd – Swimming microorganisms
47.63.mf – Low-Reynolds-number motions
PACS 05.45.Xt – Synchronization; coupled oscillators
PACS

PACS

Abstract – A model of an autonomous three-sphere microswimmer is proposed by implementing
a coupling effect between the two natural lengths of an elastic microswimmer. Such a couplingmechanism is motivated by the previous models for synchronization phenomena in coupled oscillator systems. We numerically show that a microswimmer can acquire a nonzero steady state
velocity and a finite phase difference between the oscillations in the natural lengths. These velocity
and phase differences are almost independent of the initial phase difference. There is a finite rangeof the coupling parameter for which a microswimmer can have an autonomous directed motion.
The stability of the phase difference is investigated both numerically and analytically in order todetermine its bifurcation structure.
c 2021 EPLA
Copyright !

Introduction. – Microswimmers are small machinesthat swim in a fluid and they are expected to be usedin microfluidics and microsystems [1]. Over the lengthscale of microswimmers, the fluid forces acting on them
are dominated by the frictional viscous forces. By transforming chemical energy into mechanical energy, however,microswimmers change their shape and move efficiently
in viscous environments. According to Purcell’s scallop
theorem, reciprocal body motion cannot be used for locomotion in a Newtonian fluid [2,3]. As one of the simplest
models exhibiting nonreciprocal body motion, Najafi andGolestanian proposed a three-sphere swimmer (NG swimmer) [4,5], in which three in-line spheres are linked by twoarms of varying length. In recent years, such a swimmer
has been experimentally realized by using colloidal beads
manipulated by optical tweezers [6], or ferromagnetic particles at an air-water interface [7,8].
Recently, some of the present authors have proposeda generalized three-sphere microswimmer model in which
the spheres are connected by two harmonic springs, i.e., anelastic microswimmer [9]. Compared with the NG swimmer, the main difference is that the natural length of eachspring (rather than the arm length) is assumed to undergo
a prescribed cyclic motion. A similar model was proposedby other people [10–12]. We have analytically obtained theaverage swimming velocity as a function of the frequency

of cyclic change in the natural length [9]. Using this model,
we have also discussed the hydrodynamic interaction between two elastic swimmers [13] and a thermally driven
elastic microswimmer [14–16].
In the above three-sphere microswimmer models, either
the arm lengths (NG swimmer) or the natural lengths ofthe springs (elastic swimmer) are assumed to undergo aprescribed cyclic motion. Such active motions can leadto a net locomotion if the swimming strokes are nonreciprocal. In these models, the average swimming velocity ispurely determined by the frequency and the phase difference of the prescribed motions [4,5,9,13].
On the other hand, it is beneficial for a microswimmer
if the swimming velocity is autonomously determined byitself rather than being imposed externally. Moreover, asophisticated microswimmer requires a feedback control
system in order to regulate the switching between thestatic and swimming states by tuning the system parameters. For a macroscopic quadruped robot (not a swimmer),it was demonstrated that the communication between legs
during movements is essential for interlimb coordinationin quadruped walking [17,18]. A similar mechanism is also
useful for the locomotion of a microswimmer.
In this letter, extending the mechanism of an elastic swimmer [9], we propose a new type of three-sphereswimmer which can autonomously determine its velocity.

34001-p1

Katsutomo Era et al.
the equations of motion of each sphere as [9,13]
ẋi =

3
!

Mij fj ,(2)

j=1

where the three forces fi are given by
Fig. 1: An autonomous elastic microswimmer in a viscous fluid
characterized by the shear viscosity η. Three identical spheresof radius a are connected by two harmonic springs characterized by the spring constant K. The time-dependent positionsof the spheres are denoted by x1 , x2 , and x3 which evolve intime according to eq. (2). The time-dependent natural lengthsof the springs are denoted by "A and "B whose dynamics is
described by eqs. (5) and (6), respectively, whereas the corresponding phases θA and θB obey eqs. (7) and (8), respectively.

fi = −

∂E
.
∂xi

(3)

Here the details of the hydrodynamic interactions aretaken into account through the mobility coefficients Mij .
Within Oseen’s approximation, which is justified when
the spheres are considerably far from each other (a "|xi − xj |), the expressions for the mobility coefficients Mij
can be written as
1
i = j,
 6πηa ,In order to implement such a control mechanism, we in(4)Mij =troduce a coupling between the two natural lengths of an
1



,
i
=
#
j.
elastic microswimmer by using the interaction adopted in
4πη|xi − xj |the Kuramoto model for coupled oscillators [19–22]. ImThe force-free condition, f1 + f2 + f3 = 0, is automaticallyportantly, the proposed microswimmer acquires a steadysatisfiedin the present model [9,13]. We define the centerstate velocity and a finite phase difference in the long-timeof-mass
positionof a microswimmer by X(t) = (x1 + x2 +limit without any external control. The steady state velocx3 )/3 and the swimming velocity of the whole object by
ity can be mainly tuned by changing the coupling parameter in the model. Moreover, we investigate the condition V (t) = Ẋ(t).
Next, we consider that the two natural lengths of thethat a microswimmer can attain an autonomous locomosprings
undergo the following cyclic changes in time [9,13]:tion, and further perform a linear stability analysis of thesteady state.
!A (t) = ! + d cos θA (t),(5)
Synchronization phenomena are widely observed in ac!B (t) = ! + d cos θB (t),(6)
tive biological systems such as flagella and cilia [23,24].
In particular, synchronization of a pair of flagella in where ! is the constant natural length, d is the oscillation
Chlamydomonas was observed experimentally [25,26]. For amplitude, θA (t) and θB (t) are the time-dependent phases.
a three-sphere model of Chlamydomonas in which the The most important aspect of our model is that θA (t)spheres representing the flagella move on circular trajec- and θB (t) are affected by the relative positions and thetories relative to the body sphere [27,28], the two flagella velocities of the three spheres. We employ the following
can synchronize due to the local hydrodynamic friction time evolution equations for θA (t) and θB (t) which areforces [27]. On the other hand, hydrodynamic interaction often used to describe synchronization phenomena [19]:between the flagella is indispensable for the net swimming.
θ̇ = Ω + α sin[θ (t) − φ (t)],
(7)
A

A

A

θ̇B = Ω + α sin[θB (t) − φB (t)],(8)
Model of an autonomous microswimmer. – As
schematically shown in fig. 1, the present model consists of
where Ω is the constant frequency, α is the coupling pathree hard spheres of the same radius a connected by two
rameter describing the strength of synchronization, andharmonic springs characterized by the spring constant K.
φA and φB are the mechanical phases as explained below.
The total elastic energy is given by
To define the above mechanical phases for a three-spheremodel, it is convenient to introduce the following springKK(1) lengths uA and uB with respect to !:E = (x2 − x1 − !A )2 + (x3 − x2 − !B )2 ,22
where xi (t) (i = 1, 2, 3) are the positions of the threespheres in a one-dimensional coordinate system and weassume x1 < x2 < x3 without loss of generality. In
the above, !A (t) and !B (t) are the natural lengths of thesprings and their dynamics will be explained later (seeeqs. (7) and (8)). Each sphere exerts a force on the viscous fluid of shear viscosity η and experiences an opposite
force from it.
Denoting the velocity of each sphere by ẋi = dxi /dt
and the force acting on each sphere by fi , we can writeuA (t) = x2 (t) − x1 (t) − !,(9)uB (t) = x3 (t) − x2 (t) − !.

(10)Obviously, these quantities are related to the sphere velocities as u̇A = ẋ2 − ẋ1 and u̇B = ẋ3 − ẋ2 . Then thetime-dependent mechanical phases φA and φB are introduced by the relative positions and the velocities of thespheres as34001-p2

cos φA = uA /DA ,sin φA = −u̇A /(ΩDA ),(11)

cos φB = uB /DB ,sin φB = −u̇B /(ΩDB ),(12)

Autonomous three-sphere microswimmer

Fig. 2: Dynamics of θA describing the phase of the naturallength (see eq. (5)) and φA describing the mechanical phase(see eq. (11)). When θA > φA at t, as shown in the left figure,and when α > 0 in eq. (7), the velocity θ̇A becomes larger at alater time t + ∆t, as shown in the right figure. As a result, thedifference between θA and φA also increases at t+∆t. A similar
dynamics occurs also for θB and φB .

where DA(B) = [u2A(B) + (u̇A(B) /Ω)2 ]1/2 . Physically, themechanical phase φA(B) specifies the position in the phasespace of a micromachine spanned by uA(B) and u̇A(B) as
shown in fig. 2.
The above equations complete our model for an autonomous three-sphere microswimmer. In this letter, weshall consider the case of α ≥ 0. Then the physical meaning of eqs. (7) and (8) is that the phase θA (θB ) for thenatural length and the mechanical phase φA (φB ) tend to
be different due to the coupling term, as schematicallyshown in fig. 2. Since the middle sphere is connectedto the other two spheres, our model contains a feedback
mechanism that regulates the dynamics of the two natural lengths !A and !B . Such a coupling effect in the spring
motions gives rise to a non-reciprocal body motion and results in an autonomous locomotion of the microswimmer.
Although α in eqs. (7) and (8) can be different, we shallfirst stick to the symmetric case for the sake of simplicity.
In general, the other quantities such as K, !, and Ω canalso be asymmetric.
Let us define the time-dependent phase difference between the oscillations in the natural lengths by δ(t) =θB (t) − θA (t). When α = 0, the present model reduces
to that of the original elastic microswimmer [9,13]. In
this limit, we have θA (t) = Ωt and θB (t) = Ωt + δ0 , whereδ0 = δ(0) is the initial phase difference. According to Purcell’s scallop theorem [2,3], an elastic microswimmer can
exhibit a directed motion when δ0 #= 0, ±π, i.e., a nonreciprocal motion. Hence the initial phase difference δ0 andthe frequency Ω fully determines the average velocity oflocomotion when α = 0 [9,13].
When the coupling effect is present, however, we showthat a stable phase difference δ controls the dynamics ofa micromachine irrespective of its initial value δ0 . Moreover, the transition to a nonreciprocal motion as well as
the average velocity can be precisely tuned by the coupling parameter α and the velocity is not solely fixedby the externally given frequency Ω as in the previous
models [4,5,9,13].

Fig. 3: The plots of dimensionless center-of-mass position X̂
of an autonomous three-sphere microswimmer as a functionof dimensionless time t̂ for Ω̂ = 0.1 when the initial phasedifferences are (a) δ0 = −π/2 and (b) δ0 = −39π/40. In bothplots, the dimensionless coupling parameter is chosen as α̂ =0.1 (black), 0.5 (red), and 2 (green).

For numerical simulations, it is convenient to introduce
a characteristic time scale defined by
τ=

6πηa
,
K

(13)

which represents the spring relaxation time. Then we use! to scale all the relevant lengths (such as xi , a, and d) andemploy τ to scale the quantities related to time (such asΩ and α). All the dimensionless variables and parameters
are written with a hat such as x̂i = xi /!, Ω̂ = Ωτ , andα̂ = ατ .
Simulation results. – First we have performed computer simulations by numerically solving eq. (2) together
with eqs. (5)–(12) with the use of Euler’s method. The
parameters to characterize the swimmer size are chosen asâ = 0.01 and dˆ = 0.1, satisfying the conditions a, d " !.
Concerning the initial conditions, we put the three spheresat x̂1 (0) = −1, x̂2 (0) = 0, and x̂3 (0) = 1, whereas
the initial phase difference, δ0 , is varied within the range−π ≤ δ0 ≤ π. In the present work, we focus on the lowfrequency regime, Ω̂ < 1. The following simulation results
do not depend on the initial positions of the three spheres.
In fig. 3, we plot the dimensionless center-of-mass position X̂ as a function of time t̂ for different values of
the coupling parameter α̂ when (a) δ0 = −π/2 and(b) δ0 = −39π/40, whereas the frequency is fixed toΩ̂ = 0.1. Although X̂ also oscillates in time at muchsmaller time scales, as shown later in fig. 4(a), one can extract an average steady state velocity V∞ in the long-time
limit by fitting with a straight line. We estimate such an
average velocity V∞ for each curve in fig. 3, and regard
it as an autonomously determined steady state velocity.
For α̂ = 0.1 (black), V∞ vanishes both in figs. 3(a) and(b). For α̂ = 0.5 (red), on the other hand, V∞ is finite
in fig. 3(a) but vanishes in (b). In this case, the steadystate velocity depends on δ0 . For α̂ = 2 (green), V∞ is thesame between figs. 3(a) and (b), showing that V∞ does
not depend on the initial phase difference δ0 although thesign of V∞ can change as we show later in fig. 5(a).

34001-p3

Katsutomo Era et al.

Fig. 4: The plots of (a) dimensionless center-of-mass positiondifference ∆X̂ = X̂(t̂ + ∆t̂) − X̂(t̂) and (b) the phase differenceδ = θB − θA between the oscillations in the natural lengths asa function of dimensionless time difference ∆t̂ measured fromt̂ = 199,900 when Ω̂ = 0.1 and δ0 = −π/2. In both plots, thedimensionless coupling parameter is chosen as α̂ = 0.1 (black),
0.5 (red), and 2 (green). The average steady state velocity V∞
is obtained by fitting with a straight line, whereas the steadystate phase difference δ∞ is obtained by averaging over a cyclein the oscillations of δ.

Fig. 5: The plots of (a) dimensionless stationary velocity V̂∞and (b) stationary phase difference δ∞ as a function of theinitial phase difference δ0 (−π ≤ δ0 ≤ π) when Ω̂ = 0.1. In
both plots, the dimensionless coupling parameter is chosen asα̂ = 0.1 (black), 0.5 (red), and 2 (green). Notice that, for eachcolor, there are multiple data points close to δ0 = ±π.

In figs. 4(a) and (b), the behaviors of the center-ofmass position difference ∆X̂ = X̂(t̂ + ∆t̂) − X̂(t̂) andthe phase difference δ at much smaller time scales areplotted, respectively, as a function of the time difference0 ≤ ∆t̂ ≤ 100 measured after t̂ = 199,900. Since the otherparameters are δ0 = −π/2 and Ω̂ = 0.1, fig. 4(a) is themagnification of fig. 3(a) in the long-time limit after thesteady state has been reached. It is important to notethat both ∆X̂ and δ exhibit oscillatory behaviors whose
period becomes smaller as α is increased. Such a change
in the period is consistent with a perturbation expansionof eqs. (7) and (8) in terms of α, as we shall explain later.
In fig. 4(b), the phase difference δ oscillates around a constant value that can be regarded as the steady state phasedifference δ∞ . Here, we define δ∞ as the average over acycle in the oscillations of δ. Since there is always a welldefined steady state for a given set of parameters, further

Fig. 6: The plots of (a) dimensionless stationary velocity |V̂∞ |and (b) stationary phase difference |δ∞ | as a function of thedimensionless coupling parameter α̂. In both plots, the dimensionless frequency is chosen as Ω̂ = 0.05 (black), 0.1 (orange),
and 0.2 (blue), while δ0 = −π/2 is fixed. There is a lower critical value αc above which both |V∞ | and |δ∞ | become nonzero.
|V∞ | and |δ∞ | take maximum values at αm > αc , and they
vanish for large α.

simulations have been performed for different values of δ0
to investigate the behaviors of V∞ and δ∞ systematically.
Fixing the frequency to Ω̂ = 0.1, we plot in figs. 5(a)and (b) the steady state velocity V̂∞ and the phase difference δ∞ , respectively, as a function of the initial phasedifference δ0 for the range −π ≤ δ0 ≤ π. Different colors indicate different α̂ values as we have used in fig. 3. In
fig. 5(a), we see that V∞ either vanishes or takes a nonzero
constant value within a certain range of δ0 . This means
that, under certain conditions, the proposed microswimmer can autonomously determine its steady state velocity as well as the phase difference. We also see that V̂∞changes its sign at δ0 = 0 although the absolute value isthe same. The sign of V̂∞ and δ∞ also changes for α̂ = 0.5(red) and 2 (green) when δ0 becomes close to ±π.
For α̂ = 0.5 (red) in fig. 5(a), the velocity V∞ tends
to vanish when the initial phase difference δ0 is close to±π. In this situation, we see in fig. 5(b) that the steadystate phase difference approaches δ∞ = ±π, i.e., a reciprocal motion. When V∞ is finite in fig. 5(a) for α̂ = 0.5(red) and 2 (green), on the other hand, the corresponding
phase difference is δ∞ #= 0, ±π, i.e., a nonreciprocal motion. These results are in accordance with Purcell’s scallop
theorem [2,3]. A more detailed discussion concerning thestability of the phase difference will be given later in fig. 7.
When α̂ = 0.5, there are two stable fixed points; one withfinite V∞ and the other with vanishing V∞ .
In figs. 6(a) and (b), we plot |V̂∞ | and |δ∞ |, respectively,
as a function of α̂ for different frequencies Ω̂ = 0.05, 0.1,
and 0.2. To make these plots, we have used δ0 = −π/2.
When Ω̂ = 0.1 (orange), for example, there is a finite critical value of α̂c ≈ 0.2 above which |V̂∞ | and |δ∞ | become
nonzero. For α̂ < α̂c , on the other hand, both |V̂∞ | and|δ∞ | vanish. The existence of such a finite critical value
α̂c is a nontrivial outcome of the present model. When α̂
is very large, such as α̂ ≥ 12.5 for Ω̂ = 0.1, both |V̂∞ | and|δ∞ | vanish again. Hence autonomous locomotion can be34001-p4

Autonomous three-sphere microswimmer
fixed points at δ = ±π appear. These new fixed points result in a reciprocal motion, prohibiting the locomotion ofa microswimmer. For larger coupling parameter α̂ ≥ 12.5,
δ = 0 becomes stable again. Although such a stability
diagram depends on Ω̂, the general structure of the bifurcation diagram remains the same.

Fig. 7: The numerically obtained stability diagram in the planeof the phase difference δ and the dimensionless coupling parameter α̂ when Ω̂ = 0.1. The orange circles indicate stable fixed
points where δ̇ = 0 holds. For α̂ > α̂c ≈ 0.2, the stable pointδ = 0 bifurcates into two stable fixed points. These nonzero
fixed points correspond to nonreciprocal motions (δ $= 0, ±π),leading to a nonzero velocity V∞ . There are other stable fixed
points at δ = ±π when 0.39 ≤ α̂ ≤ 1.58. For larger α̂, δ = 0
becomes stable again for α̂ ≥ 12.5. The dashed lines indicate
the numerically determined separatrices which are obtained bycomparing the initial phase difference δ0 and the stationary
phase difference δ∞ .

Linear stability analysis in the weak coupling
limit. – Although the equations of motion of our modelare highly nonlinear, we can analytically investigate thelinear stability of the phase difference δ when α̂ is smallenough. In other words, we consider the case α̂ < α̂c infig. 7, for which δ = 0 is the only stable point and a microswimmer exhibits a reciprocal motion. Here, we shallexpress δ̇ in terms of δ and perform a stability analysis.
Starting from eq. (2), we first neglect hydrodynamic interactions by considering the case a " !. Then the equationsof motion for uA and uB (see eqs. (9) and (10)) are approximated as1[−2(uA − d cos θA ) + (uB − d cos θB )], (14)τ
1u̇B ≈ [(uA − d cos θA ) − 2(uB − d cos θB )] , (15)
τ

u̇A ≈

where τ is the spring relaxation time introduced ineq. (13).
When the coupling parameter α is small enough ineqs. (7) and (8), one can assume that δ is almost constantand the phases of the two natural lengths can be approxachieved for a finite range of the coupling parameter α̂.
imated as θA (t) ≈ Ωt and θB (t) ≈ Ωt + δ. According to
Such a behavior is common for other frequencies Ω̂.
ref. [9], the coupled linear equations in eqs. (14) and (15)Moreover, it is interesting to note that |V∞ | takes maxcan be solved in the frequency domain. By performing theimum values such as at α̂m ≈ 5.0 when Ω̂ = 0.1. Hence
inverse Fourier transform, we havethe present autonomous microswimmer can maximize its
velocity by tuning the coupling parameter α. Notice thatduA (t) ≈both α̂c and α̂m depend on the frequency Ω̂, and they are
9 + 10Ω̂2 + Ω̂4
not universal quantities. However, it is worth mentioning
×[(9 + 5Ω̂2 ) cos(Ωt) − 4Ω̂2 cos(Ωt + δ)that we find the relation α̂c /Ω̂ ≈ 2 for all Ω̂ chosen in our
simulations.
+(6Ω̂ + 2Ω̂3 ) sin(Ωt)
From these simulation results, one can discuss the sta(16)+(3Ω̂ − Ω̂3 ) sin(Ωt + δ)],bility of the relative phase difference δ. Comparing itsduB (t) ≈initial value δ0 and the steady state value δ∞ , we can
9 + 10Ω̂2 + Ω̂4identify the stable fixed points. Such a stability diagram
×[−4Ω̂2 cos(Ωt) + (9 + 5Ω̂2 ) cos(Ωt + δ)
for Ω̂ = 0.1 is presented in fig. 7 in the plane of δ and α̂,describing the bifurcation structure of the present model.
+(3Ω̂ − Ω̂3 ) sin(Ωt)
The orange circles indicate the stable fixed points where(17)+(6Ω̂ + 2Ω̂3 ) sin(Ωt + δ)],δ̇ = 0 holds. In fig. 7 we have also plotted the numericallydetermined separatrices by the dashed lines. These points where Ω̂ = Ωτ as before and δ here is a constant.
are obtained by comparing the initial phase difference δ0
The above results can be inserted into eqs. (11) and (12)and the stationary phase difference δ∞ . Hence they areto obtain φA (t) and φB (t), respectively. Considering the
not mathematically obtained rigorous unstable points.
low-frequency limit, Ω̂ " 1, we expand φA and φB up to
For α̂ < α̂c ≈ 0.2, the stable fixed points exist only at
the second order in Ω̂:δ = 0. As α̂ is increased, the stable point at δ = 0 bifurcates into two stable fixed points for α̂ > α̂c , whereas δ = 0(2 − cos δ) sin δ 2
2 + cos δ
Ω̂ −
Ω̂ ,
(18)
φA (t) ≈ Ωt −becomes unstable. The stable fixed points at nonzero δ39correspond to nonreciprocal motions, leading to finite V∞(2 − cos δ) sin δ 2
2 + cos δφB (t) ≈ Ωt + δ −Ω̂ +Ω̂ . (19)and δ∞ . When α̂ satisfies 0.39 ≤ α̂ ≤ 1.58, the two stable
3934001-p5

Katsutomo Era et al.
fixed points at nonzero δ for α > αc . The difficultyarises because eqs. (16) and (17) are correct only for small
α. Moreover, hydrodynamic interactions, which are neglected in the above analysis, need to be further taken into
account to fully discuss the bifurcation structure of themodel.
It is worth mentioning, however, that the stability of thephase difference δ can be determined even in the absenceof hydrodynamic interactions as we have discussed in thissection. For a three-sphere model of Chlamydomonas, it
was shown that hydrodynamic interactions contribute little to synchronization [27]. In the present model as well as
in the previous models [4,5,9,13], hydrodynamic interactions play an essential role for the locomotion of a threesphere microswimmer.
Fig. 8: The plot of δ̇/(2α) (see eq. (22)) as a function of δfor different dimensionless frequencies Ω̂ = 0.05 (black), 0.1(orange), and 0.2 (blue). Here δ = 0 corresponds to the stable fixed point, while δ = ±π are unstable ones. This result is in agreement with the stability diagram in fig. 7 when
α < αc .

Then we substitute these expressions into eqs. (7) and (8)to obtain&'
2 + cos δ(2 − cos δ) sin δ 2Ω̂ +Ω̂ , (20)θ̇A ≈ Ω + α sin
39&'2 + cos δ(2 − cos δ) sin δ 2θ̇B ≈ Ω + α sinΩ̂ −Ω̂ , (21)39where we have used the approximations θA (t) ≈ Ωt andθB (t) ≈ Ωt + δ. The stability of δ can be discussed interms of δ̇ = θ̇B − θ̇A that is given by&''&(−2 + cos δ) sin δ 2δ̇
2 + cos δ≈ sinΩ̂ cosΩ̂ . (22)2α93Summary and discussion. – In this letter, we haveproposed a model of an autonomous three-sphere microswimmer by considering a coupling effect between thetwo natural lengths of an elastic microswimmer [9]. Ourmodel is motivated by the previous models for synchronization phenomena in coupled oscillator systems [19–22].
Performing numerical simulations, we have shown that a
microswimmer can acquire a nonzero steady state velocity
V∞ that is almost independent of the initial phase difference δ0 (see fig. 5(a)). The corresponding phase difference δ∞ between the oscillations in the natural lengths
becomes also finite (see fig. 5(b)), which is consistent withPurcell’s scallop theorem for microswimmers in a viscous
fluid [2,3].
We have explored in detail the dependencies of V∞ andδ∞ on the coupling parameter α and the frequency Ω.
We find that both |V∞ | and |δ∞ | take nonzero values for
α > αc , and they also show maximum values at αm (fig. 6).
There is a finite range of α for which a microswimmer canhave an autonomous directed motion. We have also analyzed the stability of the phase difference δ by constructing
a stability diagram (fig. 7). This result has been analytically confirmed in the limit of small α (fig. 8).
In the present work, we have discussed the case when
the frequency Ω is small enough, i.e., Ω̂ < 1. When Ωis made larger, the difference between θA(B) (phases of
the natural lengths) and φA(B) (mechanical phases) becomes also larger. In the original elastic microswimmer
without any coupling effect, it was shown that the average velocity decreases with increasing frequency in thehigh-frequency limit due to the intrinsic spring relaxation
dynamics [9,13]. Such a reduction of the velocity in thehigh-frequency regime also occurs for a three-sphere microswimmer moving in a viscoelastic medium [29–31]. Our
future numerical and analytical studies include not only
the high-frequency behavior of the model but also the caseof α < 0.

In fig. 8, we plot δ̇/(2α) as a function of δ for Ω̂ = 0.05,
0.1, and 0.2. We first note that eq. (22) is an odd functionof δ. Since δ̇ > 0 for δ < 0 and δ̇ < 0 for δ > 0, we
find that δ = 0 is a stable fixed point. This result isin accordance with the stability diagram in fig. 7 when
α < αc . The slope at δ = 0 becomes steeper as Ω̂ isincreased. We also see that δ = ±π are the unstable fixed
points.
Up to the first order in Ω̂, we see in eqs. (18) and (19)that the mechanical phases φA and φB are delayed withrespect to those of the natural length θA and θB , respectively. Thus, the time evolutions of θA and θB are accelerated by the second terms in eqs. (20) and (21) that arecontrolled by α. Hence, as α is increased, the oscillation
frequencies of ∆X̂ and δ become larger than the original spring frequency Ω. Such a change of the oscillation∗∗∗
frequency was shown in fig. 4(b).
Within the present approximation, however, we cannot
We thank T. Kato for useful discussions. YK acknowlanalytically predict the critical value αc nor the stable edges support by a Grant-in-Aid for JSPS Fellows (Grant
34001-p6

Autonomous three-sphere microswimmer
No. JP19J00365) from the Japan Society for the Promotion of Science (JSPS). YH acknowledges support by aGrant-in-Aid for JSPS Fellows (Grant No. 19J20271) fromthe JSPS. KY acknowledges support by a Grant-in-Aid
for JSPS Fellows (Grant No. 18J21231) from the JSPS.
YK, HK and SK acknowledge support by a Grant-in-Aid
for Scientific Research (C) (Grant No. 19K03765) fromthe JSPS. SK further acknowledges support by a Grantin-Aid for Scientific Research (C) (Grant No. 18K03567)
from the JSPS, and support by a Grant-in-Aid for Scientific Research on Innovative Areas “Information Physics of
Living Matters” (Grant No. 20H05538) from the Ministryof Education, Culture, Sports, Science and Technology ofJapan.
REFERENCES[1] Lauga E. and Powers T. R., Rep. Prog. Phys., 72(2009) 096601.
[2] Purcell E. M., Am. J. Phys., 45 (1977) 3.
[3] Lauga E., Soft Matter, 7 (2011) 3060.
[4] Najafi A. and Golestanian R., Phys. Rev. E, 69 (2004)
062901.
[5] Golestanian R. and Ajdari A., Phys. Rev. E, 77 (2008)
036308.
[6] Leoni M., Kotar J., Bassetti B., Cicuta P. and
Lagomarsino M. C., Soft Matter, 5 (2009) 472.
[7] Grosjean G., Hubert M., Lagubeau G. and Vandewalle N., Phys. Rev. E, 94 (2016) 021101(R).
[8] Grosjean G., Hubert M. and Vandewalle N., Adv.
Colloid Interface Sci., 255 (2018) 84.
[9] Yasuda K., Hosaka Y., Kuroda M., Okamoto R. andKomura S., J. Phys. Soc. Jpn., 86 (2017) 093801.
[10] Dunkel J. and Zaid I. M., Phys. Rev. E, 80 (2009)021903.
[11] Pande J. and Smith A.-S., Soft Matter, 11 (2015) 2364.
[12] Pande J., Merchant L., Krüger T., Harting J. andSmith A.-S., New J. Phys., 19 (2017) 053024.

[13] Kuroda M., Yasuda K. and Komura S., J. Phys. Soc.
Jpn., 88 (2019) 054804.
[14] Hosaka Y., Yasuda K., Sou I., Okamoto R. and Komura S., J. Phys. Soc. Jpn., 86 (2017) 113801.
[15] Sou I., Hosaka Y., Yasuda K. and Komura S., Phys.
Rev. E, 100 (2019) 022607.
[16] Sou I., Hosaka Y., Yasuda K. and Komura S., Physica
A, 562 (2021) 125277.
[17] Owaki D., Kano T., Nagasawa K., Tero A. and
Ishiguro A., J. R. Soc. Interface, 10 (2012) 20120669.
[18] Fukuhara A., Owaki D., Kano T., Kobayashi R. and
Ishiguro A., Adv. Robot., 32 (2018) 794.
[19] Kuramoto Y., Chemical Oscillations, Waves, and Turbulence (Springer-Verlag, New York) 1984.
[20] Pikovsky A., Rosenblum M. and Kurths J., Synchronization: A Universal Concept in Nonlinear Sciences(Cambridge University Press, Cambridge) 2001.
[21] Strogatz S. H., Sync: The Emerging Science of Spontaneous Order (Hyperion, New York) 2003.
[22] Acebrón J. A., Bonilla L. L., Pérez Vicente C. J.,
Ritort F. and Spigler R., Rev. Mod. Phys., 77 (2005)
137.
[23] Golestanian R., Yeomans J. M. and Uchida N., Soft
Matter, 7 (2011) 3074.
[24] Uchida N., Golestanian R. and Bennett R. R.,
J. Phys. Soc. Jpn., 86 (2017) 101007.
[25] Polin M., Tuval I., Drescher K., Gollub J. P. andGoldstein R. E., Science, 325 (2009) 487.
[26] Goldstein R. E., Polin M. and Tuval I., Phys. Rev.
Lett., 103 (2009) 168103.
[27] Friedrich B. M. and Jülicher F., Phys. Rev. Lett., 109(2012) 138102.
[28] Bennett R. R. and Golestanian R., Phys. Rev. Lett.,
110 (2013) 148102.
[29] Yasuda K., Okamoto R. and Komura S., J. Phys. Soc.
Jpn., 86 (2017) 043801.
[30] Yasuda K., Okamoto R. and Komura S., EPL, 123(2018) 34002.
[31] Yasuda K., Kuroda M. and Komura S., Phys. Fluids,
32 (2020) 093102.

34001-p7