On the reaction–diffusion type modelling of the self-propelled object motion

In this study, we propose a mathematical model of self-propelled objects based on the Allen–Cahn type phase-field equation. We combine it with the equation for the concentration of surfactant used in previous studies to construct a model that can handle self-propelled object motion with shape change. A distinctive feature of our mathematical model is that it can represent both deformable self-propelled objects, such as droplets, and solid objects, such as camphor disks, by controlling a single parameter. Furthermore, we demonstrate that, by taking the singular limit, this phase-field based model can be reduced to a free boundary model, which is equivalent to the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2$$\end{document}L2-gradient flow model of self-propelled objects derived by the variational principle from the interfacial energy, which gives a physical interpretation to the phase-field model.


On the reaction-diffusion type modelling of the self-propelled object motion
Masaharu Nagayama 1,7* , Harunori Monobe 2,7 , Koya Sakakibara 3,4,7 , Ken-Ichi Nakamura 5 , Yasuaki Kobayashi 1 & Hiroyuki Kitahata 6 In this study, we propose a mathematical model of self-propelled objects based on the Allen-Cahn type phase-field equation.We combine it with the equation for the concentration of surfactant used in previous studies to construct a model that can handle self-propelled object motion with shape change.A distinctive feature of our mathematical model is that it can represent both deformable selfpropelled objects, such as droplets, and solid objects, such as camphor disks, by controlling a single parameter.Furthermore, we demonstrate that, by taking the singular limit, this phase-field based model can be reduced to a free boundary model, which is equivalent to the L 2 -gradient flow model of self-propelled objects derived by the variational principle from the interfacial energy, which gives a physical interpretation to the phase-field model.
A self-propelled system is composed of a particle or droplet that exhibits self-propelled motion by consuming the free energy under non-equilibrium conditions.Such a particle or droplet can move by the internal mechanism; some can move by deforming themselves, and others move by changing the characters in the neighbouring field.Motions of living organisms, such as birds, insects, bacteria, and cells [1][2][3][4] and those in non-living materials, such as camphor disks, swimming droplets, running droplets, and Janus particles [5][6][7][8][9] , are regarded as the motion in self-propelled systems.In recent years, various types of cell motility have been analysed through mathematical models, including keratocyte motility 10 , cell population motility [11][12][13] , and cell division 14 .Mathematical models for the motility of non-living materials are also investigated by constructing the mathematical model, for example, camphor particles 15,16 , pentanol droplets 17 , running oil droplets 18 , and blebbing oil droplets 19 .In this study, we focus on the spatio-temporal behaviours in self-propelled systems, such as the motion of camphor disks and pentanol droplets, where the mechanism for the self-propelled motion is believed to be primarily governed by surface tension gradient.Such surface-tension-driven systems are classified into two cases; the systems without shape deformation, like camphor disks, and those with shape deformation, like pentanol droplets.Formulations of two-dimensional motion models for the former case have already been made, and some comparisons of experimental results with those of mathematical models and mathematical analyses have been performed 20,21 .
On the other hand, as for mathematical modelling for the self-propelled system with shape deformation, the mathematical modelling approach has not yet been satisfactorily performed.A mathematical model based on the membrane motion model has been proposed for the droplet motion; however, it is not easy to handle as a model of droplet dynamics due to its substantial computational cost 22 .The other approach for the droplet motion is to adopt the lubrication approximation.Thiele and coauthors reported the two-dimensional model for the self-propelled droplet 23 , but no systematic study has not been demonstrated.Therefore, we aimed to construct a mathematical model suitable for numerical computation at a low cost that can be connected to the natural model from the physics viewpoint.For the purpose above, we adopted the Allen-Cahn equation 24 to describe the droplet shape.By combining it with the equation for surfactant concentration used in previous studies, we constructed a mathematical model that can handle self-propelled object motion with shape change.The formulation of this mathematical model is a modelling method called the phase-field method, which has been used for crystal growth models in supercooled liquids 25 and crystal interface motion models 26 .In recent years, there have been models 10,13,14 to represent cell motility in living systems using the Allen-Cahn and the Cahn-Hilliard 27 equations.They are also used in many other fields of materials science 28 .Although the phasefield model has been widely used and has the advantage that the cost for numerical computation is low, the model also has a disadvantage in that the correspondence between the actual physical quantity and the order parameter ϕ remains unclear, which is because the order parameter is introduced artificially to connect the regions with different phases smoothly.Therefore, in this paper, we first construct the reaction-diffusion model for a self-propelled droplet using the volume-preserving Allen-Cahn type phase-field equation.Then, we derive the singular-limit model from the reaction-diffusion system model and confirm that the singular-limit model matches the self-propelled object motion model derived as the L 2 -gradient flow based on the variational principle.From these results, we can demonstrate that the reaction-diffusion system model is an ε-approximation of the L 2 -gradient flow model, which supports the physical meaning of the constructed mathematical model based on the Allen-Cahn type phase-field equation.
This study proposes the following mathematical model for a self-propelled object using the volume-preserving Allen-Cahn equation: where and Here, ϕ(x, t) represents the position and shape of the self-propelled droplet by defining 1/2 ≤ ϕ(x, t) ≤ 1 as the inside of the droplet and 0 ≤ ϕ(x, t) < 1/2 as the water surface.u(x, t) represents the surface density of sur- factant molecules on the water surface.d u , k 1 , k 2 , and k 3 are the diffusion coefficient, dissolution rate, sublimation rate, and supply rate of surfactant molecules, respectively; s 0 is the density of self-propelled surfactant molecules inside the droplet; τ , σ , and ε are positive constants such that 0 < ε ≪ 1 .is a bounded region in R 2 .The first equation for the time evolution of ϕ represents the time evolution of the droplet shape.Since this equation is obtained based on the Allen-Cahn equation, the profile ϕ should have the regions with ϕ ≃ 0 and ϕ ≃ 1 and the smoothly connecting boundary region, which correspond to the inside of the droplet, water surface and the periphery of the droplet, respectively.Function a in Eq. ( 2) represents the driving force by the surface tension at the object's periphery and also keeps the object's volume almost constant.
The equation for a is constructed so that greater γ (u) and the smaller volume decrease a, which drives the droplet periphery outward.As for the time-evolution equation for u, the first term of the righthand side shows the diffusion of the surfactant molecules at the water surface.Since the surface tension gradient induces the Marangoni convection 29 , the transport of the surfactant molecules may be partly caused by the convection.However, this transport can be incorporated into the effective diffusion coefficient, which describes the overall surfactant molecule transformation 30 .The second term corresponds to the decrease in concentration by dissolution to the aqueous bulk phase and sublimation to the gas phase.We assume that the dissolution and sublimation occur proportionally to the surface concentration u.The last term shows the supply of surfactant molecules from the object.Here we assume that the supply is proportional to ϕ since the supply only occurs in the region where the object covers.Such a quasi-conservative reaction-diffusion system was first proposed by Krischer et al. 31 .Previous results for the volume-preserving Allen-Cahn equation proved the existence and uniqueness of the classical solution and the boundedness of the solution 32 .Moreover, the mean curvature flow is derived using the singular limit 33 .Furthermore, the model (1) mentioned above, except the second equation, is also considered in terms of the geometric measure theory 34,35 , which shows that the singular limit becomes the volume-preserving mean curvature flow.
First, we present a mathematical model for a self-propelled object using the volume-preserving Allen-Cahn type phase-field equation and demonstrate by numerical simulations that the model can describe the motion of self-propelled objects with varying deformability.Then, we present a free boundary problem in which a closed curve represents the self-propelled object by performing a singular perturbation expansion.Numerical computations for this model show how the self-propelled motion can be reproduced.We also confirm numerically that the parameter σ can control the interface change.Although free boundary problems obtained from the singular limit have been reported for reaction-diffusion systems with two variables [36][37][38] , ones with integral terms, such as Eq.(1), have not been known until now.We then define the interfacial energy of the water surface, the length energy of the self-propelled object, and the area-preserving energy.Based on the variational principle, we derive an L 2 -gradient flow model representing the self-propelled object as a closed curve, showing that this mathemati- cal model is consistent with the free boundary problem.We observe that the mathematical model using the volume-preserving Allen-Cahn equation approximates the L 2 -gradient flow model, which corresponds to the free energy form.Finally, we summarise this study and discuss future work. (1) tion, we propose a phase-field type reaction-diffusion model that describes the motion of a self-propelled object.
Here, the dimensionless model is shown, whereas the non-dimensionalisation from the dimensional model is shown in "Methods".In our model, we set the order parameter ϕ(x, t) , which defines 1/2 ≤ ϕ(x, t) ≤ 1 as the self-propelled object and 0 ≤ ϕ(x, t) < 1/2 as the water surface.u(x, t) is the concentration of surfactant mol- ecules supplied by the self-propelled object.Then, we obtain the following non-dimensionalised mathematical model from ( 1)-(2) (see "Methods"): where The initial values for the mathematical model ( 3) are defined as follows: where ∂� represents the boundary of , ϕ 0 is a function with compact support, for example, given as Numerical computation by the reaction-diffusion model.Numerical computations are performed on the reaction-diffusion model ( 3)-( 4) under the initial condition (5) to investigate the properties of the mathematical model.Since the width of the interface of ϕ is O(εσ ) , ε and σ are given so that εσ = 0.025.First, consider the case σ = 0.05 (i.e., ε = 0.5 ) with τ and the integral S 0 of the initial function as free param- eters.When S 0 = 1 , a disk-shaped standing spot appears stably, as shown in Fig. 1a for large τ .As τ gradually decreases, the disk-shaped standing spot becomes unstable, and a dumbbell-shaped standing spot appears, as shown in Fig. 1b.As τ is further decreased, a banana-shaped travelling spot solution appears (Fig. 1c).Then, as τ is decreased, the velocity increases, and the banana-shaped travelling spot deforms into a rice-ball-shaped travelling spot (Fig. 1d).When S 0 = 0.5 , as τ gradually decreases, a travelling spot appears, which is almost disk- shaped, from the disk-shaped standing spot (Fig. 1e).The above result shows that a travelling spot bifurcates from the disk-shaped standing spot supercritically.
As τ is further decreased, a faster rice-ball-shaped travelling spot appears.Next, we consider the case σ = 1 (i.e., ε = 0.025 ).For S 0 = 1 , a disk-shaped standing spot appears for large τ , but as τ is gradually decreased, the disk-shaped standing spot becomes unstable, and a travelling spot close to the disk appears (Fig. 1f).Furthermore, even when S 0 = 2 , although the travelling velocity is high, the deformation is slightly elliptical, but the convexity is maintained (Fig. 1g).To confirm that the disk shape is maintained for larger values of σ , we set σ = 5 (i.e., ε = 0.005 ), and indeed, a travelling spot close to a disk appears for small values of τ .The results show that an almost disk-shaped travelling spot emerges even at small τ .Furthermore, these spot solutions are found to asymp- tote to a constant velocity from a suitable compact support function ϕ 0 (x) , as shown in Supplemental Fig. 1a-g.
(3) The banana-shaped travelling spot, as shown in Fig. 1c, is similar to a pentanol droplet exhibiting the translational motion (Fig. 1(c-1) in the Reference 17 , Fig. 1 in the Reference 39 ); the elliptic travelling spot, as shown in Fig. 1d, appears in the translational motion of a relatively small pentanol droplet (Fig. 2a in the Reference 22 ); and the travelling spot close to a disk, as shown in Fig. 1f, appears in the translational motion of an ethyl salicylate droplet (Fig. 1 in the Reference 40) .Moreover, the travelling spot near the disk shape with high velocity, as shown in Fig. 1g, approximates the non-deformable self-propelled objects such as solid camphor disks.These results suggest that the mathematical model ( 3)-( 4) can describe a wide range of motions, from deformable self-propelled objects such as droplets to non-deformable self-propelled objects such as solid camphor disks, using σ as a parameter.Reduction to the singular-limit model.In this section, we clarify that the shape of the spot pattern in (3) can be controlled by using only the parameter σ .To this end, we derive a singular limit model as ε tends to 0 in (3) and perform numerical simulations for the singular limit model.
The first equation of ( 3) is equivalent to the "Allen-Cahn equation", provided that a(S[ϕ](t), u, ε) ≡ 1/2 .It is well-known that, as ε tends to 0 under this setting, the motion of the level set at ϕ(x, t) = 1/2 is close to that of the mean curvature flow.This section investigates how the solution of (3) and the level set at ϕ(x, t) = 1/2 behave as ε tends to 0. Let Ŵ(t) be a Jordan curve depending on the time t in R 2 , that is, and � in (t) be a bounded domain with a smooth boundary Ŵ(t) .Then, the singular limit of (3) leads to the fol- lowing free boundary problem, composed of the interface equation and the reaction-diffusion equation: where The derivation of ( 6) is shown in "Methods".Numerical computation by the singular limit model.In this section, we perform simulations of the singular limit Eq. ( 6) and show that the parameter σ controls the shape of travelling spots that appear in (6).
To begin with, let σ = 0.05 .We here regard τ and S 0 defined by the initial function as free parameters.When S 0 = 1 , if τ is sufficiently large, the disk-shaped stationary solution, which is stable, appears (Fig. 2a).As τ gradu- ally decreases, as in the reaction-diffusion model, a dumbbell-shaped standing spot solution first appears (see Fig. 2b), followed by a banana-shaped travelling spot (see Fig. 2c).As τ is further reduced, a rice-ball-shaped travelling spot eventually appears (see Fig. 2d).On the other hand, as S 0 = 0.5 , the situation is slightly different.When τ decreases gradually, the travelling spot appears as well as S 0 = 1 , but the banana-shaped one does not appear.Only travelling spots, whose shapes are almost disk-shaped, are observed (Fig. 2e).Next, let σ = 1.0 .Under the condition of S 0 = 1 , as τ decreases gradually, the almost disk-shaped travelling spot appears (Fig. 2f).As S 0 = 2 , the shape of the spot changes from a disk to an ellipse and preserves the convexity (Fig. 2g).Similar to the model ( 3)-(4), the singular limit equation ( 6) also asymptotes to a constant velocity and a constant shape spot from suitable initial values, as shown in Supplemental Figs.2a-g and 3I,II.In this section, we con- struct a gradient flow model of a self-propelled motion to give a physical meaning to our proposed model (3).Let L(∂� in ) , E(� in ) and S(� in ) be the surface energy of water surface, the perimeter energy of self-propelled particles and the energy of area-preserving, respectively, which are given by where s = s(x, t) stands for the arc-length of x ∈ ∂� in , ∂� in = ∂� in (t) , � in = � in (t) and Define the total energy H( �in ) by Since we obtain from (19) (see "Methods") that assuming that the energy of ∂� in (t) decreases in time, we conclude that Combining the interface motion (10) with the reaction-diffusion model for surfactants and choosing we obtain the following L 2 -gradient flow model: which is the same as (6).The detailed calculation is shown in "Methods".This result shows that our proposed model ( 3) is an ε-approximate model of the gradient flow derived from the energy variation.The equation of motion of the interface is derived from the variational principle, where time evolution occurs in the direction of decreasing total energy.However, the variational principle is not applied to the concentration field.Note, therefore, the energy of the entire system is not conserved.Furthermore, the system is not stationary at the energy minimum.

