Autonomous elastic microswimmer

Microswimmers are small machines that swim in a fluid and they are expected to be used in microfluidics and microsystems [1]. Over the length scale 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 and Golestanian proposed a three-sphere swimmer (NG swimmer) [4,5], in which three in-line spheres are linked by two arms 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 proposed a generalized three-sphere microswimmer model in which the spheres are connected by two harmonic springs, i.e., an elastic microswimmer [9]. Compared with the NG swimmer, the main difference is that the natural length of each spring (rather than the arm length) is assumed to undergo a prescribed cyclic motion. A similar model was proposed by other people [10–12]. We have analytically obtained the average 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 of the springs (elastic swimmer) are assumed to undergo a prescribed cyclic motion. Such active motions can lead to a net locomotion if the swimming strokes are nonreciprocal. In these models, the average swimming velocity is purely 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 by itself rather than being imposed externally. Moreover, a sophisticated microswimmer requires a feedback control system in order to regulate the switching between the static 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 coordination in 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-sphere swimmer which can autonomously determine its velocity. In order to implement such a control mechanism, we introduce a coupling between the two natural lengths of an elastic microswimmer by using the interaction adopted in the Kuramoto model for coupled oscillators [19–22]. Importantly, the proposed microswimmer acquires a steady state velocity and a finite phase difference in the long-time limit without any external control. The steady state velocity can be mainly tuned by changing the coupling parameter in the model. Moreover, we investigate the condition that a microswimmer can attain an autonomous locomotion, and further perform a linear stability analysis of the steady state.

Synchronization phenomena are widely observed in active biological systems such as flagella and cilia [23,24]. In particular, synchronization of a pair of flagella in Chlamydomonas was observed experimentally [25,26]. For a three-sphere model of Chlamydomonas in which the spheres representing the flagella move on circular trajectories relative to the body sphere [27,28], the two flagella can synchronize due to the local hydrodynamic friction forces [27]. On the other hand, hydrodynamic interaction between the flagella is indispensable for the net swimming.

As schematically shown in fig. 1, the present model consists of three hard spheres of the same radius a connected by two harmonic springs characterized by the spring constant K. The total elastic energy is given by

$$\begin{eqnarray} &&E = \frac{K}{2}(x_2 - x_1 - \ell_{\rm A})^2 +\frac{K}{2}(x_3 - x_2 - \ell_{\rm B})^2, \end{eqnarray} \tag{ 1 }$$

where $x_i(t)~(i=1, 2, 3)$ are the positions of the three spheres in a one-dimensional coordinate system and we assume x₁ < x₂ < x₃ without loss of generality. In the above, $\ell_{\rm A}(t)$ and $\ell_{\rm B}(t)$ are the natural lengths of the springs and their dynamics will be explained later (see eqs. (7) and (8)). Each sphere exerts a force on the viscous fluid of shear viscosity η and experiences an opposite force from it.

Zoom In Zoom Out Reset image size

Denoting the velocity of each sphere by $\dot x_i=\mathrm{d}x_i/\mathrm{d}t$ and the force acting on each sphere by f_i , we can write the equations of motion of each sphere as [9,13]

$$\begin{equation} \dot x_i = \sum_{j=1}^{3}M_{ij}f_j, \end{equation} \tag{ 2 }$$

where the three forces f_i are given by

$$\begin{equation} f_i=-\frac{\partial E}{\partial x_i}. \end{equation} \tag{ 3 }$$

Here the details of the hydrodynamic interactions are taken into account through the mobility coefficients M_ij . Within Oseen's approximation, which is justified when the spheres are considerably far from each other $(a \ll \vert x_i-x_j \vert)$, the expressions for the mobility coefficients M_ij can be written as

$$\begin{equation} M_{ij}=\begin{cases} \displaystyle\frac{1}{6\pi\eta a}, & \quad i=j,\\ \displaystyle\frac{1}{4\pi\eta \vert x_i-x_j \vert}, & \quad i\neq j. \end{cases} \end{equation} \tag{ 4 }$$

The force-free condition, $f_1+f_2+f_3=0$, is automatically satisfied in the present model [9,13]. We define the center-of-mass position of a microswimmer by $X(t)=(x_1+x_2+x_3)/3$ and the swimming velocity of the whole object by $V(t) = \dot{X}(t)$.

Next, we consider that the two natural lengths of the springs undergo the following cyclic changes in time [9,13]:

