Spot Dynamics of a Reaction- Diffusion System on the Surface of a Torus

Quasi-stationary states consisting of localized spots in a reaction-diffusion system are considered on the surface of a torus with major radius R and minor radius r. Under the assumption that these localized spots persist stably, the evolution equation of the spot cores is derived analytically based on the higher-order matched asymptotic expansion with the analytic expression of the Green's function of the Laplace--Beltrami operator on the toroidal surface. Owing to the analytic representation, one can investigate the existence of equilibria with a single spot, two spots, and the ring configuration where N localized spots are equally spaced along a latitudinal line with mathematical rigor. We show that localized spots at the innermost/outermost locations of the torus are equilibria for any aspect ratio \alpha = R r . In addition, we find that there exists a range of the aspect ratio in which localized spots stay at a special location of the torus. The theoretical results and the linear stability of these spot equilibria are confirmed by solving the nonlinear evolution of the Brusselator reaction-diffusion model by numerical means. We also compare the spot dynamics with the point vortex dynamics, which is another model of spot structures.

1. Introduction.Self-organizing beautiful patterns of localized spot-like structures appear ubiquitously in many natural phenomena.A regular/irregular lattice of spot structures is formed in Bose--Einstein condensates (BECs) [1,13].Interaction between fluid and magnetic fields gives rise to various stationary lattice configurations of small magnetic discs on a liquid-air interface [16].We can find more examples such as the formation of lattice patterns of magnetically confined electron spots in non-neutralized plasma [11], and a ring configuration of a vortex structure of an electron [12].In chemical reaction systems, it is experimentally observed that such localized spot patterns emerge in a ferrocyanide-iodate-sulphite reaction [18], a chlorine dioxide-iodine-malonic acid reaction [10], and a gas charge system [3,4].More examples are also found in [38].
In order to understand self-organization of spot patterns theoretically, it is helpful to construct phenomenological models describing the dynamics of those localized spot structures.A well-known model is vortex dynamics, which is derived from the Euler equations for incompressible and inviscid fluid flows in two-dimensional space.Suppose that the vorticity, which is defined as the curl of the velocity field, is concentrated in discrete points like Dirac's measures.We then obtain a system of ODEs describing the evolution [24].It is efficiently utilized to understand many stationary pattern formations, called vortex crystals or vortex lattices, in superfluids, BECs, fluids, and plasmas.See the survey of the history of vortex lattice theory by Newton and Chamoun [20].Another model for interacting localized spot structures is obtained from reaction-diffusion (RD) systems, in which spatially homogeneous steady states self-organize into localized spot structures due to Turing instability [34].The evolution equation describing the interactions among those spot structures is derived from two-dimensional RD equations with the asymptotic analysis [9,17].
Dynamics and pattern formations of localized spot structures can be considered on twodimensional Riemannian manifolds.A mathematical formulation of vortex dynamics on closed surfaces is found in the survey by Turner, Vitelli, and Nelson [35].RD systems on growing surfaces are derived in [22], in which the curvature and growth effects on the stability of patterns are observed numerically.In particular, owing to the geophysical and biological relevance, there are many studies on pattern formations of localized spot structures on the surface of a sphere.For example, it is shown in [21] that point vortices become a vortex crystal when they are placed on the vertices of regular polyhedrons, and the relation between the configurations and the optimal packing problem is discussed.The linear stability analysis of a ring configuration of point vortices along the line of a latitude [25] and the nonintegrability of the system [30,31] have been investigated.On the other hand, pattern formations of localized spots in RD systems on a growing sphere are used as a model of tumor growth [8] and of evolving biological surfaces [6].Formation of Turing patterns on a growing/nongrowing sphere have been studied numerically [14,39].The spherical surface is geometrically simple since it has a constant curvature.
In the meantime, another remarkable geometric feature of compact surfaces is the existence of handle structures.Hence, it is interesting to investigate how the handles affect the dynamics and the stability of localized spot patterns.One of the simplest compact surfaces is a toroidal surface with major radius R and minor radius r.Different from the surface of a sphere, it has not only nonconstant curvature but also a handle that is measured by the aspect ratio \alpha = R/r.Towards the applications to superfluids, the evolution equation of vortex dynamics on the toroidal surface has been derived in [28], in which some vortex crystals are constructed, and the dynamics of one and two point vortices are investigated.It has also been shown in [29] that the stability of a ring configuration of N point vortices changes depending on the sign of curvature and the modulus \alpha .More vortex crystals on the toroidal surface have recently been constructed [26,27].On the other hand, in RD systems, S\' anchez-Gardu\ño et al. [32] have considered Turing--Hopf bifurcations in the FitzHugh--Nagumo RD model on a growing torus and sphere.Recently, Tzou and Tzou [36] have proposed an analytic-numerical method for computing the Green's function for Helmholtz operators on curved surfaces, which is applied to derive an ODE describing a slow dynamics of N localized spots for the Schnakenberg RD model.With this model, they numerically investigate the stability of one and two localized spots.
In the present paper, we consider an RD system of the following form on a surface \scrM : where \bigtria \scrM is the Laplace--Beltrami operator on \scrM , and F u (u, v) and F v (u, v) represent the reaction terms specified as (2.2).The parameters are A, B \in \BbbR , \tau > 0, and 0 < \epsilon \ll 1.One of the examples is the Brusselator RD (BRD) model, which is used as a mathematical model of some chemical reactions [2,18,23].It is specified by , A > 0, and B = 0 in (1.1) with a parameter 0 < f < 1 satisfying \tau = 1 f 2 .Note that the model is considered on a bounded domain of a plane [7,37] as well as on the unit sphere [23,33].Another example of this type is the Schnakenberg model [36], in which \tau > 0, A = 0, B > 0, F u (u, v) = - u + u 2 v, and Our analysis is based on the higher-order matched asymptotic expansion used in [23,33,36].In section 2, we derive an ODE describing the slow dynamics of localized spot cores in quasi-equilibrium solutions of the RD system (1.1) on a toroidal surface.In section 3, using the ODE, we investigate the existence of equilibria having one spot, two spots, and N -ring spots, and we then discuss their linear stability.In the derivation, we utilize the explicit analytic formula of the Green's function of the Laplace--Beltrami operator on the toroidal surface [28].This is different from the derivation by Tzou and Tzou [36], in which the Helmholtz Green's function is constructed numerically.Owing to the analytic formula, one can conduct a rigorous mathematical analysis of the spot dynamics.We also carry out numerical simulations of the BRD system, which are compared with our theoretical results.The last section is a summary.