Conclusion
In this study, we have demonstrated that a phase-field based reaction-diffusion model of self-propelled motion can be considered as an approximation of an energy-based gradient flow model.We have shown that a free boundary model ( 6) can be derived as a singular limit from the reaction-diffusion model (3) and that this derived model is exactly the same as the model ( 11) derived from the variational principle applied to the interfacial energy of the self-propelled object.This result gives a phase-field based model, typically regarded as a phenomenological model, a physical interpretation tied to the interfacial energy.Furthermore, we have demonstrated that, by adjusting a single control parameter σ , our reaction-diffusion model yields both deformed and near-circular traveling spots, which means that our model can represent both deformable droplet motion and solid camphor disk motion.Thus, our work offers a more comprehensive and theoretically grounded framework for the modeling of self-propelled objects.
In our model, the bifurcation from a standing to travelling spot occurs when the parameter τ is decreased.Intuitively, if a self-propelled object has a circular shape, the surface tension should be balanced due to the symmetric diffusion of surfactant from the object.However, when the object moves slightly in a certain direction, the concentration field behind it increases, and the surface tension behind it decreases, causing an imbalance of forces.When the imbalance overwhelms the effect of diffusion attempting to restore symmetry, it continues to move (see the Supplementary movies).For deformable objects such as droplets, the motion can also break the symmetry of the shape; for solid objects, like camphor particles, only the concentration field breaks the symmetry.The supplementary movies shows that, even if the self-propelled object maintains a symmetrical shape, the motion continues as the concentration field profile becomes asymmetrical.This illustrates how our model can exhibit both deformed and circular travelling spots.
Our current model cannot handle arbitrarily shaped solid objects such as an elliptic-shaped camphor particle, which requires that the steady state of the phase-field model be controlled in a desired way.Whether it is possible to represent self-propelled object motions that maintain elliptical or other shapes within our framework would be an interesting future work.We would then like to compare the results of the previous analysis 20,21 of elliptical camphor motion with the current model.www.nature.com/scientificreports/Moreover, the stable dumbbell-shaped standing spot bifurcated from the stable disk-shaped standing spot has been found by numerical calculation in the self-propelled motion of reaction-diffusion systems.Figure 1b shows that the stable dumbbell-shaped standing spot has been reported for the reaction-diffusion equation with the spatially inhomogeneous nonlinear term 41 .However, only the unstable dumbbell-shaped standing spot has been reported for spatially homogeneous reaction-diffusion systems 42,43 .In the future, we will analyse the existence and stability of stationary solutions to the reaction-diffusion type self-propelled object motion model ( 3)-( 4) and the free boundary model (6) and investigate the dynamics of the mathematical model in detail by analysing the bifurcation structure to the travelling spot.In the case of splitting phenomena such as observed in droplets, the volume-preserving reaction-diffusion model ( 3)-( 4) can be used to confirm the splitting of self-propelled objects numerically.However, this model does not conserve the volume of individual self-propelled objects.To describe the motion of multiple self-propelled objects after splitting, it is necessary to extend the mathematical model to handle the phenomenon of splitting, for example, by preparing a volume-preserving model for the number of self-propelled objects as a simple extension.

