On a numerical bifurcation analysis of a particle reaction-diffusion model for a motion of two self-propelled disks

Theoretical analysis using mathematical models is often used to understand a mechanism of collective motion in a self-propelled system. In the experimental system using camphor disks, several kinds of characteristic motions have been observed due to the interaction of two camphor disks. In this paper, we understand the emergence mechanism of the motions caused by the interaction of two self-propelled bodies by analyzing the global bifurcation structure using the numerical bifurcation method for a mathematical model. Finally, it is also shown that the irregular motion, which is one of the characteristic motions, is chaotic motion and that it arises from periodic bifurcation phenomena and quasi-periodic motions due to torus bifurcation.


Introduction
Since the work of T. Vicsek et al. many theoretical analyses using particle motion models have been conduced to understand the collective motion of biological species [1]. For example, flocks of birds [2], schools of fish [3], and insect swarms [4] have been investigated. In addition, physicists have developed many mathematical models of collective motion assuming various interactions [6][7][8][9][10][11]. Recently, there is also an analysis of collective motion that considers non-local interactions [12].
In particular, there is a macroscopic experimental system called a self-propelled system in which individual materials move by changing the state of the field. For example, camphor particles [19] and oil droplets [20,21] move in translational motion by changing the surface tension of the water.
Many researches have been performed to understand the mechanism of collective motions using such self-propelled systems, and collective motions on annular water channel using camphor boats and camphor disks [22,23] and on the water surface of petri dishes [24][25][26] have been reported. In addition, a self-assembly motion of droplets has been reported by using the chemical reaction [27]. The mathematical model has also been proposed to understand the collective motion for such macroscopic self-propelled material. For example, there is a mathematical model for the collective motion of camphor boats in the annular water channel [28] and the collective oscillatory motion of camphor disks [29], and there is also a mathematical model for the self-assembly motion of droplets with the chemical reaction [30].
One of the mathematical models to describe the above motion of self-propelled material is the particle reaction-diffusion model as follows [31]: where x i c (t) and u(x, t) are the center of an i-th disk, which is the sufactant and a surface concentration of the surfactant, respectively. , , d u and k denote an area density of the disk, a friction coefficient, a diffusion coefficient of surfactant molecules on a water surface and a combined rate of sublimation and dissolution of the surfactant, respectively. When N = 1 , (1) represents the motion of a single selfpropelled material, and when N ≥ 2 , it represents the collective motion of self-propelled material. The function G(u) is a driving force to the self-propelled material, which is given by, for example, or where the function (u) > 0 describe the surface tension on the water surface, for instance, which is the following strictly decreasing function of u > 0:

