Self-propelled particles, which can transduce free energy into kinetic energy under nonequilibrium conditions, have been investigated for last decades since they can be a model for the motion of living organisms [1–6]. The relationship between the self-propelled motion and the shape of the particle has attracted much attention, and both experimental and theoretical studies have been intensively reported [7–11]. When we discuss the mechanism of the self-propelled motion, we have to keep in mind that the conservation of momentum should hold rigourously. For example, a self-propelled liquid droplet in the other immisible liquid can be discussed based on the momentum exchange in hydrodynamics. There have been many papers on such self-propelled motions, most of which are considered in the regime of low-Reynolds number [4, 12–15]. In contrast, some systems, such as a self-propelled particle at a liquid surface, a self-propelled droplet on a solid surface, and a living cell on a substrate, can apparently break the momentum conservation [16–18]. It is sure that they exchange momentum with the environment but the momentum conservation is not treated explicitly in their model.

The self-propelled motion of a camphor disk at a water surface is one of the most studied self-propelled motion [19–26]. The camphor disk releases the camphor molecules at a water surface, and the concentration gradient of the camphor molecules induces the surface tension gradient, which drives the disk. Such a camphor disk motion can be discussed based on the hydrodynamics [27–29], but it has also been discussed from the viewpoint of the coupling between the reaction-diffusion field and the disk motion [9, 22, 30, 31]. Based on the latter approach, we have studied the self-propelled motion of a camphor particle with the shape other than a circle, and have discussed the relation between the camphor particle shape and the direction of motion [6, 9, 32–35]. For example, we show that an elliptic camphor particle moves in the dircetion of its minor axis [9, 33–35]. As for a triangular-shaped camphor particle, the direction of motion depends on the size of the particle; a smaller particle moves in the direction of a corner, while a larger particle moves in the direction of a side [35]. These studies clarify that the symmetric property of the shape is important for determining the preferable direction of the self-propelled motion.

Motivated by these previous studies, here we propose the simpler system with which the relation between the direction of self-propelled motion and symmetric property can be discussed; a camphor float with n-fold rotational symmetry, which comprises a rigid light circular plate attaching n camphor disks along a periphery with an equivalent spacing. From the viewpoint of rotational symmetry, an elliptic camphor particle corresponds to the camphor float with n = 2, while the triangular camphor particle corresponds to that with n = 3; we mainly discuss these two cases. By considering these camphor floats, we discuss the direction of the self-propelled motion affected by the symmetric properties. In the present paper, we first formulate the mathematical model for them, and then analytically discuss the preferable direction of the self-propelled motion. Finally, we check the validity of the analytical results by numerical simulations.

We construct a mathematical model that describes the self-propelled motion of the camphor float with n-fold rotational symmetry, which is floating at a water surface. The float comprises a rigid light circular plate attaching n camphor disks along a periphery with an equivalent spacing. We consider the concentration field of camphor molecules u( r , t) on a two-dimensional plane corresponding to the water surface. As for the motion of the camphor float, the position and configureation can be described by the center of mass of the camphor float r _c and the characteristic angle θ _c . For the camphor float with n-fold rotational symmetry, the center position of the k-th camphor disk (k = 0, 1, …, n − 1) is set as r k = r c + ℓ e θ c + 2 π k n , (1) where e (θ) is the unit vector defined as e θ = cos ⁡ θ e x + sin ⁡ θ e y . (2)

Here, e _x and e _y are the unit vectors in x- and y-directions, respectively. The schematic illustrations for the camphor floats with n = 2 and 3 are shown in Figure 1.

Each camphor disk has a radius of R and supplys the camphor molecules to the water surface at the rate of S ₀. The dynamics of u( r , t) is described as ∂ u ∂ t = D ∇ 2 u − a u + S 0 π R 2 ∑ k = 0 n − 1 Θ R − r − r k . (3)

The first, second, and third terms in the righthand side correspond to the diffusion, sublimation, and supply of the camphor molecules. Here, D is the effiective diffusion coefficient [36, 37], a is the sublimation rate, and Θ ⋅ is the Heaviside’s step function. It should be noted that the effective diffusion coefficient is introduced to include the effect by the hydrodynamic transport due to the Marangoni effect. The effective diffusion coefficient is estimated from the concentration profile with a resting camphor particle based on the reaction-diffusion-advection equation. Thus, this approxiation works well when the velocity is small.

