Mathematical modeling for the synchronization of two interacting active rotors

TL;DR

This paper develops a mathematical model to understand how two active rotors synchronize their rotation through chemical concentration fields, revealing conditions for in-phase and anti-phase synchronization.

Contribution

It introduces a combined numerical and theoretical framework to analyze the synchronization of active rotors influenced by surface-active chemical fields.

Findings

In-phase and anti-phase synchronization depend on rotor distance.

Theoretical predictions match numerical simulation results.

Stability of synchronization modes is analyzed via phase reduction.

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.

Equations

$${\mathbf{r}}_1={\mathbf{\ell }}_1+a\mathbf{e}({\varphi }_1)=a\mathbf{e}({\varphi }_1).$$

$${\mathbf{r}}_1={\mathbf{\ell }}_1+a\mathbf{e}({\varphi }_1)=-\frac{L}{2}{\mathbf{e}}_x+a\mathbf{e}({\varphi }_1),$$

$${\mathbf{r}}_2={\mathbf{\ell }}_2-a\mathbf{e}(\pm {\varphi }_2)=\frac{L}{2}{\mathbf{e}}_x-a\mathbf{e}(\pm {\varphi }_2).$$

$${\eta }_{i}A\frac{d{\mathbf{r}}_i}{dt}={\eta }_{i}A{\mathbf{v}}_i={\mathbf{F}}_{u,i}+{\mathbf{F}}_{c,i},$$

$${\mathbf{v}}_1=a\frac{d{\varphi }_1}{dt}\mathbf{e}\left({\varphi }_1+\frac{\pi }{2}\right),$$

$${\mathbf{v}}_2=a\frac{d{\varphi }_2}{dt}\mathbf{e}\left({\varphi }_2\pm \frac{\pi }{2}\right),$$

$${\mathbf{F}}_{u,i}=$$

$$=$$

$$\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.$$

$${\eta }_{1}A{a}^2\frac{d{\varphi }_1}{dt}=$$

$$=$$

$${\eta }_{2}A{a}^2\frac{d{\varphi }_2}{dt}=$$

$$=$$

$$\frac{\partial u}{\partial t}={\nabla }^{2}u-u+\frac{1}{A}\sum _{i=1}^{N}S(\mathbf{r}-{\mathbf{r}}_i),$$

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

$$\frac{\partial u}{\partial t}={\nabla }^{2}u-u+\sum _{i=1}^N\delta (\mathbf{r}-{\mathbf{r}}_i),$$

$$u=u_1+u_2,$$

$$\frac{\partial u_i}{\partial t}={\nabla }^2{u}_i-u_i+\delta (\mathbf{r}-{\mathbf{r}}_i),$$

$$\eta a^2\frac{d{\varphi }_1}{dt}=$$

$$\eta a^2\frac{d{\varphi }_2}{dt}=$$

$${\mathbf{F}}_{u,i}={\mathbf{F}}_{u,i,1}+{\mathbf{F}}_{u,i,2},$$

$$\frac{1}{A}{\mathbf{F}}_{u,i,j}=$$

$$\to$$

$$\frac{1}{\eta Aa}\mathbf{e}\left({\varphi }_1+\frac{\pi }{2}\right)\cdot {\mathbf{F}}_{u,1,1}=\omega ,$$

$$\frac{1}{\eta Aa}\mathbf{e}\left(\pm \left({\varphi }_2-\frac{\pi }{2}\right)\right)\cdot {\mathbf{F}}_{u,2,2}=\omega ,$$

$$-\omega \frac{\partial \tilde{u}}{\partial \tilde{\theta }}={\tilde{\nabla }}^2\tilde{u}-\tilde{u}+\delta \left(\tilde{\mathbf{r}}-a{\tilde{\mathbf{e}}}_x\right),$$

$$\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.$$

$${\tilde{u}}_{\mathrm{in}}(\tilde{r},\tilde{\theta })$$

$$=\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 }},$$

$${\tilde{u}}_{\mathrm{out}}(\tilde{r},\tilde{\theta })$$

$$=\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 }}.$$

Full text

††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.

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

[TABLE]

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

[TABLE]

[TABLE]

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

[TABLE]

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

[TABLE]

[TABLE]

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

[TABLE]

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.,

[TABLE]

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

[TABLE]

and

[TABLE]

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

[TABLE]

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 [math] 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$.

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].

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

[TABLE]

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.

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

[TABLE]

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,

[TABLE]

where

[TABLE]

for $i=1,2$.

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

[TABLE]

and

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

for rotor 1, and

[TABLE]

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

[TABLE]

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

[TABLE]

where

[TABLE]

and

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

where

[TABLE]

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.

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

[TABLE]

which leads to

[TABLE]

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

[TABLE]

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

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}$,

[TABLE]

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

[TABLE]

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$.

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

In other words, we have the relation as

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

The homogeneous equation for (68) is expressed as

[TABLE]

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

[TABLE]

we obtain

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

then we obtain

[TABLE]

[TABLE]

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:

[TABLE]

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

for $r<a$, and

[TABLE]

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

[TABLE]

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

[TABLE]

and

[TABLE]

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

[TABLE]

and

[TABLE]

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

[TABLE]

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

[TABLE]

which corresponds to the considered situation.

Considering that

[TABLE]

[TABLE]

for any nonzero integer $n$, and that

[TABLE]

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

[TABLE]

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

[TABLE]

for $r<a$, and

[TABLE]

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]

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

[TABLE]

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

[TABLE]

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,

[TABLE]

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

[TABLE]

where

[TABLE]

and

[TABLE]

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

[TABLE]

where we set

[TABLE]

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

[TABLE]

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

[TABLE]

and

[TABLE]

Therefore, we obtain

[TABLE]

[TABLE]

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

[TABLE]

and

[TABLE]

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

[TABLE]

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

[TABLE]

Considering geometric symmetry, we obtain

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

Equations (36) and (106) give

[TABLE]

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

[TABLE]

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

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

[TABLE]

[TABLE]

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

[TABLE]

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

[TABLE]

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.