3
On a numerical bifurcation analysis of a particle reaction-diffusion model where a > 0 , m ∈ ℕ , and 0 , 1 > 0 represent the surface tension of pure water and that of the critical micelle concentration of the surfactant, respectively. The function F(x) represents the supply of the surfactant molecules and a simple example of it is the supply from solid surfactant, which is described by where k u and u 0 are the supply rate and the density of the solid surfactant, respectively, and r is the radius of solid surfactant. (x) is the Dirac delta function, which corresponds the point source. In order to theoretically understand the mechanism of self-propelled material motion, much computer-aided and formal analysis using mathematical models have been investigated [32][33][34][35][36]. Furthermore, as rigorous mathematical analysis of the mathematical model (1), the global existence of solutions in the whole space and the existence of unique traveling wave solutions [37], rotational motion solutions for one-dimensional periodic boundary conditions in space have been shown for the piecewise constant function supply term [38]. In addition, the existence of weak solutions is shown for the delta function supply term [39]. To understand the collective motion of camphor disks in the annular water channel [40], in 2015 Nishi et al. analyzed the interaction of two camphor disks using the following dimensionless mathematical model without the inertia term, which assumes that camphor disks are very light [41]: where > 0 and L > 2Nr > 0 , and the function (u) and F(x) are given by where a, 0 , 1 are positive constants. L is a map from ℝ to [−L∕2, L∕2) that is equivalent to a canonical projection ℝ → ℝ∕ℤ ≅ S 1 to which the annular water channel is considered to be homeomorphic. When N = 2 , from the computer-aided analysis of (3) under the periodic boundary condition and the additional experimental results, the model system (3) was shown to qualitatively reproduce well the motion of two camphor disks [41]. The results also suggest that the traffic jam motion and billiard-like motions reported in the experiments can be reproduced. They are reproduced qualitatively, as shown in Fig. 1. In this study, the existence and stability of the motionless solution, the symmetric rotating solution (Fig. 2a), and the asymmetric rotating solution (Fig. 2b) have also been investigated by computer-aided analysis. In addition, characteristic numerical results such as an asymmetric oscillatory solution (Fig. 3a), an oscillatory rotating solution (Fig.3b, c) and an irregular motion solution( Fig. 3d) have been shown, but thebifurcation mechanism by which these solutions emerge is not known. On a numerical bifurcation analysis of a particle reaction-diffusion model In our study, we theoretically show the mechanism of the emergence of complex motions by analyzing the global bifurcation structure of (3) on the annular water channel when N = 2 and L = 10 , using the numerical bifurcation method. Furthermore, the bifurcation mechanism in which irregular motions appear is explained, and the chaotic nature of irregular motions is clarified by numerical calculation of the maximum Lyapunov exponent. This paper is organized as follows: the solution to be obtained by the numerical bifurcation is defined in Chapter 2. In Chapter 3, based on the definition of the solution, numerical solutions are obtained by the bifurcation method, and the linearization stability of the solutions is also shown numerically. In Chapter 4, we discuss the results of the bifurcation structure in Chapter 3 and focus on irregular motions to investigate the mechanism of their appearance. Finally, we summarize the entire study and discuss future work.

Mathematical definitions of disk motions
To analyze motions of disks on annular water channel, we introduce mathematical equations of which solutions correspond to actual motions. Firstly, we change the coordinate system of (3) from stationary to moving coordinate system z = x − ct where c represents a system velocity, which allows us to treat rotating motions such where . We define four types of motion as solutions of mathematical equations using (7). Note that we suppose that z 1

Definition 1 (Motionless solution)
A stationary solution of (7) with c = 0 is called a motionless solution.

Definition 2 (Rotating solution)
A stationary solution of (7) with c ≠ 0 is called a rotating solution.

Definition 3 (Oscillatory solution)
A time-periodic solution of (7) with c = 0 and period T > 0 is called an oscillatory solution.

Definition 4 (Oscillatory rotating solution)
A time-periodic solution of (7) with c ≠ 0 and period T > 0 is called an oscillatory rotating solution.

Numerical results
Based on the definition explained in Chapter 3, we numerically analyzed the global bifurcation structure of (7).

Newfound motions and their definitions
In the previous researches [40,41], only eight types of solution: symmetric motionless solution, symmetric/asymmetric rotating solution, symmetric/asymmetric oscillatory solution, symmetric/asymmetric oscillatory rotating solution and irregular motion, have been reported, where the symmetries are about positions and amplitudes of disks for motionless/rotating solutions and oscillatory/oscillatory rotating solutions, respectively. As the result of our bifurcation analysis, we found another four types of motion: round-trip-asymmetric motion ( Fig. 5(a-2)), complete in-phase motion (Fig. 5b), period-doubled motion (Fig. 5(c-2)) and quasi-periodic motion ). Round-trip-asymmetric motion is an oscillatory motion with different trajectories between the first and second half of its round-trip. In-phase motion is one of the periodic motions. Unlike the other motions, the disks oscillate in the same directions as shown in Fig. 5b. Period-doubled motion is a motion derived from a periodic motion. Although its trajectory is similar to the original one, there is a difference between the two successive periodic orbits. Its period is doubled from the original period because of the difference. Quasi-periodic motion has two different periods in its motion. In some figures, we negate constant rotating velocity by giving non-zero c for visibility. From the observation of these motions, we introduce additional mathematical definitions of motion to classify them. In the following definitions, we fix N = 2 for   The definition above means that the trajectory of the first half of their round-trip is equal to the inverted one of the second half. Figure 5(a-1) and (a-2) show example motions which correspond to a round-trip-symmetric and asymmetric solution, respectively. Note that we do not consider round-trip-symmetry on oscillatory rotating solutions because they are usually round-trip-asymmetric due to their constant rotating velocity.