As for the motion of the camphor float, we consider the equation of motion for the k-th camphor disk as m d 2 r k d t 2 = − η d r k d t + F k c o n c + F k i n t , (4) where m is the mass, and η is the viscous friction coefficient of the camphor disk. We assume that they are both proportional to the area of the camphor disk as m = πρR ² and η = πκR ². F k ( c o n c ) is the force originating from the surface tension gradient due to the camphor concentration field described as F k c o n c = ∮ ∂ Ω k γ u r ′ n r ′ d ℓ ′ (5) = ∬ Ω k ∇ ′ γ u r ′ d A ′ . (6)

Here, Ω _k is the region of the k-th camphor disk, which is defined as Ω k = r | r − r k ≤ R , (7) and ∂Ω _k is the periphery of Ω _k . dℓ′ is the line element along the periphery, n ( r ′) is the outward unit normal vector at the periphery, ∇′ is the nabla operator with respect to r ′, and dA′ is the area element in the region Ω _k . Here, the relation between the surface tension and camphor concentration is assumed to be γ = Γ 0 − Γ u , (8) where Γ₀ is the surface tension of pure water, and Γ is a positive constant connecting the camphor concentration and surface tension [37–39]. Using this relationship, the force F k ( c o n c ) is described as F k c o n c = − Γ ∮ ∂ Ω k u r ′ n r ′ d ℓ ′ (9) = − Γ ∬ Ω k ∇ ′ u r ′ d A ′ . (10)

F k ( i n t ) is the internal force that maintains the configuration of the camphor float. It can be described as the summation of the internal force, f _k,j , from the j-th disk as F k i n t = ∑ j = 0 n − 1 f k , j , (11) where f k , j = − f j , k , (12) and f k , j ∝ r k − r j . (13)

Here, it should be noted that f _k,k = 0.

By summing up the equations of motion for n disks, we obtain the evolution equation for r _c as ρ π R 2 n d 2 r c d t 2 = − κ π R 2 n d r c d t + ∑ k = 0 n − 1 F k c o n c = − κ π R 2 n d r c d t + π R 2 F n , (14) where F ( n ) is the force per unit area exerting on the camphor float with n-fold rotational symmetry. By calculating the vector product of the equations of motion with r _k − r _c and then summing up for n disks, we obtain the evolution equation for θ _c as ρ π R 2 ℓ 2 n d 2 θ c d t 2 = − κ π R 2 ℓ 2 n d θ c d t + ∑ k = 0 n − 1 r k − r c × F k c o n c = − κ π R 2 ℓ 2 n d θ c d t + π R 2 N n , (15) where N ( n ) is the torque per unit area exerting on the camphor float with n-fold rotational symmetry. Here, the operator “×” is defined so that a × b = a _x b _y − a _y b _x for a = a _x e _x + a _y e _y and b = b _x e _x + b _y e _y . It should be noted that the obtained time evolution equations for r _c and θ _c in Eqs 14, 15 does not explicitly include the terms originating from the internal force.

Finally, we derive the dimensionless version of our model. The dimensionless variables are defined as r ̃ = r D / a , r ̃ k = r k D / a , r ̃ c = r c D / a , R ̃ = R D / a , ℓ ̃ = ℓ D / a , t ̃ = a t , u ̃ = u S 0 / D , ρ ̃ = ρ D a Γ S 0 , κ ̃ = η D Γ S 0 , F ̃ k c o n c = F k c o n c Γ S 0 / D a , F ̃ n = F n Γ S 0 a / D D , N ̃ n = N n Γ S 0 / D , (16) and our model with dimensionless variables is summarized as ∂ u ̃ ∂ t ̃ = ∇ ̃ 2 u ̃ − u ̃ + ∑ k = 0 n − 1 1 π R ̃ 2 Θ R ̃ − r ̃ − r k ̃ , (17) ρ ̃ d 2 r ̃ c d t ̃ 2 = − κ ̃ d r ̃ c d t ̃ + 1 n π R ̃ 2 ∑ k = 0 n − 1 F ̃ k c o n c = − κ ̃ d r ̃ c d t ̃ + 1 n F ̃ n , (18) ρ ̃ d 2 θ c d t ̃ 2 = − κ ̃ d θ c d t ̃ + 1 n π R ̃ 2 ℓ ̃ 2 ∑ k = 0 n − 1 r ̃ k − r ̃ c × F ̃ k c o n c = − κ ̃ d θ c d t ̃ + 1 n ℓ ̃ 2 N ̃ n , (19) r ̃ k = r ̃ c + ℓ ̃ e θ c + 2 π k n , (20) F ̃ k c o n c = − ∮ ∂ Ω ̃ k u ̃ r ̃ ′ n r ̃ ′ d ℓ ̃ ′ (21) = − ∬ Ω ̃ k ∇ ̃ ′ u ̃ r ̃ ′ d A ̃ ′ . (22) Ω ̃ k = r ̃ | r ̃ − r ̃ k ≤ R ̃ . (23)