$$\begin{equation} \ell_{\rm A}(t) = \ell+d \cos \theta_{\rm A}(t), \end{equation} \tag{ 5 }$$

$$\begin{equation} \ell_{\rm B}(t) = \ell+d \cos \theta_{\rm B}(t), \end{equation} \tag{ 6 }$$

where ℓ is the constant natural length, d is the oscillation amplitude, $\theta_{\rm A}(t)$ and $\theta_{\rm B}(t)$ are the time-dependent phases. The most important aspect of our model is that $\theta_{\rm A}(t)$ and $\theta_{\rm B}(t)$ are affected by the relative positions and the velocities of the three spheres. We employ the following time evolution equations for $\theta_{\rm A}(t)$ and $\theta_{\rm B}(t)$ which are often used to describe synchronization phenomena [19]:

$$\begin{equation} \dot{\theta}_{\rm A} = \Omega+\alpha\sin[\theta_{\rm A}(t)-\phi_{\rm A}(t)], \end{equation} \tag{ 7 }$$

$$\begin{equation} \dot{\theta}_{\rm B} = \Omega+\alpha\sin[\theta_{\rm B}(t)-\phi_{\rm B}(t)], \end{equation} \tag{ 8 }$$

where Ω is the constant frequency, α is the coupling parameter describing the strength of synchronization, and $\phi_{\rm A}$ and $\phi_{\rm B}$ are the mechanical phases as explained below.

To define the above mechanical phases for a three-sphere model, it is convenient to introduce the following spring lengths $u_{\rm A}$ and $u_{\rm B}$ with respect to ℓ:

$$\begin{equation} u_{\rm A}(t) = x_2(t)-x_1(t)-\ell, \end{equation} \tag{ 9 }$$

$$\begin{equation} u_{\rm B}(t) = x_3(t)-x_2(t)-\ell. \end{equation} \tag{ 10 }$$

Obviously, these quantities are related to the sphere velocities as $\dot u_{\rm A} = \dot x_2- \dot x_1$ and $\dot u_{\rm B} = \dot x_3- \dot x_2$. Then the time-dependent mechanical phases $\phi_{\rm A}$ and $\phi_{\rm B}$ are introduced by the relative positions and the velocities of the spheres as

$$\begin{equation} \cos \phi_{\rm A} = u_{\rm A}/D_{\rm A},\quad \sin \phi_{\rm A}=-\dot{u}_{\rm A}/(\Omega D_{\rm A}), \end{equation} \tag{ 11 }$$

$$\begin{equation} \cos \phi_{\rm B} = u_{\rm B}/D_{\rm B},\quad \sin \phi_{\rm B}=-\dot{u}_{\rm B}/(\Omega D_{\rm B}), \end{equation} \tag{ 12 }$$

where $D_{\rm A(B)}=[u_{\rm A(B)}^{2}+(\dot{u}_{\rm A(B)}/\Omega)^2]^{1/2}$. Physically, the mechanical phase $\phi_{\rm A(B)}$ specifies the position in the phase space of a micromachine spanned by $u_{\rm A(B)}$ and $\dot{u}_{\rm A(B)}$ as shown in fig. 2.

Zoom In Zoom Out Reset image size