Methods
and the periodic boundary condition is imposed.The parameters α and u 1 are set to α = 1000 , u 1 = 0.005 for all settings.Numerical computations of the reaction-diffusion model (3)-(4) were performed using the alternating direction implicit method 44 .
Formal derivation of the singular-limit model.In what follows, we demonstrate the interface equation by the formal argument.Outer expansion.Assume that ϕ and u have the following expansions away from the interface Ŵ(t) : Substituting the expansion ( 12) into (3) and collecting the term ε 0 , we know that ϕ 0 satisfies ϕ 0 (1 − ϕ 0 )(ϕ 0 − 1/2) = 0 , and hence ϕ 0 = 1 in � in (t) and ϕ 0 = 0 in � \ � in (t) .Similarly, u 0 satisfies Inner expansion.We consider the formal expansion near Ŵ(t) .Let p be an arc length parameter of Ŵ(t) (coun- ter-clockwise) and q be a distance parameter along the outward normal direction at the point x 0 (p) ∈ Ŵ(t) .Any point x in the neighborhood of Ŵ(t) is represented by x = x(p, q) = x 0 (p) + qν(p) , where ν(p) is the outer nor- mal unit vector at x 0 (p) .From this, two inverse functions p = P(x, t) and q = Q(x, t) for x ∈ Ŵ(t) are defined.Denote and U by Substituting these two functions into (3), we have Here △ p stands for Laplace-Beltrami operator.We collect the ε 1 term and obtain the relation where f , κ is the mean-curvature on Ŵ(t) , V is the normal velocity of Ŵ(t) .Assume that and U have the following expansions near the interface Ŵ(t) : Substituting ( 15) into ( 14), we know that 0 and 1 satisfy where S[ϕ 0 ](t) = α( � ϕ 0 (x, t) dx − � ϕ 0 (x, 0) dx ).Meanwhile, by repeating the same process for U, we obtain that By the matched conditions of the inner and outer solutions, U and need to satisfy some boundary conditions.Then it holds U 0 (p, q, t) ≡ u 0 (x(p, 0), t) , and then U 0 is independent of q.Moreover, since the func- tion 0 satisfies ( 16) and the matched conditions, it holds that lim q→−∞ 0 = 0 , lim q→+∞ 0 = 1 and lim q→±∞ U 0 = lim q→±0 u 0 (x(p, q), t) and hence ( 17) is rewritten by Differentiating both sides of ( 16) in q, we have σ , and then L is a self-adjoint operator (see Lemma 2.2 in the reference 36 ).Hence and the following interface equation is obtained : Combining ( 13) and (18), we obtain the free boundary problem (6).
Method for the numerical computation by the singular-limit equation.The computational area is a rectangular area � = (−L x , L x ) × (−L y , L y ) , the boundary condition is periodic and the parameters L x , L y , α and u 0 are the same as in Fig. 1.The boundaries are approximated by polygonal curves in the numerical compu- tations of the interface model (6).See the supplementary information for these details.(13) ∂u 0 ∂t = △u 0 − ku 0 + ϕ 0 , x ∈ � \ Ŵ(t), t > 0.
�(p, q, t) := �(P(x, y, t), Q(x, y, t)/ε, t) = ϕ(x, t), U(p, q, t) := U(P(x, y, t), Q(x, y, t)/ε, t) = u(x, t) In what follows, we construct the gradient flow model ( 11) of a self-propelled motion.Recall that the definition of L(∂� in ),E(� in ) , S(� in ) and H( �in ) are given by ( 8) and (9).Assume that a small perturbation of ∂� in is given by where s is the arc length of ∂� in , ε is a small positive constant, p(s) is the distance function and N(s) is the normal unit vector.Each first variation of ( 8) is represented by respectively, where κ is the curvature of ∂� in .For more details, see the supplementary information.Here we remark that, in general, δH( �in ) satisfies Using p(x) instead of u(x), we have It follows from (19) that Assume that the energy of ∂� in (t) decreases in time, and then Since ∂ ∂t p(s, t) is the normal velocity of ∂� in (t) , we conclude that Combining the interface motion (20) with the dynamics model for surfactants, we obtain the following selfpropelled motion model: where F(x, � in (t)) is given by (7).By choosing (21) is the same as (6).For the computation of the first variations of the energies L(∂� in ) , E(� in ) , and S(� in ) , see the supplementary information.

Non-dimensionalisation of
the reaction-diffusion model with dimension.We perform the following non-dimensionalisation on the reaction-diffusion model (1): Then, the model Eqs.(1) and (2) become and respectively, where We obtain the following by re-writing the non-dimensionalized mathematical model with the original variables: where Method for the numerical calculation by the reaction-diffusion model.The region is fixed as