Hereafter, we adopt this dimensionless model with the tildes omitted.

In this section, we analyze our model on the camphor float with n-fold rotational symmetry, derived in the previous section. To enable us to analyze our model, we assume that the camphor disk radius R is sufficiently small and that the concentration field of camphor molecules are described as the function of the positions and velocities of the comprising n camphor disks. In this section, we first obtain the expression of the concentration field under the above assumptions, and then calculate the force F ( n ) and torque N ( n ) originating from the concentration field. Based on the ordinary differential equations obtained under the assumptions, we perform the linear stability anaysis. Finally, we discuss the motion of the camphor float by calculating the force and torque in the case that the camphor float is moving at a constant velocity with a constant characteristic angle θ _c .

In order to confirm the analytical results, we perform the numerical calculation. We used the dimensionless model in Eqs 17–23. For the numerical calculation, we used the smoothed function; the source part and the force are changed as ∂ u ∂ t = ∇ 2 u − u + ∑ k = 0 n − 1 1 π R 2 G r − r k , (61) F k c o n c = − ∬ R 2 ∇ ′ u r ′ G r ′ − r k d A ′ . (62) where G r = 1 2 1 + tanh R − r δ . (63)

It should be noted that G(r) converges to Θ(R − r) for the limit of δ → + 0.

Numerical calculation was performed with the explicit method with the time step Δt = 0.000 5, and the spatial mesh Δx = 0.05. The calculation area for the concentration field is set to 20 × 20 with periodic boundary conditions. The parameters are fixed as ρ = 0.01, R = 0.2, and δ = 0.05, and κ and ℓ were varied as parameters. The initial condition for the concentration was u = 0 in the whole region. The initial conditions for the velocity, characterisitic angle, and angular velocity were set as v _c = 0.5 e _x + 0.01 e _y , θ _c = 0.5, and ω _c = 0.01. Since the initial condition was set so that the camphor float has a high velocity and low angular velocity, the float was easy to converge to the translational motion.

In Figure 2A, we first plot the stationary translational velocity against κ for a single camphor disk. In this case, we could only consider the translational motion but not the rotational motion. Thus, the initial condition for the velocity was set as v _c = 0.5 e _x + 0.01 e _y . We confirmed that the motion of the single camprhor disk converged to a translational motion with a constant velocity at t = 100, even when the initial condition was changed. Thus, we plotted v c at t = 100 for the stationary translational velocity. The stationary translational velocity was finite for the smaller κ, while it was almost zero for the larger κ. The bifurcation point with respect to the bifurcation parameter κ for the translational motion was slightly greater than 0.12, where the bifurcation was supposed to be classified into a supercritical pitchfork bifurcation.

We also calculated the cases with n = 2 and 3. In Figures 2B,C, we plotted the stationary translational velocity for various ℓ. Also in these cases, the velocity at t = 100 was adopted as the stationary translational velocity, at which we confirmed that the velocity was converged. For the camphor float with n = 2, the bifurcation point was greater for the smaller ℓ as shown in Figure 2B. This agrees with the theoretical results in Eq. 40. For the camphor float with n = 3, the bifurcation point had a minimum value at a certain ℓ around 1.5 as shown in Figure 2C. This also corresponds to the theoretical result in Eq. 45, which indicates that the bifuracation point has the minimum value at ℓ = 1.377 98… .

