Bifurcation in the angular velocity of a circular disk propelled by symmetrically distributed camphor pills
Yuki Koyano, Hiroyuki Kitahata, Marian Gryciuk, Nadejda Akulich, Agnieszka Gorecka, Maciej Malecki, Jerzy Gorecki
aa r X i v : . [ n li n . AO ] J a n Bifurcation in the angular velocity of a circular disk propelled bysymmetrically distributed camphor pills
Yuki Koyano, a) Hiroyuki Kitahata, Marian Gryciuk, Nadejda Akulich, Agnieszka Gorecka, Maciej Malecki, and Jerzy Gorecki Department of Physics, Chiba University, Chiba 263-8522, Japan Institute of Physical Chemistry, Polish Academy of Sciences, Warsaw 01-224,Poland Department of Chemistry, Technology of Electrochemical Production and Electronic Engineering Materials,Belarusian State Technological University, Minsk 220006, Belarus School of Physics and Astronomy, Monash University, Clayton, Victoria 3800,Australia (Dated: 29 January 2019)
We studied rotation of a disk propelled by a number of camphor pills symmetrically distributed at its edge.The disk was put on a water surface so that it could rotate around a vertical axis located at the disk center.In such a system, the driving torque originates from surface tension difference resulting from inhomogeneoussurface concentration of camphor molecules released from the pills. Here we investigated the dependence ofthe stationary angular velocity on the disk radius and on the number of pills. The work extends our previousstudy on a linear rotor propelled by two camphor pills [Phys. Rev. E, 96, 012609 (2017)]. It was observedthat the angular velocity dropped to zero after a critical number of pills was exceeded. Such behavior wasconfirmed by a numerical model of time evolution of the rotor. The model predicts that, for a fixed frictioncoefficient, the speed of pills can be accurately represented by a function of the linear number density ofpills. We also present bifurcation analysis of the conditions at which the transition between a standing anda rotating disk appears.
Camphor is one of many substances that forma layer on the water surface and modify the sur-face characteristics. The presence of camphormolecules at water surface reduces its surface ten-sion. The camphor surface concentration profileresults from the balance between camphor releasefrom the source, its transport, evaporation anddissolution. Inhomogeneities in camphor concen-tration around a floating object activate the mo-tion because the object is propelled towards theregion characterized by the lowest concentration.In this paper, we study rotation of a disk pro-pelled by a number of camphor pills. The disk an-gular velocity nonlinearly depends on the numberof pills and falls to zero when the number of pillsexceeds the critical value. The developed modeltreats the time evolution of camphor surface con-centration as a reaction-diffusion process. It re-produces qualitatively experimental observations,which confirms model usefulness for simulationsof systems with surface interactions. Moreover,the model can be reduced and allows for an an-alytical investigation of bifurcation between thestanding and the rotating states of a disk. We be-lieve that our results are important because theydescribe a realistic complex system for which abifurcation can be investigated analytically. a) Electronic mail: [email protected].
I. INTRODUCTION
Studies on self-propelled objects have become popularin the recent years because the behavior of many suchsystems shows similar characters of motion to that ex-pressed by living organisms. Self-propelled motion canbe observed in systems with embedded asymmetry ofsystem structure and interactions. For example, Janusparticles, characterized by different rates of reactions atdifferent parts of their surface, can move in the direc-tion determined by the chemical activity . There arealso objects in which the boundaries direct a jet of reac-tion products and force the motion . The self-propelledmotion can be also observed for symmetric objects inwhich the symmetry is broken by processes that gener-ate the motion. Such systems include droplets whereBelousov-Zhabotinsky (BZ) reaction proceeds . The in-terfacial tension between a droplet and the surroundingoil phase is related to the level of catalyst oxydization .If a droplet is sufficiently large then homogeneous oscilla-tions are observed and for yet larger droplets, a propagat-ing excitation pulses can appear . The related changesin interfacial tension generate a jump of the droplet inthe direction of pulse propagation . However, since thedirection of an excitation pulse is random, a symmetricBZ droplet can be shifted in a stochastic direction andthere are no factors that can stabilize the direction ofmotion.In this paper, we are concerned with self-propelledmotion induced by interfacial phenomena related to dy-namically changing surface concentration of camphormolecules. It is known that if a piece of camphor is placedon the water surface then camphor molecules hardly dis-solve in water, but the majority of them forms a layeron the water surface . In typical experimental condi-tions, this layer is unstable, because camphor moleculescontinuously evaporate. The water surface tension de-creases as camphor surface concentration increases .As a consequence, the force acting on a camphor pieceis directed towards the neighboring region with the low-est camphor surface concentration. One of the simplestand most known camphor-propelled objects is a cam-phor boat i.e., a boat-shaped piece of plastic with a bitof camphor glued at its stern . Such configurationof camphor-propelled objects breaks system symmetry.The surface concentration of released camphor moleculesaround the stern is higher than that around the bow,which decreases the surface tension in the stern area. Asa result, the boat moves forward.Geometrically symmetric objects can be also propelledby camphor pieces, because there is a positive feedbackbetween the generated force (or torque) and the directionof object motion. Let us consider a camphor disk reclin-ing on the water surface. It releases camphor moleculesaround, but both formation of a camphor layer and theevaporation of camphor molecules are subject to fluctu-ations. If an area characterized by a low surface concen-tration of camphor appears close to the disk, then thedisk is shifted towards the area because it is attracted bythe region with higher interfacial tension. When the diskis shifted from the original position, the surface camphorconcentration in front of the disk is lower than in theregion behind the disk, because the area in front of thedisk has been more distant from the camphor source thanthe region behind the disk. Therefore, the disk motioncontinues up to the moment the disk hits the boundaryor it is repelled by water meniscus near the boundary.Studies on self-propelled rotational motion are inter-esting since such motion occurs in a confined space, thuseffect of boundaries can be neglected. There have beenseveral reports on systems that show spontaneous rota-tion including systems with broken chiral symmetry and systems in which rotation occurs through the spon-taneous breaking of chiral symmetry . For camphordriven systems, it has been demonstrated that fluc-tuations of surface camphor concentration can induce ini-tial rotation that is supported by the positive feedbackbetween the direction of motion and the concentrationgradient, like for the translational motion of a camphordisk mentioned above.However, there have not been too many studies inwhich mathematical modeling of self-propelled rotationalmotion has been compared with experimental results.In this respect, systems that are propelled by camphorpieces are worth considering because they can be ana-lyzed using a simplified model of their time evolution.This model is based on a reaction-diffusion equation forcamphor surface concentration coupled with the Newto-nian equation of motion for the camphor pieces .In our previous paper , we considered a camphor ro- tor with two camphor pills at the ends of a plastic stripe.The pills were floating on the water surface, whereas thestripe was elevated above the surface. The system wasallowed to rotate around a vertical axis at the center ofthe stripe. We observed that such rotor can move only af-ter a distance between the camphor pills was larger thanthe critical one. For this system, the rotor radius can beconsidered as a bifurcation parameter. The mathemati-cal model of the spontaneous symmetry breaking can beformulated in terms of pitchfork bifurcation in dynamicalsystems. Here, as a generalization of previously studiedproblem, we consider disk-shaped rotors powered by anumber of camphor pills. The pills are symmetricallydistributed at the disk edge. There are two parametersthat describe the system: the disk radius and the num-ber of pills. It can be expected that a disk propelled bygreater number of pills rotates at a higher angular ve-locity. On the other hand, a disk with many camphorpills seems equivalent to a disk with a continuous cam-phor source along its edge, which clearly does not rotate.Therefore, there is a question on how the angular velocityof a disk depends on the number of camphor pills. Wehave performed experiments and analyzed the stationaryangular velocity of the disk as a function of the both pa-rameters. The results are reported in Section II. It turnsout that, for small disk radius, the angular velocity dropsto zero when the number of pills is large. If the disk islarge, then the angular velocity weakly depends on thenumber of pills attached in the same range of pill num-bers. In Section III, we present a mathematical modeldescribing the disk rotation. Results of numerical simu-lations presented in Section IV allowed us to determinethe values of model parameters for which the qualitativeagreement with experimental results is obtained. Sec-tion V is concerned with the analytical methods used tostudy disk evolution and with the analysis of bifurcationbetween the rotating and the still disk. Finally in theSection VI, we present numerical arguments that, for afixed friction coefficient, the speed of pills can be accu-rately represented by a function of a single argument: thelinear number density of pills. We demonstrate that suchbehavior can be found in experimental results. II. EXPERIMENTS
We study the angular velocity of a disk propelled by anumber of camphor pills located on water surface. Thesystem is illustrated in Fig. 1. The disk could rotatearound a vertical axis located at its center. The pillswere symmetrically distributed close to the disk edge.They were glued to the columns located below the diskssuch that the pills were in contact with water, whereasthe disk was elevated over the water surface to avoid gen-eration of an extensive hydrodynamic flows by the mov-ing disk . The camphor surface concentration on wateris determined by a number of physicochemical processessuch as the inflow of camphor molecules from the pills, (a)(b)(c) FIG. 1. Disk with pills glued at the bottom of supportingcolumns studied in experiments: (a) the side view, (b) thetop view and (c) the slanted view of the disk propelled by4 pills (the pills are marked yellow). The disk could rotatearound the vertical axis fixed at the tank center. The columnswere made of bended plastic stripe. diffusion at the surface, evaporation into the air, and dis-solution into the water . The randomness embedded inthese processes can lead to a nonuniform surface tensiongenerating the driving torque of the disk. Like in thecamphor disk motion discussed in Introduction, there isa positive feedback between the disk angular shift andthe force generating the torque.We investigated the stationary angular velocity de-pending on the number of pills N and on the distancebetween the axis and the pill center ℓ .The pills of mass m and radius ρ were tangent to thedisk, so the disk radius R = ℓ + ρ . In the experiments, weused 150 ml of water poured into a square tank (tank side12 cm) so that the water level was 1 cm. Water was puri-fied using a Millipore system (Elix 5) and its temperaturewas 22 ± ◦ C. In the experiments, we used commerciallyavailable camphor (99% purity, Sigma-Aldrich) withoutfurther purification. The pills were made by pressingcamphor in a pill maker. The camphor pill radius was ρ = 0 .
15 cm and the height was 0.1 cm.Figure 2 illustrates the analysis of experimental data.A marker was attached on the disk surface. Its posi-tion was recorded on a movie and analyzed using ImageJprogram . We applied two methods to obtain the an-gular velocity. In one of them, we calculated the angu-lar velocity as a function of time using positions of the (a) C o l o r i n t e n s it y A ngu l a r v e l o c it y (r a d / s ) FIG. 2. Illustration of the experimental data processing. (a)A disk ( R = 2 cm, N = 16) seen from above. A black markerglued on the disk was used to measure the angular velocity.(b) Time dependent color intensity at the center of red line (cf.(a)) during the disk rotation read from the frames of filmedexperiment. The minima correspond to moments when themarker crossed the line. (c) Angular velocity of the disk as afunction of time. marker in the consecutive frames of the movie. In thesecond method, we analyzed time-dependent color inten-sity in the region marked as a red line in Fig. 2 (a). Weobserved a minimum in color intensity every time themarker passed this region. The angular velocity aver-aged over a single rotation was calculated from the timedifference corresponding to the successive minima. Bothmethods gave similar results. In a typical experiment,the angular velocity was quite stable (cf. Fig. 2 (c)) andweakly depended on time at the time scale of a few min-utes, which corresponded to over 100 rotations in thepresented case. A small decrease in angular velocity withtime can be related to the increase in camphor concentra-tion in bulk aqueous phase. There is a slow dissolutionof camphor into the bulk aqueous phase. The dissolved N A ngu l a r v e l o c it y ω (r a d / s ) FIG. 3. Experimental results on the angular velocity ω ( ℓ, N )as a function of the number of camphor pills N for 3 selectedvalues of ℓ – the distance between disk axis and the dot center(0.85, 1.35 and 1.85 cm). Dots represent experimental dataand curves show their fit using a quadratic polynomial withzeros at N = 14, 22 and 25, respectively. Blue, green, andorange points correspond to ℓ = 0 .
85, 1 .
35, and 1 .
85 cm,respectively. molecules migrate to the surface and contribute to thesurface concentration of camphor. The contribution ishomogeneous and reduces gradient of camphor surfaceconcentration resulting from local release and evapora-tion. As a result, the inhomogeneities in surface tensionare decreased and so is the torque.Figure 3 summarizes the experimental results obtainedfor 3 different disk radii. It shows the average angularvelocity ω ( ℓ, N ) measured within first 6 minutes of rota-tion. For a small number of pills ( N ≤ N was a decreasing function of ℓ . Con-sidering the dependence of ω ( ℓ, N ) on N , we observedthat for small disk radii the angular velocity rapidly de-creased with the increasing number of pills. The disk of ℓ = 0 .
85 cm powered by N = 14 pills randomly moved inboth directions, but it did not show unidirectional rota-tion lasting more than a second. For a larger disk radii( ℓ = 1 .
35 cm and ℓ = 1 .
85 cm), the angular velocity wasslowly decreasing with N in the range of 5 < N <
20. For ℓ = 1 .
35 cm and N = 21, the angular velocity was aroundthe half of the value observed at small N . The disk didnot rotate when N ≥
22. In the case of ℓ = 1 .
85 cm,we observed a stable rotation with a very long period of115 s for N = 24 and the disk stopped at N = 25. Inall experimental results ω ( N = 3) < ω ( N = 4), so wecannot exclude a maximum of angular velocity at small N ∈ { , } . The lines in Figure 3 show a fit using the sec-ond order polynomials with zeros at the smallest numberof pills for which the disk was not rotating. III. MODEL
In this section, we introduce a mathematical model fora disk propelled with N camphor pills attached. We setthe coordinates so that the center of the disk, i.e., therotation axis, is located at the origin. The motion of thedisk is described by the characteristic angle φ ( t ). All N camphor pills are located at the distance of ℓ from theorigin and have equal spacing. Thus the position of j -thcamphor pill, ℓ j ( t ), can be described as: ℓ j ( t ) = ℓ (cid:20) cos (cid:18) φ ( t ) + 2 πjN (cid:19) e x + sin (cid:18) φ ( t ) + 2 πjN (cid:19) e y (cid:21) , (1)for j = 0 , · · · , N −
1, where e x and e y denote the unitvectors in x - and y -directions, respectively.The camphor pills release camphor molecules to thewater surface, and the camphor molecules diffuse at thewater surface. Some camphor molecules sublimate to theair and some dissolve to the water bulk phase. All theseprocesses are taken into account in the equation for thetime evolution of the camphor surface concentration: ∂c ( r , t ) ∂t = ∇ c ( r , t ) − c ( r , t ) + N − X j =0 f ( r ; ℓ j ) , (2)where c ( r , t ) is the surface concentration of camphormolecules at the position r and time t . The first termcorresponds to the diffusion, and the second one to subli-mation to the air and dissolution to the water bulk phase.The last term f ( r ; ℓ j ( t )) denotes the release of camphormolecules from the j -th pill, which is represented as f ( r , ℓ j ( t )) = 1 πρ Θ ( ρ − | r − ℓ j ( t ) | ) , (3)where ρ corresponds to the camphor pill radius and Θ( ξ )is the Heaviside’s step function, i.e., Θ( ξ ) = 1 for ξ ≥ ξ ) = 0 for ξ < I d φdt = − η r dφdt + T . (4)Here, I is the momentum of inertia of the disk, and it isapproximated using the mass of one camphor pill and itssupporting column m as I = N mℓ = πN σρ ℓ , (5)where σ is the average surface density. − η r dφ/dt is thetorque originating from the friction force working on thecamphor pills and η r is described using the friction coef-ficient per unit area, κ , η r = πN κρ ℓ , (6)as is derived in the previous paper . T is the torqueexerting on the disk, which is represented as T = N − X j =0 ℓ j × F j , (7)where F j is the driving force of the j -th camphor pillinduced by the surface tension gradient. It should benoted that constraint force should work on each pill tomaintain the composition of the disk, but the direction ofthe constraint force working on the j -th pill is the sameas ℓ j , and therefore it does not affect the torque. F j isdescribed as F j = Z π γ ( c ( ℓ j + ρ e ( θ ))) e ( θ ) ρdθ, (8)where γ ( c ) is the surface tension depending on the cam-phor surface concentration, and e ( θ ) is a unit vector inthe direction of θ , i.e., e ( θ ) = cos θ e x + sin θ e y . For sim-plicity, we set γ ( c ( r , t )) = γ − kc ( r , t ) , (9)where γ is water surface tension, and k is a positiveconstant. Hereafter, we set k = 1.Taken in all, we obtain πN σρ ℓ d φdt = − πN κρ ℓ dφdt + T , (10)or σ d φdt = − κ dφdt + 1 πN ρ ℓ T . (11)All variables in the equations above are dimensionless.Let us assume that D is the diffusion constant of camphormolecules, a is the combined rate of sublimation and dis-solution, and f is the release rate of camphor moleculesfrom one camphor pill. The dimensionless variables aredefined such that the time unit is the characteristic timeof combined sublimation and dissolution 1 /a , the lengthunit is the diffusion length p D/a , and the concentrationunit is the ratio f /a . Our approach approximately treatsthe hydrodynamic effects. Due to surface tension gradi-ents, the Marangoni flow should appear in the system.In our approach, D is an “effective” diffusion constantthat allows to include the camphor transport related tothe Marangoni effect . IV. NUMERICAL SIMULATIONS
Based on the model introduced in the previous section,we performed numerical calculation. The release of cam-phor molecules, described by Eq. (3), was approximatedusing the expression: f ( r , ℓ j ( t )) = 12 πρ (cid:20) − ( | r − ℓ j | − ρ ) δ (cid:21) , (12) Number of camphor pills N S t a ti on a r y a ngu l a r v e l o c it y ω ℓ = 1.35 ℓ = 0.85 ℓ = 1.85(a) Number of camphor pills N S t a ti on a r y a ngu l a r v e l o c it y ω ℓ = 1.35 ℓ = 0.85 ℓ = 1.85(b) Number of camphor pills N S t a ti on a r y a ngu l a r v e l o c it y ω ℓ = 1.35 ℓ = 0.85 ℓ = 1.85(c) 110 202020 FIG. 4. Results of numerical calculation on the stationary an-gular velocity, ω ( ℓ, N ), expressed in the dimensionless units,depending on the number of camphor pills, N , for each ℓ .Cyan, green, and red plots correspond to ℓ = 0 .
85, 1 . .
85, respectively. (a) κ = 0 .
1, (b) κ = 0 .
01, and (c) κ = 0 . in order to reduce the effect of discretization, where δ isa positive constant for smoothing.The parameters were set to be σ = 0 . ρ = 0 . δ = 0 . ℓ , the number of camphor pills N , and the coefficient of the friction κ were varied asparameters. The time evolution was calculated with theEuler algorithm, and the diffusion was calculated withthe explicit method. Time step was 10 − and the spatialmesh was 0 . φ = 1 and dφ/dt = 0 .
1. The terminal angular velocity ω is set tobe dφ/dt at t = 100, since we confirmed that the systemwas close to the stable stationary state at t = 100.The numerical results are shown in Fig. 4, in whichwe simultaneously plotted ω against N for ℓ = 0 . , . .
85. It should be noted that the ratios betweenthe camphor pill radius and the disk radius are the samein experiments and numerical simulation. In Fig. 4 (a),(b), and (c), we plotted the results on ω against N fordifferent values of κ . The intersections of the plots withdifferent ℓ changed depending on κ , and therefore wehope we can estimate κ from the experimental results. Inthis case, the plot in (b) is most close to the experimentalresults, and thus we can estimate κ ≃ . V. THEORETICAL ANALYSIS
In this section, we derive the reduced evolution equa-tion for the spinning of a disk with N camphor pills. Weadopted Eq. (2) for the time evolution of concentrationfield, assuming that the support of function describingthe release of camphor molecules is infinitesimally small: f ( r , ℓ j ( t )) = δ ( r − ℓ j ( t )) , (13)instead of a finite size (cf. Eq. (3)). Here δ ( · ) is Dirac’sdelta function in a two-dimensional space. The sourceterm given by Eq. (13) delivers the same amount of cam-phor as those described by Eq. (3).In our previous study , we derived the explicit formof the concentration field expanded with respect to thevelocity, acceleration, jerk (i.e. the rate of change of ac-celeration), and so on of a camphor pill. The concentra-tion field at the position r originating from a camphorpill whose position is ℓ ( t ) is described as: c s ( r ; ℓ ) = 12 π K ( d ) − π K ( d ) h d · ˙ ℓ i + 116 π d K ( d ) h d · ¨ ℓ i − π d K ( d ) | ˙ ℓ | + 116 π K ( d ) h d · ˙ ℓ i + 132 π d K ( d ) | ˙ ℓ | h d · ˙ ℓ i − π K ( d ) h d · ˙ ℓ i + 132 π d K ( d ) h ˙ ℓ · ¨ ℓ i − π d K ( d ) h d · ˙ ℓ i h d · ¨ ℓ i − π d K ( d ) [ d · ... ℓ ] + (higher order terms) , (14) N A ( ℓ , N ) ℓ = 2 ℓ = 3 ℓ = 5 ℓ = 40 20100-0.1-0.3 Number of camphor pills N C ( ℓ , N ) ℓ = 2 ℓ = 3 ℓ = 5 ℓ = 40 201010.50 Number of camphor pills N ω ( ℓ , N ) ℓ = 2 ℓ = 3 ℓ = 5 ℓ = 4(a)(b)(c) -0.21.2 FIG. 5. Plots of (a) A ( ℓ, N ), (b) C ( ℓ, N ), and (c) ω ( ℓ, N )against N . The distance ℓ is set to be ℓ = 2, 3, 4, and 5. Theradius of camphor disks ρ is fixed to be 0.1. The value of κ isset to be 0 . where d = r − ℓ and d = | d | . K n is the second-kindmodified Bessel function of n -th order, and a dot (˙) meansthe time derivative. In this expansion, we assume thatthe camphor pill speed is sufficiently small.In the case of a disk with N camphor pills, the con-centration field is expressed by summing up the concen-tration field originating from each camphor pill since theevolution equation for concentration field is linear. Thus,the concentration field made by the disk whose center isat the origin is given by c ( r ) = N − X j =0 c s ( r ; ℓ j ) . (15)Considering that the driving force originates from theimbalance of surface tension, the driving force workingon the j -th pill F j can be calculated as follows:1 πρ F j = − N − X k =0 lim ρ → +0 πρ Z π c s ( ℓ j + ρ e ( θ ); ℓ k ) e ( θ ) ρdθ = 14 π (cid:18) − γ + log 2 ρ (cid:19) ˙ ℓ j − π ¨ ℓ j − π (cid:12)(cid:12)(cid:12) ˙ ℓ j (cid:12)(cid:12)(cid:12) ˙ ℓ j + 148 π ... ℓ j − X k = j ∇ c s ( ℓ j ; ℓ k ) + O ( ρ ) , (16)where ρ is considered to be an infinitesimally small pa-rameter corresponding to the radius of a camphor pill. Itshould be noted that ∇ c s ( ℓ j ; ℓ k ) = ∇ c s ( r ; ℓ k ) | r = ℓ j .By explicitly calculating ∇ c s ( ℓ j ; ℓ k ) and substitutingthe results into Eq. (16), the torque working on the j -thpill, T j , is obtained by taking the vector product of theradial vector and the force, T j = ℓ j × F j . (17)Since the torques acting on pills are identical due to thegeometric symmetry, we finally obtain the total torque T working on the disk T = N − X j =0 T j = N T . (18)Therefore the reduced equation for the time evolution ofthe angle describing disk position φ ( t ) reads:( σ + B ( ℓ, N )) d φdt = ( A ( ℓ, N ) − κ ) dφdt + C ( ℓ, N ) (cid:18) dφdt (cid:19) , (19)where A ( ℓ, N ), B ( ℓ, N ), and C ( ℓ, N ) are given as A ( ℓ, N ) = 14 π (cid:18) − γ Euler + log 2 ρ + N − X j =1 " − ℓ K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πjN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) sin (cid:0) πjN (cid:1) (cid:12)(cid:12) sin (cid:0) πjN (cid:1)(cid:12)(cid:12) + K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πjN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) cos (cid:18) πjN (cid:19)(cid:21)(cid:19) , (20) B ( ℓ, N )= 116 π − N − X j =1 (cid:20) K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πjN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ sin (cid:18) πjN (cid:19) − K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πjN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πjN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) cos (cid:18) πjN (cid:19)(cid:21)(cid:19) , (21) C ( ℓ, N )= 1192 π (cid:0) − ℓ − N − X j =1 (cid:20) K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πjN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ sin (cid:18) πjN (cid:19) cos (cid:18) πjN (cid:19) −K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πjN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ sin (cid:0) πjN (cid:1)(cid:12)(cid:12) sin (cid:0) πjN (cid:1)(cid:12)(cid:12) +4 K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πjN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πjN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) × (cid:26) − (cid:18) πjN (cid:19) + 4 sin (cid:18) πjN (cid:19)(cid:27) +8 K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πjN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ sin (cid:18) πjN (cid:19) cos (cid:18) πjN (cid:19)(cid:21)(cid:19) , (22)where γ Euler is Euler-Mascheroni constant ( ≃ . d φ/dt . The detailed process of calcu-lation is shown in Appendix A. It should be noted that B ( ℓ, N ) is always positive, while A ( ℓ, N ) and C ( ℓ, N )change their signs depending on the parameters. Thepositive value of B ( ℓ, N ) means that the stationary so-lution and its stability does not change in the viscositylimit when σ in Eq. (19) goes to zero.At the bifurcation point, the coefficient of dφ/dt inthe right side of Eq. (19) should be 0, i.e., κ = A ( ℓ, N ).The friction coefficient κ depends on the water level, butin our experiments the water level was fixed so it canbe regarded as a constant. The coefficient A ( ℓ, N ) isplotted against N in Fig. 5(a) with ℓ = 2, 3, 4, and 5.The coefficient A ( ℓ, N ) monotonically decreases with anincrease in N for every ℓ .If A ( ℓ, N ) is smaller than κ , then the disk stops. Thus,a monotonical decrease in A ( ℓ, N ) with an increase in N means that there is a critical number N c such that thedisk stops for N ≥ N c .The stable angular velocity ω ( ℓ, N ) is given by ω ( ℓ, N ) = ± p ( A ( ℓ, N ) − κ ) / ( − C ( ℓ, N )), if A ( ℓ, N ) − κ is positive and C ( ℓ, N ) is negative. Figure 5(c) showsthe plot of ω against N with ℓ = 2, 3, 4, and 5. Whenthe coefficient C ( ℓ, N ) is positive, then the bifurcationis subcritical, and a higher order term (e.g. fifth-orderterm) is supposed to suppress the divergence of the angu-lar velocity, which is out of scope of the present analysis. VI. DISCUSSION AND CONCLUSION
In the paper, we analyzed bifurcation between standingand rotating states of a disk propelled by a number ofcamphor pills.Experiments on the stationary rotation of a disk lo-cated on the water surface were motivation for thepresent study. The experiments revealed highly nonlin-ear behavior of the angular velocity ω ( ℓ, N ) as a func-tion of the disk size ℓ and of the number of pills pro-pelling it N . We observed that for a moderate numberof pills (5 ≤ N ≤
12) the angular velocity was a decreas-ing function of the disk size: larger disks rotated slowerthan smaller ones. For a fixed disk radius and a moderatenumber of pills, the angular velocity ω ( ℓ, N ) was a slowlydecreasing function of N . However, if the number ofpills approached a critical value N c then ω ( ℓ, N ) rapidlydropped, and the disk did not rotate for all N ≥ N c . Thevalue of N c depended on the disk radius, and it was 14,22 and 25 for ℓ = 0 .
85, 1 .
35 and 1 .
85 cm, respectively.We investigated the nonlinear behavior of the angularvelocity of a disk driven by multiple camphor pills fo-cusing on the bifurcation between rotation and rest. Itwas reproduced by a mathematical model describing thedynamics of camphor surface concentration, its influenceon the surface tension, and driving torque exerting onthe disk through concentration-dependent surface ten-sion. The model we have used describes hydrodynamiceffects and associated camphor transport via effective dif-fusion constant, that should be optimized for a particularsystem. It uses scaled variables selected in such way thatthe effective diffusion constant D as well as the combinedrates of camphor evaporation and dissolution in water areequal to 1. One of the adjustable model parameters isthe friction coefficient κ .The numerical results for the speed of pills v ℓ ( ℓ, N ) = ℓω ( ℓ, N ) as a function of the linear number density of pills N/ (2 πℓ ) are shown in Fig. 6. It can be seen that, outsidethe bifurcation region, v ℓ ( ℓ, N ) was the same for differentpairs ( ℓ , N ) that lead to the same linear number densityof pills. The linear number density of pills is a localquantity. It can be expected that the disk behavior solelydepends on the density when the diameter of the diskis greater than the characteristic length of the diffusion(in our case it is equal to 1) because the camphor pillshardly interact with the pills at the opposite side of thedisk. Such universal behavior of v ℓ ( ℓ, N ) observed forlarge disks, especially the density at which the transitionbetween rotation and rest occurs, can be used to estimatethe value of κ that matches experimental results.Figure 7 illustrates the phase diagram of disk behaviorin the variables κ and the linear number density of pills N/ (2 πℓ ) calculated for ρ = 0 .
15. Dots mark pairs ( κ , N c / (2 πℓ )) corresponding to the smallest number of pillsfor which v ℓ ( ℓ, N ) = 0. They separate the phase spaceinto regions where the disk rotates (below the dots) andwhere it does not (above the dots).The combined experimental results for v ℓ ( N/ (2 πℓ )) as (a) Linear number density of pills N / (2 π ℓ )0 41 2 34023 S p ee d v ℓ ℓ = 1.35 ℓ = 0.85 ℓ = 1.851(b) 400203010 Linear number density of pills N / (2 π ℓ )0 41 2 3 S p ee d v ℓ ℓ = 1.35 ℓ = 0.85 ℓ = 1.85(c) 150050100 Linear number density of pills N / (2 π ℓ )0 41 2 3 S p ee d v ℓ ℓ = 1.35 ℓ = 0.85 ℓ = 1.85 FIG. 6. Results of numerical calculation on the stationaryspeed, v ℓ = ℓω ( ℓ, N ), as a function of the linear number den-sity of pills, N/ (2 πℓ ). Cyan, green, and red plots correspondto ℓ = 0 .
85, 1 .
35, and 1 .
85, respectively. (a) κ = 0 .
1, (b) κ = 0 .
01, and (c) κ = 0 . Friction coefficient per unit area 0.001 0.10.014012 L i n ea r nu m b e r d e n s it y o f p ill s N / ( π ℓ ) κ Upper density limit
FIG. 7. Phase diagram illustrating different disk behaviordepending on κ and the linear number density of pills. Dotsmark pairs ( κ , N c / (2 πℓ )). The dashed line represents thehighest linear disk density for ρ = 0 .
15. In all cases repre-sented by points ( κ , N/ (2 πℓ )) below the dots, the disk ro-tates. For the systems characterized by ( κ , N/ (2 πℓ )) abovethe dots, the disk does not move. The numerical results for ℓ = 1 .
85 were used for the present plot. a function of the linear number density of pills are shownin Fig. 8. It can be seen that results for ℓ = 0 .
85 cm and ℓ = 1 .
35 cm nicely overlap, confirming universal behaviorof v ℓ ( N/ (2 πℓ )). The location of a critical point can beapproximated by selecting κ ∼ .
01. The majority of v ℓ ( N/ (2 πℓ )) values for ℓ = 1 .
85 cm was smaller thanthose for disks of smaller radii. Also the critical linearnumber density of pills for ℓ = 1 .
85 cm was smaller thanthose for ℓ = 0 .
85 cm or ℓ = 1 .
35 cm. Results shown inFig. 6 suggest that the value of κ corresponding to ℓ =1 .
85 cm should be larger than 0 .
01. The discrepancy canbe explained by the fact that the largest disk was heavierand the pills were immersed deeper in water, and thus alarger friction coefficient should be applied to describe itsmotion.We also discussed the model of disk rotation inthe limit where camphor pills are infinitesimally small.Within this assumption, we derived an analytical expres-sion for camphor concentration profile and reduced equa-tion of motion to the form that allows to find the bifurca-tion point and to estimate the stationary value of ω ( ℓ, N ).The results of the analytical model correctly reflect thebasic experimental results: the value of N c increases withthe disk radius and, outside the neighborhood of bifur-cation, ω ( ℓ, N ) is a decreasing function of ℓ . Also theshapes of ω ( ℓ, N ) as a function of N are similar to thoseseen in experiments (cf. Fig. 3 and Fig. 5(c)). ACKNOWLEDGMENTS
The authors are grateful to Professor S. Nakata forhis helpful comments. This work was supported by N / (2 π ℓ ) (1/cm)0 31 2 S p ee d v ℓ ( c m / s ) FIG. 8. Experimental results for v ℓ as a function of N/ (2 πℓ ).Blue, green, and orange points correspond to ℓ = 0 .
85, 1 . .
85 cm, respectively.
JSPS-PAN Bilateral Joint Research Program “ Spatio-temporal patterns of elements driven by self-generated,geometrically constrained flows” between Japan and thePolish Academy of Sciences and by the European UnionsHorizon 2020 research and innovation programme underthe Marie Skodowska-Curie grant agreement No 734276(N.A.) with additional support from the Ministry of Sci-ence and Higher Education of Poland, agreement no3854/H2020/17/2018/2. Another author (Y.K.) is grate-ful for the support within JSPS KAKENHI Grant Num-ber JP17J05270 and the Cooperative Research Programof “Network Joint Research Center for Materials and De-vices” No. 20181023.
Appendix A: Detailed calculation
In this section, we show the detailed derivation ofEq. (19) with Eqs. (20) to (22).In Eq. (16), the last listed term represents the forceoriginating from the camphor concentration releasedfrom the other pills, while the first four terms describethe force originating from the camphor concentration re-leased from the considered pill itself. Here we set F j = F (self) j + X k = j F (other) j,k , (A1)where1 πρ F (self) j = 14 π (cid:18) − γ + log 2 ρ (cid:19) ˙ ℓ j − π ¨ ℓ j − π (cid:12)(cid:12)(cid:12) ˙ ℓ j (cid:12)(cid:12)(cid:12) ˙ ℓ j + 148 π ... ℓ j + O ( ρ ) , (A2)and 1 πρ F (other) j,k = −∇ c s ( ℓ j ; ℓ k ) . (A3)0By explicitly calculating the gradient of c s ( r ; ℓ ) as ∇ c s ( r ; ℓ ) = ∂∂ r c s ( r ; ℓ ) (A4)= 12 π K ′ ( | r − ℓ | ) r − ℓ | r − ℓ | − π K ′ ( | r − ℓ | ) h ( r − ℓ ) · ˙ ℓ i r − ℓ | r − ℓ | − π K ( | r − ℓ | ) ˙ ℓ − π K ( | r − ℓ | ) h ( r − ℓ ) · ¨ ℓ i ( r − ℓ ) + 116 π | r − ℓ |K ( | r − ℓ | ) ¨ ℓ + 116 π K ( | r − ℓ | ) | ˙ ℓ | ( r − ℓ )+ 116 π K ′ ( | r − ℓ | ) h ( r − ℓ ) · ˙ ℓ i r − ℓ | r − ℓ | + 18 π K ( | r − ℓ | ) h ( r − ℓ ) · ˙ ℓ i ˙ ℓ − π K ( | r − ℓ | ) | ˙ ℓ | h ( r − ℓ ) · ˙ ℓ i ( r − ℓ ) + 132 π | r − ℓ |K ( | r − ℓ | ) | ˙ ℓ | ˙ ℓ − π K ′ ( | r − ℓ | ) h ( r − ℓ ) · ˙ ℓ i r − ℓ | r − ℓ | − π K ( | r − ℓ | ) h ( r − ℓ ) · ˙ ℓ i ˙ ℓ − π | r − ℓ |K ( | r − ℓ | ) (cid:16) ˙ ℓ · ¨ ℓ (cid:17) ( r − ℓ ) + 132 π K ( | r − ℓ | ) h ( r − ℓ ) · ˙ ℓ i h ( r − ℓ ) · ¨ ℓ i ( r − ℓ ) − π | r − ℓ |K ( | r − ℓ | ) h ( r − ℓ ) · ¨ ℓ i ˙ ℓ − π | r − ℓ |K ( | r − ℓ | ) h ( r − ℓ ) · ˙ ℓ i ¨ ℓ + 196 π | r − ℓ |K ( | r − ℓ | ) [( r − ℓ ) · ... ℓ ] ( r − ℓ ) − π | r − ℓ | K ( | r − ℓ | ) ... ℓ , (A5)we can obtain the explicit form of F (other) j,k . Here we usedthe property of modified Bessel function : z K ν ′ ( z ) + ν K ν ( z ) = − z K ν − ( z ).The torque T j,j working on the j -th camphor pill orig-inating from F (self) j , is given by1 πρ T j,j = lim r → ℓ (cid:20) π K ( | r − ℓ | ) (cid:16) ℓ × ˙ ℓ (cid:17) − π | r − ℓ |K ( | r − ℓ | ) (cid:16) ℓ × ¨ ℓ (cid:17) − π | r − ℓ |K ( | r − ℓ | ) | ˙ ℓ | (cid:16) ℓ × ˙ ℓ (cid:17) + 196 π | r − ℓ | K ( | r − ℓ | ) ( ℓ × ... ℓ ) (cid:21) = 14 π (cid:18) − γ Euler + log 2 ǫ (cid:19) ℓ ˙ φ − π ℓ ¨ φ − π ℓ ˙ φ + 148 π ℓ (... φ − ˙ φ ) . (A6)Here we used lim x → +0 K ( x ) = − γ Euler + log(2 /x ),lim x → +0 x K ( x ) = 1, lim x → +0 x K ( x ) = 2.From Eq. (18), we only have to obtain T consideringthe system symmetry, and it is calculated as: T = ℓ × F = ℓ × F (self)0 + ℓ × N − X k =1 F (other)0 ,k = T , + N − X k =1 T ,k . (A7) Therefore, T ,k is obtained from Eq. (A5) as11 πρ T ,k = − [ ℓ × ∇ c s ( ℓ ; ℓ k )]= 12 π K ′ ( | ℓ − ℓ k | ) ℓ × ℓ k | ℓ − ℓ k | − π K ′ ( | ℓ − ℓ k | ) (cid:16) ℓ · ˙ ℓ k (cid:17) ℓ × ℓ k | ℓ − ℓ k | + 14 π K ( | ℓ − ℓ k | ) (cid:16) ℓ × ˙ ℓ k (cid:17) − π K ( | ℓ − ℓ k | ) h ( ℓ − ℓ k ) · ¨ ℓ k i ( ℓ × ℓ k ) − π | ℓ − ℓ k |K ( | ℓ − ℓ k | ) (cid:16) ℓ × ¨ ℓ k (cid:17) + 116 π K ( | ℓ − ℓ k | ) | ˙ ℓ k | ( ℓ × ℓ k ) + 116 π K ′ ( | ℓ − ℓ k | ) (cid:16) ℓ · ˙ ℓ k (cid:17) ℓ × ℓ k | ℓ − ℓ k |− π K ( | ℓ − ℓ k | ) (cid:16) ℓ · ˙ ℓ k (cid:17) (cid:16) ℓ × ˙ ℓ k (cid:17) − π K ( | ℓ − ℓ k | ) | ˙ ℓ k | (cid:16) ℓ · ˙ ℓ k (cid:17) ( ℓ × ℓ k ) − π | ℓ − ℓ k |K ( | ℓ − ℓ k | ) | ˙ ℓ k | ( ℓ × ˙ ℓ k ) − π K ′ ( | ℓ − ℓ k | ) (cid:16) ℓ · ˙ ℓ k (cid:17) ℓ × ℓ k | ℓ − ℓ k | + 132 π K ( | ℓ − ℓ k | ) (cid:16) ℓ · ˙ ℓ k (cid:17) (cid:16) ℓ × ˙ ℓ k (cid:17) − π | ℓ − ℓ k |K ( | ℓ − ℓ k | ) (cid:16) ˙ ℓ k · ¨ ℓ k (cid:17) ( ℓ × ℓ k )+ 132 π K ( | ℓ − ℓ k | ) (cid:16) ℓ · ˙ ℓ k (cid:17) h ( ℓ − ℓ k ) · ¨ ℓ k i ( ℓ × ℓ k )+ 132 π | ℓ − ℓ k |K ( | ℓ − ℓ k | ) h ( ℓ − ℓ k ) · ¨ ℓ k i (cid:16) ℓ × ˙ ℓ k (cid:17) + 132 π | ℓ − ℓ k |K ( | ℓ − ℓ k | ) (cid:16) ℓ · ˙ ℓ k (cid:17) (cid:16) ℓ × ¨ ℓ k (cid:17) + 196 π | ℓ − ℓ k |K ( | ℓ − ℓ k | ) [( ℓ − ℓ k ) · ... ℓ k ] ( ℓ × ℓ k ) + 196 π | ℓ − ℓ k | K ( | ℓ − ℓ k | ) ( ℓ × ... ℓ k )= 14 π K ′ (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ sin (cid:0) πkN (cid:1)(cid:12)(cid:12) sin (cid:0) πkN (cid:1)(cid:12)(cid:12) + 18 π K ′ (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ sin (cid:0) πkN (cid:1)(cid:12)(cid:12) sin (cid:0) πkN (cid:1)(cid:12)(cid:12) ˙ φ + 14 π K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ cos (cid:18) πkN (cid:19) ˙ φ + 116 π K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ sin (cid:18) πkN (cid:19) (cid:18) ˙ φ (cid:18) cos 2 πkN − (cid:19) + ¨ φ sin 2 πkN (cid:19) − π K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) (cid:18) ¨ φ cos (cid:18) πkN (cid:19) − ˙ φ sin (cid:18) πkN (cid:19)(cid:19) + 116 π K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ sin (cid:18) πkN (cid:19) ˙ φ + 132 π K ′ (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ sin (cid:0) πkN (cid:1)(cid:12)(cid:12) sin (cid:0) πkN (cid:1)(cid:12)(cid:12) ˙ φ + 18 π K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ sin (cid:18) πkN (cid:19) cos (cid:18) πkN (cid:19) ˙ φ + 132 π K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ sin (cid:18) πkN (cid:19) ˙ φ − π K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) cos (cid:18) πkN (cid:19) ˙ φ + 1192 π K ′ (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ sin (cid:0) πkN (cid:1)(cid:12)(cid:12) sin (cid:0) πkN (cid:1)(cid:12)(cid:12) ˙ φ + 132 π K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ sin (cid:18) πkN (cid:19) cos (cid:18) πkN (cid:19) ˙ φ − π K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19) ˙ φ ¨ φ + 132 π K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ sin (cid:18) πkN (cid:19) (cid:20) ˙ φ (cid:18) cos (cid:18) πkN (cid:19) − (cid:19) + ¨ φ sin (cid:18) πkN (cid:19)(cid:21) ˙ φ − π K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) cos (cid:18) πkN (cid:19) (cid:20) ˙ φ (cid:18) cos (cid:18) πkN (cid:19) − (cid:19) + ¨ φ sin (cid:18) πkN (cid:19)(cid:21) ˙ φ + 116 π K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19) (cid:20) ˙ φ sin (cid:18) πkN (cid:19) − ¨ φ cos (cid:18) πkN (cid:19)(cid:21) ˙ φ − π K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19) (cid:20) (... φ − ˙ φ ) sin (cid:18) πkN (cid:19) + 3 ˙ φ ¨ φ (cid:18) cos (cid:18) πkN (cid:19) − (cid:19)(cid:21) + 124 π K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ sin (cid:18) πkN (cid:19) (cid:20) (... φ − ˙ φ ) cos (cid:18) πkN (cid:19) − φ ¨ φ sin (cid:18) πkN (cid:19)(cid:21) . (A8)2Here, we used the time derivative of Eq. (1). We thus have the expression for the torque T :1 πρ T = 1 πρ T , + N − X k =1 πρ T ,k = ℓ π − γ + log 2 ρ + N − X k =1 " − ℓ K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) sin (cid:0) πkN (cid:1) (cid:12)(cid:12) sin (cid:0) πkN (cid:1)(cid:12)(cid:12) + K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) cos (cid:18) πkN (cid:19) ˙ φ + ℓ π − N − X k =1 (cid:20) K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ sin (cid:18) πkN (cid:19) − K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) cos (cid:18) πkN (cid:19)(cid:21)(cid:19) ¨ φ + ℓ π − ℓ − N − X k =1 (cid:20) K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ sin (cid:18) πkN (cid:19) cos (cid:18) πkN (cid:19) −K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ sin (cid:0) πkN (cid:1)(cid:12)(cid:12) sin (cid:0) πkN (cid:1)(cid:12)(cid:12) + 4 K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) (cid:26) − (cid:18) πkN (cid:19) + 4 sin (cid:18) πkN (cid:19)(cid:27) +8 K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ sin (cid:18) πkN (cid:19) cos (cid:18) πkN (cid:19)(cid:21)(cid:19) ˙ φ + ℓ π N − X k =1 (cid:20) −K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19) +2 K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ sin (cid:18) πkN (cid:19) cos (cid:18) πkN (cid:19)(cid:21)(cid:19) ... φ . (A9)Reducing the number of terms we used: N − X k =1 f (cid:18)(cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) sin (cid:18) πkN (cid:19) = 0 , (A10) N − X k =1 f (cid:18)(cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) sin (cid:18) πkN (cid:19) = 0 , (A11) N − X k =1 f (cid:18)(cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) sin (cid:18) πkN (cid:19) = 0 . (A12) Thus, the equation for the position of the disk reads:3 N σℓ ¨ φ = − κN ℓ ˙ φ + N ℓ ˙ φ π − γ + log 2 ρ + N − X k =1 " − ℓ K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) sin (cid:0) πkN (cid:1) (cid:12)(cid:12) sin (cid:0) πkN (cid:1)(cid:12)(cid:12) + K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) cos (cid:18) πkN (cid:19) + N ℓ π − N − X k =1 (cid:20) K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ sin (cid:18) πkN (cid:19) − K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) cos (cid:18) πkN (cid:19)(cid:21)! ¨ φ + N ℓ π − ℓ − N − X k =1 " K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ sin (cid:18) πkN (cid:19) cos (cid:18) πkN (cid:19) − K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ sin (cid:0) πkN (cid:1)(cid:12)(cid:12) sin (cid:0) πkN (cid:1)(cid:12)(cid:12) + 4 K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) (cid:26) − (cid:18) πkN (cid:19) + 4 sin (cid:18) πkN (cid:19)(cid:27) +8 K (cid:18) ℓ (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:18) πkN (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:19) ℓ sin (cid:18) πkN (cid:19) cos (cid:18) πkN (cid:19)(cid:21)(cid:19) ˙ φ . (A13)Here we neglect the term proportional to ... φ . By dividingthe both sides of Eq. (A13) with N ℓ , we obtain Eqs. (19)to (22). W. F. Paxton, K. C. Kistler, C. C. Olmeda, A. Sen, S. K. St. An-gelo, Y. Cao, T. E. Mallouk, P. E. Lammert, and V. H. Crespi,
J. Am. Chem. Soc. , 13424 (2004). H. Ke, S. Ye, R. L. Carroll, and K. Showalter,
J. Phys. Chem.A , 5462 (2010). M. Nijemeisland, L. K. E. A. Abdelmohsen, W. T. S. Huck,D. A. Wilson, and J. C. M. van Hest,
ACS Central Sci. , 843(2016). H. Kitahata, R. Aihara, N. Magome, and K. Yoshikawa,
J. Chem.Phys. , 5666 (2002). T. Kusumi, K. Yoshikawa, S. Nakata, Oscillating Chemical Re-action in Oil/Water Systems: Generation of Macroscopic Os-cillatory Force. In Proceedings of the 3rd International Sympo-sium Far-From-Equilibrium Dynamics of Chemical Systems; edsJ. Gorecki, A. S. Cukrowski, A. L. Kawczynski, B. Nowakowski(World Scientific: River Edge, NJ, 1994) pp 87–93. J. Guzowski, K. Gizynski, J. Gorecki and P. Garstecki,
Lab Chip , 764 (2016). H. Kitahata, N. Yoshinaga, K. H. Nagai, and Y. Sumino,
Phys.Rev. E , 015101 (2011). H. Kitahata, N. Yoshinaga, K. H. Nagai, and Y. Sumino,
Chem.Lett. , 1052 (2012). S. Nakata, M. Nagayama, H. Kitahata, N. J. Suematsu, andT. Hasegawa,
Phys. Chem. Chem. Phys. , 10326 (2015). S. Nakata, Y. Iguchi, S. Ose, M. Kuboyama, T. Ishii,K. Yoshikawa,
Langmuir , 4454 (1997). Many chapters describing self-motion of camphor particles at dif-ferent conditions can be found in the recent e-book, S. Nakata,V. Pimienta, I. Lagzi, H. Kitahata and N. J. Suematsu, Eds. Self-organized motion: Physicochemical design based on nonlineardynamics. (Royal Society of Chemistry, Cambridge, UK, 2018). C. Tomlinson,
Proc. R. Soc. London , 575 (1862). L. Rayleigh,
Proc. R. Soc. London , , 364 (1889). S. Soh, K. J. M. Bishop, and B. A. Grzybowski,
J. Phys. Chem.B , 10848 (2008). N. J. Suematsu, T. Sasaki, S. Nakata, and H. Kitahata,
Langmuir , 8101 (2014). Y. Karasawa, S. Oshima, T. Nomoto, T. Toyota, and M. Fuji-nami,
Chem. Lett. , 1002 (2014). S. Nakata, M. I. Kohira and Y. Hayashima,
Chem. Phys. Lett. , 419 (2000). M. Shimokawa, M. Oho, K. Tokuda, and H. Kitahata,
Phys. Rev.E , 022606 (2018). R. Di Leonardo, L. Angelani, D. Dell’Arciprete, G. Ruocco,V. Iebba, S. Schippa, M. P. Conte, F. Mecarini, F. De Ange-lis, and E. Di Fabrizio,
Proc. Natl. Acad. Sci. USA , 9541(2010). M. Hayakawa, H. Onoe, K. H. Nagai, and M. Takinoue,
Sci. Rep. , 20793 (2016). M. Hayakawa, H. Onoe, K. H. Nagai, and M. Takinoue,
Micro-machines , 229 (2016). F. K¨ummel, B. ten Hagen, R. Wittkowski, I. Buttinoni, R. Eich-horn, G. Volpe, H. L¨owen, and C. Bechinger,
Phys. Rev. Lett. , 198302 (2013). T. Yamamoto, M. Kuroda, and M. Sano,
Europhys. Lett. ,46001 (2015). M. Frenkel, G. Whyman, E. Shulzinger, A. Starostin, and E. Bor-mashenko,
Appl. Phys. Lett. , 131604 (2017). T. M. Squires and M. Z. Bazant,
J. Fluid. Mech. , 65 (2006). T. Mitsumata, J. P. Gong, and Y. Osada,
Polym. Adv. Technol. , 136 (2001). V. Pimienta, M. Brost, N. Kovalchuk, S. Bresch, and O. Stein-bock,
Angew. Chem. Int. Ed. , 10728 (2011). V. Pimienta and C. Antonie,
Curr. Opinion Colloid InterfaceSci.
290 (2014). F. Takabatake, N. Magome, M. Ichikawa, and K. Yoshikawa,
J.Chem. Phys. , 114704 (2011). K. H. Nagai, F. Takabatake, Y. Sumino, H. Kitahata,M. Ichikawa and N. Yoshinaga,
Phys. Rev. E , 013009 (2013). F. Takabatake, K. Yoshikawa, and M. Ichikawa,
J. Chem. Phys. , 051103 (2014). N. Bassik, B. T. Abebe, and D. H. Gracias,
Langmuir , 12158(2008). H. Ebata and M. Sano,
Sci. Rep. , 8546 (2015). S.-I. Ei, H. Kitahata, Y. Koyano, and M. Nagayama,
Physica D , , 10 (2018). S. Nakata, K. Kayahara, H. Yamamoto, P. Skrobanska,J. Gorecki, A. Awazu, H. Nishimori, and H. Kitahata,
J. Phys.Chem. C , 3482 (2018). Y. Koyano, M. Gryciuk, P. Skrobanska, M. Malecki, Y. Sumino,H. Kitahata, and J. Gorecki,
Phys. Rev. E , 012609 (2017). K. Iida, H. Kitahata, and M. Nagayama,
Physica D , 39(2014). M. Nagayama, S. Nakata, Y. Doi, and Y. Hayashima,
Physica D , 151 (2004). H. Kitahata and K. Yoshikawa,
Physica D , 283 (2005). Y. Koyano, T. Sakurai, and H. Kitahata,
Phys. Rev. E , Y. Koyano, N. J. Suematsu, and H. Kitahata, arXiv:1806.04961(2018). W. S. Rasband, ImageJ, U. S. National Institutes of Health,Bethesda, Maryland, USA, https://imagej.nih.gov/ij/, (1997-2018). The frictional force working on the j -th camphor pill is repre-sented as − πρ κ ( d ℓ j /dt ). Therefore the torque originating from this force is − πρ κ ℓ j × ( d ℓ j /dt ) = − πρ κℓ ( dφ/dt ). By summingup the torque on N particles, we obtain Eq. (6). H. Kitahata and N. Yoshinaga,
J. Chem. Phys. , 134906(2018). G. N. Watson,