The above equations complete our model for an autonomous three-sphere microswimmer. In this letter, we shall consider the case of $\alpha \ge 0$. Then the physical meaning of eqs. (7) and (8) is that the phase $\theta_{\rm A}~(\theta_{\rm B})$ for the natural length and the mechanical phase $\phi_{\rm A}~(\phi_{\rm B})$ tend to be different due to the coupling term, as schematically shown in fig. 2. Since the middle sphere is connected to the other two spheres, our model contains a feedback mechanism that regulates the dynamics of the two natural lengths $\ell_{\rm A}$ and $\ell_{\rm 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 shall first stick to the symmetric case for the sake of simplicity. In general, the other quantities such as K, ℓ, and Ω can also be asymmetric.

Let us define the time-dependent phase difference between the oscillations in the natural lengths by $\delta (t)=\theta_{\rm B}(t)-\theta_{\rm A}(t)$. When $\alpha=0$, the present model reduces to that of the original elastic microswimmer [9,13]. In this limit, we have $\theta_{\rm A}(t) = \Omega t$ and $\theta_{\rm B}(t) = \Omega t + \delta_0$, where $\delta_0=\delta(0)$ is the initial phase difference. According to Purcell's scallop theorem [2,3], an elastic microswimmer can exhibit a directed motion when $\delta_0 \neq 0, \pm \pi$, i.e., a nonreciprocal motion. Hence the initial phase difference $\delta_0$ and the frequency Ω fully determines the average velocity of locomotion when $\alpha=0$ [9,13].

When the coupling effect is present, however, we show that a stable phase difference δ controls the dynamics of a micromachine irrespective of its initial value $\delta_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 fixed by the externally given frequency Ω as in the previous models [4,5,9,13].

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

$$\begin{equation} \tau=\frac{6\pi\eta a}{K}, \end{equation} \tag{ 13 }$$

which represents the spring relaxation time. Then we use ℓ to scale all the relevant lengths (such as x_i , a, and d) and employ τ to scale the quantities related to time (such as Ω and α). All the dimensionless variables and parameters are written with a hat such as $\hat{x}_i=x_i/\ell$, $\hat{\Omega}=\Omega \tau$, and $\hat{\alpha}=\alpha \tau$.

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 as $\hat{a}=0.01$ and $\hat{d}=0.1$, satisfying the conditions $a, d \ll \ell$. Concerning the initial conditions, we put the three spheres at $\hat{x}_1(0)=-1$, $\hat{x}_2(0)=0$, and $\hat{x}_3(0)=1$, whereas the initial phase difference, $\delta_0$, is varied within the range $-\pi \leq \delta_0 \leq \pi$. In the present work, we focus on the low-frequency regime, $\hat{\Omega}<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 $\hat{X}$ as a function of time $\hat{t}$ for different values of the coupling parameter $\hat{\alpha}$ when (a) $\delta_0=-\pi/2$ and (b) $\delta_0=-39\pi/40$, whereas the frequency is fixed to $\hat{\Omega}=0.1$. Although $\hat{X}$ also oscillates in time at much smaller time scales, as shown later in fig. 4(a), one can extract an average steady state velocity $V_\infty$ in the long-time limit by fitting with a straight line. We estimate such an average velocity $V_\infty$ for each curve in fig. 3, and regard it as an autonomously determined steady state velocity. For $\hat{\alpha}=0.1$ (black), $V_\infty$ vanishes both in figs. 3(a) and (b). For $\hat{\alpha}=0.5$ (red), on the other hand, $V_\infty$ is finite in fig. 3(a) but vanishes in (b). In this case, the steady state velocity depends on $\delta_0$. For $\hat{\alpha}=2$ (green), $V_\infty$ is the same between figs. 3(a) and (b), showing that $V_\infty$ does not depend on the initial phase difference $\delta_0$ although the sign of $V_\infty$ can change as we show later in fig. 5(a).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In figs. 4(a) and (b), the behaviors of the center-of-mass position difference $\Delta \hat{X}=\hat{X}(\hat{t}+\Delta\hat{t})-\hat{X}(\hat{t})$ and the phase difference δ at much smaller time scales are plotted, respectively, as a function of the time difference $0 \le \Delta \hat{t} \le 100$ measured after $\hat{t}=199{,}900$. Since the other parameters are $\delta_0=-\pi/2$ and $\hat{\Omega}=0.1$, fig. 4(a) is the magnification of fig. 3(a) in the long-time limit after the steady state has been reached. It is important to note that both $\Delta \hat{X}$ and δ exhibit oscillatory behaviors whose period becomes smaller as α is increased. Such a change in the period is consistent with a perturbation expansion of 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 phase difference $\delta_\infty$. Here, we define $\delta_\infty$ as the average over a cycle in the oscillations of δ. Since there is always a well-defined steady state for a given set of parameters, further simulations have been performed for different values of $\delta_0$ to investigate the behaviors of $V_\infty$ and $\delta_\infty$ systematically.

Fixing the frequency to $\hat{\Omega}=0.1$, we plot in figs. 5(a) and (b) the steady state velocity $\hat{V}_\infty$ and the phase difference $\delta_\infty$, respectively, as a function of the initial phase difference $\delta_0$ for the range $-\pi \le \delta_0 \le \pi$. Different colors indicate different $\hat{\alpha}$ values as we have used in fig. 3. In fig. 5(a), we see that $V_\infty$ either vanishes or takes a nonzero constant value within a certain range of $\delta_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 $\hat{V}_\infty$ changes its sign at $\delta_0=0$ although the absolute value is the same. The sign of $\hat{V}_\infty$ and $\delta_\infty$ also changes for $\hat{\alpha}=0.5$ (red) and 2 (green) when $\delta_0$ becomes close to $\pm \pi$.

For $\hat{\alpha}=0.5$ (red) in fig. 5(a), the velocity $V_\infty$ tends to vanish when the initial phase difference $\delta_0$ is close to $\pm \pi$. In this situation, we see in fig. 5(b) that the steady state phase difference approaches $\delta_\infty = \pm \pi$, i.e., a reciprocal motion. When $V_\infty$ is finite in fig. 5(a) for $\hat{\alpha}=0.5$ (red) and 2 (green), on the other hand, the corresponding phase difference is $\delta_\infty \neq 0, \pm \pi$, i.e., a nonreciprocal motion. These results are in accordance with Purcell's scallop theorem [2,3]. A more detailed discussion concerning the stability of the phase difference will be given later in fig. 7. When $\hat{\alpha}=0.5$, there are two stable fixed points; one with finite $V_\infty$ and the other with vanishing $V_\infty$.

In figs. 6(a) and (b), we plot $\vert \hat{V}_\infty \vert$ and $\vert \delta_\infty \vert$, respectively, as a function of $\hat{\alpha}$ for different frequencies $\hat{\Omega}=0.05$, 0.1, and 0.2. To make these plots, we have used $\delta_0=-\pi/2$. When $\hat{\Omega}=0.1$ (orange), for example, there is a finite critical value of $\hat{\alpha}_{\rm c} \approx 0.2$ above which $\vert \hat{V}_\infty \vert$ and $\vert \delta_\infty \vert$ become nonzero. For $\hat{\alpha} < \hat{\alpha}_{\rm c}$, on the other hand, both $\vert \hat{V}_\infty \vert$ and $\vert \delta_\infty \vert$ vanish. The existence of such a finite critical value $\hat{\alpha}_{\rm c}$ is a nontrivial outcome of the present model. When $\hat{\alpha}$ is very large, such as $\hat{\alpha}\ge 12.5$ for $\hat{\Omega}=0.1$, both $\vert \hat{V}_\infty \vert$ and $\vert \delta_\infty \vert$ vanish again. Hence autonomous locomotion can be achieved for a finite range of the coupling parameter $\hat{\alpha}$. Such a behavior is common for other frequencies $\hat{\Omega}$.

Zoom In Zoom Out Reset image size

Moreover, it is interesting to note that $\vert V_\infty \vert$ takes maximum values such as at $\hat{\alpha}_{\rm m} \approx 5.0$ when $\hat{\Omega}=0.1$. Hence the present autonomous microswimmer can maximize its velocity by tuning the coupling parameter α. Notice that both $\hat{\alpha}_{\rm c}$ and $\hat{\alpha}_{\rm m}$ depend on the frequency $\hat{\Omega}$, and they are not universal quantities. However, it is worth mentioning that we find the relation $\hat{\alpha}_{\rm c}/\hat{\Omega} \approx 2$ for all $\hat{\Omega}$ chosen in our simulations.

From these simulation results, one can discuss the stability of the relative phase difference δ. Comparing its initial value $\delta_0$ and the steady state value $\delta_\infty$, we can identify the stable fixed points. Such a stability diagram for $\hat{\Omega}=0.1$ is presented in fig. 7 in the plane of δ and $\hat{\alpha}$, describing the bifurcation structure of the present model. The orange circles indicate the stable fixed points where $\dot{\delta}=0$ holds. In fig. 7 we have also plotted the numerically determined separatrices by the dashed lines. These points are obtained by comparing the initial phase difference $\delta_{0}$ and the stationary phase difference $\delta_{\infty}$. Hence they are not mathematically obtained rigorous unstable points.

Zoom In Zoom Out Reset image size

For $\hat{\alpha} < \hat{\alpha}_{\rm c} \approx 0.2$, the stable fixed points exist only at $\delta =0$. As $\hat{\alpha}$ is increased, the stable point at $\delta =0$ bifurcates into two stable fixed points for $\hat{\alpha} > \hat{\alpha}_{\rm c}$, whereas $\delta =0$ becomes unstable. The stable fixed points at nonzero δ correspond to nonreciprocal motions, leading to finite $V_\infty$ and $\delta_\infty$. When $\hat{\alpha}$ satisfies $0.39 \le \hat{\alpha} \le 1.58$, the two stable fixed points at $\delta = \pm \pi$ appear. These new fixed points result in a reciprocal motion, prohibiting the locomotion of a microswimmer. For larger coupling parameter $\hat{\alpha}\ge 12.5$, $\delta =0$ becomes stable again. Although such a stability diagram depends on $\hat{\Omega}$, the general structure of the bifurcation diagram remains the same.

Although the equations of motion of our model are highly nonlinear, we can analytically investigate the linear stability of the phase difference δ when $\hat{\alpha}$ is small enough. In other words, we consider the case $\hat{\alpha} < \hat{\alpha}_{\rm c}$ in fig. 7, for which $\delta=0$ is the only stable point and a microswimmer exhibits a reciprocal motion. Here, we shall express $\dot{\delta}$ in terms of δ and perform a stability analysis. Starting from eq. (2), we first neglect hydrodynamic interactions by considering the case $a \ll \ell$. Then the equations of motion for $u_{\rm A}$ and $u_{\rm B}$ (see eqs. (9) and (10)) are approximated as

$$\begin{equation} \dot{u}_{\rm A} \approx \frac{1}{\tau}[ -2(u_{\rm A}-d\cos\theta_{\rm A}) +(u_{\rm B}-d\cos\theta_{\rm B})], \end{equation} \tag{ 14 }$$

$$\begin{equation} \dot{u}_{\rm B} \approx \frac{1}{\tau}\left[(u_{\rm A}-d\cos\theta_{\rm A}) -2(u_{\rm B}-d\cos\theta_{\rm B})\right], \end{equation} \tag{ 15 }$$

where τ is the spring relaxation time introduced in eq. (13).

When the coupling parameter α is small enough in eqs. (7) and (8), one can assume that δ is almost constant and the phases of the two natural lengths can be approximated as $\theta_{\rm A}(t) \approx \Omega t$ and $\theta_{\rm B}(t) \approx \Omega t +\delta$. According to ref. [9], the coupled linear equations in eqs. (14) and (15) can be solved in the frequency domain. By performing the inverse Fourier transform, we have

$$\begin{eqnarray} u_{\rm A}(t)&\approx \frac{d}{9+10\hat{\Omega}^2+\hat{\Omega}^4}\nonumber\\ & \quad\times [(9+5\hat{\Omega}^2)\cos(\Omega t)-4\hat{\Omega}^2\cos(\Omega t+\delta) \nonumber\\ & \quad+(6\hat{\Omega}+2\hat{\Omega}^3)\sin(\Omega t)\nonumber\\ & \quad+(3\hat{\Omega}-\hat{\Omega}^3)\sin(\Omega t+\delta)],\qquad \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} u_{\rm B}(t)&\approx \frac{d}{9+10\hat{\Omega}^2+\hat{\Omega}^4}\nonumber\\ & \quad\times [-4\hat{\Omega}^2\cos(\Omega t)+(9+5\hat{\Omega}^2) \cos(\Omega t+\delta) \nonumber\\ & \quad+(3\hat{\Omega}-\hat{\Omega}^3)\sin(\Omega t)\nonumber\\ & \quad+(6\hat{\Omega}+2\hat{\Omega}^3)\sin(\Omega t+\delta)],\qquad \end{eqnarray} \tag{ 17 }$$