In Figure 3, the time series of the velocity v c , angular velocity ω _c , and the angle difference ψ between the direction of velocity and the characteristic angle θ _c are plotted for the camphor float with n = 2. The angle difference ψ is defined as ψ = ϕ − θ c , (64) where ϕ is defined as v c = v c cos ⁡ ϕ e x + sin ⁡ ϕ e y . We set κ = 0.12 since we wanted to discuss the behavior close to the bifurcation point for the translational motion of a single camphor disk. The convergence to the stationary translational motion was faster for small ℓ, whereas the system almost reached the stationary states around t = 30 even in the case of ℓ = 2. The angular velocity decayed to zero and the angle difference ψ approached π/2, which corresponds to the translational motion in the direction perpendicular to the line connecting the two camphor disks.

As for the camphor float with n = 3, the time series are shown in Figure 4. It should be noted that the time range in Figure 4 is ten times as long as that in Figure 3, since the convergence of the angle difference ψ was much slower than the case with n = 2, though the convergence of the velocity and angular velocity was as fast as that in the case with n = 2. The angle difference ψ converged to zero for ℓ = 0.25, 0.5, 0.75, and 1.0, while it converged to − π/3 for ℓ = 2.0. This agrees with the analytical results in Eqs 59, 60. These results show that the smaller camphor float with n = 3 moved with one disk at its front, while the larger camphor float moved with one disk at its rear, which agrees with the analytical results. In the case with ℓ = 1.5, the convergence of ψ was very slow, which can be understood from the fact that ℓ = 1.5 is close to the threshold value ℓ _3,t ≃ 1.464 58… . Considering that 1.5 > ℓ _3,t , the angle difference should converge to ± π/3 but it appeared to approach 0. We consider that the difference between the theoretical and analytical results owes the slight difference in the setting; e.g. the finite size of the camphor disk, etc.

In Figure 5, the snapshots of the camphor concentration u together with the camphor disk positions are shown for the cases of the camphor float with n = 2 and 3 for ℓ = 0.5. Considering that the camphor concentration profile has a tail in the opposite direction of the traveling direction, the obtained concentration field also indicated that the angle difference ψ between the velocity and the characteristic angle was π/2 for the case with n = 2, and 0 for the case with n = 3.

In the analysis in Section 3, we assume that the system is close to the bifurcation point. In such a situation, our analytical results match the results by numerical calculation; i.e., the camphor float with n = 2 moves in the perpendicular direction to the line connecting the two camphor disks. For the case with n = 3, the smaller camphor float moves with one disk at its front, while the larger camphor float moves with one disk at its rear. However, in the condition far from the bifurcation point, the system often becomes multistable and one of several modes of motions is selected depending on the initial condition. Actually, the numerical results that exhibit multistability are shown in Figures 6, 7 for n = 2 and 3, respectively. The parameters were set as κ = 0.08 and ℓ = 1.5. In the case with n = 2 shown in Figure 6, the translational motion and spinning motion were observed depending on the initial conditions. The direction of the translational motion was always perpendicular to the line connecting the two disks. In the case with n = 3, the two types of translational motions, i.e., the motion with one disk at its front and that with one disk at its rear, and the spinning motion were observed depending on the initial conditions. Such multistability could be observed in the wider parameter regions for the larger ℓ, but the system was often monostable for the smaller ℓ.

In our model, the concentration field around one camphor disk has the leading term described by the second-kind modified Bessel function as far as the traveling velocity is low. This means that the concentration field decays exponentially. Taking this into consideration, the coupling between the direction of motion and the camphor float configuration should be stronger for the smaller camphor float. Thus, the smaller camphor float can easily change the moving direction and quickly takes the preferable direction. In contrast, for the larger camphor float, the interaction is so small that it takes much time for the camphor float to take the preferable direction. Moreover, the system can often show the multistability for the larger camphor float as discussed in the previous paragraph. Threfore, the present discussion can be adequately adopted for the smaller camphor float.