2.
Quasi-stationary spot solution on the surface of a torus.
It follows from that we obtain where \bigtria \bfity = \partial

Derivation of evolution equation for spot cores.
Based on the asymptotic analysis in [33], the evolution equation of N spot centers is derived from the second-order inner core problem (2.9) with the operator \scrL containing the temporal derivative in terms of \sigma .The where .

.43)
A Self-archived copy in Kyoto University Research Information Repository https://repository.kulib.kyoto-u.ac.jp Let us notice that \bfitalp is determined from the Green's function only.As shown in [36], w e j1 = \mathrm{ \mathrm{ \mathrm{ \theta j R - r \mathrm{ \mathrm{ \mathrm{ \theta j ( - \partialy 2 ) is the solution of the first equation.Since it contains no temporal derivative term, it has nothing to do with the spot dynamics.Hence, we construct the evolution equation for the jth spot by solving the second equation of (2.41) for \bfitw d j1 .By differentiating (2.8), we obtain \scrP \partial\bfitw j0 \partialy i = 0 for j = 1, 2, which means the dimension of the null-space of the adjoint operator \scrP \ast = (\bigtria \bfity + \scrM T j ) is at least two.Let us consider the homogeneous adjoint problem \scrP \ast \Psi = 0.This is solved by the separation of variables in terms of the local coordinates \bfity = (\rho cos \omega , \rho sin \omega ) T , \Psi (\rho , \omega ) = \bfitP (\rho )T (\omega ), \bfitP (\rho ) = where T (\omega ) = cos \omega or sin \omega .Substituting (2.44) into the equation, we obtain the following equation for \bfitP (\rho ): The boundary condition of \bfitP as \rho \rightar \infty is obtained as follows.Owing to (2.9) with u j0 \rightar 0 and u j0 v j0 \rightar 0 as \rho \rightar \infty , \scrM T j should satisfy This yields \bigtria \rho P 2 -\rho - 2 P 2 = 0 as \rho \rightar \infty , and we thus have P 2 = \scrO (\rho - 1 ) as \rho \rightar \infty .Normalizing \bfitP so that P 2 \sim 1 \rho as \rho \rightar \infty , we have P 1 \sim - On the other hand, since \scrP \ast \Psi = 0, substituting (2.41) into the left-hand side of (2.47) and using (2.51) Here, the constant \scrC j is defined by We note that since the solution u j0 of (2.10) depends on the strength S j and the reaction terms F u , F v , so does \scrC j .Equating (2.50) and (2.51) for T (\omega ) = cos \omega and T (\omega ) = sin \omega , we obtain the equation of the jth spot, (2.53) The evolution equation is valid as long as the localized spots of RD model (1.1) with the strengths \bfitS persist stably for a long time, and the constant \scrC j has a fixed sign independently of S j .These conditions are validated numerically for BRD model (1.2) in the next section.
2.4.Validation of the theory for Brusselator reaction-diffusion system.We construct quasi-stationary solutions u \mathrm{ \mathrm{ and v \mathrm{ \mathrm{ for BRD model (1.2) by numerical means to validate the existence of stable localized spots.That is to say, we determine the source strength \bfitS \in \BbbR N , \bfitchi \in \BbbR N , and v \in \BbbR so that they satisfy (2.10), (2.29), and (2.30) and check their stability.Let us first consider the following boundary value problem on 0 \leq \rho \leq \rho 0 for \rho 0 \gg 1 for a given scalar S: Taking \rho 0 = 20, we solve this equation with the COLNEW method [5] in the bvpSolve R library [19].We then set \chi (S) = \widehat v(\rho 0 ) -S log \rho 0 .This defines a map \chi : S \in \BbbR \mapsto \rightar \chi (S) \in \BbbR .Then, for the jth component S j of \bfitS , we obtain the approximation \chi (S j ) \approx \chi (S j ).Consequently, \widehat u and \widehat v are the approximate solutions u j0 and v j0 of (2.10) with S j .In addition, it is important to observe that the shape of the solution depends on the parameters f and S. Figure 1(a) shows that the radial solution \widehat u(\rho ) is localized when S = 2, but it tends to be volcano-shaped as S increases for f = 0.7.As a matter of fact, it is numerically confirmed that the radial solution remains localized for S \leq 3.44.Since the solution is assumed to be localized in the present asymptotic analysis, we need to restrict our attention to small S. The algorithm solving \bfitg (\bfitS ) = 0 is described in Appendix B. The plot of \chi (S) for various f is shown in Figure 1(b).Note that Figure 1(a) and (b) are the same as those in [23], although the chosen parameters are different.
Next we confirm the stability of the localized spot solutions of BRD model ( The boundary condition is given by \widehat \psi \prime (0) = \widehat \phi \prime (0) = 0, \widehat \psi \rightar 0, and \widehat \phi \rightar 0 as \rho \rightar \infty .For the approximate solutions u j0 and v j0 of the core problem (2.10) and given m, we solve (2.55) by using the finite central differences on 0 < \rho < \rho 0 = 20, which gives rise to a generalized matrix eigenvalue problem.We pay attention to the eigenvalue of (2.55) having the largest real part, say the principal eigenvalue \lambda max .Figure 2(a) shows the real part of \lambda max for fixed f = 0.7 and m = 2, 3, 4, which is the same plot as that in [23].It indicates that \lambda max is negative for small S and gets larger as S increases monotonically, and it finally becomes positive for large S. Hence, there exists a unique threshold, denoted by \Sigma m (f ), where the principal eigenvalue becomes zero.If S > \Sigma m (f ), since the real part of the principal eigenvalue is positive, the spot becomes unstable, while it is stable for S < \Sigma m (f ).Since \Sigma 2 (f ) < \Sigma 3 (f ) < \Sigma 4 (f ) for f = 0.7, the spot is stable for the modes of perturbations with m \geq 2 if S < \Sigma 2 (f ).It is important to notice that the stability of the localized spot depends not on the locations but on the strength \bfitS , the parameter f , and the mode m.
Finally, the value of C j is computed.We solve the following boundary value problem on 0 \leq \rho \leq \rho 0 with \rho 0 \gg 1 to approximate (P 1 , P 2 ) satisfying (2.45):We thus have \scrC j \approx \scrC (S j ) for given S j .Figure 2(b) shows the plot of \scrC (S) of BRD model (1.2) with f = 0.7, which is the same plot as that in [33].Let us note that \scrC j is independent of the location of the jth spot by construction, and it is negative for 0 < S < \Sigma 2 (0.7).Consequently, we conclude that the stable localized spots with 0 < S < \Sigma 2 (0.7) with a negative \scrC exist, where the equation (2.53) of the spot cores remains valid.
3. Dynamics of quasi-stationary localized spots.Suppose that localized N spots persist stably and \scrC j < 0, j = 1, . . ., N .Based on (2.53), we then find equilibrium states, meaning that N spots of RD model (1.1) are in a quasi-equilibrium state moving very slowly with \scrO (\epsilon - 2 ) time scale.The stationary N localized spots at (\theta j , \varphi j ) exist if and only if \alpha j,1 = \alpha j,2 = 0, j = 1, . . ., N .It is important to remember that \alpha j,1 and \alpha j,2 are independent of the choice of the reaction terms F u and F v , and so is the existence of the stationary spot cores.On the other hand, we need to specify the reaction terms to discuss the linear stability, since the matrix generally depends on \partialS j \partial\theta i and \partialS j \partial\varphi i , i, j = 1, 2, . . ., N , except the one-spot case.The theoretical results are compared with the nonlinear evolutions of BRD model (1.2) that are obtained numerically.
To confirm the linear stability of the one-spot case, we solve BRD model (1.2) numerically for the initial condition (2.31) and (2.32) having one spot on the torus of (R, r) = (1.1,1.0) and (R, r) = (1.3,1.0).The numerical parameters are given by \varepsi = 0.05, f = 0.7, S 1 = 3 < \Sigma 2 (0.7), and A = S 1 2\pi Rr .After computing the solution up to t = 100 when the localized spot is formed, we add a 2\% random perturbation to the solution.For \alpha = 1.1 < \alpha s , the present theory expects that the spot at \vargamm s (1.1) \approx 0.64295 is stable, whereas that at \theta = 0 and \pi are unstable.Figure 3 shows that the spot centered at \theta 1 = 0 is moving toward the stable one spot at \theta 1 = \vargamm s (1.1) after the perturbation.When \alpha = 1.3 > \alpha s , the spot at \theta 1 = 0 is stable and that at \theta 1 = \pi is unstable.Figure 4 confirms that the spot centered at \theta 1 = \pi is moving toward \theta 1 = 0 after a long-time evolution.

The \bfitN -ring configuration.
Let us consider the ring configuration of N spots located at \theta j = \vargamm N and \varphi j = (2j -1)\pi /N for j = 1, . . ., N on the toroidal surface, which we call the N -ring at \vargamm N .Then the strengths of the N spots become identical according to (2.30) and are set as S j = S c = 2\pi rRE N .This means that the existence of the N -ring is independent of the choice of the reaction terms F u and F v .It follows from (2.43) with (A.4) and (A.6) that we have From (2.43), we have .  .
Substituting \theta j = \vargamm , we have It is easy to see that if \vargamm = 0 or \pi , \alpha j,1 = 0 for j = 1, 2, . . ., N .Hence, the N -ring at the innermost/outermost location of the torus becomes an equilibrium state for any \alpha > 1.
For \vargamm \not = 0, \pi , it is sufficient to consider the existence of stationary N -rings at \vargamm \in (0, \pi ) by symmetry.We have the following theorem.
We observe the linear stability of the N -ring configuration of BRD model (1.2) for N = 2, . . ., 6 on the torus of (R, r) = ( \alpha 2 , 1 2 ) with \alpha = [1.01,10] by numerical means.The parameters are \varepsi = 0.05, f = 0.7, S c = 1.5 < \Sigma 2 (0.7), and A = N Sc 2\pi Rr .We compute the eigenvalues of the linearized matrix of \partial\theta j \partial\sigma and \partial\varphi j \partial\sigma (2.53) for j = 1, 2, . . ., N at the equilibria, thereby observing the principal eigenvalue, say \lambda max .We note that 0 is always an eigenvalue of this equilibrium originated from the invariance of the infinitesimal translation of the torus in the \varphi direction.Figure 5(a) shows the real part of the principal eigenvalue, indicating that there exists \alpha s (N ) such that the N -ring at \vargamm = 0 is neutrally stable for \alpha > \alpha s (N ), and it is unstable otherwise.Figure 5(b) shows that the N -ring at \vargamm = \pi is always unstable.The real part of the principal eigenvalue \lambda max (\alpha ) for the N -ring at \vargamm N (\alpha ) \in (0, \pi ) with N = 2, . . ., 6 in the range of \alpha \in (\alpha m (N ), \alpha M (N )) is shown in Figure 5(c).This indicates that it is unstable.Let us compare the result with that of the one-spot case in the previous section, which is equivalent to the 1-ring.According to Theorem 3.1, we find that the stable 1-ring at 1 (\alpha ) = \vargamm s (\alpha ) exists for 1 < \alpha M (1) = \alpha s (1) = \alpha s , although \alpha m (1) is not defined.On the other hand, Figure 5 indicates that \alpha m (N ) < \alpha M (N ) < \alpha s (N ) for N \geq 2.Moreover, the stability of the 1-ring at \vargamm 1 (\alpha ) is stable, whereas the N -ring at \vargamm N (\alpha ) for N \geq 2 is unstable.
The N -rings (N \geq 2) at the outermost (\theta = \pi ) and the innermost (\theta = 0) latitudinal lines are equilibria for \alpha > 1.We also obtain a range of the aspect ratio \alpha \in (\alpha m (N ), \alpha M (N )) where there exists \vargamm N (\alpha ) \in (0, \pi ) such that the N -spots at \vargamm N (\alpha ) and 2\pi -\vargamm N (\alpha ) are equilibria.We observe the linear stability of these N -ring configurations of the BRD model.The outermost N -ring is always unstable, while there exists an aspect ratio \alpha s (N ) such that the innermost one is unstable (resp., neutrally stable) for 1 < \alpha < \alpha s (N ) (resp., \alpha > \alpha s (N )).The Nrings at \vargamm N (\alpha ) and 2\pi -\vargamm N (\alpha ) are unstable equilibria.This is in contrast to the fact that the N -ring configuration of point vortices at any location is a relative equilibrium, whose linear stability is stable (resp., unstable) in the innermost (resp., outermost) region of the torus for a sufficiently large aspect ratio [29].Quasi-stationary solutions of the BRD model consisting of the unstable N -ring are numerically investigated.The unstable N -ring spots are moving toward stable quasi-stationary states having nonsymmetric configuration of N spots, indicating the existence of more nontrivial spot equilibria that are globally stable.

Figure 4 . 1 -
Figure 4. Evolution of BRD model (1.2) from a one-spot initial condition (2.31) and (2.32) centered at \theta 1 = \pi and \varphi 1 = \pi on the torus of R = 1.3 and r = 1.0, i.e., \alpha = 1.3.The numerical parameters are the same as those for Figure 3.The unstable spot starts moving toward the stable spot at \theta = 0 as expected.