# Mathematical modeling for the synchronization of two interacting active
rotors

Hiroyuki Kitahata<EMAIL_ADDRESS>Department of Physics, Graduate School
of Science, Chiba University, Chiba 263-8522, Japan  Yuki Koyano
<EMAIL_ADDRESS>Department of Human Environmental Science, Graduate
School of Human Development and Environment, Kobe University, Kobe 657-0011,
Japan

###### Abstract

We investigate the synchronization of active rotors. A rotor is composed of a
free-rotating arm with a particle that releases a surface-active chemical
compound. It exhibits self-rotation due to the surface tension gradient
originating from the concentration field of the surface-active compound
released from the rotor. In a system with two active rotors, they should
interact through the concentration field. Thus, the interaction between them
does not depend only on the instantaneous positions but also on the dynamics
of the concentration field. By numerical simulations, we show that in-phase
and anti-phase synchronizations occur depending on the distance between the
two rotors. The stability of the synchronization mode is analyzed based on
phase reduction theorem through the calculation of the concentration field in
the co-rotating frame with the active rotor. We also confirm that the
numerical results meet the prediction by theoretical analyses.

††preprint: APS/123-QED

##  I Introduction

Self-propelled particles have been intensively investigated for decades both
experimentally and theoretically [1, 2]. Motions of living organisms such as
cells, bacteria, fish, birds, and insects attract much interest as examples of
the self-propulsion. For such motions, several mechanisms on the motion are
suggested, e.g., hydrodynamic interaction due to the surface deformation [3,
4], momentum exchange with the substrate [5, 6], and the tactic motions [7, 8,
9]. As for the tactic motions, they are classified into several types such as
chemotaxis, phototaxis, mechanotaxis, geotaxis, and so on. Here, we focus on
the chemotactic motion, in which the direction of motion is determined by the
concentration gradient. For positive and negative chemotaxes, the object moves
in the positive and negative directions of the gradient of the concentration
field, respectively. If the object releases a chemical compound around itself
and exhibits negative chemotaxis, then the rest state where the object stands
still can become unstable since the object is likely to move away from the
original position with higher concentration. The motion can be sustained since
the particle motion can keep the anisotropy in the concentration field around
the particle. This is one of the mechanism for the self-propulsion with the
taxis [10].

An experimental example for such self-propulsion with negative chemotaxis is a
camphor particle floating at a water surface. The camphor particle releases
the camphor molecules at the water surface, and the molecules reduce the
surface tension. The object is pulled toward the region with higher surface
tension reflecting lower camphor concentration, which can be understood as
negative chemotaxis [11, 12, 13, 14, 15, 16]. Several types of active rotors,
or self-propelled rotors, were recently reported using camphor or some other
chemicals with surface activity [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27,
28, 29]. For example, an elliptic camphor paper can rotate around an axis
penetrating a hole at the center of the paper [17, 18]. We also reported a
camphor rotor, in which a plastic plate with two or more camphor particles can
freely rotate around an axis penetrating at the center of the plate [19, 20].
Some other types of rotors using camphor or other surface-active compounds
have also been reported.

As theoretical approaches for the camphor particle motion, the reaction-
diffusion equation for camphor concentration coupled with the Newtonian
equation for a camphor particle motion has often been adopted [30, 31, 14,
15]. We previously discussed the bifurcation of a camphor rotor by considering
the time evolution for the camphor concentration coupled with the Newtonian
equation for the rotation of the arm attached with camphor particles [19, 20].
We theoretically derived the simplified ordinary differential equation on the
angle of the camphor rotor using the perturbation method and showed that the
self-propelled rotation emerges through supercritical pitchfork bifurcation by
changing the friction coefficient as a bifurcation parameter. Sharma and co-
authors recently reported the results of experiments and numerical simulation
based on the simple mathematical model. In their experimental system, the two
rotors made of rectangular camphor papers were closely located, and they
observed the desynchronization and anti-phase synchronization in our
definition, i.e., the particles alternately come close to the other rotor. In
their numerical simulation, they assume that a camphor rotor has an intrinsic
angular velocity and interacts with the other rotor by a Yukawa-type potential
depending on their relative position. They succeeded in reproducing their
experimental results [21]. They also investigated the interaction between the
multiple rotors and reported many interesting states such as synchronous,
quasiperiodic, and chaotic states [22, 23, 24, 25].

As mentioned above, a single rotor composed of the particle with a surface-
active compound rotates due to the surface tension gradient at the water
surface, and the surface tension is a function of the concentration of the
compound. In the typical experiments with rotors with camphor particles, the
rotation period of the rotors is of the order of 1 s and the effective
diffusion length of the concentration field is several tens of millimeters.
Thus, the characteristic time for the chemical compound to diffuse to the
other rotor is comparable to the rotation periods. Therefore, we consider that
the interaction through the concentration field is important when the distance
between rotors is large. Thus, here we discuss the dynamics of camphor rotors
including the time evolution of the concentration field. We found that the in-
phase or anti-phase synchronization mode occurs depending on the distance
between the rotors. In addition, we succeeded in the mathematical analysis
based on the phase description.

Our manuscript is constructed as follows: We first describe the mathematical
model for the active rotors in Sec. II and then show the results of numerical
simulation based on the model in Sec. III. Then, we perform the theoretical
analyses using the phase reduction. The procedure and results of the analyses
are described in Sec. IV. Then, we discuss the validity of the phase
description by directly calculating the phase coupling function using
numerical simulation in Sec. V. Finally, we summarize the results and show a
possible extension of our model in Sec. VI.

##  II Mathematical Model

We construct a mathematical model for a system with symmetric active rotors,
which can rotate in either clockwise or counterclockwise direction, based on
the previous studies [31, 14, 15]. We mainly consider a two-rotor system, but
also consider a single-rotor one to clarify the characteristics of a composing
rotor. Our model comprises the time-evolution equation for the configurations
of the rotors with particles (camphor particles) that release a surface-active
compound (camphor molecules), and that for the concentration of a surface-
active compound.

Figure 1:  Schematic illustration for the active rotors and definition of
phases. (a) A single rotor. (b) Two coupled rotors which rotate in the same
direction (counterclockwise). (c) Two coupled rotors, which rotate in
different directions. The first (left) and second (right) ones rotate
counterclockwise and clockwise, respectively.

The $i$th particle can move along a circle with a radius of $a$ and the center
position of $\bm{\ell}_{i}$. Therefore, the particle position can be described
only by using one variable $\phi_{i}$, which is called the phase of the $i$th
particle. We define the origin and positive direction of each phase as
schematically shown in Fig. 1. For a single-rotor system as in Fig. 1(a), we
set $\bm{\ell}_{1}=\bm{0}$ and thus we can express the position of the
particle position $\bm{r}_{1}$ as

$\displaystyle\bm{r}_{1}=\bm{\ell}_{1}+a\bm{e}(\phi_{1})=a\bm{e}(\phi_{1}).$
(1)

For a two-rotor system as in Figs. 1(b) and (c), we set
$\bm{\ell}_{1}=-(L/2)\bm{e}_{x}$ and $\bm{\ell}_{2}=(L/2)\bm{e}_{x}$. For the
symmetric expression between the first and second rotors on the time evolution
of each phase, the origins of the phases are set so that $\phi_{i}=0$
corresponds to the direction toward the center of the other rotor, and
positive directions of the phase are set as the rotation direction of each
rotor. That is to say, the positions of the two particles $\bm{r}_{1}$ and
$\bm{r}_{2}$ are expressed using the phases $\phi_{1}$ and $\phi_{2}$ as

$\displaystyle\bm{r}_{1}=\bm{\ell}_{1}+a\bm{e}(\phi_{1})=-\dfrac{L}{2}\bm{e}_{x}+a\bm{e}(\phi_{1}),$
(2)
$\displaystyle\bm{r}_{2}=\bm{\ell}_{2}-a\bm{e}(\pm\phi_{2})=\dfrac{L}{2}\bm{e}_{x}-a\bm{e}(\pm\phi_{2}).$
(3)

In Eq. (3), the positive and negative signs correspond to the rotation in the
same direction (Fig. 1(b)) and in the opposite direction (Fig. 1(c)),
respectively. Here, we set the Cartesian coordinates so that the origin meets
the midpoint of the centers of the two rotors and the line connecting the
centers of the two rotors meets the $x$ axis. The unit vector in the $x$ and
$y$ axes is set as $\bm{e}_{x}$ and $\bm{e}_{y}$, respectively, and
$\bm{e}(\theta)$ is a unit vector in the direction of $\theta$, i.e.,
$\bm{e}(\theta)=\cos\theta\bm{e}_{x}+\sin\theta\bm{e}_{y}$.

The time evolution equation for the rotor is obtained based on the Newtonian
equation with the overdamped scheme. That is to say, the equation is written
as

$\displaystyle\eta_{i}A\frac{d\bm{r}_{i}}{dt}=\eta_{i}A\bm{v}_{i}=\bm{F}_{u,i}+\bm{F}_{c,i},$
(4)

where $\bm{v}_{i}$ is the velocity of the particle composed of the $i$th
rotor, $\eta_{i}$ is the friction coefficient per area for the $i$th rotor,
and $A$ is the area of the particle. The velocity $\bm{v}_{i}$ is expressed
using the phase $\phi_{i}$ as

$\displaystyle\bm{v}_{1}=a\frac{d\phi_{1}}{dt}\bm{e}\left(\phi_{1}+\frac{\pi}{2}\right),$
(5)
$\displaystyle\bm{v}_{2}=a\frac{d\phi_{2}}{dt}\bm{e}\left(\phi_{2}\pm\frac{\pi}{2}\right),$
(6)

where the positive and negative signs correspond to the cases with the same
and opposite rotation directions, respectively. $\bm{F}_{u,i}$ and
$\bm{F}_{c,i}$ are the force exerted on the $i$th particle due to the surface
tension gradient and the constraint force in the direction of
$\bm{e}(\phi_{i})$. The force due to the surface tension gradient is expressed
using area integration as

$\displaystyle\bm{F}_{u,i}=$
$\displaystyle\oint_{\partial\Omega_{i}}\left[-\Gamma u\bm{n}\right]d\ell$
$\displaystyle=$ $\displaystyle\iint_{\Omega_{i}}\left[-\Gamma\nabla
u\right]dA$ $\displaystyle=$
$\displaystyle-\Gamma\iint_{\mathbb{R}^{2}}\left(\nabla
u\right)S\left(\bm{r}-\bm{r}_{i}\right)dA,$ (7)

where $u$ is the concentration field of the surface-active chemical compound,
$\Omega_{i}$ is the region of the particle composing the $i$th rotor,
$\partial\Omega_{i}$ is the periphery of $\Omega_{i}$, and $\bm{n}$ is the
outward normal unit vector at the particle periphery. $dA$ and $d\ell$ are the
area and line elements. The surface tension should be a decreasing function of
the concentration of the surface-active compound. Here, we assume a linearity
between surface tension and the concentration, where the proportionality
constant is $-\Gamma$. $S\left(\cdot\right)$ is a level function. In the
numerical simulation, we adopt a smoothed step function which has values close
to unity inside the particle and values close to zero outside of it, i.e.,