where $\hat{\Omega}=\Omega \tau$ as before and δ here is a constant.

The above results can be inserted into eqs. (11) and (12) to obtain $\phi_{\rm A}(t)$ and $\phi_{\rm B}(t)$, respectively. Considering the low-frequency limit, $\hat{\Omega} \ll 1$, we expand $\phi_{\rm A}$ and $\phi_{\rm B}$ up to the second order in $\hat{\Omega}$:

$$\begin{equation} \phi_{\rm A}(t) \approx \Omega t -\frac{2+\cos \delta}{3}\hat{\Omega}-\frac{(2-\cos \delta)\sin \delta}{9} \hat{\Omega}^2, \end{equation} \tag{ 18 }$$

$$\begin{equation} \phi_{\rm B}(t) \approx \Omega t +\delta -\frac{2+\cos \delta}{3}\hat{\Omega} +\frac{(2-\cos \delta)\sin \delta}{9} \hat{\Omega}^2. \end{equation} \tag{ 19 }$$

Then we substitute these expressions into eqs. (7) and (8) to obtain

$$\begin{equation} \dot{\theta}_{\mathrm A} \approx \Omega + \alpha \sin \left [ \frac{2+\cos \delta}{3}\hat{\Omega} + \frac{(2-\cos \delta)\sin \delta}{9} \hat{\Omega}^2 \right], \end{equation} \tag{ 20 }$$

