Quantitative Analysis of the Complex Time Evolution of a Camphor Boat

: The motion of a camphor boat on the water’s surface is a long-studied example of the direct transformation of chemical energy into a mechanical one. Recent experimental papers have reported a complex character of boat motion depending on the location of the camphor source. If the source is close to the stern, the boat moves at a constant speed. When it is shifted towards the boat center, oscillations of speed are observed. When the source is close to the boat center, pulses of speed followed by oscillations appear. Here, we focus on numerical simulations of camphor boat motion. We discuss approximations that allow us to reduce the numerical complexity of the problem and formulate a model in which the equation for boat velocity is coupled with a one-dimensional reaction–diffusion equation for camphor surface concentration. We scanned the phase space of model parameters and found the values that give qualitative agreement with the experiments. The model predicts all types of boat motion (continuous, oscillating, and pulsating) observed in experiments. Moreover, the model with selected parameter values shows that for specific locations of the camphor source, a spike in speed is followed by transient oscillations, which are an inherent part of speed relaxation.


Introduction
Self-motion of a camphor piece on the water surface is a spectacular example of a system in which the chemical energy of surface active molecules is transformed into kinetic energy of a floating object [1][2][3][4].The phenomenon of camphor "self-motion" has been attracting scientific attention for almost 200 years [3][4][5][6].In normal conditions, camphor evaporation is fast, and a piece of camphor located at a water surface can be regarded as an open system with a continuous flux of molecules from the source to the air through a surface layer.Camphor is amphiphilic, and the water surface tension is a decreasing function of the camphor surface concentration that generates Marangoni flows around the piece.An asymmetric profile of surface tension around a camphor piece generated, for example, by an inhomogeneous camphor concentration profile can produce a force on the order of 10 µN [7], which is sufficient to drive the piece.A camphor boat is made of a boat-shaped float with a piece of camphor located in the stern area.As a result of asymmetric camphor location, the directed motion appears [3,[7][8][9][10][11][12].It is hard to tell where and when the first camphor boats were made, but they were sold as toys in Edo-period Japan.Typical speeds of toy camphor boats are in the order of a few cm/s.
For a long time, researchers who investigated camphor boats were just satisfied with the observation of their motion.In the paper [13], the authors proposed a very simple mathematical model of the boat motion.A pioneering study on the complexity of the boat's motion was published in 2010 by Suematsu and collaborators [14], who recorded the boat speed as a function of time for different distances d between the location of a camphor piece driving the boat and the stern (cf. Figure 1).Their results can be summarized as follows: when the camphor piece was located close to the stern, a continuous motion with a constant speed was observed.If the piece was distant from the stern, then the boat motion was pulsating, and it was characterized by narrow intervals of high-speed motion separated by long time intervals when the boat did not move.The experiments reported in [14] were performed with boats freely moving inside a Petri dish (155 mm), which means that the boat collided with dish walls a few seconds after it was placed on the water surface.In our recent experiments [15], we considered a boat cut as a piece of U-shaped angle aluminum profile.The profile width was 4 mm, and the height of the profile was 2 mm.We used such a profile to localize camphor surface dissipation to the boat ends only.In experiments, we used boats of lengths 10 and 15 mm.The boats were propelled by a camphor cylinder or cuboid, which were made in a pill maker.The camphor source was located at the distance d from the stern, where d was the experiment's control parameter.An example of 10 mm long and 4 mm wide boat with a 4 mm × 4 mm square pill below is illustrated in Figure 1b.We used the boats that could rotate along a circle with the radius of 70 mm, whose center was set so that it met the center of the Petri dish (cf. Figure 1).The boat trajectory was one-dimensional, and we were able to evaluate the motion at stationary conditions for more than 1 h without collisions on the wall.Our experiments confirmed both types of previously reported motion: continuous when the camphor source was located at the boat edge and pulsating if the pill was close to the boat center.Moreover, at the intermediate distances between the pill and the boat edge, we observed a new type of motion characterized by an oscillating speed always larger than zero [15].Videos illustrating continuous, oscillating, and pulsating motion have been included as supplementary material to this paper.We used a one-dimensional model of boat motion introduced in [14] and adjusted its dimensionless parameters such that all three types of motion mentioned above were reproduced in numerical simulations.For the source close to the boat center, we observed in experiments a number of speed oscillations that follow after a peak.In the discussion presented in [15], we speculated that this type of behavior is related to the metastability of motion types for a range of d values.If so, then the relaxation of a speed pulse towards a standing boat occurs through metastable oscillations.However, we did not observe speed oscillations following a speed pulse in numerical simulations within the model considered in [15].This paper is concerned with modeling the complexity of camphor boat motion as a function of distance d.Unlike in [15], we focused our attention on a quantitative description of boat speed.At the beginning (Section 2), we describe the reduction of a two-dimensional model to a one-dimensional form, allowing for faster numerical simulations and discuss the approximations involved.In the next Section (Section 3), we consider the values of model parameters.The parameters introduced in [14] and adjusted in [15] allow for a rough qualitative description of the changes in the character of boat speed as the function of d.Here, we focus our attention on realistic values of parameters describing the boat and the dissipation of camphor to obtain results matching experimental measurements.We critically verify the values of parameters given in [7] and propose some modifications.The analysis of simulation results allows us to understand the mechanism leading to different types of boat motion.We obtain oscillations following the speed pulse for the optimal set of model parameters.In the final section, we discuss the obtained results and include suggestions for further study on the boat motion.

Reduction of the Two-Dimensional Reaction-Diffusion Model to the One-Dimensional One
In this section, we describe approximations that allow us to reduce the model of the boat motion to a set of differential equations in which the most complex part is a reactiondiffusion equation for a one-dimensional spatial variable.Table 1 lists symbols used and model parameters.
Rigorous numerical simulations of time evolution for a self-propelled object on the water surface are a hard numerical problem because the model should describe a few coupled phenomena.The driving force from the boat motion comes from the continuous flux of camphor molecules from the source through the surface layer to the environment.Camphor molecules decrease water surface tension depending on their local concentration.The local differences in surface tension generate Marangoni flows in the water phase and can produce forces and torques acting on the object and changing the position and orientation of the camphor source.The appearing flows influence camphor transport on the water surface.As a result, the description of object motion requires coupled equations for flows in all water volume (a vector function in three-dimensional space), the surface concentration of camphor represented by a scalar function of two spatial variables and a few scalar functions that describe object location and orientation [17,18].Of course, solving these equations is computationally complex.As a workable approximation to numerical simulations of the time evolution for a camphor-propelled object, the reaction-diffusion model was proposed in [8,10].Within this approach, the flow is not directly included in the model, and the surface transport of camphor is approximated by a reaction-diffusion equation with an "effective" diffusion coefficient [7,19].Therefore, the system is described by the surface concentration of camphor c(x, y, t) and all parameters that define the location of the considered object p i (t).Here, x and y are coordinates defining a point on the surface.
The evolution equation for c(x, y, t) has the form: The first term of the equation denotes the diffusive transport of the camphor molecules on the water surface [7,19].Here, D(x, y, t) is the effective diffusion coefficient of the molecules at the water surface selected and introduced to approximate the flow-induced transport.As we discuss later, D(x, y, t) can depend both on the position and time because the transport of camphor molecules below the object can be different from those on the free surface.Let us assume that the time-dependent shape of the object (in the following a camphor boat) is given by the function Ω B (x, y, t).The diffusion coefficient D can be represented by a two-valued function: where D 1 and D 2 are the diffusion coefficients below the boat and on the free water surface off the boat, respectively.We assume that the values of D 1 and D 2 (D 1 ≪ D 2 ) can differ by a few orders of magnitude [14].
The dissipation term combines camphor sublimation from the water surface with its dissolution in the bulk water.This rate of camphor dissipation from the surface at the point (x, y) is proportional to c(x, y, t) thus: The proportionality constant k can depend on the position and time because dissipation from the free water surface is different from that in the area covered by the object.It is common to assume that: where k 1 is the constant outflow rate to the water phase under the boat, and k 2 is the constant sublimation rate beyond the boat to the air phase.In the numerical experiments in this paper, we followed the assumption k 1 = 0 in accordance with the paper [14].
The source term describes the inflow of camphor and it also depends on position and time via the location of camphor source.Let us assume that the function Ω S (x, y, t) describes the location of the source as the function of time.Thus: COMMENT {camphor inflow from the source} = S(x, y, t) where c max is the maximum surface concentration of camphor at the specific conditions, and the factor (c max − c(x, y, t)) takes into account the decrease in camphor dissipation if its concentration at the source is close to c max .
Combining the terms discussed above, we obtain: where X Ω S (x, y, t) is the characteristic function of camphor source location.Equations describing the time evolution of parameters defining the location of the object (p i , including cartesian coordinates and angles) are of the Newton type and include forces and torques generated by gradients of the surface tension γ(x, y).They should also include friction resulting from the water phase.Formally, they can be formulated as: d 2 p i dt 2 ∼ {forces and torques depending on γ(x, y)} − {friction terms}.(7) The force acting on the object is generated by a difference in the surface tension around its boundary δΩ B (x, y, t) and can be calculated as: where ds is a normal outward vector element of the object boundary.
Similarly, a torque acting on the object is calculated with respect to an arbitrarily selected axis O (cf. Figure 2) as: where r = (x, y) is a vector pointing towards the surface element.Friction terms are standard and they can include translational and angular velocities and nonlinear combinations.Now let us formulate the evolution equations discussed above to the problem of the rotating boat illustrated in Figure 1.The boat considered in experiments has the form of a rectangle with the length l and width 2a, and it is marked with a black line in Figure 3.The camphor pill propelling the boat is located at a distance d from the stern, with d being the control parameter in our analysis.Let us assume that the center of the rectangular boat is at the point (x 0 , y 0 ).Therefore, the corners of the rectangle are at the points: (x 0 − l/2, y 0 − a), (x 0 − l/2, y 0 + a), (x 0 + l/2, y 0 + a) and (x 0 + l/2, y 0 − a).As shown below, the rotating boat's time evolution can be simplified if the shape is approximated as a rectangle in the polar coordinates.The center of the rectangle is at the point (r 0 , θ 0 ) and here The time evolution of the boat is described by a single variable, the time-dependent angle θ 0 (t), because the radius describing the location of the boat center r 0 remains unchanged.The evolution equation reads: The momentum I of the boat and arm complex can be approximated as: with the boat mass m 0 .The argument for such approximation is that the boat is narrow so the radius can be well approximated by r 0 and the arm of the boat is cut from thin plastic transparent foil and it is light if compared with the boat.The last term in Equation ( 10) describes the momentum of friction acting on the boat and it is approximated to be the function of the speed of the boat center.Such approximation is justified because the boat is narrow.In this term, µ 0 is the friction coefficient related to water viscosity and η 0 is the coefficient contributing to the quadratic dependence of friction on the velocity.Let us also note that the momentum of force acting on the boat edges can be approximated by the surface tension corresponding to camphor concentration averaged over the boat edge: If we introduce the new time-dependent variable z(θ, t) corresponding to camphor concentration averaged over the range of radii covering the boat width, then: Then, the Newton Equation (10) can be formulated as: For the value of I approximated according to Equation ( 11) and defining v(t) = r 0 dθ 0 /dt we obtain: To solve the above Newton equation, we need to couple it with the equation for the time evolution of z(θ, t).Such an equation can be derived directly from the reactiondiffusion equation in Equation ( 6) transformed into polar coordinates: Now, let us analyze the terms appearing on the right side of Equation ( 16).If the boat is approximated by a rectangle in polar coordinates, then the diffusion coefficient D(s, θ) for s ∈ [r 0 − a, r 0 + a] depends only on θ.
We can approximate all terms that appear on the right side of Equation ( 16).Let us note that the boat is narrow if compared to the arm length.Therefore, for the range of s integration [r 0 − a, r 0 + a] the values of s and 1/s are not much different from r 0 and 1/r 0 , respectively.Using this simplification, we obtain: The terms (∂c/∂s)| r 0 +a and (∂c/∂s)| r 0 −a are equal to 0 below the boat because, due to the U-shaped angle of the boat, there is no outflow of camphor in the radial direction.On the other hand, there is a camphor outflow through the dashed lines plotted in Figure 2. It can be expected that the outflow rate increases with an increase in the concentration of camphor between the dashed lines so, as the approximation, we can use: where L ≡ 0 below the boat and L = const.> 0 on the free surface.Within such approximation this term can be combined with the dissipation term discussed below.
The second term related to diffusion operator in Equation ( 16) can be approximated as follows: If we assume that the camphor dissipation is a two-value function given by Equation ( 4), then k(s, θ, t) does not depend on s in for s ∈ [r 0 − a, r 0 + a].As a result, we obtain: Finally, the source term includes information on the shape of the camphor piece used to propel the boat.If a piece is rectangular in polar coordinates and if it extends from one boat side to another (like a square pill illustrated in Figure 1b or Figure 3), then we can assume that S 0 (s, θ, t) does not depend on s for s ∈ [r 0 − a, r 0 + a].Therefore, Combining the terms listed above, we arrive at the one-dimensional reaction-diffusion equation for z(θ, t) approximating the two-dimensional reaction-diffusion equation for c(s, θ, t): This equation is solved for θ ∈ [0, 2π] as the spatial argument, or alternatively, we can use q = r 0 θ ∈ [0, 2πr 0 ] as the distance variable and transform all functions of θ to functions of q using θ = q/r 0 .Finally, the reaction-diffusion equation for z(q, t) has the form: As the following simplification step, we can introduce scaled surface concentration z(q, t) = z(q, t)/c max .Equation ( 23) is linear in z(q, t) and thus its structure remains unchanged if it is transformed into the equation for the time evolution of z(q, t): It is convenient to transform the evolution equations from the laboratory reference system to the moving coordinate system tied to the boat because the functions D(q, t), k(q, t), L(q, t) and S 0 (q, t) do not depend on time in the reference frame related to the moving boat.It can be carried out by introducing the space coordinate p such that: In such a variable, Equations ( 15) and ( 24) have the following form: and The set of equations Equations ( 26) and ( 27) was the basis of numerical simulations we discuss it the next Section.
The main difference between Equations ( 26) and ( 4) in [14] is the presence of the convection term v(t)∂ p z(p, t) that appears naturally when a diffusion equation is transformed to the moving coordinate system tied to the boat.This step is convenient computationally because the position of the boat is fixed and, consequently, the positions of the nodes in the corresponding mesh for Ω = {p ∈ R : 0 ≤ p ≤ Λ 0 } used to solve Equation (26) numerically are fixed with respect to the boat.Friction coefficient related to water viscosity η 0 Coefficient contributing to the quadratic dependence of friction on velocity L Outflow rate of camphor through the dashed lines (see Figure 2)

Results
In this section, we present numerical results.They were obtained considering the boat located close to the center of container Ω = {p ∈ R : 0 ≤ p ≤ Λ 0 } where Λ 0 = 750 mm.For the container boundaries, the non-flux condition for concentration was used: where t max denotes the simulation time.
In the majority of our calculations, the initial conditions corresponded to a standing boat and zero camphor concentration in Ω, which means z(p, t = 0) = 0, p ∈ Ω v(t = 0) = 0 for the reaction-diffusion-convection equation (Equation ( 26)), and v(t = 0) = 0 for the Newton equation (Equation ( 27)).The boat ends were at p 0 + l/2 = 382.5 mm and p 0 − l/2 = 367.5 mm.Equations were solved on spatial mesh ∆ p = 0.15 mm with the time step dt = 0.005 s.
For the numerical space discretization of Equation ( 26), we chose the finite volume method [20].We covered the region Ω with so-called control volumes, and then we integrated Equation (26) over each one.Thus, we received partial discretization with the values at the ends of the control volumes as auxiliary unknowns as well as integrals of the reaction and source terms.Then, for the reaction term ∂ ∂p D(p) ∂ z ∂p after integration, we further approximated D(p) ∂ z ∂p on the boundaries of the control volumes with finite differences in discrete values in the middle of them, which were the unknowns that we looked for.For the convection term −v(t) ∂ z ∂p , after integration, we replaced the values at the boundaries with the ones corresponding to the middle of the control volumes.For the reaction and source terms, we approximated integrals with values in the middle of the volumes multiplied by the volumes themselves.In the case of the control volumes of the same size (uniform mesh of the sought variables), this procedure is essentially equivalent to the finite difference scheme, but we decided to use it to manage the discontinuity of the diffusion coefficient D(p) at the ends of the camphor boat, which would be impossible for the approximation of the reaction term with finite differences, in which case, the term at the jump of D(p) generally approaches to infinity for the mesh size (control volumes) tending to zero.To cope with this obstacle, we used the fact that the term D(p) ∂ z ∂p is continuous at the jump (for details see [20]), and the numerical equivalent of the term at the jump involves the harmonic average of the values of D(p) on both sides of discontinuity.
The time discretization of Equation ( 26) was based on the implicit first-order finite difference scheme.The implicit scheme was chosen to have unconditional stability with respect to the time step size ∆t.In the case of Equation ( 27), the Euler (explicit first order) scheme was applied with the same time step.Throughout the calculations, the time step remained unchanged.Both camphor concentration z(p, t) and velocity v(t) at a given time t were calculated simultaneously using the values of these quantities calculated at t − ∆t.

Selection of Model Parameters
The boat motion is a complex phenomenon, and its model involves many parameters.The selection of the majority of their values can be determined in independent experiments.The droplet mass can be measured directly.The diffusion coefficient D 2 on the free water surface was estimated by observing tracer particles spreading out after a camphor disk was placed on the surface [7].The source term S 0 was obtained by measuring the mass of the camphor pill as the function of time.Boat friction can also be estimated in a direct experiment [7].We hope experimental information on the other quantities needed for the model will be available shortly.For example, the maximum surface concentration of camphor c max and the functional relationship between the surface camphor concentration and air-water interfacial tension and parameters describing it can be obtained from the quasi-elastic laser scattering technique [21].However, there are also parameters such as the ratio between surface diffusion below the boat and on the free water surface ∆ = D 1 /D 2 , the outflow rate of camphor to the water phase under the boat k 1 or the nonlinear contribution to friction η 0 , which seems to be hard to measure directly.Here, we applied the simplest approximation to η 0 and k 1 , assuming that both quantities were small and could be regarded as equal to zero.To observe boat motion, the value of ∆ has to be positive, and here, we fixed its value by the agreement between boat speed observed in experiments and in simulations.The same approach was used to determine γ(c).
The results of numerical simulations of boat evolution equations presented in the first paper, where the pulsating mode of boat motion was reported [14], were obtained only for the dimensionless parameters because they were focused on qualitative agreement with experimental results.In [14], as well as in the other papers cited below, the onedimensional model of the boat motion based on Equations ( 26) and ( 27) was used.The results confirmed continuous motion for small d and pulsating motion when d is large.The same model but with modified values of model parameters has also been used in our recent publication [15].After the modification of model parameters, the speed oscillations appeared for intermediate values of d.
A step towards the quantitative description of camphor boat motion was reported in [7] and followed by [12], where the values of some model parameters include: D 2 = 3940.0mm 2 s −1 , k 2 = 0.018 s −1 , µ 0 = 0.0175 g s −1 mm −1 and the rate of supply of camphor to the system S 0 • c max was experimentally measured.In [12], the camphor outflow from disk-shaped pill with 3 mm radius was approximated as 4 nmol s −1 .The boat mass m 0 = 0.064 g corresponds to the 15 mm boat used in our experiments.For the value of ∆ = 0.0001, we follow [14].These values were used in qualitative calculations [15] to match types of boat motion with experimental results.
We used two methods to determine the parameter values for a quantitative model.One of them aimed to transform dimensionless parameters used in [15] to the dimensional form.Such a procedure gives: In another approach, we used the experimental results of [7], the values of ∆ and η 0 mentioned above and adjusted the hard-to-measure parameter S 0 to obtain the best match with v(t) measured in experiments.Now, the parameter values are: The function γ : R + → R + appearing in Equation ( 27) relates the camphor concentration z to the surface tension.Following the previous papers [14,15], we use the following formula to define it: where γ 0 = 72 µN/mm is the surface tension of the pure water and γ 1 is the smallest surface tension possible induced by the presence of camphor.The function γ( z) is decreasing, which reflects the drop of surface tension with the increasing camphor concentration, and the parameters c 1 , c 2 , c 3 determine its shape.We used γ 1 = 0.1 × 72 µN/mm and γ 1 = 0.5 × 72 µN/mm for the parameter sets ( 28) and ( 29), respectively.The first value is suggested in [14] and appears as the ratio of dimensionless counterparts of γ 1 and γ 0 .The second seems to be in good agreement with experiments.The different values of the thresholds c 1 , c 2 and c 3 were used for the considered parameter sets: c 1 = 0.001, c 2 = 0.2005, c 3 = 0.4 for the parameter set (28) and c 1 = 0.015, c 2 = 0.1174, c 3 = 0.2198 for the case (29).By the definition of γ, the camphor concentrations below c 1 do not reduce surface tension, and this value is substantially greater for the experimental parameter set (29).Moreover, for the parameter set (29), the function γ is approximately two times narrower than in the case (28).Figure 4 shows the plot of γ( z) for the case corresponding to the experimentally adjusted parameter set (29).

Simulation Results
Let us start the discussion of camphor boat motion with the presentation of results for boat speed obtained for the selected set of model parameters.Simulations show that if we start with the initial condition described at the beginning of the section, then the time of relaxation towards the long-term behavior of v(t) increases with d.For d = 2.7 mm, the stable speed is observed after 200 s.In the case of complex pulsating motion, the transient evolution can be longer than 600 s.The results for the speed, the range of distances where a particular type of motion appears, and the periods of oscillatory and pulsating motion are presented in Table 2.It is clear that the set of model parameters adjusted for experiment (Equation (29); columns 2 and 3 in Table 2) give more realistic values of continuous motion speed, the correct range of d where the oscillatory motion is observed, and there is a better estimation of period P for d ∼ l/2 than the set of parameters in Equation (28) (cf.[15]).Therefore, in the following, we focus our attention on results obtained for the parameter set Equation (29).
The time-dependent speed of a boat for a few selected values of the distance d is illustrated in Figure 5. Here, the colors are coded as follows: blue, green, red, and magenta curves show v(t) for d = 2.7 mm, d = 3.0 mm, d = 5.4 mm and d = 5.7 mm, respectively.They illustrate v(t) for continuous, oscillatory low-and high-amplitude motion and for the pulsating mode.Next, we use concentration profiles from the selected intervals of time to explain the appearance of a particular motion type.
Table 2. Velocity v and period P of the continuous, oscillatory and pulsating motion of the camphor boat for various source distances d and parameters from Equation (29) (second and third columns) and from Equation (28) (fourth and fifth columns).For the continuous motion, we use P = 0.0; for the pulsating motion, the velocity range starts from zero; for the oscillatory motion, the lower limit of the velocity range is greater than zero.Before discussing the time evolution of camphor surface concentration, we compare the results of simulations for the parameter set (29) with experiments and numerical simulations for a dimensionless model (cf.Figures 7 and 8 in [15]).To carry this out, we introduce dimensionless variables d 0 , v and T 0 related to the distance measured in half-widths of the camphor pill w, speed in units of the boat speed when the pill located at the boat stern (v 0 = v(d∼ 0)) and time in units defined by w/v 0 [15].The graphical presentation of these results is given in Figure 6.The comparison with experimental results (cf. Figure 7 in [15]) shows much better agreement than the qualitative model presented in the above-mentioned paper.In particular, in the presented simulation results, we do not see a rapid increase in the interspike time when the camphor source approaches the boat center.Moreover, the model with parameters (29) predicts that the period of oscillations following a speed spike decreases with d what was observed in experiments (cf.Figures 6b and 7b in [15]).
Figures 7-10 illustrate the time evolution of the camphor concentration around the boat.In all figures, the flat, green rectangle shows the boat's position.A tall yellow bar indicates the location of the camphor pill.In all cases, the camphor was closer to the right end of the boat; thus, the convection in the reference frame related to the boat is always directed to the right.
It can be noticed that for the selected model parameters, the camphor concentration at the boat bow (left end) is always small, which can be explained by a small diffusion below the boat and a large diffusion outside it.Therefore, the value of concentration at the boat stern (right boat end) decides if the boat is propelled.If it is larger than c 1 , the force propelling the boat appears.Another factor that influences the boat's motion is the convection.A slow diffusion below the boat helps to accumulate camphor molecules close to the camphor pill.The rate at which they are transported to the stern region is decided by the boat velocity that defines convection in the considered reference frame.Let us start the discussion of the camphor concentration profile with the case when the boat moves at a constant speed.The stable concentration profile for a pill located at d = 2.7 mm from the stern is illustrated in Figure 7.The stationary speed of such a boat is reached in simulations after 200 s starting with the initial condition defined at the beginning of the section, and it is equal to 8.96 mm/s (cf. Figure 5a).The camphor concentration profile within the pill is determined by convection that moves camphor to the area between the pill and the stern.The drop in concentration at the right boat end comes from a high diffusion outside the boat.The value of camphor concentration at the stern decreases the surface tension by ∼0.1 µN/mm, and such a difference between it and the surface tension at the bow is high enough to support the continuous motion with the speed given above.The surface tension around continuously moving camphor boats has been recently investigated by light scattering technique [21].Our relationship between the surface tension difference and the boat speed agrees with these measurements.Figure 8 illustrates the time evolution of the camphor concentration profile for lowamplitude oscillatory motion (cf. Figure 5b).The dots in the inserts illustrate the times for which the results are presented.Concentration profiles for the accelerating boat are shown in Figure 8a, and those for the decelerating boat (right branch) are in Figure 8b.In all cases, the changes in concentration outside the boat are marginal.Moreover, the concentration profile in the pill area seems to be mainly related to convection and changes a little with boat speed.What matters is the concentration profile below the boat in the area between the pill and the stern.We observe an increase in camphor concentration when the boat is slow (blue vs. black curve in Figure 8a), a slight decrease at high speed (black vs. red curve in Figure 8a), a decrease in concentration during the fast motion to the area where there is no camphor (the black curve in Figure 8b) and, finally, an increase in concentration for a slow boat (the red curve in Figure 8b).Therefore, the instability of camphor concentration for speeds exceeding 10 mm/s is responsible for oscillatory boat motion.7.
For the high-amplitude oscillatory motion (cf. Figure 5c), the time evolution of camphor concentration is more complex, as shown in Figure 9. Unlike the low-amplitude oscillatory motion, the maxima of concentration observed for accelerating and decelerating speeds differ a few times.When the boat is slow (the blue curve in Figure 9a), camphor accumulates in the pill area.The effect is similar to that shown in Figure 8a, but due to slower speed, the maximum concentration is twice as high.This translates into a larger accelerating force and high speed at the peak.During the fast motion, the boat moves to the area where the camphor concentration is low (the red curve in Figure 9b); the propelling force is small, and the boat is slowed down by friction.When the boat moves at a low speed, camphor concentration increases (the black curve in Figure 9b) and the cycle repeats.We believe that the key factor allowing high-amplitude oscillations of speed is the area between the pill and the stern.If it is large, then more camphor can accumulate there and produce a stronger speed spike.If the distance traveled with a large speed is long, then the boat moves to the area with a low camphor concentration and stops there.As a result, the pulsating motion appears.7.
The time evolution of camphor concentration before and during the initial peak of speed in a pulsating motion is illustrated in Figure 10. Figure 10a illustrates the concentration when the boat slows down (the blue curve), when the boat does not move (the black curve), and when it starts to move again (the red curve).The scenario is similar to that describing slow motion during the large-amplitude oscillations shown in Figure 9b.During the spike of speed, the boat moves to an area where camphor concentration is low, and the motion is that fast if compared to camphor release that it has not cumulated below the boat.The blue curve in Figure 10a illustrates concentration for a slowing boat.As seen, it is below c 1 , which defines the threshold for non-zero propelling force acting on the boat.If there is no force propelling the boat, the boat stops, and camphor concentration increases symmetrically around the pill, as shown by the black curve.Asymmetry of concentration develops as the result of non-symmetrical pill location with respect to the bow and stern.It takes over 50 s to develop the force that can move the boat.During that time, a high concentration of camphor cumulates below the boat, and it is shifted towards the stern when the boat starts to move (cf.the blue curve in Figure 10b).As a result, the driving force increases and leads to fast boat acceleration.And again, after a speed spike, the concentration of camphor below the boat drops (cf. the red curve in Figure 10b and the blue curve in Figure 10c); thus, the boat slows down due to the friction.When the speed is low, the process of camphor cumulation below the boat starts again, as indicated by the red curve in Figure 10c.

Discussion
In this paper, we formulated a model of low numerical complexity that describes the motion of a boat propelled by camphor dissipation from an asymmetrically located source.The model combines the Newton equation for boat speed with a one-dimensional reaction-diffusion equation for surface camphor concentration.It characterizes the time evolution of the boat by its velocity v(t) and the profile of camphor concentration along the boat's trajectory.The model was applied to successfully explain the complexity of motion types for boats with different values of d-the distance between the position of the camphor source and the boat stern.The experiments with camphor boats of the considered geometry indicate continuous motion for small d, speed oscillations for their intermediate value, and pulsating motion d corresponding to camphor sources located near the boat center.These observations suggest a general tendency of changes in motion type: a steady supply of camphor produces continuous motion; when it is reduced, oscillations appear, and for a small supply, one can observe pulsating motion.Other physical factors can also modify camphor release, diffusion, and dissipation.One of them is the water temperature [12], which is expected to speed up all these processes.As may be anticipated, the boat motion changes from pulsating to continuous at a critical temperature.The experiments reported in [12] were conducted for a free-moving boat.Thus, the motion was affected by collisions with the dish walls.We believe an experiment performed in carefully controlled stationary conditions will also reveal boat oscillations.Actually, a trace of them with a period of 2 s can be seen in Figure 2(c-2) of [12].Time is another factor that can modify the character of motion.Camphor slowly dissolves in water, and the increase in dissolved camphor concentration reduces the gradients of interfacial tension.There are experiments showing the time transition between a continuous motion to a pulsating motion for a droplet of paraffin solution of camphor located on the water surface [22,23].It is interesting to note that similar motion types were also observed during the time evolution of a gel float containing camphoric acid as the surface active agent [24], and we expect the presented model can be also used to explain this phenomenon.
Using the experimental values of some model parameters and adjusting those that are hard to measure directly, we were able to reproduce the time-dependent speed observed for boats characterized by different values of d.Our results go beyond the qualitative agreement with experimental results because they qualitatively characterize velocity profiles.In the previously published experimental paper [15], we speculated that oscillations following a peak of speed can be explained by the coexistence of two metastable states of the boat motion: one corresponding to separated spikes of speed and another representing oscillatory motion.The presented simulation results suggest that for some values of d, speed oscillations are an inherent part of spike relaxation and do not prove the existence of a separated metastable oscillatory mode.
In parallel to v(t), the model predicts the time-dependent camphor concentration profile.This important information has not attracted much scientific attention so far because the non-invasive measurements of surface tension were not popular.The recent studies ( [21] and references therein) have changed it.In this respect, our model, which provides information about the motion and state of the surface around the object, is important for the verification of theoretical results and for a better understanding of phenomena observed in experiments.Moreover, as demonstrated in Section 3, the model allows us to analyze the origin of transitions between different motion modes.For example, the transition between the continuous motion and small amplitude oscillations can be explained by the instability of the camphor concentration profile in the region between the edge of the pill and the boat stern.Therefore, the presented model can be useful for the future detailed study of bifurcations between the motion types of self-propelled objects.
The stability of numerical results concerning parameter values depends on the motion mode.The continuous motion and the pulsating motion are the most stable (and the same comment applies to experiments).They can be observed for an extensive range of parameters used.The oscillation window is narrow in the parameter space.Complex pulsations (cf. Figure 5d) are the most unstable and therefore most challenging to reproduce.Among the parameters, the numerical results are the most sensitive with respect to S 0 , µ 0 , and ∆, whereas they are pretty stable for D and k 2 .
It is interesting to point out mathematical similarities between camphor boat and the problem of object motion on a spatially distributed viscoelastic foundation [25].From the applied mechanics viewpoint, the stability of motion is essential, whereas nonlinear chemistry is more interested in locating bifurcations in motion type in the parameter space.Nevertheless, mathematically, the problems seem compatible, and the methods for investigating system stability described in the abovementioned reference can be adaptable to the boat motion problem.

Figure 1 .
Figure 1.(a) An example of the setup for our experiments with a rotating camphor boat.The diameter of the Petri dish was 20 cm.The camphor boat was located at the end of a 7 cm-long arm and can rotate around the vertical axis located at the dish center.A typical water level was 5 mm.(b) The boat was 10 mm long and 4 mm wide, as seen below.It was cut from an aluminum U-shaped angle profile.The mass of the cuboid camphor pill (4 mm × 4 mm × 1 mm) was 0.008 g.The mass of the aluminum 15 mm boat profile was 0.056 g.The mass of the arm was 0.020 g, but we neglected it in calculations of inertia.The square marks a pill of camphor at the distance d from the stern [16].

Figure 2 .Figure 3 .
Figure 2.A section of the Petri dish (the black arc) illustrating the geometry of the experimental system.The camphor boat (blue rectangle, cf.Figure3) is located at the end of an arm (green line, radius r 0 = 70 mm) and can rotate around the vertical axis at point O.The trajectory of the boat center is illustrated with a red arc.The boat width is 2a = 4 mm.The dashed blue rings show trajectories of two points located below the arm at the opposite edges of the boat.

Figure 5 .
Figure 5. Boat speed of as the function of time for different values of the distance d illustrating changes in the character of motion from stationary to complex pulsations.Sub figures (a), (b), (c) and (d) correspond to d = 2.7 mm, d = 3.0 mm, d = 5.4 mm and d = 5.7 mm, respectively.

Figure 6 .Fi gure 7 .
Figure 6.(a) The range of boat speeds observed for different modes of motion as a function of the distance d.Red disks, blue squares, and green circles correspond to continuous, oscillating, and pulsating motion.(b) The period of oscillating motion and the mean time between speed bursts in the pulsating mode as a function of d.

Table 1 .
Symbols used in the paper and their explanation.
= 402.195s (v(t 2 ) ∼ 10.1 mm/s) and red-t 3 = 402.545s (v(t 3 ) ∼ 8.3 mm/s).The meaning of the bar and rectangle is the same as in Figure Concentration of camphor as a function of the space variable for the boat with d = 5.4 mm when oscillations characterized by a high amplitude are observed (cf.Figure5c).The color of the plots corresponds to the points plotted in the inset, describing the time dependence of the boat speed: (a) blue-t 1 = 703.36s (v(t 1 ) ∼ 2.7 mm/s), black-t 2 = 704.37s (v(t 2 ) ∼ 24.7 mm/s) and red-t 3 = 704.42s (v(t 3 ) ∼ 35.7 mm/s), (b) blue-t 1 = 704.5 s (v(t 1 ) ∼ 45.9 mm/s), black-t 2 = 704.76s (v(t 2 ) ∼ 35.1 mm/s) and red-t 3 = 705.63s (v(t 3 ) ∼ 13.5 mm/s).The meaning of the bar and rectangle is the same as in Figure Concentration of camphor as a function of the space variable for the boat with d = 5.7 mm when oscillations characterized by a small amplitude are observed (cf.Figure5d).The color of the plots corresponds to the points plotted in the inset, describing the time dependence of the boat speed: (a) blue-t 1 = 626 s (v(t 1 ) ∼ 1.09 mm/s), black-t 2 = 671.73s(v(t 2 ) ∼ 0.0 mm/s) and red-t 3 = 683.5 s (v(t 3 ) ∼ 1.12 mm/s), (b) blue-t 1 = 683.78s(v(t 1 ) ∼ 159 mm/s) and red-t 2 = 683.84s(v(t 2 ) ∼ 250 mm/s), (c) blue-t 1 = 684.31s(v(t 1 ) ∼ 149 mm/s) and red-t 2 = 688.3s (v(t 2 ) ∼ 3.0 mm/s).The meaning of the bar and rectangle is the same as in Figure7.