$\displaystyle
S(\bm{\ell})=\left\\{\begin{array}[]{ll}1,&\left|\bm{\ell}\right|\leq
R-\varepsilon,\\\
1-\left(\left|\bm{\ell}\right|-R+\varepsilon\right)^{2}/\left(2\varepsilon^{2}\right),&R-\varepsilon<\left|\bm{\ell}\right|\leq
R,\\\
\left(\left|\bm{\ell}\right|-R-\varepsilon\right)^{2}/\left(2\varepsilon^{2}\right),&R<\left|\bm{\ell}\right|<R+\varepsilon,\\\
0,&\left|\bm{\ell}\right|\geq R+\varepsilon,\end{array}\right.$ (12)

where $R$ is the radius of a particle of the surface-active compound and
$\varepsilon$ is a smoothing factor. Then, we obtain the evolution equation by
calculating the vector product of Eq. (4) with $a\bm{e}(\phi_{i})$ as

$\displaystyle\eta_{1}Aa^{2}\frac{d\phi_{1}}{dt}=$
$\displaystyle\left(\bm{r}_{1}-\bm{\ell}_{1}\right)\times\bm{F}_{1}$
$\displaystyle=$ $\displaystyle
a\bm{e}(\phi_{1})\times\left(\bm{F}_{u,1}+\bm{F}_{c,1}\right)$
$\displaystyle=$ $\displaystyle
a\bm{e}\left(\phi_{1}+\frac{\pi}{2}\right)\cdot\bm{F}_{u,1},$ (13)

and

$\displaystyle\eta_{2}Aa^{2}\frac{d\phi_{2}}{dt}=$
$\displaystyle\pm\left(\bm{r}_{2}-\bm{\ell}_{2}\right)\times\bm{F}_{2}$
$\displaystyle=$ $\displaystyle\mp
a\bm{e}(\pm\phi_{2})\times\left(\bm{F}_{u,2}+\bm{F}_{c,2}\right)$
$\displaystyle=$ $\displaystyle\mp
a\bm{e}\left(\pm\phi_{2}+\frac{\pi}{2}\right)\cdot\bm{F}_{u,2}$
$\displaystyle=$ $\displaystyle
a\bm{e}\left(\pm\left(\phi_{2}-\frac{\pi}{2}\right)\right)\cdot\bm{F}_{u,2}.$
(14)

Here, the upper and lower signs correspond to the rotation in the same
rotation direction and the opposite rotation direction, respectively. The
operator “$\times$” denotes the vector product in two dimensions, i.e.,
$\bm{\alpha}\times\bm{\beta}=\alpha_{x}\beta_{y}-\alpha_{y}\beta_{x}$ for
$\bm{\alpha}=\alpha_{x}\bm{e}_{x}+\alpha_{y}\bm{e}_{y}$ and
$\bm{\beta}=\beta_{x}\bm{e}_{x}+\beta_{y}\bm{e}_{y}$. In the calculation, we
used that the constraint force $\bm{F}_{c,1}$ is parallel to
$\bm{e}(\phi_{1})$ and thus $\bm{F}_{c,1}\times\bm{e}(\phi_{1})=0$. We also
used that $\bm{F}_{c,2}$ is parallel to $\bm{e}(\pm\phi_{2})$ and thus
$\bm{F}_{c,2}\times\bm{e}(\pm\phi_{2})=0$.

The dynamics of the concentration field is described as

$\displaystyle\frac{\partial u}{\partial
t}=\nabla^{2}u-u+\frac{1}{A}\sum_{i=1}^{N}S(\bm{r}-\bm{r}_{i}),$ (15)

where the first, second, and third terms on the right side correspond to the
diffusion, evaporation, and supply of the surface-active chemical compound.
$S(\bm{r}-\bm{r}_{i})/A$ denotes the supply of the surface-active compound
from the $i$th particle located at $\bm{r}_{i}$, and $N=1$ for a single-rotor
system and $N=2$ for a two-coupled rotor system.

It should be noted that we used the equations with dimensionless variables.
The length, time, and concentration are normalized with the diffusion length,
$\sqrt{D/\kappa}$, the characteristic time of sublimation, $1/\kappa$, and the
ratio between the supply rate and sublimation rate, $f/\kappa$. Here, $D$ is
the effective diffusion constant of the surface-active chemical compound [32,
33], $\kappa$ is the sublimation rate, and $f$ is the supply rate of the
compound for each disk.

##  III Numerical simulation

Numerical simulation was performed based on the model in Sec. II. For the
numerical simulation, we adopted the Euler method for the dynamics of the
rotors in Eqs. (13) and (14) and with an explicit method for the dynamics of
the concentration in Eq. (15). The program code was prepared by ourselves in
the C language. The calculation was performed in the region with $-X/2\leq
x\leq X/2$ and $-Y/2\leq y\leq Y/2$. The Robin boundary condition $\nabla
u\cdot\bm{n}_{b}+u=0$ was adopted [34], where $\bm{n}_{b}$ is the outward
normal unit vector at the boundary. The initial condition was set as $u=0$ in
all of the calculation region. In order to stabilize the rotation direction,
we fixed $d\phi_{i}/dt=20$ for $0\leq t\leq 1$. Due to the asymmetry in the
concentration field around the particle maintained by the rotational motion,
we succeeded in obtaining the stably rotating rotors. $\phi_{1}$ and
$\phi_{2}$ at $t=0$ were set to be $0$ and $\pi/2$, respectively, so that the
convergence to the in-phase and anti-phase synchronization modes could be
easily obtained. The parameters were fixed as $\Gamma=1$, $R=0.1$, $a=0.2$,
and $\varepsilon=0.025$. The spatial mesh and time step were $\Delta x=0.025$
and $\Delta t=0.0001$. The calculation region was $X=16$ and $Y=10$.

Figure 2:  Numerical results for a single rotor. (a) Time series of
$d\phi_{1}/dt$ for $\eta_{1}=0.05$ (red), $0.1$ (green), $0.15$ (cyan) and
$0.2$ (blue). (b) The stable angular velocity $\omega$ depending on the
friction coefficient $\eta_{1}$. A single rotor exhibited a rotation with a
finite angular velocity for $\eta_{1}<\eta_{c}\simeq 0.17$, while it stopped
for $\eta_{1}>\eta_{c}$.

First, we performed a numerical simulation for a single rotor. In Fig. 2(a),
we show the time series of $d\phi_{1}/dt$ for $\eta_{1}=0.05$, $0.1$, $0.15$,
and $0.2$. For $\eta_{1}=0.05$, $0.1$, and $0.15$, $d\phi_{1}/dt$ converged to
a finite positive value, while $d\phi_{1}/dt$ decayed to zero for
$\eta_{1}=0.2$. We also confirmed that at $t=10$, $d\phi_{1}/dt$ reached a
steady value. Therefore, we defined a stable angular velocity $\omega$ as
$d\phi_{1}/dt$ at $t=100$. In Fig. 2(b), the plot of $\omega$ against the
friction coefficient $\eta_{1}$ is shown. The results suggest that a single
rotor exhibited a stable rotation at a finite constant angular velocity for
$\eta_{1}<\eta_{c}$, while it stopped for $\eta_{1}>\eta_{c}$, where
$\eta_{c}\simeq 0.17$. This can be understood as a supercritical pitchfork
bifurcation just as shown in the previous study [19].

Figure 3: Snapshots representing the particle position and camphor
concentration at $t=1000$, at which the coupled rotors for $L=2$ and $3$ reach
the stable synchronization mode. (a) $L=2$ and (b) $L=3$ in the case with the
same rotation direction. (c) $L=2$ and (d)$L=3$ in the opposite rotation
direction. The yellow arrows show the rotation direction. The cross points of
the yellow dotted lines correspond to the centers of the rotors,
$\bm{\ell}_{1}$ and $\bm{\ell}_{2}$. Figure 4:  Numerical results for coupled
rotors. Time series of $\Delta\phi$ are plotted for (a) the case with the same
rotation direction and (b) the case with the opposite rotation direction. The
distance between the two rotor centers $L$ was (a) $L=2$ (red), (b) $L=3$
(green), and (c) $L=4$ (cyan). In both cases, in-phase synchronization was
observed for $L=2$ and $4$, while anti-phase synchronization was observed for
$L=3$.

Then, we fixed $\eta_{1}=\eta_{2}=0.1$ and calculated the behavior of the
coupled system. We demonstrated the two cases: (i) the two rotors rotate in
the same direction (cf. Fig. 1(b)) and (ii) they do so in the opposite
directions (cf. Fig. 1(c)). In order to clarify the mode of synchronization,
we detected the time $\tau^{(i)}_{\mu}$ at which the rotor $i$ passes through
the line segment connecting $\bm{\ell}_{1}$ and $\bm{\ell}_{2}$ for the
$\mu$th time. Then, the phase difference is defined as

$\displaystyle\Delta\phi=2\pi\frac{\tau^{(2)}_{\mu}-\tau^{(1)}_{\nu}}{\tau^{(1)}_{\nu+1}-\tau^{(1)}_{\nu}},$
(16)

where $\mu$ and $\nu$ holds
$\tau^{(1)}_{\nu}\leq\tau^{(2)}_{\mu}<\tau^{(1)}_{\nu+1}$.

We changed $L$ and calculated the dynamics for the coupled rotors until
$t=10000$. The snapshots after the synchronized state becomes stable are shown
for $L=2$ and $3$ in each rotation direction in Fig. 3. The time evolution of
$\Delta\phi$ is shown in Fig. 4. In both cases with the same and opposite
rotation directions, the in-phase synchronization ($\Delta\phi=0$) was
observed for $L=2$ and $4$, while the anti-phase synchronization
($\Delta\phi=\pi$) was observed for $L=3$. The relaxation time to the
synchronized state was intensively dependent on the distance $L$ between the
two rotors.

In order to clearly show the dependence of the synchronization mode on the
distance $L$, we plotted $\Delta\phi$ at $t=5000$, $10000$, $15000$, and
$20000$ against $L$ in Fig. 5. It clearly exhibits that the in-phase and anti-
phase synchronization alternates with an increase in $L$. The preferable
synchronization mode was almost the same in both cases. In Fig. 5, the phase
difference was not converged for the region greater than 4 and those close to
the transition points between the in-phase and anti-phase synchronization.
This should be because in these regions the interaction between the two rotors
was so small that it takes much time to reach the stable synchronized state of
the rotors. In such a case, synchronization cannot be observed in actual
systems due to the inevitable fluctuation.

Figure 5:  Numerical results for coupled rotors. The phase differences
$\Delta\phi$ at $t=5000$ (light green), $10000$ (yellow), $15000$ (orange),
and $20000$ (red) are plotted in (a) the case with the same rotation direction
and (b) the case with the opposite rotation direction. In-phase
synchronization ($\Delta\phi=0$) and anti-phase synchronization
($\Delta\phi=\pi$) alternate with an increase in $L$. Figure 6:  Numerical
results on the synchronization for two coupled rotors with different intrinsic
angular velocities by setting $\eta_{1}\neq\eta_{2}$ with $L=2$. Here, we
fixed $\eta_{1}=0.1$ and varied $\eta_{2}$. (a), (b) The difference
$\Delta\omega$ in the averaged angular velocities depending on
$\eta_{2}-\eta_{1}$ in the cases with (a) the same and (b) the opposite
rotation directions. (c), (d) Phase difference for $\Delta\phi$ when the two
rotors synchronized ($\left|\Delta\omega\right|<0.01$) depending on
$\eta_{2}-\eta_{1}$ in the cases with (c) the same and (d) the opposite
rotation directions.

In the above paragraphs, we only showed the numerical simulation results on
the synchronization of the two identical rotors. In the synchronization of
nonlinear oscillators, it is known that the synchronization can be observed
for the two oscillators with slightly different intrinsic angular velocities
[35, 36]. Thus, we considered the two slightly different rotors by fixing
$\eta_{1}=0.1$ and varying $\eta_{2}$ when $L=2$. We measured the averaged
angular velocities $\omega_{1}$ and $\omega_{2}$ over $100<t<200$ and
calculated the difference between them $\Delta\omega=\omega_{2}-\omega_{1}$.
The plots of $\Delta\omega$ against $\eta_{2}-\eta_{1}$ are shown in Figs. 6
(a) and (b) in the case with the same and opposite rotation directions,
respectively. We observed the synchronization ranges in
$-0.0005\lesssim\eta_{2}-\eta_{1}\lesssim 0.0005$, where $\Delta\omega=0$. We
also measured the phase difference $\Delta\phi=\phi_{2}-\phi_{1}$ when the two
rotors synchronized. As shown in Figs. 6(c) and (d), the phase differences
$\Delta\phi$ decreased with an increase in $\eta_{2}-\eta_{1}$, and they are 0
when $\eta_{2}=\eta_{1}$. This means that perfect in-phase synchronization was
realized for the two identical rotors, while in-phase synchronization with a
slight phase shift was realized for the two rotors with slightly different
intrinsic angular velocities. These results suggest that the behavior of the
coupled two-rotor system is translated in terms of the synchronization of the
two nonlinear oscillators.

##  IV Theoretical analysis

In order to discuss the mechanism of the alternation of the stable
synchronization modes depending on $L$, we perform the theoretical analysis to
discuss the synchronization of the two active rotors based on the phase
reduction method. Hereafter, we only consider the case with two identical
rotors, i.e., $\eta_{1}=\eta_{2}=\eta$. The model in Eqs. (1)–(15) is used,
but a point source is adopted, i.e. $R\to+0$ as

$\displaystyle\frac{\partial u}{\partial
t}=\nabla^{2}u-u+\sum_{i=1}^{N}\delta(\bm{r}-\bm{r}_{i}),$ (17)

in the place of Eq. (15), where $\delta(\cdot)$ is the Dirac’s delta function.

Since the time-evolution equation for the concentration field $u$ in Eq. (15)
is linear, $u$ is described as the summation of $u_{1}$ and $u_{2}$, which
originate from the supply from the particles 1 and 2, respectively. That is to
say,

$\displaystyle u=u_{1}+u_{2},$ (18)

where

$\displaystyle\frac{\partial u_{i}}{\partial
t}=\nabla^{2}u_{i}-u_{i}+\delta(\bm{r}-\bm{r}_{i}),$ (19)

for $i=1,2$.

As for the time evolution of $\phi_{i}$, Eqs. (13) and (14) give

$\displaystyle\eta a^{2}\frac{d\phi_{1}}{dt}=$ $\displaystyle
a\bm{e}\left(\phi_{1}+\frac{\pi}{2}\right)\cdot\frac{1}{A}\bm{F}_{u,1},$ (20)

and

$\displaystyle\eta a^{2}\frac{d\phi_{2}}{dt}=$ $\displaystyle
a\bm{e}\left(\pm\left(\phi_{2}-\frac{\pi}{2}\right)\right)\cdot\frac{1}{A}\bm{F}_{u,2}.$
(21)

The force originating from the surface tension $\bm{F}_{u,i}$ is also
decomposed into two terms:

$\displaystyle\bm{F}_{u,i}=\bm{F}_{u,i,1}+\bm{F}_{u,i,2},$ (22)

in the same way as in the concentration $u$. Under the point-source
approximation,

$\displaystyle\frac{1}{A}\bm{F}_{u,i,j}=$
$\displaystyle-\frac{\Gamma}{A}\iint_{\mathbb{R}^{2}}\left(\nabla
u_{j}\right)\delta\left(\bm{r}-\bm{r}_{i}\right)dA,$ $\displaystyle\to$
$\displaystyle-\Gamma\left.\nabla u_{j}\right|_{\bm{r}=\bm{r}_{i}},$ (23)

for $i\neq j$. It should be noted that the expression in Eq. (23) for $i=j$
does not hold since the force $\bm{F}_{u,i,i}/A$ shows the logarithmic
divergence and we should introduce a small positive value corresponding to the
particle radius [19, 37]. Nevertheless, from the physical point of view, the
$i$th rotor should rotate with a constant angular velocity for $t\to\infty$
when it is driven only by $\bm{F}_{u,i,i}$. We define the terminal angular
velocity to be $\omega$ and then the following equations hold:

$\displaystyle\frac{1}{\eta
Aa}\bm{e}\left(\phi_{1}+\frac{\pi}{2}\right)\cdot\bm{F}_{u,1,1}=\omega,$ (24)

for rotor 1, and

$\displaystyle\frac{1}{\eta
Aa}\bm{e}\left(\pm\left(\phi_{2}-\frac{\pi}{2}\right)\right)\cdot\bm{F}_{u,2,2}=\omega,$
(25)

for rotor 2. The positive and negative signs correspond to the same and
opposite rotation directions, respectively.

Hereinafter, the effect of the concentration field of the surface-active
compound released from the other rotor is treated as a perturbation. We first
calculate the concentration field of the chemical compound released from one
rotor and then the force originating from the concentration field exerting on
the other particle.

For the construction of the concentration field generated by one camphor
rotor, we consider that a single rotor is rotating at a constant angular
velocity $\omega$, that is, the position is described as
$\bm{r}=a\bm{e}(\phi_{1})=a\bm{e}(\omega t+\phi_{0})$. We introduce a co-
rotating frame with the rotor, where the variables in the frame are expressed
with tildes. The single point source is located at
$\tilde{\bm{r}}=a\tilde{\bm{e}}_{x}$ and the concentration field $\tilde{u}$
in the co-rotating frame should be in a steady state. Here,
$\tilde{\bm{e}}_{x}$ is a unit vector in the $\tilde{x}$ direction in the co-
rotating frame. Then, the steady-state concentration field should hold

$\displaystyle-\omega\frac{\partial\tilde{u}}{\partial\tilde{\theta}}=\tilde{\nabla}^{2}\tilde{u}-\tilde{u}+\delta\left(\tilde{\bm{r}}-a\tilde{\bm{e}}_{x}\right),$
(26)

where $\tilde{\nabla}$ is the nabla operator in the co-rotating frame. After
lengthy calculation, we obtain

$\displaystyle\tilde{u}(\tilde{r},\tilde{\theta})=\left\\{\begin{array}[]{ll}\tilde{u}_{\mathrm{in}}(\tilde{r},\tilde{\theta}),&\tilde{r}<a,\\\
\tilde{u}_{\mathrm{out}}(\tilde{r},\tilde{\theta}),&\tilde{r}\geq
a,\end{array}\right.$ (29)

where

$\displaystyle\tilde{u}_{\rm in}(\tilde{r},\tilde{\theta})$
$\displaystyle\qquad=\frac{1}{2\pi}\sum_{n=-\infty}^{\infty}\mathcal{K}_{n}\left(a\sqrt{1-in\omega}\right)\mathcal{I}_{n}\left(\tilde{r}\sqrt{1-in\omega}\right)e^{in\tilde{\theta}},$
(30)

and

$\displaystyle\tilde{u}_{\rm out}(\tilde{r},\tilde{\theta})$
$\displaystyle\qquad=\frac{1}{2\pi}\sum_{n=-\infty}^{\infty}\mathcal{I}_{n}\left(a\sqrt{1-in\omega}\right)\mathcal{K}_{n}\left(\tilde{r}\sqrt{1-in\omega}\right)e^{in\tilde{\theta}}.$
(31)

Here, $\mathcal{I}_{n}(\cdot)$ and $\mathcal{K}_{n}(\cdot)$ are the modified
Bessel functions of the first and second kinds with the degree of $n$,
respectively. The detailed derivation and notes for the Bessel function with
complex parameters are shown in Appendix A.

Based on Eq. (31), we obtain the asymptotic form of the concentration field
far from the rotor. Here we consider $\tilde{r}\gg 1\gg a$. That is to say, we
consider the case that the distance between two rotors is much greater than
the diffusion length and the arm of the rotor is much less than the diffusion
length. We take into account the order up to $\mathcal{O}(a)$. Using Eq. (31),
we obtain the asymptotic form as

$\displaystyle\tilde{u}(\tilde{r},\tilde{\theta})\simeq$
$\displaystyle\frac{1}{2\sqrt{2\pi\tilde{r}}}e^{-\tilde{r}}+\frac{a}{2\sqrt{2\pi\tilde{r}}}\rho^{1/4}e^{-\tilde{r}\sqrt{\rho}\cos(\chi/2)}$
$\displaystyle\quad\times\cos\left(\tilde{\theta}-\frac{\chi}{4}+\tilde{r}\sqrt{\rho}\sin\frac{\chi}{2}\right)+\mathcal{O}\left(a^{2}\right).$
(32)

Here, we set $1-i\omega=\rho e^{-i\chi}$, that is, $\rho=\sqrt{1+\omega^{2}}$
and $\chi=\arctan\omega$. The detailed calculation is found in Appendix B.

In the laboratory frame, the asymptotic form of the concentration field
generated by a rotor, which is located at the origin, whose phase is $\phi$,
and which rotates in the counterclockwise rotation, is described as

$\displaystyle u(r,\theta,\phi)$ $\displaystyle\simeq\frac{1}{2\sqrt{2\pi
r}}e^{-r}+\frac{a}{2\sqrt{2\pi r}}\rho^{1/4}e^{-r\sqrt{\rho}\cos(\chi/2)}$
$\displaystyle\qquad\times\cos\left(\theta-\phi-\frac{\chi}{4}+r\sqrt{\rho}\sin\frac{\chi}{2}\right)+\mathcal{O}\left(a^{2}\right),$
(33)

for $r\gg 1$ and $t\to\infty$. It should be noted that the first term with
$\mathcal{O}(1)$ does not depend on time and that the second term with
$\mathcal{O}(a)$ depends on time which can induce the synchronization between
multiple rotors.

We adopt the asymptotic form in Eq. (33) for $u_{j}$ to calculate
$\bm{F}_{u,i,j}$ $(i\neq j)$ and consider the interaction between two rotors
in Eq. (23). Considering that the asymptotic form of $u_{j}$ is described as a
function of $\phi_{j}$ and that the position $\bm{r}_{i}$ is a function of
$\phi_{i}$, the force $\bm{F}_{u,i,j}$ is a function of $\phi_{i}$ and
$\phi_{j}$, i.e., $\bm{F}_{u,i,j}(\phi_{i},\phi_{j})$. We assume that the
interaction is so weak that the phase difference hardly changes in one period
$2\pi/\omega$ and we can adopt the averaging method in the phase description
[38].

First, we calculate the time evolution of $\phi_{1}$ from Eqs. (13) and (24)
as

$\displaystyle\frac{d\phi_{1}}{dt}=$ $\displaystyle\frac{1}{\eta
Aa}\bm{e}\left(\phi_{1}+\frac{\pi}{2}\right)\cdot\bm{F}_{u,1}$
$\displaystyle=$ $\displaystyle\omega+\frac{1}{\eta
Aa}\bm{e}\left(\phi_{1}+\frac{\pi}{2}\right)\cdot\bm{F}_{u,1,2}(\phi_{1},\phi_{2}),$
$\displaystyle\simeq$ $\displaystyle\omega+\frac{\omega}{2\pi\eta
Aa}\int_{0}^{2\pi/\omega}\bm{e}\left(\phi_{1}+\frac{\pi}{2}\right)\cdot\bm{F}_{u,1,2}(\phi_{1},\phi_{2})dt$
$\displaystyle=$ $\displaystyle\omega+\frac{1}{2\pi\eta Aa}$
$\displaystyle\times\int_{0}^{2\pi}\bm{e}\left(\phi_{1}+\frac{\pi}{2}\right)\cdot\bm{F}_{u,1,2}(\phi_{1},\phi_{1}+\Delta\phi)d\phi_{1}$
$\displaystyle\equiv$
$\displaystyle\left\\{\begin{array}[]{ll}\omega+G_{s}(\Delta\phi),&(\mathrm{same\,rotation\,direction}),\\\
\omega+G_{o}(\Delta\phi),&(\mathrm{opposite\,rotation\,direction}).\end{array}\right.$
(36)

Here, we set $\Delta\phi=\phi_{2}-\phi_{1}$ and calculate the integral under
the assumption that $\Delta\phi$ is constant. It should be noted that
$G_{s}(\Delta\phi)$ and $G_{o}(\Delta\phi)$ are different since
$\bm{F}_{u,1,2}$ depends on the rotation direction of rotor 2. In the same
manner, we obtain from Eqs. (14) and (25) as

$\displaystyle\frac{d\phi_{2}}{dt}\simeq$
$\displaystyle\omega+\frac{1}{2\pi\eta Aa}$
$\displaystyle\times\int_{0}^{2\pi}\bm{e}\left(\pm\left(\phi_{2}-\frac{\pi}{2}\right)\right)\cdot\bm{F}_{u,2,1}(\phi_{2}-\Delta\phi,\phi_{2})d\phi_{2}$
$\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}\omega+G_{s}(-\Delta\phi),&(\mathrm{same\,rotation\,direction}),\\\
\omega+G_{o}(-\Delta\phi),&(\mathrm{opposite\,rotation\,direction}),\end{array}\right.$
(39)

where the positive and negative signs in the second term on the right side
correspond to the cases with the same and opposite rotation directions,
respectively. $G_{s}(\Delta\phi)$ and $G_{o}(\Delta\phi)$ are the so-called
phase coupling functions in terms of coupled oscillators [36].

In the case with the same rotation direction, the second term on the right
side in Eq. (39) is calculated as

$\displaystyle G_{s}(\Delta\phi)$
$\displaystyle\simeq\frac{\Gamma\rho^{1/4}e^{-\sqrt{\rho}L\cos(\chi/2)}}{4\eta\sqrt{2\pi
L}}$
$\displaystyle\quad\times\left\\{\frac{1}{2L}\sin\left[\Delta\phi+\frac{\chi}{4}-L\sqrt{\rho}\sin\left(\frac{\chi}{2}\right)\right]\right.$
$\displaystyle\quad\left.\quad-\sqrt{\rho}\sin\left[\Delta\phi+\frac{3\chi}{4}-L\sqrt{\rho}\sin\left(\frac{\chi}{2}\right)\right]\right\\}$
$\displaystyle\equiv g_{s}(\Delta\phi),$ (40)

which is plotted in Fig. 7(a). The detailed calculation is shown in Appendix
C. Then, we have

$\displaystyle\frac{d\phi_{1}}{dt}\simeq$
$\displaystyle\omega+g_{s}(\Delta\phi),$ (41)
$\displaystyle\frac{d\phi_{2}}{dt}\simeq$
$\displaystyle\omega+g_{s}(-\Delta\phi).$ (42)

By calculating the difference between two equations, we obtain the time-
evolution equations for the slow dynamics of $\Delta\phi$ as

$\displaystyle\frac{d\Delta\phi}{dt}\simeq$ $\displaystyle
g_{s}(-\Delta\phi)-g_{s}(\Delta\phi)\equiv h_{s}(\Delta\phi),$ (43)

where

$\displaystyle h_{s}(\Delta\phi)=$
$\displaystyle\frac{\Gamma\rho^{1/4}e^{-\sqrt{\rho}L\cos(\chi/2)}}{2\eta\sqrt{2\pi
L}}\sin\Delta\phi$
$\displaystyle\times\left\\{-\frac{1}{2L}\cos\left[-\frac{\chi}{4}+L\sqrt{\rho}\sin\left(\frac{\chi}{2}\right)\right]\right.$
$\displaystyle\left.\quad+\sqrt{\rho}\cos\left[-\frac{3\chi}{4}+L\sqrt{\rho}\sin\left(\frac{\chi}{2}\right)\right]\right\\}$
$\displaystyle\equiv$
$\displaystyle\frac{\Gamma\rho^{1/4}e^{-\sqrt{\rho}L\cos(\chi/2)}}{2\eta\sqrt{2\pi
L}}C_{s}\sin\Delta\phi$ $\displaystyle\equiv$ $\displaystyle
h_{s}^{(1)}\sin\Delta\phi.$ (44)

As shown in Fig. 7(c), $\Delta\phi=0$ and $\pi$ are fixed points of Eq. (43)
and their stability is determined by the sign of $C_{s}$ since the factors
other than $C_{s}$ are positive. That is to say, the fixed point at
$\Delta\phi=0$ is stable when $C_{s}<0$ and it is unstable when $C_{s}>0$. On
the other hand, the fixed point at $\Delta\phi=\pi$ is unstable when $C_{s}<0$
and it is stable when $C_{s}>0$. Thus, when $C_{s}<0$ and $C_{s}>0$, in-phase
and anti-phase synchronization should occur, respectively. It should be noted
that $C_{s}$ depends only on $\omega$ and $L$ since $\rho$ and $\chi$ are
functions of $\omega$, as shown just below Eq. (32). $C_{s}$ and $h_{s}^{(1)}$
for $\eta=0.1$ are plotted as a function of $L$ in Fig. 8, where $\omega$ is
set to be constant at $13.23$ from the numerical results in Fig. 2.

Figure 7:  Plots of (a) $g_{s}$, (b) $g_{o}$, (c) $h_{s}$, and (d) $h_{o}$
against $\Delta\phi$ for $L=2$, $2.5$, $3$, and $3.5$. We adopt
$\omega=13.23$. Figure 8:  Plots of (a) $C_{s}$, (b) $C_{o}$, (c)
$h_{s}^{(1)}$, and (d) $h_{o}^{(1)}$ against $L$. Here, we adopt
$\omega=13.23$. Positive (red) and negative (cyan) signs represent the
preference of the anti-phase and in-phase synchronization, respectively.

In the case with the opposite rotation direction, we obtain in the same manner
as in the case with the same rotation direction,

$\displaystyle G_{o}(\Delta\phi)$
$\displaystyle\simeq-\frac{\Gamma\rho^{1/4}e^{-\sqrt{\rho}L\cos(\chi/2)}}{4\eta\sqrt{2\pi
L}}$
$\displaystyle\quad\times\left\\{\frac{3}{2L}\sin\left[\Delta\phi+\frac{\chi}{4}-L\sqrt{\rho}\sin\left(\frac{\chi}{2}\right)\right]\right.$
$\displaystyle\left.\quad+\sqrt{\rho}\sin\left[\Delta\phi+\frac{3\chi}{4}-L\sqrt{\rho}\sin\left(\frac{\chi}{2}\right)\right]\right\\}$
$\displaystyle\equiv g_{o}(\Delta\phi),$ (45)

which leads to

$\displaystyle\frac{d\Delta\phi}{dt}\simeq
g_{o}(-\Delta\phi)-g_{o}(\Delta\phi)\equiv h_{o}(\Delta\phi).$ (46)

Here, we have calculated $h_{o}(\Delta\phi)$ as

$\displaystyle h_{o}(\Delta\phi)=$
$\displaystyle\frac{\Gamma\rho^{1/4}e^{-\sqrt{\rho}L\cos(\chi/2)}}{2\eta\sqrt{2\pi
L}}\sin\Delta\phi$
$\displaystyle\times\left\\{\frac{3}{2L}\cos\left[-\frac{\chi}{4}+L\sqrt{\rho}\sin\left(\frac{\chi}{2}\right)\right]\right.$
$\displaystyle\left.\quad+\sqrt{\rho}\cos\left[-\frac{3\chi}{4}+L\sqrt{\rho}\sin\left(\frac{\chi}{2}\right)\right]\right\\}$
$\displaystyle\equiv$
$\displaystyle\frac{\Gamma\rho^{1/4}e^{-\sqrt{\rho}L\cos(\chi/2)}}{2\eta\sqrt{2\pi
L}}C_{o}\sin\Delta\phi$ $\displaystyle\equiv$ $\displaystyle
h_{o}^{(1)}\sin\Delta\phi.$ (47)

The plots of $g_{o}(\Delta\phi)$ and $h_{o}(\Delta\phi)$ are shown in Figs.
7(b) and (d). In this case, $\Delta\phi=0$ and $\pi$ are also the fixed
points. We can discuss the stability of the synchronization mode in the
parallel manner, and thus the sign of the coefficient $C_{o}$ plays an
important role. To exemplify the stable synchronization mode, $C_{o}$ are also
plotted against $L$ in Fig. 8. The signs of $C_{s}$ and $C_{o}$ almost
coincide for fixed $L$, which means that the stable synchronization mode is
the same in the cases with the same and opposite rotation directions for each
$L$.

##  V Discussion

Figure 9:  (a) Snapshot of the concentration field for a single rotor rotating
at a constant angular velocity at $t=t_{0}\simeq 100.1436$, when $\phi_{1}=0$
first holds after $t=100$. (b) Averaged concentration field $\bar{u}$ over a
period. (c) Difference $\Delta u$ between the concentration field $u$ in (a)
from the averaged field $\bar{u}$ in (b). (d) Enhanced profile of (c), where
the color range for concentration is magnified. The region with the size
$8\times 8$ is shown.

Here, we discuss the mechanism on the synchronization of the coupled rotors
based on the theoretical results. The concentration field of the surface-
active compound that one rotor releases is expressed in Eq. (33). The first
term on the right side does not depend on the phase, but the second term does.
The effect of the dynamics of the other rotor is approximately obtained by
averaging the effect over one period. In such an averaging process, the effect
of the first term is canceled out and only the second term matters. This means
that the periodically changing concentration field only affects the stability
of the synchronization mode. To visualize the time-dependent component of the
concentration field, we numerically calculated the averaged concentration
field $\bar{u}$,

$\displaystyle\bar{u}(x,y)=\frac{1}{T}\int_{t_{0}}^{t_{0}+T}u(x,y,t)dt,$ (48)

where $t_{0}$ was the time when the rotor motion reached the stationary state
and corresponded to $\phi_{1}\simeq 0$. Figure 9 shows the plot of $\Delta
u(x,y,t_{0})$, where $\Delta u(x,y,t)=u(x,y,t)-\bar{u}(x,y)$. The time-
dependent component $\Delta u$ has a spiral structure, which is consistent
with the second term on the right side of Eq. (33). The pitch of the spiral
$L_{0}$ is almost double of the length $L$ for which the stable
synchronization mode changes. Considering that such spiral structure mainly
comes from the second term on the right-hand side in Eq. (33), the pitch
$L_{0}$ is estimated from

$\displaystyle L_{0}\sqrt{\rho}\sin\frac{\chi}{2}=2\pi.$ (49)

Actually, $L_{0}$ is calculated as 2.54 with $\omega=13.23$. This value for
$L_{0}$ well corresponds to both the results by numerical simulations and
theoretical analyses, where the stability of the synchronization mode changes
every $\simeq 1.2$ in $L$.

Figure 10:  Plots of (a) $G_{s}$, (b) $G_{o}$, (c) $H_{s}$, and (d) $H_{o}$
against $\Delta\phi$. The results with $L=2$ (red), $L=2.5$ (orange), $L=3$
(green), and $L=3.5$ (cyan) are shown in each panel.  Figure 11:  Plots of
$\hat{H}_{s}^{(1)}$ and $\hat{H}_{o}^{(1)}$ against $L$. The positive and
negative values are indicated with red and cyan points, respectively. The
positive $\hat{H}_{s}^{(1)}$ or $\hat{H}_{o}^{(1)}$ means that the anti-phase
synchronization state is stable, while the negative $\hat{H}_{s}^{(1)}$ or
$\hat{H}_{o}^{(1)}$ means that the in-phase synchronization state is stable.

Next, we directly calculated the phase coupling function using the numerical
simulation in order to justify the approximation adopted in the theoretical
analysis. We directly calculated the functions $G_{s}(\Delta\phi)$ and
$G_{o}(\Delta\phi)$ in Eqs. (36) and (39). $H_{s}(\Delta\phi)$ and
$H_{o}(\Delta\phi)$ are defined as

$\displaystyle H_{i}(\Delta\phi)=G_{i}(-\Delta\phi)-G_{i}(\Delta\phi)$ (50)

where $i$ denotes $s$ or $o$. The plots for $G_{s}(\Delta\phi)$,
$G_{o}(\Delta\phi)$, $H_{s}(\Delta\phi)$, and $H_{o}(\Delta\phi)$ obtained
from the numerical calculation are shown in Fig. 10. In the calculation, we
first calculated the dynamics of the $i$th rotor only considering the
concentration field released from itself until it reached a stationary angular
velocity, and then calculated the force $\bm{F}_{u,j,i}$ working on the $j$th
particle rotating at the given angular velocity, which was the same as the
$i$th rotor’s. Using the obtained force, $G_{s}(\Delta\phi)$ and
$G_{o}(\Delta\phi)$ were calculated by averaging over a period. The numerical
simulation was performed in the same procedure as in Sec. II. We detected
$\tau^{(1)}_{\nu}$ just after $t=100$ and calculated the average from
$t=\tau^{(1)}_{\nu}$ to $\tau^{(1)}_{\nu+1}$. $G_{i}$ and $H_{i}$ were
calculated for $\Delta\phi=2\pi\lambda/360$ where $\lambda=0,\cdots,359$, and
the averaging was performed for each time step during one period. The
functions $H_{s}(\Delta\phi)$ and $H_{o}(\Delta\phi)$ are odd functions, which
take the value of zero at $\Delta\phi=0,\pi$ and have one positive and one
negative peak. In order to discuss the stability in the synchronization mode,
the slopes at $\Delta\phi=0,\pi$ are important. Therefore, we consider the
Fourier sine expansion of the functions as

$\displaystyle H_{i}(\Delta\phi)=\sum_{k=1}^{\infty}\hat{H}_{i}^{(k)}\sin
k\Delta\phi.$ (51)

where $i$ denotes $s$ or $o$. The first-mode coefficient determines the
stability of the synchronization mode; the in-phase and anti-phase
synchronization is stable for $\hat{H}_{i}^{(1)}<0$ and $\hat{H}_{i}^{(1)}>0$,
respectively. In Fig. 11, $\hat{H}_{s}^{(1)}$ and $\hat{H}_{o}^{(1)}$ are
plotted against $L$, which is close to the plot in Figs. 8(c) and (d) obtained
by theoretical analysis.

In the calculation of the phase coupling function shown in the previous
paragraph, we assume that the two rotors rotate at a constant angular
velocity, i.e., intrinsic angular velocity; however, the angular velocity of
the rotor should be affected by the concentration field originating from the
other rotor. Therefore, we also calculated the phase coupling function
including the effect of the change in the angular velocity, and found that the
time evolution of the phase difference is qualitatively the same as the one
shown in Figs. 10 and 11. The details are shown in Appendix D.

Here, we discuss the parameter set used in the numerical simulation. From the
previous reports on experiments [39], the diffusion length is considered to be
several tens of millimeters. Considering that the unity in the length scale is
set to be the diffusion length, the arm length of the rotor and the distance
between the rotors are estimated to be several millimeters and several
centimeters, respectively. These scales are in the same scale adopted in the
experimental setup [23]. The time scale is normalized by the characteristic
evaporation time, which is around several seconds. The angular velocity of the
rotors and the characteristic time for the relaxation time to the
synchronization mode should be affected by the combination of the particle
size, supply rate of the chemicals, and proportionality coefficient between
the surface tension and concentration. Since these parameters are difficult to
be directly measured from the experiments, the more precise correspondence to
the experimental system should be done as a future study.

In the present system, the stable synchronization mode changes depending on
the distance between the two rotors. There have been several studies of
similar behaviors in the other systems such as the cell thickness pattern in
slime mold [40] and the coupled system of a flickering candle flame [41, 42].
The time delay in the interaction plays an important role in the former case,
while the nonlinear coupling manner is dominant in the latter case. In the
coupling between the camphor rotor discussed in the present paper, the
interaction between two rotors is through the concentration field that obeys
the linear equation, and thus the time delay seems to play an important role
in the present system. Actually, the spiral structure shown in Fig. 9 is the
result of the supply from the rotating rotor and diffusion. Due to the
linearity of the equation, we succeeded to write the time-delay effect
directly through the concentration field, and such time-delay effect is
reduced to the interaction term in the phase dynamics.

##  VI Conclusion

We investigated the coupled active rotors, which spontaneously exhibit
rotation due to the surface tension gradient originating from the surface-
active chemicals released from itself. It has been reported that by
experiments such a coupled rotor system show both in-phase and anti-phase
synchronization, but the mechanism was not fully clarified. We consider the
mathematical model which includes the time evolution of the concentration
field and the motion of the rotors, and obtain the results that the stable
synchronization mode alternates between in-phase and anti-phase
synchronization with an increase in the distance between the two rotors. By
adopting the phase description, which has often been used for the coupled
oscillator systems, we derive the time evolution equation for the phase
difference between the two rotors. The theoretical results suggest the
alternate stable synchronization mode depending on $L$, which well corresponds
to the numerical calculation results. We also directly evaluated the phase
coupling function from the numerical calculation and confirmed that our
theoretical approach works well. We hope that our results will be reproduced
in the experimental systems with two active rotors with well-controlled
intrinsic angular velocities. As for the extension of the present study,
three-or-more rotor systems should be interesting since the system has many
possible stable modes and may exhibit chaotic behaviors. As another extension,
the theoretical description for the case with strong interaction should be
interesting. In such a case, the rotation of each rotor should change to the
reciprocal motion along an arc and such reciprocal motion of the two rotors
can synchronize.

###### Acknowledgements.

This work was supported by JSPS KAKENHI Grants No. JP19K03765, No. JP20K14370,
No. JP20H02712, and No. JP21H01004, and also the Cooperative Research Program
of “Network Joint Research Center for Materials and Devices” (Grants No.
20224003 and No. 20221173). This work was also supported by JSPS and PAN under
the Japan-Poland Research Cooperative Program (Grant No. JPJSBP120204602).

##  Appendix A Derivation of concentration field in a co-rotating frame with a
single rotor

In this section, we derive a steady concentration field in a co-rotating frame
with a single rotor with an angular velocity $\omega$. The positions in the
original system and in the co-rotating system are denoted as
$\bm{r}={}^{t}(r\cos\theta,r\sin\theta)$ and
$\tilde{\bm{r}}={}^{t}(\tilde{r}\cos\tilde{\theta},\tilde{r}\sin\tilde{\theta})$,
respectively. Then, we have

$\displaystyle\tilde{\bm{r}}=\mathcal{R}(-\omega t){\bm{r}},$ (52)

where $\mathcal{R}(\psi)$ is a matrix for rotation in a two-dimensional
system, i.e.,

$\displaystyle\mathcal{R}(\psi)=\left(\begin{array}[]{cc}\cos\psi&-\sin\psi\\\
\sin\psi&\cos\psi\end{array}\right).$ (55)

In other words, we have the relation as

$\displaystyle\tilde{\theta}=\theta-\omega t.$ (56)

Then the operator for the time derivative is rewritten in the $\tilde{\bm{r}}$
system as

$\displaystyle\frac{\partial}{\partial
t}+\omega\left(\begin{array}[]{cc}-\sin\omega t&\cos\omega t\\\ -\cos\omega
t&-\sin\omega t\end{array}\right)\bm{r}\cdot\tilde{\nabla}$ (59)
$\displaystyle=$ $\displaystyle\frac{\partial}{\partial
t}+r\omega\left(\begin{array}[]{cc}\sin(\theta-\omega t)\\\
-\cos(\theta-\omega t)\end{array}\right)\cdot\tilde{\nabla}$ (62)
$\displaystyle=$ $\displaystyle\frac{\partial}{\partial
t}+\tilde{r}\omega\left(\begin{array}[]{cc}\sin\tilde{\theta}\\\
-\cos\tilde{\theta}\end{array}\right)\cdot\tilde{\nabla}$ (65)
$\displaystyle=$ $\displaystyle\frac{\partial}{\partial
t}-\omega\tilde{r}\tilde{\bm{e}}_{\theta}\cdot\tilde{\nabla}.$ (66)

Here $\tilde{\nabla}$ is the nabla operator in the $\tilde{\bm{r}}$ system and
$\tilde{\bm{e}}_{\theta}$ is a unit vector in the co-rotating system,
$\tilde{\bm{e}}_{\theta}={}^{t}(-\sin\tilde{\theta},\cos\tilde{\theta})$.

The position of the point source $\tilde{\bm{a}}$ in the co-rotating frame is
no longer dependent on time and is written as

$\displaystyle\tilde{\bm{a}}=a\tilde{\bm{e}}_{x},$ (67)

where $\tilde{\bm{e}}_{x}$ is a unit vector along the $\tilde{x}$ axis in the
$\tilde{\bm{r}}$ system.

Hereafter, we omit the tilde ($\,\tilde{}\,$) for the variables in the co-
moving frame. In order to consider the steady-state concentration field
$u(r,\theta)$ in the co-rotating system, we set the time derivative to be 0.
Therefore, the equation to be considered is explicitly written as

$\displaystyle-\omega\frac{\partial
u}{\partial\theta}=\nabla^{2}u-u+\delta\left(\bm{r}-\bm{a}\right).$ (68)

The homogeneous equation for (68) is expressed as

$\displaystyle-\omega\frac{\partial u}{\partial\theta}=\nabla^{2}u-u.$ (69)

By assuming that Eq. (69) has a solution in the form of

$\displaystyle u(r,\theta)=f_{n}(r)e^{in\theta},$ (70)

we obtain

$\displaystyle\frac{d^{2}f_{n}}{dr^{2}}+\frac{1}{r}\frac{df_{n}}{dr}-\frac{n^{2}}{r^{2}}f_{n}-(1-in\omega)f_{n}=0.$
(71)

By setting $\hat{r}=r\sqrt{1-in\omega}$, we get

$\displaystyle\frac{d^{2}f_{n}}{d{\hat{r}}^{2}}+\frac{1}{\hat{r}}\frac{df_{n}}{d\hat{r}}-\frac{n^{2}}{{\hat{r}}^{2}}f_{n}-f_{n}=0,$
(72)

which is the so-called modified Bessel equation and the solution is given as

$\displaystyle
f_{n}(\hat{r})=A_{n}\mathcal{I}_{n}\left(\hat{r}\right)+B_{n}\mathcal{K}_{n}\left(\hat{r}\right).$
(73)

Thus, $f_{n}(r)$ is described as

$\displaystyle f_{n}(r)=$ $\displaystyle
A_{n}\mathcal{I}_{n}\left(r\sqrt{1-in\omega}\right)+B_{n}\mathcal{K}_{n}\left(r\sqrt{1-in\omega}\right).$
(74)

Then, the general solution of Eq. (69) is given as

$\displaystyle u(r,\theta)=$
$\displaystyle\sum_{n=-\infty}^{\infty}\left[A_{n}\mathcal{I}_{n}\left(r\sqrt{1-in\omega}\right)\right.$
$\displaystyle\left.\qquad+B_{n}\mathcal{K}_{n}\left(r\sqrt{1-in\omega}\right)\right]e^{in\theta},$
(75)

where $A_{n}$ and $B_{n}$ are complex constants. Considering that $u$ is real,
$A_{n}=A^{*}_{-n}$ and $B_{n}=B^{*}_{-n}$ should hold, where the superscript
“∗” indicates the complex conjugate.

It is notable that the solution of Eq. (71) should be complex. If we set it to
be

$\displaystyle f_{n}(r)=p_{n}(r)+iq_{n}(r),$ (76)

then we obtain

$\displaystyle\frac{d^{2}p_{n}}{dr^{2}}+\frac{1}{r}\frac{dp_{n}}{dr}-\frac{n^{2}}{r^{2}}p_{n}-p_{n}-n\omega
q_{n}=0,$ (77)
$\displaystyle\frac{d^{2}q_{n}}{dr^{2}}+\frac{1}{r}\frac{dq_{n}}{dr}-\frac{n^{2}}{r^{2}}q_{n}-q_{n}+n\omega
p_{n}=0.$ (78)

Considering that the modified Bessel function of the first kind,
$\mathcal{I}_{n}(\cdot)$, is analytic in $\mathbb{C}$, and the modified Bessel
function of the second kind, $\mathcal{K}_{n}(\cdot)$, is analytic in
$\mathbb{C}$ except negative real numbers, the following relations hold:

$\displaystyle\mathcal{I}_{n}\left(z^{*}\right)=\left[\mathcal{I}_{n}(z)\right]^{*},$
(79)
$\displaystyle\mathcal{K}_{n}\left(z^{*}\right)=\left[\mathcal{K}_{n}(z)\right]^{*}.$
(80)

These expressions are derived from Eq. (75) considering that
$\mathcal{I}_{n}(z)$ does not diverge at $\left|z\right|\to 0$ and that
$\mathcal{K}_{n}(z)$ does not diverge at $\left|z\right|\to\infty$ for
$\Re(z)>0$. Now we set

$\displaystyle 1-in\omega=\rho e^{-i\chi_{n}},$ (81)

where $\rho>0$ and $-\pi/2<\chi_{n}<\pi/2$ holds, and define

$\displaystyle\sqrt{1-in\omega}=\sqrt{\rho}e^{-i\chi_{n}/2}.$ (82)

Then, we can show that $\mathcal{K}_{n}\left(r\sqrt{1-in\omega}\right)$ does
not diverge for $r\to\infty$.

Thus, we can describe

$\displaystyle u_{\rm in}(r,\theta)=$ $\displaystyle A_{0}\mathcal{I}_{0}(r)$
$\displaystyle+2\sum_{n=1}^{\infty}\Re\left[A_{n}\mathcal{I}_{n}\left(r\sqrt{1-in\omega}\right)e^{in\theta}\right],$
(83)

for $r<a$, and

$\displaystyle u_{\rm out}(r,\theta)=$ $\displaystyle B_{0}\mathcal{K}_{0}(r)$
$\displaystyle+2\sum_{n=1}^{\infty}\Re\left[B_{n}\mathcal{K}_{n}\left(r\sqrt{1-in\omega}\right)e^{in\theta}\right],$
(84)

for $r>a$.

The coefficients $A_{n}$ and $B_{n}$ are determined by the condition that $u$
and $\nabla u$ are continuous at $r=a$ and that $u$ satisfies the
inhomogeneous equation in Eq. (68). From the continuity condition for $u$, we
obtain

$\displaystyle
A_{n}\mathcal{I}_{n}\left(a\sqrt{1-in\omega}\right)=B_{n}\mathcal{K}_{n}\left(a\sqrt{1-in\omega}\right),$
(85)

and thus we newly set $C_{n}$ so it holds that

$\displaystyle A_{n}=C_{n}\mathcal{K}_{n}\left(a\sqrt{1-in\omega}\right),$
(86)

and

$\displaystyle B_{n}=C_{n}\mathcal{I}_{n}\left(a\sqrt{1-in\omega}\right).$
(87)

Then, we calculate the difference in the derivative in the $r$ direction. Note
that

$\displaystyle\lim_{r\to a-0}\frac{\partial u_{\rm in}}{\partial
r}=\hat{C}_{0}\mathcal{K}_{0}\left(\hat{a}_{0}\right)\mathcal{I}_{1}\left(\hat{a}_{0}\right)$
$\displaystyle\qquad+2\sum_{n=1}^{\infty}\Re\left[\hat{C}_{n}\mathcal{K}_{n}\left(\hat{a}_{n}\right)\frac{\mathcal{I}_{n-1}\left(\hat{a}_{n}\right)+\mathcal{I}_{n+1}\left(\hat{a}_{n}\right)}{2}e^{in\theta}\right],$
(88)

and

$\displaystyle\lim_{r\to a+0}\frac{\partial u_{\rm out}}{\partial
r}=-\hat{C}_{0}\mathcal{I}_{0}\left(\hat{a}_{n}\right)\mathcal{K}_{1}\left(\hat{a}_{n}\right)$
$\displaystyle\qquad-2\sum_{n=1}^{\infty}\Re\left[\hat{C}_{n}\mathcal{I}_{n}\left(\hat{a}_{n}\right)\frac{\mathcal{K}_{n-1}\left(\hat{a}_{n}\right)+\mathcal{K}_{n+1}\left(\hat{a}_{n}\right)}{2}e^{in\theta}\right],$
(89)

where $\hat{a}_{n}=a\sqrt{1-in\omega}$ and
$\hat{C}_{n}=C_{n}\sqrt{1-in\omega}$. Considering that $\hat{C}_{0}=C_{0}$ and
$\hat{a}_{0}=a$, we obtain

$\displaystyle\lim_{r\to a-0}\frac{\partial u_{\rm in}}{\partial r}-\lim_{r\to
a+0}\frac{\partial u_{\rm out}}{\partial r}$
$\displaystyle=C_{0}\left[\mathcal{K}_{0}\left(a\right)\mathcal{I}_{1}\left(a\right)+\mathcal{K}_{1}\left(a\right)\mathcal{I}_{0}\left(a\right)\right]$
$\displaystyle\quad+2\sum_{n=1}^{\infty}\Re\biggl{[}\hat{C}_{n}\left[\mathcal{K}_{n}\left(\hat{a}_{n}\right)\mathcal{I}_{n+1}\left(\hat{a}_{n}\right)+\mathcal{K}_{n+1}\left(\hat{a}_{n}\right)\mathcal{I}_{n}\left(\hat{a}_{n}\right)\right.$
$\displaystyle\qquad\left.+\mathcal{K}_{n}\left(\hat{a}_{n}\right)\mathcal{I}_{n-1}\left(\hat{a}_{n}\right)+\mathcal{K}_{n-1}\left(\hat{a}_{n}\right)\mathcal{I}_{n}\left(\hat{a}_{n}\right)\right]e^{in\theta}\biggr{]}$
$\displaystyle\quad=\frac{C_{0}}{a}+2\sum_{n=1}^{\infty}\Re\left[\frac{\hat{C}_{n}}{\hat{a}_{n}}e^{in\theta}\right]$
$\displaystyle\quad=\frac{C_{0}}{a}+2\sum_{n=1}^{\infty}\Re\left[\frac{C_{n}}{a}e^{in\theta}\right].$
(90)

Therefore, we set $C_{n}=1/(2\pi)$ for all $n$, and we obtain

$\displaystyle\lim_{r\to a-0}\frac{\partial u_{\rm in}}{\partial r}-\lim_{r\to
a+0}\frac{\partial u_{\rm out}}{\partial r}=$ $\displaystyle\frac{1}{2\pi
a}\sum_{n=-\infty}^{\infty}e^{in\theta}$ $\displaystyle=$
$\displaystyle\frac{1}{a}\delta(\theta),$ (91)

which corresponds to the considered situation.

Considering that

$\displaystyle\int_{0}^{2\pi}\mathcal{I}_{n}\left(r\sqrt{1-in\omega}\right)e^{in\theta}d\theta=0,$
(92)
$\displaystyle\int_{0}^{2\pi}\mathcal{K}_{n}\left(r\sqrt{1-in\omega}\right)e^{in\theta}d\theta=0,$
(93)

for any nonzero integer $n$, and that

$\displaystyle\int\left[\omega\frac{\partial
u}{\partial\theta}+\nabla^{2}u-u+\delta\left(\bm{r}-a\bm{e}_{x}\right)\right]d\bm{r}$
$\displaystyle\quad=-\int u(\bm{r})d\bm{r}+1$ $\displaystyle\quad=0,$ (94)

we can explicitly execute the integration of $u(\bm{r})$ as

$\displaystyle\int_{0}^{2\pi}\int_{0}^{\infty}u(r,\theta)d\bm{r}$
$\displaystyle\qquad=2\pi\left[\int_{0}^{a}\frac{1}{2\pi}\mathcal{K}_{0}(a)\mathcal{I}_{0}(r)rdr\right.$
$\displaystyle\qquad\qquad\left.+\int_{a}^{\infty}\frac{1}{2\pi}\mathcal{I}_{0}(a)\mathcal{K}_{0}(r)rdr\right]$
$\displaystyle\qquad=2\pi\left[\frac{1}{2\pi}\mathcal{K}_{0}(a)a\mathcal{I}_{1}(a)+\frac{1}{2\pi}\mathcal{I}_{0}(a)a\mathcal{K}_{1}(a)\right]$
$\displaystyle\qquad=a\left[\mathcal{K}_{0}(a)\mathcal{I}_{1}(a)+\mathcal{I}_{0}(a)\mathcal{K}_{1}(a)\right]$
$\displaystyle\qquad=1,$ (95)

which satisfies Eq. (94).

Therefore, the steady-state concentration field for a single point source at
$\tilde{\bm{r}}=a\tilde{\bm{e}}_{x}$ in the co-rotating frame is written as

$\displaystyle u_{\rm in}(r,\theta)$
$\displaystyle=\frac{1}{2\pi}\mathcal{K}_{0}(a)\mathcal{I}_{0}(r)$
$\displaystyle\qquad+\sum_{n=1}^{\infty}\Re\left[\frac{1}{\pi}\mathcal{K}_{n}\left(a\sqrt{1-in\omega}\right)\mathcal{I}_{n}\left(r\sqrt{1-in\omega}\right)e^{in\theta}\right]$
$\displaystyle=\frac{1}{2\pi}\sum_{n=-\infty}^{\infty}\mathcal{K}_{n}\left(a\sqrt{1-in\omega}\right)\mathcal{I}_{n}\left(r\sqrt{1-in\omega}\right)e^{in\theta},$
(96)

for $r<a$, and

$\displaystyle u_{\rm out}(r,\theta)$
$\displaystyle=\frac{1}{2\pi}\mathcal{I}_{0}(a)\mathcal{K}_{0}(r)$
$\displaystyle\qquad+\sum_{n=1}^{\infty}\Re\left[\frac{1}{\pi}\mathcal{I}_{n}\left(a\sqrt{1-in\omega}\right)\mathcal{K}_{n}\left(r\sqrt{1-in\omega}\right)e^{in\theta}\right]$
$\displaystyle=\frac{1}{2\pi}\sum_{n=-\infty}^{\infty}\mathcal{I}_{n}\left(a\sqrt{1-in\omega}\right)\mathcal{K}_{n}\left(r\sqrt{1-in\omega}\right)e^{in\theta},$
(97)

for $r>a$, which are Eqs. (30) and (31) in the main text. In the calculation,
we used the equalities given in Ref. [43].

##  Appendix B Derivation of the asymptotic form

Here, we derive the asymptotic form of the concentration field far from the
source ($r\gg 1$) for $a\ll 1$.

Considering the Maclaurin expansion of $\mathcal{I}_{n}(z)$ is given as [43]

$\displaystyle\mathcal{I}_{n}(z)=\sum_{k=0}^{\infty}\frac{1}{k!\Gamma(n+k+1)}\left(\frac{z}{2}\right)^{n+2k},$
(98)

the leading term of $I_{n}(z)$ is

$\displaystyle\mathcal{I}_{n}(z)=\frac{1}{2^{n}n!}z^{n}.$ (99)

Therefore, we only need to consider $n=0$ and $n=\pm 1$ as far as we consider
the first order of $a$. Consider that the asymptotic form of
$\mathcal{K}_{n}(z)$ is

$\displaystyle\mathcal{K}_{n}(z)\sim\sqrt{\frac{\pi}{2z}}e^{-z}\sum_{k=0}^{\infty}\frac{\Gamma(n+k+1/2)}{k!\Gamma(n-k+1/2)}\frac{1}{(2z)^{k}},$
(100)

for $\left|\arg z\right|<3\pi/2$ [43]. Here, $\Gamma(\cdot)$ is the gamma
function. For $n=0$ and $n=1$, we obtain

$\displaystyle\mathcal{K}_{0}(z)=\sqrt{\frac{\pi}{2z}}e^{-z}\left(1-\frac{1}{8z}+\mathcal{O}\left(\frac{1}{z^{2}}\right)\right),$
(101)
$\displaystyle\mathcal{K}_{1}(z)=\mathcal{K}_{-1}(z)=\sqrt{\frac{\pi}{2z}}e^{-z}\left(1+\frac{3}{8z}+\mathcal{O}\left(\frac{1}{z^{2}}\right)\right),$
(102)

Then, we only need to consider the terms with $n=0,\pm 1$. By setting
$\chi=\chi_{1}=-\chi_{-1}$, we obtain

$\displaystyle u_{\rm out}(r,\theta)$
$\displaystyle\quad=\frac{1}{2\pi}\left[\mathcal{K}_{0}\left(r\right)+\frac{a}{2}\sqrt{\rho}e^{-i\chi/2}\mathcal{K}_{1}\left(r\sqrt{\rho}e^{-i\chi/2}\right)e^{i\theta}\right.$
$\displaystyle\qquad\left.+\frac{a}{2}\sqrt{\rho}e^{i\chi/2}\mathcal{K}_{1}\left(r\sqrt{\rho}e^{i\chi/2}\right)e^{-i\theta}\right]+\mathcal{O}(a^{2})$
$\displaystyle\quad\simeq\frac{1}{2\pi}\left[\sqrt{\frac{\pi}{2r}}e^{-r}\right.$
(103)
$\displaystyle\qquad\left.+\frac{a}{2}\sqrt{\rho}e^{-i\chi/2}\frac{\sqrt{\pi}}{\sqrt{2r}\rho^{1/4}e^{-i\chi/4}}e^{-r\sqrt{\rho}e^{-i\chi/2}}e^{i\theta}\right.$
$\displaystyle\qquad\left.+\frac{a}{2}\sqrt{\rho}e^{i\chi/2}\frac{\sqrt{\pi}}{\sqrt{2r}\rho^{1/4}e^{i\chi/4}}e^{-r\sqrt{\rho}e^{i\chi/2}}e^{-i\theta}\right]+\mathcal{O}\left(a^{2}\right)$
$\displaystyle\quad=\frac{1}{2\sqrt{2\pi
r}}e^{-r}\left[1+\frac{a}{2}\rho^{1/4}e^{r(1-\sqrt{\rho}\cos(\chi/2))}\right.$
$\displaystyle\qquad\times\left.\left(e^{i(\theta-\chi/4+r\sqrt{\rho}\sin(\chi/2))}+e^{-i(\theta-\chi/4+r\sqrt{\rho}\sin(\chi/2))}\right)\right]$
$\displaystyle\qquad+\mathcal{O}\left(a^{2}\right)$
$\displaystyle\quad=\frac{1}{2\sqrt{2\pi r}}e^{-r}+\frac{a}{2\sqrt{2\pi
r}}\rho^{1/4}e^{-r\sqrt{\rho}\cos(\chi/2)}$
$\displaystyle\qquad\times\cos\left(\theta-\frac{\chi}{4}+r\sqrt{\rho}\sin\frac{\chi}{2}\right)+\mathcal{O}\left(a^{2}\right),$
(104)

and thus we obtain Eq. (32) in the main text. Here, we only considered the
leading terms in Eqs. (101) and (102), and used Eq. (82). In the case that the
phase of the rotor is $\phi_{1}$, we obtain Eq. (33) by replacing $\theta$
with $\theta-\phi_{1}$.

##  Appendix C Calculation on the phase description

Here, we show the detailed calculation of the time evolution using the
averaging method. From Eqs. (23) and (39), we need to obtain $\left.-\nabla
u\cdot\bm{e}(\phi_{2}-\pi/2)\right|_{\bm{r}=\bm{r}_{2}}$. Thus, we first
calculate the gradient of the concentration field in the polar coordinates,

$\displaystyle\nabla u=\frac{\partial u}{\partial
r}\bm{e}(\theta)+\frac{1}{r}\frac{\partial
u}{\partial\theta}\bm{e}\left(\theta+\frac{\pi}{2}\right).$ (105)

The asymptotic form in Eq. (33) is separated into two parts,

$\displaystyle
u(r,\theta,\phi_{1})=u^{(0)}(r)+u^{(1)}(r,\theta,\phi_{1})a+\mathcal{O}(a^{2}),$
(106)

where

$\displaystyle u^{(0)}(r)=\frac{1}{2\sqrt{2\pi r}}e^{-r},$ (107)

and

$\displaystyle u^{(1)}(r,\theta,\phi_{1})=$ $\displaystyle\frac{1}{2\sqrt{2\pi
r}}\rho^{1/4}e^{-r\sqrt{\rho}\cos(\chi/2)}$
$\displaystyle\times\cos\left(\theta-\phi_{1}-\frac{\chi}{4}+r\sqrt{\rho}\sin\frac{\chi}{2}\right).$
(108)

Considering that $u^{(0)}$ does not depend on $\phi_{1}$ or $\phi_{2}$,
$-\left.\nabla u^{(0)}\cdot\bm{e}(\phi_{2}-\pi/2)\right|_{\bm{r}=\bm{r}_{2}}$
is a function of only $\phi_{2}$ and not $\phi_{1}$. As a result of the
averaging, the dependence on $\phi_{2}$ should be omitted, and it gives only a
constant value. Therefore, $u^{(0)}$ only secondarily affects the stability of
the synchronization mode.

Therefore, we consider the effect by $u^{(1)}$. First we calculate the
gradient of $u^{(1)}$ as

$\displaystyle\nabla u^{(1)}=$
$\displaystyle\left[-\frac{1}{2r}u^{(1)}-\sqrt{\rho}u^{(1)}\cos\frac{\chi}{2}-\sqrt{\rho}\hat{u}^{(1)}\sin\frac{\chi}{2}\right]\bm{e}(\theta)$
$\displaystyle-\frac{1}{r}\hat{u}^{(1)}\bm{e}\left(\theta+\frac{\pi}{2}\right),$
(109)

where we set

$\displaystyle\hat{u}^{(1)}=$ $\displaystyle\frac{1}{2\sqrt{2\pi
r}}\rho^{1/4}e^{-\sqrt{\rho}\cos(\chi/2)r}$
$\displaystyle\quad\times\sin\left[\theta-\phi_{1}-\frac{\chi}{4}+\sqrt{\rho}\sin\left(\frac{\chi}{2}\right)r\right].$
(110)

Hereafter, we separately calculate the two cases, i.e., the case with the same
rotation direction and that with the opposite rotation direction. First, we
consider the case with the same rotation. We calculate $\nabla
u_{1}\cdot\bm{e}(\phi_{2}-\pi/2)$ at
$\bm{r}=\bm{r}_{2}=L\bm{e}_{x}-a\bm{e}(\phi_{2})$. Considering that
$r=\sqrt{L^{2}+a^{2}-2La\cos\phi_{2}}$ and
$\tan\theta=-a\sin\phi_{2}/(L-a\cos\phi_{2})$ in the polar coordinates, we
obtain

$\displaystyle\left.-\nabla
u^{(1)}\cdot\bm{e}\left(\phi_{2}-\frac{\pi}{2}\right)\right|_{\bm{r}=\bm{r}_{2}}$
$\displaystyle\quad=\left.-\frac{\partial u^{(1)}}{\partial
r}\bm{e}(\theta)\cdot\bm{e}\left(\phi_{2}-\frac{\pi}{2}\right)\right|_{\bm{r}=\bm{r}_{2}}$
$\displaystyle\qquad\left.-\frac{1}{r}\frac{\partial
u^{(1)}}{\partial\theta}\bm{e}\left(\theta+\frac{\pi}{2}\right)\cdot\bm{e}\left(\phi_{2}-\frac{\pi}{2}\right)\right|_{\bm{r}=\bm{r}_{2}}$
$\displaystyle\quad=\left.\left(\frac{1}{2r}+\sqrt{\rho}\cos\frac{\chi}{2}\right)u^{(1)}\sin(\phi_{2}-\theta)\right|_{\bm{r}=\bm{r}_{2}}$
$\displaystyle\qquad\left.-\left(\frac{1}{r}\cos(\phi_{2}-\theta)-\sqrt{\rho}\sin\frac{\chi}{2}\sin(\phi_{2}-\theta)\right)\hat{u}^{(1)}\right|_{\bm{r}=\bm{r}_{2}}.$
(111)

By considering the Maclaurin expansion of $r$ and $\theta$ with respect to
$a$, we obtain

$\displaystyle r=$ $\displaystyle L+\mathcal{O}(a).$ (112)

and

$\displaystyle\theta=\mathcal{O}(a).$ (113)

Therefore, we obtain

$\displaystyle\sin(\phi_{2}-\theta)=$
$\displaystyle\sin\phi_{2}+\mathcal{O}(a),$ (114)
$\displaystyle\cos(\phi_{2}-\theta)=$
$\displaystyle\cos\phi_{2}+\mathcal{O}(a).$ (115)

We also calculate $u^{(1)}$ and $\hat{u}^{(1)}$ at $\bm{r}=\bm{r}_{2}$ as

$\displaystyle\left.u^{(1)}\right|_{\bm{r}=\bm{r}_{2}}$
$\displaystyle=\frac{\rho^{1/4}e^{-\sqrt{\rho}(L-a\cos\phi)\cos(\chi/2)}}{2\sqrt{2\pi(L-a\cos\phi_{2})}}\cos\left[-\frac{a}{L}\sin\phi_{2}-\phi_{1}-\frac{\chi}{4}\right.$
$\displaystyle\qquad\left.+\sqrt{\rho}(L-a\cos\phi_{2})\sin\left(\frac{\chi}{2}\right)\right]$
$\displaystyle=\frac{\rho^{1/4}e^{-\sqrt{\rho}L\cos(\chi/2)}}{2\sqrt{2\pi
L}}\cos\left[-\phi_{1}-\frac{\chi}{4}+\sqrt{\rho}\sin\left(\frac{\chi}{2}\right)L\right]$
$\displaystyle\quad+\mathcal{O}(a),$ (116)

and

$\displaystyle\left.\hat{u}^{(1)}\right|_{\bm{r}=\bm{r}_{2}}$
$\displaystyle=\frac{\rho^{1/4}e^{-\sqrt{\rho}(L-a\cos\phi_{2})\cos(\chi/2)}}{2\sqrt{2\pi(L-a\cos\phi)}}\sin\left[-\frac{a}{L}\sin\phi_{2}-\phi_{1}-\frac{\chi}{4}\right.$
$\displaystyle\qquad\left.+\sqrt{\rho}(L-a\cos\phi_{2})\sin\left(\frac{\chi}{2}\right)\right]$
$\displaystyle=\frac{\rho^{1/4}e^{-\sqrt{\rho}L\cos(\chi/2)}}{2\sqrt{2\pi
L}}\sin\left[-\phi_{1}-\frac{\chi}{4}+\sqrt{\rho}\sin\left(\frac{\chi}{2}\right)L\right]$
$\displaystyle\quad+\mathcal{O}(a).$ (117)

Equation (111) with Eqs. (114)–(117) leads to

$\displaystyle\left.-\nabla
u^{(1)}\cdot\bm{e}\left(\phi_{2}-\frac{\pi}{2}\right)\right|_{\bm{r}=\bm{r}_{2}}$
$\displaystyle\quad=\frac{\rho^{1/4}e^{-\sqrt{\rho}L\cos(\chi/2)}}{2\sqrt{2\pi
L}}$
$\displaystyle\qquad\times\left\\{-\frac{1}{2L}\sin\left[-\phi_{2}-\phi_{1}-\frac{\chi}{4}+L\sqrt{\rho}\sin\left(\frac{\chi}{2}\right)\right]\right.$
$\displaystyle\qquad\left.-\frac{1}{2L}\cos\phi_{2}\sin\left[-\phi_{1}-\frac{\chi}{4}+L\sqrt{\rho}\sin\left(\frac{\chi}{2}\right)\right]\right.$
$\displaystyle\qquad\left.+\sqrt{\rho}\sin\phi_{2}\cos\left[-\phi_{1}-\frac{3\chi}{4}+L\sqrt{\rho}\sin\left(\frac{\chi}{2}\right)\right]\right\\}$
$\displaystyle\qquad+\mathcal{O}(a),$ (118)

and thus Eqs. (39) and (106) give

$\displaystyle\frac{d\phi_{2}}{dt}\simeq\omega+\frac{\Gamma\rho^{1/4}e^{-\sqrt{\rho}L\cos(\chi/2)}}{4\eta\sqrt{2\pi
L}}$
$\displaystyle\times\left\\{-\frac{1}{2L}\sin\left[\Delta\phi-\frac{\chi}{4}+L\sqrt{\rho}\sin\left(\frac{\chi}{2}\right)\right]\right.$
$\displaystyle\left.\quad+\sqrt{\rho}\sin\left[\Delta\phi-\frac{3\chi}{4}+L\sqrt{\rho}\sin\left(\frac{\chi}{2}\right)\right]\right\\}.$
(119)

Considering geometric symmetry, we obtain

$\displaystyle\frac{d\phi_{1}}{dt}\simeq\omega+\frac{\Gamma\rho^{1/4}e^{-\sqrt{\rho}L\cos(\chi/2)}}{4\eta\sqrt{2\pi
L}}$
$\displaystyle\times\left\\{-\frac{1}{2L}\sin\left[-\Delta\phi-\frac{\chi}{4}+L\sqrt{\rho}\sin\left(\frac{\chi}{2}\right)\right]\right.$
$\displaystyle\left.\quad+\sqrt{\rho}\sin\left[-\Delta\phi-\frac{3\chi}{4}+L\sqrt{\rho}\sin\left(\frac{\chi}{2}\right)\right]\right\\}.$
(120)

From these equations, we obtain the time evolution of
$\Delta\phi=\phi_{2}-\phi_{1}$ as

$\displaystyle\frac{d\Delta\phi}{dt}=$
$\displaystyle\frac{\Gamma\rho^{1/4}e^{-\sqrt{\rho}L\cos(\chi/2)}}{2\eta\sqrt{2\pi
L}}\sin\Delta\phi$
$\displaystyle\times\left\\{-\frac{1}{2L}\cos\left[-\frac{\chi}{4}+L\sqrt{\rho}\sin\left(\frac{\chi}{2}\right)\right]\right.$
$\displaystyle\left.\quad+\sqrt{\rho}\cos\left[-\frac{3\chi}{4}+L\sqrt{\rho}\sin\left(\frac{\chi}{2}\right)\right]\right\\},$
(121)

which corresponds to Eq. (44). In the calculation, we use

$\displaystyle\cos(\Delta\phi+\Xi)-\cos(-\Delta\phi+\Xi)=2\cos\Xi\sin\Delta\phi$
(122)

Next, we consider the case with the opposite rotation direction. In this case,
we have to obtain $\left.-\nabla
u\cdot\bm{e}(-\phi_{2}+\pi/2)\right|_{\bm{r}=\bm{r}_{2}}$. The polar
coordinates corresponding to $\bm{r}_{2}$ change as
$r=\sqrt{L^{2}+a^{2}-2La\cos\phi_{2}}$ and
$\tan\theta=a\sin\phi_{2}/(L-a\cos\phi_{2})$. Then, we obtain

$\displaystyle\left.-\nabla
u^{(1)}\cdot\bm{e}\left(-\phi_{2}+\frac{\pi}{2}\right)\right|_{\bm{r}=\bm{r}_{2}}$
$\displaystyle\quad=\left.-\frac{\partial u^{(1)}}{\partial
r}\bm{e}(\theta)\cdot\bm{e}\left(-\phi_{2}+\frac{\pi}{2}\right)\right|_{\bm{r}=\bm{r}_{2}}$
$\displaystyle\qquad\left.-\frac{1}{r}\frac{\partial
u^{(1)}}{\partial\theta}\bm{e}\left(\theta+\frac{\pi}{2}\right)\cdot\bm{e}\left(-\phi_{2}+\frac{\pi}{2}\right)\right|_{\bm{r}=\bm{r}_{2}}$
$\displaystyle\quad=\left.\left(\frac{1}{2r}+\sqrt{\rho}\cos\frac{\chi}{2}\right)u^{(1)}\sin(\phi_{2}+\theta)\right|_{\bm{r}=\bm{r}_{2}}$
$\displaystyle\qquad\left.+\left(\frac{1}{r}\cos(\phi_{2}+\theta)+\sqrt{\rho}\sin\frac{\chi}{2}\sin(\phi_{2}+\theta)\right)\hat{u}^{(1)}\right|_{\bm{r}=\bm{r}_{2}}.$
(123)

Equations (112)–(115) do not change irrespective of the rotation direction
within the order of $\mathcal{O}(1)$, and thus we can adopt the expression in
Eqs. (116) and (117). Therefore, we obtain

$\displaystyle\left.-\nabla
u^{(1)}\cdot\bm{e}\left(-\phi_{2}+\frac{\pi}{2}\right)\right|_{\bm{r}=\bm{r}_{2}}$
$\displaystyle\quad=\frac{\rho^{1/4}e^{-\sqrt{\rho}L\cos(\chi/2)}}{2\sqrt{2\pi
L}}$
$\displaystyle\qquad\times\left\\{-\frac{1}{2L}\sin\left[-\phi_{2}-\phi_{1}-\frac{\chi}{4}+L\sqrt{\rho}\sin\left(\frac{\chi}{2}\right)\right]\right.$
$\displaystyle\qquad\left.+\frac{3}{2L}\cos\phi_{2}\sin\left[-\phi_{1}-\frac{\chi}{4}+L\sqrt{\rho}\sin\left(\frac{\chi}{2}\right)\right]\right.$
$\displaystyle\qquad\left.+\sqrt{\rho}\sin\phi_{2}\cos\left[-\phi_{1}-\frac{3\chi}{4}+L\sqrt{\rho}\sin\left(\frac{\chi}{2}\right)\right]\right\\}$
$\displaystyle\qquad+\mathcal{O}(a),$ (124)

Equations (36) and (106) give

$\displaystyle\frac{d\phi_{2}}{dt}\simeq\omega+\frac{\Gamma\rho^{1/4}e^{-\sqrt{\rho}L\cos(\chi/2)}}{4\eta\sqrt{2\pi
L}}$
$\displaystyle\times\left\\{\frac{3}{2L}\sin\left[\Delta\phi-\frac{\chi}{4}+L\sqrt{\rho}\sin\left(\frac{\chi}{2}\right)\right]\right.$
$\displaystyle\left.\quad+\sqrt{\rho}\sin\left[\Delta\phi-\frac{3\chi}{4}+L\sqrt{\rho}\sin\left(\frac{\chi}{2}\right)\right]\right\\}.$
(125)

In the same manner as that with the same rotation direction, the time-
evolution equation is obtained by considering the geometric symmetry as

$\displaystyle\frac{d\phi_{1}}{dt}\simeq\omega+\frac{\Gamma\rho^{1/4}e^{-\sqrt{\rho}L\cos(\chi/2)}}{4\eta\sqrt{2\pi
L}}$
$\displaystyle\times\left\\{\frac{3}{2L}\sin\left[-\Delta\phi-\frac{\chi}{4}+L\sqrt{\rho}\sin\left(\frac{\chi}{2}\right)\right]\right.$
$\displaystyle\left.\quad+\sqrt{\rho}\sin\left[-\Delta\phi-\frac{3\chi}{4}+L\sqrt{\rho}\sin\left(\frac{\chi}{2}\right)\right]\right\\}.$
(126)

From these equations, the time-evolution equation for $\Delta\phi$ is obtained
as in Eq. (47).

##  Appendix D Phase coupling function including the time change in angular
velocity

Figure 12:  Plots of (a) $\mathcal{G}_{s}$, (b) $\mathcal{G}_{o}$, (c)
$\mathcal{H}_{s}$, and (d) $\mathcal{H}_{o}$ against $\Delta\varphi$. The
results with $L=2$ (red), $L=2.5$ (orange), $L=3$ (green), and $L=3.5$ (cyan)
are shown in each panel. Figure 13:  Plots of $\hat{\mathcal{H}}_{s}^{(1)}$
and $\hat{\mathcal{H}}_{o}^{(1)}$ against $L$. The positive and negative
values are indicated with red and cyan points, respectively. The positive
$\hat{\mathcal{H}}_{s}^{(1)}$ or $\hat{\mathcal{H}}_{o}^{(1)}$ means that the
anti-phase synchronization state is stable, while the negative
$\hat{\mathcal{H}}_{s}^{(1)}$ or $\hat{\mathcal{H}}_{o}^{(1)}$ means that the
in-phase synchronization state is stable.

In this section, we discuss the phase coupling function considering the time
change in the angular velocity. We calculated the dynamics of one rotor (rotor
1) taking into consideration the concentration field of chemicals generated by
the other rotor (rotor 2), though the position of the particle composed of
rotor 2 is approximated to be located at the center of it, i.e.,
$\bm{r}_{2}(t)=\bm{\ell}_{2}$, in order to neglect the dependence of
$\phi_{2}$. In this case, the angular velocity of the first rotor
$d\phi_{1}/dt$ is no longer constant but depends on the phase of the rotor 1,
$\phi_{1}$. Thus we have to define a new phase $\varphi_{i}$, which is defined
by $d\varphi_{i}/dt=2\pi/T$, where $T$ is a period. $\varphi_{i}$ is described
as a monotonous increasing function of $\phi_{i}$. The functions
$\mathcal{G}_{s}(\Delta\varphi)$ and $\mathcal{G}_{o}(\Delta\varphi)$ are
defined as

$\displaystyle\mathcal{G}_{s}(\Delta\varphi)$ $\displaystyle=\frac{1}{2\pi\eta
a}\int_{0}^{2\pi}\bm{e}\left(\varphi_{2}-\frac{\pi}{2}\right)\cdot\bm{F}_{u,1,2}(\phi_{2}+\Delta\varphi,\varphi_{2})d\varphi_{2},$
(127) $\displaystyle\mathcal{G}_{o}(\Delta\varphi)$
$\displaystyle=\frac{1}{2\pi\eta
a}\int_{0}^{2\pi}\bm{e}\left(-\left(\varphi_{2}-\frac{\pi}{2}\right)\right)\cdot\bm{F}_{u,1,2}(\phi_{2}+\Delta\varphi,\varphi_{2})d\varphi_{2}.$
(128)

Here, we calculated the force $\bm{F}_{u,1,2}$ by taking the rotation of rotor
2 into consideration. Then, the time-evolution equation of $\Delta\varphi$ is
expressed as

$\displaystyle\frac{d\Delta\varphi}{dt}=\mathcal{G}_{i}(-\Delta\varphi)-\mathcal{G}_{i}(\Delta\varphi)\equiv\mathcal{H}_{i}(\Delta\varphi)$
(129)

where $i$ denotes $s$ or $o$. The parameters for the simulation were the same
as in Sec. V. The plots of $\mathcal{G}_{s}(\Delta\varphi)$,
$\mathcal{G}_{o}(\Delta\varphi)$, $\mathcal{H}_{s}(\Delta\varphi)$, and
$\mathcal{H}_{o}(\Delta\varphi)$ obtained by the numerical simulation are
shown in Fig. 12. From Figs. 12(a) and (b), $\mathcal{G}_{s}(\Delta\varphi)$
and $\mathcal{G}_{o}(\Delta\varphi)$ are less than those in Figs. 10(a) and
(b). This is due to the interaction through the concentration field, and shows
that the effect by the concentration field from the other rotor reduces the
averaged angular velocity. Despite the difference in
$\mathcal{G}_{s}(\Delta\varphi)$ and $\mathcal{G}_{o}(\Delta\varphi)$,
$\mathcal{H}_{s}(\Delta\varphi)$ and $\mathcal{H}_{o}(\Delta\varphi)$ in Figs.
12(c) and (d) are almost the same as those in Figs. 10(c) and (d). In the same
manner as in the previous paragraph, we consider the Fourier expansion of
$\mathcal{H}_{s}(\Delta\varphi)$ and $\mathcal{H}_{o}(\Delta\varphi)$ as

$\displaystyle\mathcal{H}_{i}(\Delta\varphi)=\sum_{k=1}^{\infty}\hat{\mathcal{H}}_{i}^{(k)}\sin
k\Delta\varphi,$ (130)

where $i$ denotes $s$ or $o$. In Fig. 13,
$\hat{\mathcal{H}}^{(1)}_{s}(\Delta\varphi)$ and
$\hat{\mathcal{H}}^{(1)}_{o}(\Delta\varphi)$ are plotted against $L$, which
are almost the same as those in Figs. 8 and 11. This indicates that the
decrease in the averaged angular velocity does not matter for the stable
synchronization mode, but the interaction between the rotor position and time-
dependent concentration field, which is shown in Figs. 9 (c) and (d), plays an
important role.

## References

  * Michelin  _et al._ [2013] S. Michelin, E. Lauga, and D. Bartolo, Spontaneous autophoretic motion of isotropic particles, Phys. Fluids 25, 061701 (2013).
  * Ohta [2017] T. Ohta, Dynamics of deformable active particles, J. Phys. Soc. Jpn. 86, 072001 (2017).
  * Klindt  _et al._ [2017] G. S. Klindt, C. Ruloff, C. Wagner, and B. M. Friedrich, In-phase and anti-phase flagellar synchronization by waveform compliance and basal coupling, New J. Phys. 19, 113052 (2017).
  * Uchida  _et al._ [2017] N. Uchida, R. Golestanian, and R. R. Bennett, Synchronization and collective dynamics of flagella and cilia as hydrodynamically coupled oscillators, J. Phys. Soc. Jpn. 86, 101007 (2017).
  * Sens [2020] P. Sens, Stick-slip model for actin-driven cell protrusions, cell polarization, and crawling, Proc. Natl. Acad. Sci. 117, 24670 (2020).
  * Tarama and Yamamoto [2018] M. Tarama and R. Yamamoto, Mechanics of cell crawling by means of force-free cyclic motion, J. Phys. Soc. Jpn. 87, 044803 (2018).
  * Budrene and Berg [1995] E. O. Budrene and H. C. Berg, Dynamics of formation of symmetrical patterns by chemotactic bacteria, Nature 376, 49 (1995).
  * Shoji  _et al._ [2014] E. Shoji, H. Nishimori, A. Awazu, S. Izumi, and M. Iima, Localized bioconvection patterns and their initial state dependency in euglena gracilis suspensions in an annular container, J. Phys. Soc. Jpn. 83, 043001 (2014).
  * Nakajima  _et al._ [2014] A. Nakajima, S. Ishihara, D. Imoto, and S. Sawai, Rectified directional sensing in long-range cell migration, Nat. Commun. 5 (2014).
  * Mikhailov and Calenbuhr [2002] A. Mikhailov and V. Calenbuhr,  _From Cells to Societies: Models of Complex Coherent Action_ (Springer, Berlin, 2002).
  * Tomlinson and Miller [1862] C. Tomlinson and W. A. Miller, On the motions of camphor on the surface of water, Proc. R. Soc. Lond. 11, 575 (1862).
  * Strutt [1890] R. J. Strutt, Measurements of the amount of oil necessary in order to check the motions of camphor upon water, Proc. R. Soc. Lond. 47, 364 (1890).
  * Nakata  _et al._ [1997] S. Nakata, Y. Iguchi, S. Ose, M. Kuboyama, T. Ishii, and K. Yoshikawa, Self-rotation of a camphor scraping on water: New insight into the old problem, Langmuir 13, 4454 (1997).
  * Nakata  _et al._ [2015] S. Nakata, M. Nagayama, H. Kitahata, N. J. Suematsu, and T. Hasegawa, Physicochemical design and analysis of self-propelled objects that are characteristically sensitive to environments, Phys. Chem. Chem. Phys. 17, 10326 (2015).
  * Nakata  _et al._ [2019] S. Nakata, V. Pimienta, I. Lagzi, H. Kitahata, and N. J. Suematsu, eds.,  _Self-organized Motion_, Theoretical and Computational Chemistry Series (The Royal Society of Chemistry, London, 2019).
  * Boniface  _et al._ [2019] D. Boniface, C. Cottin-Bizonne, R. Kervil, C. Ybert, and F. Detcheverry, Self-propulsion of symmetric chemically active particles: Point-source model and experiments on camphor disks, Phys. Rev. E 99, 062605 (2019).
  * Ei  _et al._ [2018] S.-I. Ei, H. Kitahata, Y. Koyano, and M. Nagayama, Interaction of non-radially symmetric camphor particles, Physica D 366, 10 (2018).
  * Shimokawa and Sakaguchi [2022] M. Shimokawa and H. Sakaguchi, Hula-hoop-like motion and subcritical bifurcation of the rotation of an elliptical paper disk coated with camphor, J. Phys. Soc. Jpn. 91, 074002 (2022).
  * Koyano  _et al._ [2017] Y. Koyano, M. Gryciuk, P. Skrobanska, M. Malecki, Y. Sumino, H. Kitahata, and J. Gorecki, Relationship between the size of a camphor-driven rotor and its angular velocity, Phys. Rev. E 96, 012609 (2017).
  * Koyano  _et al._ [2019a] Y. Koyano, H. Kitahata, M. Gryciuk, N. Akulich, A. Gorecka, M. Malecki, and J. Gorecki, Bifurcation in the angular velocity of a circular disk propelled by symmetrically distributed camphor pills, Chaos 29, 013125 (2019a).
  * Sharma  _et al._ [2019] J. Sharma, I. Tiwari, D. Das, P. Parmananda, V. S. Akella, and V. Pimienta, Rotational synchronization of camphor ribbons, Phys. Rev. E 99, 012204 (2019).
  * Sharma  _et al._ [2020] J. Sharma, I. Tiwari, D. Das, P. Parmananda, and V. Pimienta, Rotational synchronization of camphor ribbons in different geometries, Phys. Rev. E 101, 052202 (2020).
  * Sharma  _et al._ [2021] J. Sharma, I. Tiwari, D. Das, and P. Parmananda, Chimeralike states in a minimal network of active camphor ribbons, Phys. Rev. E 103, 012214 (2021).
  * Sharma  _et al._ [2022] J. Sharma, I. Tiwari, P. Parmananda, and M. Rivera, Aperiodic bursting dynamics of active rotors, Phys. Rev. E 105, 014216 (2022).
  * Jain  _et al._ [2022] R. Jain, J. Sharma, I. Tiwari, S. D. Gadre, S. Kumarasamy, P. Parmananda, and A. Prasad, Generation of aperiodic motion due to sporadic collisions of camphor ribbons, Phys. Rev. E 106, 024201 (2022).
  * Morohashi  _et al._ [2019] H. Morohashi, M. Imai, and T. Toyota, Construction of a chemical motor-movable frame assembly based on camphor grains using water-floating 3d-printed models, Chem. Phys. Lett. 721, 104 (2019).
  * Frenkel  _et al._ [2017] M. Frenkel, G. Whyman, E. Shulzinger, A. Starostin, and E. Bormashenko, Self-propelling rotator driven by soluto-capillary Marangoni flows, Appl. Phys. Lett. 110, 131604 (2017).
  * Frenkel  _et al._ [2019] M. Frenkel, A. Vilk, I. Legchenkova, S. Shoval, and E. Bormashenko, Mini-generator of electrical power exploiting the Marangoni flow inspired self-propulsion, ACS Omega 4, 15265 (2019).
  * Koyano and Kitahata [2021] Y. Koyano and H. Kitahata, Imperfect bifurcation in the rotation of a propeller-shaped camphor rotor, Phys. Rev. E 103, 012202 (2021).
  * Hayashima  _et al._ [2001] Y. Hayashima, M. Nagayama, and S. Nakata, A camphor grain oscillates while breaking symmetry, J. Phys. Chem. B 105, 5353 (2001).
  * Nagayama  _et al._ [2004] M. Nagayama, S. Nakata, Y. Doi, and Y. Hayashima, A theoretical and experimental study on the unidirectional motion of a camphor disk, Physica D 194, 151 (2004).
  * Suematsu  _et al._ [2014] N. J. Suematsu, T. Sasaki, S. Nakata, and H. Kitahata, Quantitative estimation of the parameters for self-motion driven by difference in surface tension, Langmuir 30, 8101 (2014).
  * Kitahata and Yoshinaga [2018] H. Kitahata and N. Yoshinaga, Effective diffusion coefficient including the Marangoni effect, J. Chem. Phys. 148, 134906 (2018).
  * Gustafson and Abe [1998] K. Gustafson and T. Abe, The third boundary condition—was it Robin’s?, Math. Intell. 20, 63 (1998).
  * Pikovsky  _et al._ [2001] A. Pikovsky, M. G. Rosenblum, and J. Kurths,  _Synchronization, A Universal Concept in Nonlinear Sciences_ (Cambridge University Press, Cambridge, 2001).
  * Stankovski  _et al._ [2017] T. Stankovski, T. Pereira, P. V. E. McClintock, and A. Stefanovska, Coupling functions: Universal insights into dynamical interaction mechanisms, Rev. Mod. Phys. 89, 045001 (2017).
  * Koyano  _et al._ [2019b] Y. Koyano, N. J. Suematsu, and H. Kitahata, Rotational motion of a camphor disk in a circular region, Phys. Rev. E 99, 022211 (2019b).
  * Mori and Kuramoto [1998] H. Mori and Y. Kuramoto,  _Dissipative structures and chaos_ , 2nd ed. (Springer, Berlin, 1998).
  * Kitahata  _et al._ [2022] H. Kitahata, Y. Koyano, R. J. Löffler, and J. Górecki, Complexity and bifurcations in the motion of a self-propelled rectangle confined in a circular water chamber, Phys. Chem. Chem. Phys. 24, 20326 (2022).
  * Takamatsu  _et al._ [2000] A. Takamatsu, T. Fujii, and I. Endo, Time delay effect in a living coupled oscillator system with the plasmodium of physarum polycephalum, Phys. Rev. Lett. 85, 2026 (2000).
  * Kitahata  _et al._ [2009] H. Kitahata, J. Taguchi, M. Nagayama, T. Sakurai, Y. S. Ikura, A. Osa, Y. Sumino, M. Tanaka, E. Yokoyama, and H. Miike, Oscillation and synchronization in the combustion of candles., J. Phys. Chem. A 113 29, 8164 (2009).
  * Manoj  _et al._ [2018] K. Manoj, S. A. Pawar, and R. I. Sujith, Experimental evidence of amplitude death and phase-flip bifurcation between in-phase and anti-phase synchronization, Sci. Rep. 8, 11626 (2018).
  * Watson [1944] G. N. Watson,  _A Treatise on the Theory of Bessel Functions_ (Cambridge University Press, Cambridge, 1944).