Definition 8 (Complete in-phase solution) Suppose
is an oscillatory or oscillatory rotating solution of (7) with N = 2 . The solution is called a complete in-phase solution if holds for all t > 0. Figure 5b shows an example motion which corresponds to a complete in-phase solution.
For clarity of explanation, we assign a symbol to each type of motion we found based on the following assignment table (Table 1). A sequence of letters combined along each row corresponds to each type of motion. We omit unnecessary letters for simplicity.

Bifurcation diagram
We numerically obtained the bifurcation diagram in the case of N = 2 and L = 10 ( Fig. 6), which is known as the parameter region where complicated solutions appear. Figure 7 shows typical examples of all the types of motion we found. For visibility of period-doubling transitions, we added z 1 c -z 2 c plots of 2 n -periodic solutions (Fig. 8).
Previously only limited parts of the bifurcation structure, which were stationary solution branches, were examined. This numerical bifurcation analysis revealed the structure among not only stationary solutions but periodic ones. Owing to that, we were able to make mainly five notable observations about periodic motions including chaotic ones: (1) there are stable branches of OSA (round-trip-asymmetric motion) which is a newfound motion. (2) There is an OrSI (complete-in-phase motion, which is also newfound) branch which is not directly connected to any stable branches. (3) There is a complicated part of the bifurcation structure hidden by the stable branch of OrSN around 0.002 ≤ ≤ 0.005 (Fig. 9), which was not discovered in previous researches even though it has some stable branches. (4) A bifurcation process that is probably a period-doubling cascade exists near = 0.016 On a numerical bifurcation analysis of a particle reaction-diffusion model ( Fig. 10). (5) A torus bifurcation point exists near = 0.004 . The last two observations can give us some hints of the mechanism whereby the chaotic motion emerges, which has never been discussed before even though it appears in a large portion of the parameter space. Thus we discuss the mechanism in the next chapter.

Analysis of the Chaotic Motion
Considering the bifurcation structures we obtained, there are two possible bifurcation routes to the chaotic motion via period-doubling bifurcation and torus one. Both have well-known mechanisms by which chaotic motions emerge: period-doubling cascade and torus breakdown [42,43], respectively. In this section, we numerically compute the maximum Lyapunov exponent [44], which is an indicator of chaotic motion [42], and investigate the relationship between those two bifurcation processes and chaotic motion in this system. Figure 11 shows that both ends of the parameter region where chaotic motions in terms of Lyapunov exponent occur, that is, the region where the maximum Lyapunov exponent is positive, coincide with the and v i (t) represent the distance to the next disk of i-th disk and velocity on the original(stationary) coordinate system at time t, respectively, is defined to be able to separate each branch curve. Solid lines indicate stable solutions, whereas dashed lines indicate unstable solutions. Five types of bifurcation point: Hopf, pitchfork, saddle-node, period-doubling and torus are represented by square, triangle, diamond, star and hexagon, respectively. Some branches disappear in the middle because the continuation of those branches was impossible due to their strong instability. 1 The branch of Qpr is not on this diagram because it does not match with Definition 1, 2, 3 or 4. However, we numerically obtained stable Qpr solutions ( Fig. 9(m-1) and (m-2)) near the torus bifurcation point by computing time evolution. Therefore, we state that there exists a branch of Qpr parameters where the period-doubling bifurcations repeatedly occur, and the saddlenode bifurcation point of the stable OrSN branch lies. The coincidence at the right side of the region (near = 0.016 ) implies that the repeated period-doubling bifurcations are most likely to be a part of period-doubling cascade, which leads to chaos. In addition, oscillatory and oscillatory rotating solutions sporadically exist around 0.014 ≤ ≤ 0.016 which are likely to be periodic windows, which are incidental to period-doubling cascade. Figure 13 shows trajectories of periodic solutions found in the periodic windows. Moreover, chaotic attractors [42] are likely to exist near the period-doubling cascade (Fig. 12). The left end of the chaotic region coincides with the saddle-node point of the OrSN branch and the transition from OrSN to