The analytical results for the cases with n = 2 and n = 3 are essentially different; the bifurcation point for the translational motion depends on the moving direction for n = 2, while it is equivalent for any direction for n = 3. Moreover, we calculated the force and torque exerting on the camphor float which is moving at a constant velocity V without rotation. For n = 2, the leading terms of the force and torque which depend on the direction of the motion are at the order of V and V ², respectively. In contrast, for n = 3, they are at the order of V ⁵ and V ³, respectively. This reflects that the linear stability analysis can lead to the nontrivial results for n = 2 but not for n = 3. Considering the orders of V in the force and torque, the configuration and the velocity of the camphor float are more strongly coupled with each other for n = 2 than for n = 3 near the bifurcation point as shown in Figures 3, 4. In addition, the force is more greatly affected to the dynamics than the torque for n = 2, and thus the preferable configuration is achieved by changing the velocity direction. In contrast, the torque is more greatly affected to the dynamics than the force for n = 3, and thus the preferable direction is achieved by changing the configuration by spinning. For the camphor float with n ≥ 4, we can perform the similar approach. However, for greater n, the effect of the configuration on the translational motion becomes smaller. Therefore, the analysis for the case with smaller n is more important, and here we especially discuss on the cases with n = 2 and n = 3.

We believe that the present theoretical results can be confirmed by the experiments. In the dimensionless model, the length and time are scaled by the diffusion length D / a and the characteristic sublimation time 1/a, respectively. From the previous results [37], they correspond to ca. 100 mm and ca. 50 s. Using these scale units, the typical size of the camphor float corresponds to ca. 100 mm in the actual system, since the dimensionless size of the camphor float adopted in this study is on the order of 1. As for the velocity, the obtained dimensionless value is on the order of 1 from Figures 2–7, which corresponds to ca. 2 mm/s. This value is rather slower than the typical velocity observed in experiments, ca. 10 mm/s, but these values are of the same order. Considering that the size of the camphor float is difficult to exceed 100 mm, the scale available in the experimental setup should be less than 1. Moreover, as discussed in the previous paragraph, the coupling between the configuration and the translational motion is stronger for the small camphor float. Therefore, we suppose that the experimental confirmation is possible only for the camphor float with small ℓ. Considering that the effect of the shape is stronger for n = 2, we suppose that the shape effect due to the 2-fold rotational symmetry can be easily confirmed, but that the effect for n = 3 or higher is more difficult to confirm.

In our previous work, we showed that an elliptic camphor particle moves in its minor-axis direction in theoretical analysis, numerical calculation, and experiments [9, 33, 34]. The symmetric property of the distribution of the mass and the camphor molecule source for an elliptic camphor particle is similar to that of the camphor float with 2-fold rotational symmetry discussed in the present study. In our result, it moves in the direction perpendicular to the line connecting the two camphors. Thus, the present results and our previous results show good correspondence. As for the camphor float with 3-fold rotational symmetry, it corresponds to the triangular camphor particle discussed in our previous study [35] from the viewpoint of the symmetric property. In our previous work, we concluded that the smaller triangular camphor particle moves in the direction of the corner for the smaller particle, while it moves in the direction of the side for the larger particle. The present results correspond to our previous ones in that the preferable direction changes depending on the size by considering the deviation of the mass and camphor molecule source. Taking the present analytical results into consideration, we suggest that the low-wave-number modulation from a symmetric circular shape most effectively affects the self-propelled motion, and the high-wave-number modulation is not so important.

In the present study, we theoretically investigated the self-propulsion of a camphor float with n-fold rotational symmetry, which comprises a rigid light circular plate attached with n camphor disks along a periphery with an equivalent spacing. We constructed a mathematical model for the camphor float, and analyzed it under the assumptions that the camphor disk radius is negligibly small and that the camphor concentration only depends on the positions and velocities of the camphor disks. From the linear stability analysis and the calculation of the force and torque for the situation that a camphor float is moving at a constant velocity with a certain characteristic angle, we concluded that the camphor float with 2-fold rotational symmetry moves in the direction perpendicular to the line connecting the two camphor disks. As for the camphor float with 3-fold rotational symmetry, the smaller float moves with one camphor disk at the front, while the larger float moves with one camphor disk at the rear.

As the extensions of the present study, we consider the particles with an asymmetric configuration. By considering such particles, we can generally discuss the manner of coupling between the motion and the symmtric featrures of self-propelled particles. We are also interested in the collective motion of the particles with n-fold rotational symmetry. In such systems, we expect some interesting structures with a nematic order and/or a hexagonal order can be realized. The symmetric properties embedded in the external condition such as the system boundary and the spatial modulation of the parameters also seem to be interesting. We hope the present results will provide the fundamental knowledge to understand the above-mentioned more complex systems, and moreover, will help understanding the universal mechanism on the shape effect on self-propelled motions driven by the self-generated concentration field.