$$\begin{equation} \dot{\theta}_{\mathrm B} \approx \Omega + \alpha \sin \left [ \frac{2+\cos \delta}{3}\hat{\Omega} - \frac{(2-\cos \delta)\sin \delta}{9} \hat{\Omega}^2 \right], \end{equation} \tag{ 21 }$$

where we have used the approximations $\theta_{\rm A}(t) \approx \Omega t$ and $\theta_{\rm B}(t) \approx \Omega t +\delta$. The stability of δ can be discussed in terms of $\dot{\delta}=\dot{\theta}_{\rm B} - \dot{\theta}_{\rm A}$ that is given by

$$\begin{align} \frac{\dot{\delta}}{2\alpha} \approx \sin\left[ \frac{(-2+\cos \delta)\sin \delta}{9} \hat{\Omega}^{2}\right] \cos\left[\frac{2+\cos \delta}{3} \hat{\Omega}\right]. \end{align} \tag{ 22 }$$

In fig. 8, we plot $\dot{\delta}/(2\alpha)$ as a function of δ for $\hat{\Omega}=0.05$, 0.1, and 0.2. We first note that eq. (22) is an odd function of δ. Since $\dot{\delta}>0$ for $\delta<0$ and $\dot{\delta}<0$ for $\delta>0$, we find that $\delta=0$ is a stable fixed point. This result is in accordance with the stability diagram in fig. 7 when $\alpha < \alpha_{\rm c}$. The slope at $\delta=0$ becomes steeper as $\hat{\Omega}$ is increased. We also see that $\delta=\pm \pi$ are the unstable fixed points.