3
On a numerical bifurcation analysis of a particle reaction-diffusion model chaos occurs suddenly according to the very sharp change of the maximum Lyapunov exponent. By our bifurcation analysis, it can be explained as follows: there should exist a torus breakdown point which emerges from the torus bifurcation point near = 0.004 . That should lead to chaos, however, the basin of the stable OrSN is so large that usual numerical calculations result in converging to OrSN solutions quickly (Fig. 14), which is the possible reason why chaotic behaviors suddenly appear after the OrSN branch disappears at the saddle-node point.

Conclusion
This study clarifies the motion of two camphor disks in an annular water channel by numerical bifurcation calculations for the particle reaction-diffusion model under periodic boundary conditions. As the result, we were able to find new stable solutions (the round-trip-asymmetric solution, N-period solution ( N = 2, 4, 8 ) and quasi-periodic solution) that could not be found in previous studies. In addition, the  On a numerical bifurcation analysis of a particle reaction-diffusion model numerical calculation of the maximum Lyapunov exponent showed that the motion shown as irregular motion in the previous result [41] is chaotic motion. The mathematical analysis of our model of the two self-propelled materials is limited to the existence of motionless, rotating, and asymmetric rotating solutions. Therefore, we want to analyze the existence and stability of periodic solutions in the future.  On a numerical bifurcation analysis of a particle reaction-diffusion model

Appendix A Discretization
In order to numerically analyze the bifurcation structure of (7), we need to discretize the system. It is preferable to employ Fourier series expansion rather than naive finite-difference method for discretization of v(z, t) because the term v( L (z i c ± r), t) may require values of v on arbitrary points in [−L∕2, L∕2).
Let (Z 1 c (t), … , Z N c (t)) be the positions of disks in the discretized system. Let us define a function V ∶ ℝ 2 → ℝ as follows: which represents v in the form of truncated Fourier series with wave number M and real coefficients (a 0 (t), a 1 (t), … , a M (t), b 1 (t), … , b M (t)) . In the following, we obtain the discretized system of (7), that is, the set of ordinary differential equations between (Z 1 c (t), … , Z N c (t)) and (a 0 (t), a 1 (t), … , a M (t), b 1 (t), … , b M (t)). First, we substitute v and z i c in (7) with V and Z i c . Then, we obtain The left-hand side of the second equation is The right-hand side is where and for k = 1, … , M . Note that here F is Next, let us integrate both sides from z = −L∕2 to L/2 then we obtain F(x) = 1, |x| ≤ r, 0, |x| > r.

3
On a numerical bifurcation analysis of a particle reaction-diffusion model Similarly, we multiply 2 cos(2 kz∕L)∕L on both sides and integrate them from z = −L∕2 to L/2 then the following is obtained: for k = 1, … , M . For b k , we multiply 2 sin(2 kz∕L)∕L on both sides and integrate them from z = −L∕2 to L/2 then the following is obtained: for k = 1, … , M . To sum up, we obtained Combining the first equation of (8), the discretized system of (8) to be obtained is Typically we used M = 128 in this paper.

Appendix B Numerical methods
For solving the discretized system (9), we employed a four-stage fourth-order L-stable Rosenbrock method called ROS4 [46] by E. Hairer and G. Wanner [47]. Pseudoarclength continuation method [48] was used for branch continuation. To find periodic solutions, we employed Newton Shooting Method [52]. To compute eigenvalues and eigenvectors, we used ARPACK [51] for sparse matrices and the routine implemented in Eigen 3.3.9 [49] which is based on EISPACK [50] for dense matrices.