Zoom In Zoom Out Reset image size

Up to the first order in $\hat{\Omega}$, we see in eqs. (18) and (19) that the mechanical phases $\phi_\mathrm{A}$ and $\phi_\mathrm{B}$ are delayed with respect to those of the natural length $\theta_\mathrm{A}$ and $\theta_\mathrm{B}$, respectively. Thus, the time evolutions of $\theta_\mathrm{A}$ and $\theta_\mathrm{B}$ are accelerated by the second terms in eqs. (20) and (21) that are controlled by α. Hence, as α is increased, the oscillation frequencies of $\Delta \hat{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 analytically predict the critical value $\alpha_{\rm c}$ nor the stable fixed points at nonzero δ for $\alpha > \alpha_{\rm c}$. The difficulty arises 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 the model.

It is worth mentioning, however, that the stability of the phase difference δ can be determined even in the absence of hydrodynamic interactions as we have discussed in this section. 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 three-sphere microswimmer.

In this letter, we have proposed a model of an autonomous three-sphere microswimmer by considering a coupling effect between the two natural lengths of an elastic microswimmer [9]. Our model 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_\infty$ that is almost independent of the initial phase difference $\delta_0$ (see fig. 5(a)). The corresponding phase difference $\delta_\infty$ between the oscillations in the natural lengths becomes also finite (see fig. 5(b)), which is consistent with Purcell's scallop theorem for microswimmers in a viscous fluid [2,3].

We have explored in detail the dependencies of $V_\infty$ and $\delta_\infty$ on the coupling parameter α and the frequency Ω. We find that both $\vert V_\infty \vert$ and $\vert \delta_\infty \vert$ take nonzero values for $\alpha > \alpha_{\rm c}$, and they also show maximum values at $\alpha_{\rm m}$ (fig. 6). There is a finite range of α for which a microswimmer can have 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., $\hat{\Omega} <1$. When Ω is made larger, the difference between $\theta_{\rm A(B)}$ (phases of the natural lengths) and $\phi_{\rm 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 the high-frequency limit due to the intrinsic spring relaxation dynamics [9,13]. Such a reduction of the velocity in the high-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 case of $\alpha < 0$.

We thank T. Kato for useful discussions. YK acknowledges support by a Grant-in-Aid for JSPS Fellows (Grant No. JP19J00365) from the Japan Society for the Promotion of Science (JSPS). YH acknowledges support by a Grant-in-Aid for JSPS Fellows (Grant No. 19J20271) from the 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) from the JSPS. SK further acknowledges support by a Grant-in-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 Ministry of Education, Culture, Sports, Science and Technology of Japan.