An Asymptotic Analysis of Localized 3-D Spot Patterns for Gierer-Meinhardt Model: Existence, Linear Stability and Slow Dynamics
AAn Asymptotic Analysis of Localized 3-D Spot Patterns for the Gierer-Meinhardt Model:Existence, Linear Stability and Slow Dynamics
Daniel Gomez ∗ , Michael J. Ward † , and Juncheng Wei ‡ Abstract.
Localized spot patterns, where one or more solution components concentrates at certain points in the domain, area common class of localized pattern for reaction-diffusion systems, and they arise in a wide range of modelingscenarios. Although there is a rather well-developed theoretical understanding for this class of localized pattern inone and two space dimensions, a theoretical study of such patterns in a 3-D setting is, largely, a new frontier. In anarbitrary bounded 3-D domain, the existence, linear stability, and slow dynamics of localized multi-spot patternsis analyzed for the well-known singularly perturbed Gierer-Meinhardt (GM) activator-inhibitor system in the limitof a small activator diffusivity ε (cid:28)
1. Our main focus is to classify the different types of multi-spot patterns, andpredict their linear stability properties, for different asymptotic ranges of the inhibitor diffusivity D . For the range D = O ( ε − ) (cid:29)
1, although both symmetric and asymmetric quasi-equilibrium spot patterns can be constructed,the asymmetric patterns are shown to be always unstable. On this range of D , it is shown that symmetric spotpatterns can undergo either competition instabilities or a Hopf bifurcation, leading to spot annihilation or temporalspot amplitude oscillations, respectively. For D = O (1), only symmetric spot quasi-equilibria exist and they arelinearly stable on O (1) time intervals. On this range, it is shown that the spot locations evolve slowly on an O ( ε − )time scale towards their equilibrium locations according to an ODE gradient flow, which is determined by a discreteenergy involving the reduced-wave Green’s function. The central role of the far-field behavior of a certain coreproblem, which characterizes the profile of a localized spot, for the construction of quasi-equilibria in the D = O (1)and D = O ( ε − ) regimes, and in establishing some of their linear stability properties, is emphasized. Finally, forthe range D = O ( ε ), it is shown that spot quasi-equilibria can undergo a peanut-splitting instability, which leadsto a cascade of spot self-replication events. Predictions of the linear stability theory are all illustrated with full PDEnumerical simulations of the GM model.
1. Introduction.
We analyze the existence, linear stability, and slow dynamics of localized N -spotpatterns for the singularly perturbed dimensionless Gierer-Meinhardt (GM) reaction-diffusion (RD) model(cf. [7])(1.1) v t = ε ∆ v − v + v u , τ u t = D ∆ u − u + ε − v , x ∈ Ω ; ∂ n v = ∂ n u = 0 , x ∈ ∂ Ω , where Ω ⊂ R is a bounded domain, ε (cid:28)
1, and v and u denote the activator and inhibitor fields,respectively. While the shadow limit in which D → ∞ has been extensively studied (cf. [20], [22], [19]),there have relatively few studies of localized RD patterns in 3-D with a finite inhibitor diffusivity D (see [2],[5], [10], [16] and some references therein). For 3-D spot patterns, the existence, stability, and slow-dynamicsof multi-spot quasi-equilibrium solutions for the singularly perturbed Schnakenberg RD model was analyzedusing asymptotic methods in [16]. Although our current study is heavily influenced by [16], our resultsfor the GM model offer some new insights into the structure of localized spot solutions for RD systemsin three-dimensions. In particular, one of our key findings is the existence of two regimes, the D = O (1)and D = O ( ε − ) regimes, for which localized patterns can be constructed in the GM-model, in contrastto the single D = O ( ε − ) regime where such patterns occur for the Schnakenberg model. Furthermore,our analysis traces this distinction back to the specific far-field behaviour of the appropriate core problem, ∗ Dept. of Mathematics, UBC, Vancouver, Canada. (corresponding author [email protected] ) † Dept. of Mathematics, UBC, Vancouver, Canada. [email protected] ‡ Dept. of Mathematics, UBC, Vancouver, Canada. [email protected] a r X i v : . [ n li n . PS ] A ug haracterizing the local behavior of a spot, for the GM-model. By numerically solving the core problem, weformulate a conjecture regarding the far-field limiting behavior of the solution to the core problem. Withthe numerically established properties of the core problem, strong localized perturbation theory (cf. [17])is used to construct N -spot quasi-equilibrium solutions to (1.1), to study their linear stability, and todetermine their slow-dynamics. We now give a more detailed outline of this paper.In the limit ε →
0, in § N -spot quasi-equilibrium solutions to (1.1). To do so, we firstformulate an appropriate core problem for a localized spot, from which we numerically compute certainkey properties of its far field behavior. Using the method of matched asymptotic expansions, we thenestablish two distinguished regimes for the inhibitor diffusivity D , the D = O (1) and D = O ( ε − ) regimes,for which N -spot quasi-equilibria exist. By formulating and analyzing a nonlinear algebraic system, wethen demonstrate that only symmetric patterns can be constructed in the D = O (1) regime, whereas bothsymmetric and asymmetric patterns can be constructed in the D = O ( ε − ) regime.In § O (1) time scale of the N -spot quasi-equilibria constructedin §
2. More specifically, we use the method of matched asymptotic expansions to reduce a linearizedeigenvalue problem to a single globally coupled eigenvalue problem. We determine that the symmetricquasi-equilibrium patterns analyzed in § D = O (1) regime but thatthey may undergo both Hopf and competition instabilities in the D = O ( ε − ) regime. Furthermore, wedemonstrate that the asymmetric patterns studied in § D = O ( ε − ) regime are always unstable.Our stability predictions are then illustrated in § § D = O ( ε ), where localized spots interactweakly through exponentially small terms. In this regime, (1.1) can be reduced to a modified core problemfrom which we numerically calculate quasi-equilibria and determine their linear stability properties. Unlikein the D = O (1) and D = O ( ε − ) regimes, we establish that spot solutions in the D = O ( ε ) regimecan undergo peanut-splitting instabilities. By performing full numerical simulations using FlexPDE6 [6],we demonstrate that these instabilities lead to a cascade of spot self-replication events in 3-D. Althoughspike self-replication for the 1-D GM model have been studied previously in the weak interaction regime D = O ( ε ) (cf. [4], [8], [11]), spot self-replication for the 3-D GM model has not previously been reported.In § § §
2. Asymptotic Construction of an N -Spot Quasi-Equilibrium Solution. In this section we asymp-totically construct an N -spot quasi-equilibrium solution where the activator is concentrated at N spec-ified points that are well-separated in the sense that x , . . . , x N ∈ Ω, | x i − x j | = O (1) for i (cid:54) = j , anddist( x i , ∂ Ω) = O (1) for i = 1 , . . . , N . In particular, we first outline the relevant core problem and de-scribe some of its properties using asymptotic and numerical calculations. Then, the method of matchedasymptotic expansions is used to derive a nonlinear algebraic system whose solution determines the quasi-equilibrium pattern. A key feature of this nonlinear system, in contrast to that derived in [16] for the3-D Schnakenberg model, is is that it supports different solutions depending on whether D = O (1) or D = O ( ε − ). More specifically, we will show that the D = O (1) regime admits only spot quasi-equilibriathat are symmetric to leading order, whereas the D = O ( ε − ) regime admits both symmetric and asym- .00 0.05 0.10 0.15 0.20 0.25 S ( S ) (a) V ( , S ) S0.010.050.1 S (b) U ( , S ) S0.010.050.1 S (c)Figure 1: Plots of numerical solutions of the core problem (2.1): (a) µ ( S ) versus S , as well as the (b) activator V and (c) inhibitor U , at a few select values of S . The value S = S (cid:63) ≈ . µ ( S ) = 0. metric N -spot quasi-equilibria. A key step in the application of the method of matched asymptotic expansionsto construct localized spot patterns is the study of the core problem∆ ρ V − V + U − V = 0 , ∆ ρ U = − V , ρ > , (2.1a) ∂ ρ V (0) = ∂ ρ U (0) = 0 ; V −→ U ∼ µ ( S ) + S/ρ , ρ → ∞ , (2.1b)where ∆ ρ ≡ ρ − ∂ ρ (cid:2) ρ ∂ ρ (cid:3) . For a given value of S >
0, (2.1) is to be solved for V = V ( ρ ; S ), U = U ( ρ ; S ),and µ = µ ( S ). Specifying the value of S > L ( R ) norm of V , as can beverified by applying the divergence theorem to the second equation in (2.1a) over an infinitely large ball,which yields the identity S = (cid:82) ∞ ρ [ V ( ρ )] dρ .When S (cid:28) V = O ( √ S ). By applying the divergence theorem tothe first equation in (2.1a) we get U = O ( √ S ), while from (2.1b) we conclude that µ = O ( √ S ). It is thenstraightforward to compute the leading order asymptotics(2.2) V ( ρ ; S ) ∼ (cid:114) Sb w c ( ρ ) , U ( ρ ; S ) ∼ (cid:114) Sb , µ ( S ) ∼ (cid:114) Sb , for S (cid:28) , where b ≡ (cid:82) ∞ ρ [ w c ( ρ )] dρ ≈ .
423 and w c > ρ w c − w c + w c = 0 , ρ > ∂ ρ w c (0) = 0 , w c → ρ → ∞ . We remark that (2.3) has been well studied, with existence being proved using a constrained variationalmethod, while its symmetry and decay properties are established by a maximum principle (see for exampleAppendix 13.2 of [22]). The limit case S (cid:28) shadow limit obtained by taking D → ∞ ,for which numerous rigorous and asymptotic results have previously been obtained (cf. [20], [22], [19]).Although the existence of solutions to (2.1) have not been rigorously established, we can use the small S asymptotics given in (2.2) as an initial guess to numerically path-follow solutions to (2.1) as S is increased.The results of our numerical computations are shown in Figure 1 where we have plotted µ ( S ), V ( ρ ; S ), nd U ( ρ ; S ) for select values of S >
0. A key feature of the plot of µ ( S ) is that it has a zero crossingat S = 0 and S = S (cid:63) ≈ . ≤ S ≤ S (cid:63) at S = S crit ≈ . µ (cid:48)(cid:48) ( S ) < < S ≤ S (cid:63) . Themajority of our subsequent analysis hinges on these numerically determined properties of µ ( S ). We leavethe task of rigorously proving the existence of solutions to (2.1) and establishing the numerically verifiedproperties of µ ( S ) as an open problem, which we summarize in the following conjecture: Conjecture 2.1.
There exists a unique value of S (cid:63) > such that (2.1) admits a ground state solutionwith the properties that V, U > in ρ > and for which µ ( S (cid:63) ) = 0 . Moreover, µ ( S ) satisfies µ ( S ) > and µ (cid:48)(cid:48) ( S ) < for all < S < S (cid:63) . We now proceed with the method ofmatched asymptotic expansions to construct quasi-equilibria for (1.1). First we seek an inner solution byintroducing local coordinates y = ε − ( x − x i ) near the i th spot and letting v ∼ DV i ( y ) and u ∼ DU i ( y ) sothat the local steady-state problem for (1.1) becomes(2.4) ∆ y V i − V i + U − i V i = 0 , ∆ y U i − ε D − U i + V i = 0 , y ∈ R . In terms of the solution to the core problem (2.1) we determine that(2.5) V i ∼ V ( ρ, S iε ) + O ( D − ε ) , U i ∼ U ( ρ, S iε ) + O ( D − ε ) , ρ ≡ | y | = ε − | x − x i | , where S iε is an unknown constant that depends weakly on ε . We remark that the derivation of the nextorder term requires that x , . . . , x N be allowed to vary on a slow time scale. This higher order analysis isdone in § S ε , . . . , S Nε we now derive a nonlinear algebraic system (NAS) by matching inner andouter solutions for the inhibitor field. As a first step, we calculate in the sense of distributions that ε − v −→ πD (cid:80) Nj =1 S jε δ ( x − x j ) + O ( ε ) as ε → + . Therefore, in the outer region the inhibitor satisfies(2.6) ∆ u − D − u = − πεD N (cid:88) j =1 S jε δ ( x − x j ) + O ( ε ) , x ∈ Ω ; ∂ n u = 0 , x ∈ ∂ Ω . To solve (2.6), we let G ( x ; ξ ) denote the reduced-wave Green’s function satisfying∆ G − D − G = − δ ( x − ξ ) , x ∈ Ω ; ∂ n G = 0 , x ∈ ∂ Ω ,G ( x ; ξ ) ∼ π | x − ξ | + R ( ξ ) + ∇ x R ( x ; ξ ) · ( x − ξ ) , as x → ξ , (2.7)where R ( ξ ) is the regular part of G . The solution to (2.6) can be written as(2.8) u ∼ πεD N (cid:88) j =1 S jε G ( x ; x j ) + O ( ε ) . Before we begin matching inner and outer expansions to determine S ε , . . . , S Nε we first motivate twodistinguished limits for the relative size of D with respect to ε . To do so, we note that when D (cid:29) G ( x, ξ ) ∼ D | Ω | − + G ( x, ξ ) + O ( D − ) , here G ( x, ξ ) is the Neumann Green’s function satisfying∆ G = 1 | Ω | − δ ( x − ξ ) , x ∈ Ω ; ∂ n G = 0 , x ∈ ∂ Ω ; (cid:90) Ω G dx = 0 , (2.10a) G ( x, ξ ) ∼ π | x − ξ | + R ( ξ ) + ∇ x R ( x ; ξ ) · ( x − ξ ) , as x → ξ , (2.10b)and R ( ξ ) is the regular part of G . In summary, for the two ranges of D we have(2.11) G ( x, ξ ) ∼ π | x − ξ | + (cid:40) R ( ξ ) + o (1) , D = O (1) ,D | Ω | − + R ( ξ ) + o (1) , D (cid:29) , as | x − ξ | → , where R ( ξ ) is the regular part of G ( x, ξ ). By matching the ρ → ∞ behaviour of U i ( ρ ) given by (2.5) withthe behaviour of u given by (2.8) as | x − x i | →
0, we obtain in the two regimes of D that(2.12) µ ( S iε ) = 4 πε (cid:40) S iε R ( x i ) + (cid:80) j (cid:54) = i S jε G ( x i , x j ) , D = O (1) ,S iε R ( x i ) + (cid:80) j (cid:54) = i S jε G ( x i , x j ) + D | Ω | − (cid:80) Nj =1 S jε , D (cid:29) . ¿From the D (cid:29) D = O ( ε − ) is a distinguished regime for which the right-hand sidehas an O (1) contribution. Defining the vectors SSS ε ≡ ( S ε , . . . , S Nε ) T , µ ( SSS ε ) ≡ ( µ ( S ε ) , . . . , µ ( S Nε )) T , and eee ≡ (1 , . . . , T , as well as the matrices E N , G , and G by(2.13) E N ≡ N eeeeee T , ( G ) ij = (cid:40) R ( x i ) , i = jG ( x i , x j ) , i (cid:54) = j , ( G ) ij = (cid:40) R ( x i ) , i = jG ( x i , x j ) , i (cid:54) = j , we obtain from (2.12) that the unknowns S ε , . . . , S Nε must satisfy the NAS µ ( SSS ε ) = 4 πε G SSS ε , for D = O (1) , (2.14a) µ ( SSS ε ) = κ E N SSS ε + 4 πε G SSS ε , for D = ε − D , where κ ≡ πN D | Ω | . (2.14b) N -Spot Quasi-Equilibrium. We now determine solutions to the NAS(2.14) in both the D = O (1) and the D = O ( ε − ) regimes. In particular, we show that it is possible toconstruct symmetric N -spot solutions to (1.1) by finding a solution to the NAS (2.14) with SSS ε = S cε eee inboth the D = O (1) and D = O ( ε − ) regimes. Moreover, when D = O ( ε − ) we will show that it is possibleto construct asymmetric quasi-equilibria to (1.1) characterized by spots each having one of two strengths.When D = O (1) the NAS (2.14a) implies that to leading order µ ( S iε ) = 0 for all i = 1 , . . . , N . Fromthe properties of µ ( S ) outlined in § µ ( S ) in Figure 1a, we deduce that S iε ∼ S (cid:63) for all i = 1 , . . . , N . Thus, to leading order, N -spot quasi-equilibria in the D = O (1) regime havespots with a common height, which we refer to as a symmetric pattern. By calculating the next orderterm using (2.14a) we readily obtain the two term result(2.15) SSS ε ∼ S (cid:63) eee + 4 πεS (cid:63) µ (cid:48) ( S (cid:63) ) G eee . We conclude that the configuration x , . . . , x N of spots only affects the spot strengths at O ( ε ) through theGreen’s matrix G . Note that if eee is an eigenvector of G with eigenvalue g then the solution to (2.14a) is SSS iε = S cε eee where S cε satisfies the scalar equation µ ( S cε ) = 4 πεg S cε . ext, we consider solutions to the NAS (2.14b) in the D = ε − D regime. Seeking a solution SSS ε ∼ SSS + εSSS + · · · we obtain the leading order problem(2.16) µ ( SSS ) = κ E N SSS . Note that the concavity of µ ( S ) (see Figure 1a) implies the existence of two values 0 < S l < S r < S (cid:63) suchthat µ ( S l ) = µ ( S r ). Thus, in addition to the symmetric solutions already encountered in the D = O (1)regime, we also have the possibility of asymmetric solutions, where the spots can have two different heights.We first consider symmetric solutions, where to leading order SSS = S c eee in which S c satisfies(2.17) µ ( S c ) = κS c . The plot of µ ( S ) in Figure 1a, together with the S (cid:28) < S c ≤ S (cid:63) for all κ >
0. In Figure 3a we illustrate graphically thatthe common spot strength S c is obtained by the intersection of µ ( S ) with the line κS . We refer to Figure4 for plots of the symmetric solution strengths as a function of κ . In addition, we readily calculate that(2.18) S c ∼ S (cid:63) (cid:18) κµ (cid:48) ( S (cid:63) ) (cid:19) + O ( κ ) , for κ (cid:28) S c ∼ bκ + O ( κ − ) , for κ (cid:29) , which provide a connection between the D = O (1) and D → ∞ (shadow limit) regimes, respectively.From (2.14b), the next order correction SSS satisfies µ (cid:48) ( S c ) SSS − κ E N SSS = 4 πS c G eee . Upon left-multiplyingthis expression by eee T we can determine eee T SSS . Then, by recalling the definition of E N ≡ N − eeeeee T we cancalculate SSS . Summarizing, a two term asymptotic expansion for the symmetric solution to (2.14b) is(2.19) SSS ε ∼ S c eee + 4 πεµ (cid:48) ( S c ) (cid:18) S c I N + µ ( S c ) µ (cid:48) ( S c ) − κ E N (cid:19) G eee , provided that µ (cid:48) ( S c ) (cid:54) = 0 (i.e. S c (cid:54) = S crit ). Note that µ (cid:48) ( S c ) − κ = 0 is impossible by the following simpleargument. First, for this equality to hold we require that 0 < S < S crit since otherwise µ (cid:48) ( S c ) < κ to get µ (cid:48) ( S c ) − κ = S − c g ( S c ) where g ( S ) ≡ Sµ (cid:48) ( S ) − µ ( S ). However,we calculate g (cid:48) ( S ) = Sµ (cid:48)(cid:48) ( S ) < S asymptotics found in (2.2) we determinethat g ( S ) ∼ − (cid:112) S/ (4 b ) < S → + . Therefore, g ( S c ) < < S c < S crit so that µ (cid:48) ( S c ) < κ holds. Finally, as for the D = O (1) case, if G eee = g eee then the common source values extends to higherorder and we have SSS ε = S cε eee where S cε is the unique solution to the scalar problem(2.20) µ ( S cε ) = ( κ + 4 πεg ) S cε . Next, we construct of asymmetric N -spot configurations. The plot of µ ( S ) indicates that for any valueof S r ∈ ( S crit , S (cid:63) ] there exists a unique value S l = S l ( S r ) ∈ [0 , S crit ) satisfying µ ( S l ) = µ ( S r ). A plot of S l ( S r ) is shown in Figure 2a. Clearly S l ( S crit ) = S crit and S l ( S (cid:63) ) = 0. We suppose that to leading orderthe N -spot configuration has n large spots of strength S r and N − n small spots of strengths S l . Morespecifically, we seek a solution of the form(2.21) SSS ε ∼ ( S r , . . . , S r , S l ( S r ) , . . . , S l ( S r )) T , so that (2.16) reduces to the single scalar nonlinear equation(2.22) µ ( S r ) = κf ( S r ; n/N ) , on S crit < S r < S (cid:63) , where f ( S ; θ ) ≡ θS + (1 − θ ) S l ( S ) . .05 0.10 0.15 0.20 S r S l ( S r ) (a) S r S l ( S r ) (b) f ( S ; ) (c)Figure 2: Plots of (a) S l ( S r ) and (b) S (cid:48) l ( S r ) for the construction of asymmetric N -spot patterns. (c) Plots of f ( S, θ )for select values of θ ≡ n/N . For 0 < θ < . f ( S, θ ) attains an interior minimum in S crit < S < S (cid:63) . Since µ ( S crit ) − κf ( S crit ; n/N ) = µ ( S crit ) − κS crit and µ ( S (cid:63) ) − κf ( S (cid:63) ; n/N ) = − κnS (cid:63) /N <
0, we obtain bythe intermediate value theorem that there exists at least one solution to (2.22) for any 0 < n ≤ N when0 < κ < κ c ≡ µ ( S crit ) /S crit ≈ . . Next, we calculate f (cid:48) ( S ; θ ) = (1 − θ ) (cid:18) θ − θ + S (cid:48) l ( S ) (cid:19) , where S (cid:48) l ( S ) is computed numerically (see Figure 2b). We observe that − ≤ S (cid:48) l ( S r ) ≤ S (cid:48) l ( S crit ) = − S (cid:48) l ( S (cid:63) ) = 0. In particular, f ( S ; n/N ) is monotone increasing if θ/ (1 − θ ) = n/ ( N − n ) >
1, while itattains a local minimum in ( S crit , S (cid:63) ) if n/ ( N − n ) <
1. A plot of f ( S ; θ ) is shown in Figure 2c. In eithercase, we deduce that the solution to (2.22) when 0 < κ < κ c is unique (see Figure 3a). On the other hand,when n/ ( N − n ) < κ c < κ < κ c for which (2.22) has two distinct solutions S crit < ˜ S r < S r < S (cid:63) . Indeed, this threshold can be found by demanding that µ ( S ) and κf ( S ; n/N ) intersect tangentially. In this way, we find that the threshold κ c can be written as(2.23a) κ c = κ c ( n/N ) ≡ µ ( S (cid:63)r ) f ( S (cid:63)r ; n/N ) , where S (cid:63)r is the unique solution to(2.23b) f ( S (cid:63)r ; n/N ) µ (cid:48) ( S (cid:63)r ) = f (cid:48) ( S (cid:63)r ; n/N ) µ ( S (cid:63)r ) . In Figure 3c we plot κ c − κ c as a functions of n/N where we observe that κ c > κ c with κ c − κ c → + and κ c − κ c → ∞ as n/N → . − and n/N → + respectively. Furthermore, in Figure 3b we graphicallyillustrate how multiple solutions to (2.22) arise as θ = n/N and κ are varied. We remark that the condition n/ ( N − n ) < n < N/
2, so that there are more small than large spots. The appearance oftwo distinct asymmetric patterns in this regime has a direct analogy to results obtained for the 1-D and2-D GM model in [18] and [21], respectively. The resulting bifurcation diagrams are shown in Figure 4 for n/N = 0 . , . , .
6. We summarize our results for quasi-equilibria in the following proposition. .00 0.05 0.10 0.15 0.20 0.25 S ( S ) and S ( S ) S (a) ( S ) and f ( S ; ) ( S )0.75 f(S,0.2)0.75 f(S,0.6)0.2 f(S,0.2)0.2 f(S,0.6) (b) n / N c ( n / N ) c (c)Figure 3: (a) Illustration of solutions to (2.17) as the intersection between µ ( S ) and κS . There is a unique solutionif κ < κ c ≡ µ ( S crit ) /S crit . (b) Illustration of solutions to (2.22) as the intersection between µ ( S ) and κf ( S, θ )where θ = n/N denotes the fraction of large spots in an asymmetric pattern. Note that when θ = 0 . < . κ > κ c ≈ . κ c − κ c versus n/N . Observe that κ c − κ c increases asthe fraction of large spots decreases. Proposition 2.1. (Quasi-Equilibria): Let ε → and x , . . . , x N ∈ Ω be well-separated. Then, the 3-DGM model (1.1) admits an N -spot quasi-equilibrium solution with inner asymptotics (2.24) v ∼ DV i ( ε − | x − x i | ) , u ∼ DU i ( ε − | x − x i | ) , as x → x i for each i = 1 , . . . , N where V i and U i are given by (2.5) . When | x − x i | = O (1) , the activator isexponentially small while the inhibitor is given by (2.8) . The spot strengths S iε for i = 1 , . . . , N completelydetermine the asymptotic solution and there are two distinguished limits. When D = O (1) the spot strengthssatisfy the NAS (2.14a) , which has the leading order asymptotics (2.15) . In particular, S iε ∼ S (cid:63) so all N -spot patterns are symmetric to leading order. When D = ε − D the spot strengths satisfy the NAS (2.14b) .A symmetric solution with asymptotics (2.19) where S c satisfies (2.17) always exists. Moreover, if < πN D | Ω | < κ c ≈ . , then an asymmetric pattern with n large spots of strength S r ∈ ( S crit , S (cid:63) ) and N − n small spots of strength S l ∈ (0 , S crit ) can be found by solving (2.22) for S r and calculating S l from µ ( S l ) = µ ( S r ) . If, in additionwe have n/ ( N − n ) < , then (2.22) admits two solutions on the range . ≈ κ c < πN D | Ω | < κ c ( n/N ) , where κ c ( n/N ) is found by solving the system (2.23) . As we have already remarked, in the D = D /ε regime, if D (cid:28) N -spot solution(2.19) coincides with the symmetric solution for the D = O (1) regime given by (2.15). The asymmetricsolutions predicted for the D = D /ε regime persist as D decreases and it is, therefore, natural to askwhat these solutions correspond to in the D = O (1) regime. From the small S asymptotics (2.2) we note .00 0.25 0.50 0.75 1.00 1.250.000.050.100.150.20 Source Strengths for n / N = 0.20 S c S r S l S r S l (a) Source Strengths for n / N = 0.40 S c S r S l S r S l (b) Source Strengths for n / N = 0.60 S c S r S l (c)Figure 4: Bifurcation diagram illustrating the dependence on κ of the common spot strength S c as well as theasymmetric spot strengths S r and S l or ˜ S r and ˜ S l . In (a) and (b) we have n/N < . S r and S l or ˜ S r and ˜ S l . In (c) the number of large spots exceeds that of small spots and onlyone type of asymmetric pattern is possible. that the NAS (2.14a) does admit an asymmetric solution, albeit one in which the source strengths of thesmall spots are of O ( ε ). Specifically, for a given integer n in 1 < n ≤ N we can construct a solution where(2.25) SSS ε ∼ ( S (cid:63) , . . . , S (cid:63) , ε S n +1 , , . . . , ε S N, ) T . By using the small S asymptotic expansion for µ ( S ) given in (2.2), we obtain from (2.14a) that(2.26) S i, = b πS (cid:63) n (cid:88) j =1 G ( x i , x j ) , i = n + 1 , . . . , N . We observe that in order to support N − n spots of strength O ( ε ), we require at least one spot of strength O (1). Setting D = D /ε , we use the large D asymptotics for G ( x, ξ ) in (2.9) to reduce (2.26) to(2.27) S i, ∼ bε − (cid:18) πD nS (cid:63) | Ω | (cid:19) , i = n + 1 , . . . , N . Alternatively, by taking κ (cid:28) D = D /ε regime, we conclude that S r ∼ S (cid:63) and S l ∼ b ( κnS (cid:63) /N ) . Since κn/N = 4 πD n/ | Ω | , as obtained from (2.14b), we confirm that the asymmetricpatterns in the D = D /ε regime lead to an asymmetric pattern consisting of spots of strength O (1) and O ( ε ) in the D = O (1) regime.
3. Linear Stability.
Let ( v qe , u qe ) be an N -spot quasi-equilibrium solution as constructed in §
2. Wewill analyze instabilities for quasi-equilibria that occur on O (1) time-scales. To do so, we substitute(3.1) v = v qe + e λt φ , u = u qe + e λt ψ , into (1.1) and, upon linearizing, we obtain the eigenvalue problem(3.2) ε ∆ φ − φ + 2 v qe u qe φ − v qe u qe ψ = λφ , D ∆ ψ − ψ + 2 ε − v qe φ = τ λψ , here ∂ n φ = ∂ n ψ = 0 on ∂ Ω. In the inner region near the j th spot, we introduce a local expansion in termsof the associated Legendre polynomials P ml (cos θ ) of degree l = 0 , , , . . . , and order m = 0 , , . . . , l (3.3) φ ∼ c j DP ml (cos θ ) e imϕ Φ j ( ρ ) , ψ ∼ c j DP ml (cos θ ) e imϕ Ψ j ( ρ ) , where ρ = ε − | x − x j | , and ( θ, ϕ ) ∈ (0 , π ) × [0 , π ). Suppressing subscripts for the moment, and assumingthat ε τ λ/D (cid:28)
1, we obtain the leading order inner problem(3.4a) ∆ ρ Φ − l ( l + 1) ρ Φ − Φ + 2 VU Φ − V U Ψ = λ Φ , ∆ ρ Ψ − l ( l + 1) ρ Ψ + 2 V Φ = 0 , ρ > , with the boundary conditions Φ (cid:48) (0) = Ψ (cid:48) (0) = 0, and Φ → ρ → ∞ . Here ( V, U ) satisfy the coreproblem (2.1). The behaviour of Ψ as ρ → ∞ depends on the parameter l . More specifically, we have that(3.4b) Ψ ∼ (cid:40) B ( λ, S ) + ρ − , for l = 0 ,ρ − (1 / γ l ) , for l > , as ρ → ∞ , where γ l ≡ (cid:113) + l ( l + 1) and B ( λ, S ) is a constant. Here we have normalized Ψ by fixing to unity themultiplicative factor in the decay rate in (3.4b). Next, we introduce the Green’s function G l ( ρ, ˜ ρ ) solving(3.5) ∆ ρ G l − l ( l + 1) ρ G l = − ρ − δ ( ρ − ˜ ρ ) , given by G l ( ρ, ˜ ρ ) = 12 γ l √ ρ ˜ ρ (cid:40) ( ρ/ ˜ ρ ) γ l , < ρ < ˜ ρ , ( ˜ ρ/ρ ) γ l , ρ > ˜ ρ , when l >
0. For l = 0 the same expression applies, but an arbitrary constant may be added. For conveniencewe fix this constant to be zero. In terms of this Green’s function we can solve for Ψ explicitly in (3.4a) as(3.6) Ψ = 2 (cid:90) ∞ G l ( ρ, ˜ ρ ) V ( ˜ ρ )Φ( ˜ ρ ) ˜ ρ d ˜ ρ + (cid:40) B ( λ, S ) , for l = 0 , , for l > . Upon substituting this expression into (3.4a) we obtain the nonlocal spectral problems(3.7a) M Φ = λ Φ + B ( λ, S ) V U , for l = 0 ; M l Φ = λ Φ , for l > . Here the integro-differential operator M l is defined for every l ≥ M l Φ ≡ ∆ ρ Φ − l ( l + 1) ρ Φ − Φ + 2 VU Φ − V U (cid:90) ∞ G l ( ρ, ˜ ρ ) V ( ˜ ρ )Φ( ˜ ρ ) ˜ ρ d ˜ ρ . A key difference between the l = 0 and l > B ( λ, S ) in the l = 0 equation. This unknown constant is determined by matching the far-fieldbehaviour of the inner inhibitor expansion with the outer solution. In this sense, we expect that B ( λ, S )will encapsulate global contributions from all spots, so that instabilities for the mode l = 0 are due tothe interactions between spots. In contrast, the absence of an unknown constant for instabilities for the l > l > M l . In Figure 5a we plot the numerically-computed .00 0.05 0.10 0.15 0.20 0.25S1.00.50.00.51.0 Real Part of Dominant Eigenvalues of l l 023 (a) (b)Figure 5: (a) Spectrum of the operator M l defined in (3.7b). The dashed blue line indicates the eigenvalue withsecond largest real part for l = 0. Notice that the dominant eigenvalue of M is zero when S = S crit ≈ . µ ( S ) (see Figure 1a). (b) Plot of B ( λ, S ). The dashed line black indicates thelargest positive eigenvalue of M ( S ) and also corresponds to the contour B ( λ, S ) = 0. We observe that B ( λ, S ) isboth continuous and negative for S > S crit ≈ . dominant eigenvalue of M l for l = 0 , , l = 0 for 0 < S < S (cid:63) .This spectrum is calculated from the discretization of M l obtained by truncating the infinite domain to0 < ρ < L , with L (cid:29)
1, and using a finite difference approximation for spatial derivatives combined witha trapezoidal rule discretization of the integral terms. The l = 1 mode always admits a zero eigenvalue,as this simply reflects the translation invariance of the inner problem. Indeed, these instabilities will bebriefly considered in Section 4 where we consider the slow dynamics of quasi-equilibrium spot patterns.From Figure 5a we observe that the dominant eigenvalues of M l for l = 2 , λ ) < l ). Therefore, since the modes l > linearly stable , forthe 3-D GM model there will be no peanut-splitting or spot self-replication instabilities such as observedfor the 3-D Schnakenberg model in [16]. In the next subsection we will focus on analyzing instabilitiesassociated with l = 0 mode, which involves a global coupling between localized spots. l = 0 Mode. ¿From (3.7a) we observe that λ is in thespectrum of M if and only if B ( λ, S ) = 0. Assuming that B ( λ, S ) (cid:54) = 0 we can then solve for Φ in (3.7a) as(3.8) Φ = B ( λ, S )( M − λ ) − ( V /U ) . Upon substituting (3.8) into the expression (3.6) for Ψ when l = 0, we let ρ → ∞ and use G ( ρ, ˜ ρ ) ∼ /ρ as ρ → ∞ , as obtained from (3.5), to deduce the far-field behavior(3.9) Ψ ∼ B + 2 Bρ (cid:90) ∞ V ( M − λ ) − ( V /U ) ρ d ρ , as ρ → ∞ . We compare this expression with the normalized decay condition on Ψ in (3.4b) for l = 0 to conclude that(3.10) B ( λ, S ) = 12 (cid:82) ∞ V ( M − λ ) − ( V /U ) ρ dρ . We now solve the outer problem and through a matching condition derive an algebraic equation for theeigenvalue λ . Since the interaction of spots will be important for analyzing instabilities for the l = 0 mode, e re-introduce the subscript j to label the spot. First, since ∂ ρ Ψ j ∼ − ρ − as ρ → ∞ , as obtained from(3.4b) for l = 0, an application of the divergence theorem to ∆ ρ Ψ j = − V j Φ j yields that (cid:82) ∞ V j Φ j ρ dρ =1 /
2. Next, by using v qe ∼ DV j ( ρ ) and φ ∼ c j D Φ j ( ρ ) for | x − x j | = O ( ε ) as obtained from (2.24) and (3.3),respectively, we calculate in the sense of distributions for ε → ε − v qe φ → πεD N (cid:88) j =1 c j (cid:18)(cid:90) ∞ V j Φ j ρ dρ (cid:19) δ ( x − x j ) = 4 πεD N (cid:88) j =1 c j δ ( x − x j ) . Therefore, by using this distributional limit in the equation for ψ in (3.2), the outer problem for ψ is(3.11) ∆ ψ − (1 + τ λ ) D ψ = − πεD N (cid:88) j =1 c j δ ( x − x j ) , x ∈ Ω ; ∂ n ψ = 0 , x ∈ ∂ Ω . The solution to (3.11) is represented as(3.12) ψ = 4 πεD N (cid:88) j =1 c j G λ ( x, x j ) , where G λ ( x, ξ ) is the eigenvalue-dependent Green’s function satisfying∆ G λ − (1 + τ λ ) D G λ = − δ ( x − ξ ) , x ∈ Ω ; ∂ n G λ = 0 , x ∈ ∂ Ω ,G λ ( x, ξ ) ∼ π | x − ξ | + R λ ( ξ ) + o (1) , as x → ξ . (3.13)By matching the limit as x → x i of ψ in (3.12) with the far-field behaviour ψ ∼ Dc i B ( λ, S i ) of the innersolution, as obtained from (3.9) and (3.3), we obtain the matching condition(3.14) B ( λ, S i ) c i = 4 πε (cid:18) c i R λ ( x i ) + N (cid:88) j (cid:54) = i c j G λ ( x i , x j ) (cid:19) . As similar to the construction of quasi-equilibria in §
2, there are two distinguished limits D = O (1) and D = D /ε to consider. The stability properties are shown to be significantly different in these two regimes.In the D = O (1) regime, we recall that S i ∼ S (cid:63) for i = 1 , . . . , N where µ ( S (cid:63) ) = 0. From (3.14), weconclude to leading order that B ( λ, S (cid:63) ) = 0, so that λ must be an eigenvalue of M when S = S (cid:63) . However,from Figure 5a we find that all eigenvalues of M when S = S (cid:63) satisfy Re( λ ) <
0. As such, from our leadingorder calculation we conclude that N -spot quasi-equilibria in the D = O (1) regime are all linearly stable.For the remainder of this section we focus exclusively on the D = D /ε regime. Assuming that ε | τ λ | /D (cid:28) G λ ( x, ξ ) ∼ ε − D / [(1 + τ λ ) | Ω | ] + G ( x, ξ ), where G is the NeumannGreen’s function satisfying (2.10). We substitute this limiting behavior into (3.14) and, after rewriting thethe resulting homogeneous linear system for ccc ≡ ( c , . . . , c N ) T in matrix form, we obtain(3.15) B ccc = κ τ λ E N ccc + 4 πε G ccc , where B ≡ diag( B ( λ, S ) , . . . , B ( λ, S N )) , E N ≡ N − eeeeee T . Here G is the Neumann Green’s matrix and κ ≡ πN D / | Ω | (see (2.14b)). Next, we separate the proceedinganalysis into the two cases: symmetric quasi-equilbrium patterns and asymmetric quasi-equilibria. .1.1. Stability of Symmetric Patterns in the D = D /ε Regime.
We suppose that the quasi-equilibrium solution is symmetric so that to leading order S = . . . = S N = S c where S c is found bysolving the nonlinear algebraic equation (2.17). Then, from (3.15), the leading order stability problem is(3.16) B ( λ, S c ) ccc = κ τ λ E N ccc . We first consider competition instabilities for N ≥ ccc T eee = 0 so that E N ccc = 0. Since B ( λ, S c ) = 0 from (3.16), it follows that λ must be an eigenvalue of M , defined in (3.7b), at S = S c .From Figure 5a we deduce that the pattern is unstable for S below some threshold where the dominanteigenvalue of M equals zero. In fact, this threshold is easily determined to correspond to S c = S crit , where µ (cid:48) ( S crit ) = 0, since by differentiating the core problem (2.1) with respect to S and comparing the resultingsystem with (3.4) when l = 0, we conclude that B (0 , S c ) = µ (cid:48) ( S c ). The dotted curve in Figure 5b shows thatthe zero level curve B ( λ, S c ) = 0 is such that λ > S c < S crit . As such, we conclude from (2.17) thatsymmetric N -spot quasi-equilibria are unstable to competition instabilities when κ > κ c ≡ µ ( S crit ) /S crit .For special spot configurations { x , . . . , x N } where eee is an eigenvector of G we can easily calculatea higher order correction to this instability threshold. Since G is symmetric, there are N − qqq , . . . , qqq N such that G qqq k = g k qqq k with qqq Tk eee = 0. Setting ccc = qqq k in (3.15), and using B (0 , S ) ∼ εµ (cid:48)(cid:48) ( S crit ) δ for S = S crit + εδ , we can determine the perturbed stability threshold where λ = 0associated with each eigenvector qqq k . By taking the minimum of such values, and by recalling the refinedapproximation (2.20), we obtain that N -spot symmetric quasi-equilibria are all unstable on the range(3.17) S cε < S crit + 4 πεµ (cid:48)(cid:48) ( S crit ) min k =2 ,...,N g k . Next we consider the case ccc = eee for which we find from (3.15) that, to leading order, λ satisfies(3.18) B ( λ, S c ) − κ τ λ = 0 . First, we note that λ = 0 is not a solution of (3.18) since, by using B (0 , S ) = µ (cid:48) ( S ), this would require that µ (cid:48) ( S c ) = κ , which the short argument following (2.19) demonstrates is impossible. Therefore, the ccc = eee mode does not admit a zero-eigenvalue crossing and any instability that arises must occur through a Hopfbifurcation. We will seek a leading order threshold τ = τ h ( κ ) beyond which a Hopf bifurcation is triggered.To motivate the existence of such a threshold we consider first the κ → ∞ limit for which the asymptotics(2.18) implies that S c = 1 / ( bκ ) (cid:28) S expansion (2.2) of the core solution wecalculate from (3.7b) that M Φ ∼ ∆ ρ Φ − Φ + 2 w c Φ + O ( κ − ). Then, by substituting this expression,together with the small S asymptotics (2.2) where S c ∼ /bκ (cid:28)
1, into (3.10) we can determine B ( λ, S c )when κ (cid:29)
1. Then, by using the resulting expression for B in (3.18), we obtain the following well-knownnonlocal eigenvalue problem (NLEP) corresponding to the shadow limit κ = 4 πN D / | Ω | → ∞ :(3.19) 1 + τ λ − (cid:82) ∞ w c (∆ ρ − w c − λ ) − w c ρ dρ (cid:82) ∞ w c ρ dρ = 0 . ¿From Table 1 in [19], this NLEP has a Hopf bifurcation at τ = τ ∞ h ≈ .
373 with corresponding criticaleigenvalue λ = iλ ∞ h with λ ∞ h ≈ . τ h ( κ ) for κ = O (1), we set λ = iλ h in (3.18) andseparate the resulting expression into real and imaginary parts to obtain(3.20) τ h = − Im ( B ( iλ h , S c )) λ h Re ( B ( iλ h , S c )) , | B ( iλ h , S c ) | Re ( B ( iλ h , S c )) − κ = 0 , Hopf Bifurcation Threshold leading orderhigher order( N = 1, = 0.01) (a) I Hopf Threshold Eigenvalue leading orderhigher order( N = 1, = 0.01) (b) Hopf Bifurcation Threshold leading orderhigher order( N = 1, = 0.01) (c)Figure 6: Leading order (a) Hopf bifurcation threshold τ h ( κ ) and (b) critical eigenvalue λ = iλ h for a symmetric N -spot pattern as calculated by solving (3.20) numerically. The leading order theory assumes ε | τ λ | /D (cid:28) N = 1 spotpattern centered at the origin of the unit ball with ε = 0 .
01 by solving (3.14) directly (note κ = 3 D ). In (c) we seethat although the leading order Hopf bifurcation threshold diverges as κ → κ c , going to higher order demonstratesthat a large but finite threshold persists. where S c depends on κ from (2.17). Starting with κ = 50 we solve the second equation for λ h using Newton’smethod with λ h = λ ∞ h as an initial guess. We then use the first equation to calculate τ h . Decreasing κ andusing the previous solution as an initial guess we obtain the curves τ h ( κ ) and λ h ( κ ) as shown in Figure 6.We conclude this section by noting that as seen in Figures 6a and 6c the leading order Hopf bifurcationthreshold diverges as κ → κ + c , where κ c = µ ( S crit ) /S crit . This is a direct consequence of the assumptionthat ε | τ λ | /D (cid:28) τ gets increasingly large. Indeed, by using the series expansionin (3.12)–(3.14) of [12] for the reduced wave Green’s function in the sphere, we can solve (3.14) directly usingNewton’s method for an N = 1 spot configuration centered at the origin of the unit ball. Fixing ε = 0 . κ ≤ κ c . Moreover, it hints at an ε dependent rescaling of τ in the region κ ≤ κ c for which acounterpart to (3.16) may be derived. While we do not undertake this rescaling in this paper we remarkthat for 2-D spot patterns this rescaling led to the discovery in [15] of an anomalous scaling law for theHopf threshold. D = D /ε Regime.
When the N -spot pattern consistsof n large spots of strength S = . . . = S n = S r and N − n small spots of strength S n +1 = . . . = S N = S l ,the leading order linear stability is characterized by the blocked matrix system(3.21) (cid:18) B ( λ, S r ) I n B ( λ, S l ) I N − n (cid:19) ccc = κ τ λ E N ccc , where I m denotes the m × m identity matrix. In particular, an asymmetric quasi-equilibrium solution islinearly unstable if this system admits any nontrivial modes, ccc , for which λ has a positive real part. Wewill show that asymmetric patterns are always unstable by explicitly constructing unstable modes. irst, we assume that 1 ≤ n < N − ccc to be a mode satisfying(3.22) c = · · · = c n = 0 , c n +1 + · · · + c N = 0 . Note that this mode describes competition among the N − n small spots of strength S l . For such a mode,(3.21) reduces to the single equation B ( λ, S l ) = 0, which implies that λ must be an eigenvalue of M at S = S l . However, since S l < S crit , we deduce from Figure 5a that there exists a real and positive λ for M at S = S l . As such, any mode ccc satisfying (3.22) is linearly unstable.We must consider the n = N − n ≥ N − n , for which n = N − ccc of the form(3.23) c = . . . = c n = c r , c n +1 = . . . = c N = c l . Then, (3.21) reduces to the system of two equations (cid:16) B ( λ, S r ) − κ τλ nN (cid:17) c r − κ τλ ( N − n ) N c l = 0 , − κ τλ nN c r + (cid:16) B ( λ, S l ) − κ τλ ( N − n ) N (cid:17) c l = 0 , which admits a nontrivial solution if and only if the determinant of this 2 × F ( λ ) ≡ B ( λ, S l ) B ( λ, S r ) − κ τ λ (cid:18) nN B ( λ, S l ) + ( N − n ) N B ( λ, S r ) (cid:19) = 0 , has a solution λ >
0. To establish this, we first differentiate µ ( S r ) = µ ( S l ) with respect to S r to obtain theidentity µ (cid:48) ( S l ) S (cid:48) l ( S r ) = µ (cid:48) ( S r ). Combining this result with B (0 , S ) = µ (cid:48) ( S ) we calculate that(3.25) F (0) = µ (cid:48) ( S l ) (cid:20) µ (cid:48) ( S r ) − κ ( N − n ) N (cid:18) n ( N − n ) + dS l dS r (cid:19)(cid:21) . Using µ (cid:48) ( S l ) > µ (cid:48) ( S r ) < S (cid:48) l ( S r ) > − n/ ( N − n ) ≥
1, we immediately deduce that F (0) <
0. Next, we let λ > M when S = S l (see Figure 5a) so that B ( λ , S l ) = 0. Then, from (3.24) we obtain(3.26) F ( λ ) = − κ τ λ ( N − n ) N B ( λ , S r ) . However, since M at S = S r > S crit has no positive eigenvalues (see Figure 5a), we deduce that B ( λ, S r )is of one sign for λ ≥ B (0 , S r ) = µ (cid:48) ( S r ) < B showing both its continuity and negativity for all λ > S > S crit ). Therefore, we have F ( λ ) > F ( λ ) = 0 has apositive solution. We summarize our leading order linear stability results in the following proposition: Proposition 3.1. (Linear Stability): Let ε (cid:28) and assume that t (cid:28) O ( ε − ) . When D = O (1) , the N -spot symmetric pattern from Proposition 2.1 is linearly stable. If D = ε − D then the symmetric N -spot pattern from Proposition 2.1 is linearly stable with respect to zero-eigenvalue crossing instabilities if κ < κ c ≡ µ ( S crit ) /S crit ≈ . and is unstable otherwise. Moreover, it is stable with respect to Hopfinstabilities on the range κ > κ c if τ < τ h ( κ ) where τ h ( κ ) is plotted in Figure 6a. Finally, every asymmetric N -spot pattern in the D = ε − D regime is always linearly unstable. . Slow Spot Dynamics. A wide variety of singularly perturbed RD systems are known to exhibit slowdynamics of multi-spot solutions in 2-D domains (cf. [9], [3], [13], [17]). In this section we derive a systemof ODE’s which characterize the motion of the spot locations x , . . . , x N for the 3-D GM model on a slowtime scale. Since the only N -spot patterns that may be stable on an O (1) time scale are (to leading order)symmetric we find that the ODE system reduces to a gradient flow. We remark that both the derivationand final ODE system are closely related to those in [16] for the 3-D Schnakenberg model.The derivation of slow spot dynamics hinges on establishing a solvability condition for higher orderterms in the asymptotic expansion in the inner region near each spot. As a result, we begin by collectinghigher order expansions of the limiting behaviour as | x − x i | → G ( x, x j ) and G ( x, x j ) that satisfy (2.7) and (2.10), respectively. In particular, we calculate that(4.1a) G ( x i + εy, x j ) ∼ (cid:40) G ( x i , x j ) + εy · ∇ G ( x i , x j ) , i (cid:54) = j , περ + R ( x i ) + εy · ∇ R ( x i ; x i ) , i = j , as | x − x i | → , where ρ = | y | and ∇ R ( x i ; x i ) ≡ ∇ x R ( x ; x ) | x = x . Likewise, for the Neumann Green’s function, we have(4.1b) G ( x i + εy, x j ) ∼ D ε | Ω | + (cid:40) G ( x i , x j ) + εy · ∇ G ( x i , x j ) , i (cid:54) = j , περ + R ( x i ) + εy · ∇ R ( x i ; x i ) , i = j , as | x − x i | → , where ∇ again denotes the gradient with respect to the first argument. We next extend the asymptoticconstruction of quasi-equilibrium patterns in § x i = x i ( σ ) where σ = ε t . For x near x i we introduce the two term inner expansion(4.2) v ∼ DV i ∼ D ( V iε ( ρ ) + ε V i ( y ) + · · · ) , u ∼ DU i ∼ D (cid:0) U iε ( ρ ) + ε U i ( y ) + · · · (cid:1) , where we note the leading order terms are V iε ( ρ ) ≡ V ( ρ, S iε ) and U iε ( ρ ) ≡ U ( ρ, S iε ). By using the chainrule we calculate ∂ t V i = − ε x (cid:48) i ( σ ) · ∇ y V i and ∂ t U i = − ε x (cid:48) i ( σ ) · ∇ y U i . In this way, upon substituting (4.2)into (1.1) we collect the O ( ε ) terms to obtain that V i and U i satisfy(4.3a) L iε WWW i ≡ ∆ y WWW i + Q iε WWW i = − fff iε , y ∈ R , where(4.3b) WWW i ≡ (cid:18) V i U i (cid:19) , fff iε ≡ (cid:18) ρ − V (cid:48) iε ( ρ ) x (cid:48) i ( σ ) · y − D − U iε (cid:19) , Q iε ≡ (cid:18) − U − iε V iε − U − iε V iε V iε (cid:19) . It remains to determine the appropriate limiting behaviour as ρ → ∞ . From the first row of Q iε , weconclude that V i → ρ → ∞ . However, the limiting behaviour of U i must be establishedby matching with the outer solution. To perform this matching, we first use the distributional limit ε − v −→ πεD N (cid:88) j =1 S jε δ ( x − x j ) + 2 ε D N (cid:88) j =1 (cid:18)(cid:90) R V jε V j dy (cid:19) δ ( x − x j ) , where the localization at each x , . . . , x N eliminates all cross terms. We then update (2.8) to include the O ( ε ) correction term. This leads to the refined approximation for the outer solution(4.4) u ∼ πεD N (cid:88) j =1 S jε G ( x ; x j ) + 2 ε D N (cid:88) j =1 (cid:18)(cid:90) R V jε V j dy (cid:19) G ( x ; x j ) . e observe that the leading order matching condition is immediately satisfied in both the D = O (1) and the D = D /ε regimes. To establish the higher order matching condition we distinguish between the D = O (1)and D = ε − D regimes and use the higher order expansions of the Green’s functions as given by (4.1a)and (4.1b). In this way, in the D = O (1) regime we obtain the far-field behaviour as | y | → ∞ given by(4.5) U i ∼ πρ (cid:90) R V iε V i dy + y · b iε , b iε π ≡ S iε ∇ R ( x i ; x i ) + (cid:88) j (cid:54) = i S jε ∇ G ( x i , x j ) . Similarly, in the D = D /ε regime we obtain the following far-field matching condition as | y | → ∞ :(4.6) U i ∼ πρ (cid:90) R V iε V i dy + 2 D | Ω | N (cid:88) j =1 (cid:90) R V jε V j dy + y · b iε , b iε π ≡ S iε ∇ R ( x i ; x i ) + (cid:88) j (cid:54) = i S jε ∇ G ( x i , x j ) . In both cases, our calculations below will show that only b iε and b iε affect the slow spot dynamics.To characterize slow spot dynamics we calculate x (cid:48) i ( σ ) by formulating an appropriate solvability con-dition. We observe for each k = 1 , , ∂ y k WWW iε where WWW iε ≡ ( V iε , U iε ) T satisfy thehomogeneous problem L iε ∂ y k WWW iε = 0. Therefore, the null-space of the adjoint operator L (cid:63)iε is at leastthree-dimensional. Assuming it is exactly three dimensional we consider the three linearly independentsolutions ΨΨΨ ik ≡ y k PPP i ( ρ ) /ρ to the homogeneous adjoint problem, where each PPP i ( ρ ) = ( P i ( ρ ) , P i ( ρ ) T solves(4.7) ∆ ρ PPP i − ρ PPP i + Q Tiε
PPP i = 0 , ρ > PPP (cid:48) i (0) = (cid:18) (cid:19) ; with Q Tiε −→ (cid:18) − (cid:19) as ρ → ∞ . Owing to this limiting far-field behavior of the matrix Q Tiε , we immediately deduce that P i = O ( ρ − )and that P i decays exponentially to zero as ρ → ∞ . Enforcing, for convenience, the point normalizationcondition P i ∼ ρ − as ρ → ∞ , we find that (4.7) admits a unique solution. We use each ΨΨΨ ik to impose asolvability condition by multiplying (4.3a) by ΨΨΨ Tik and integrating over the ball, B ρ , centered at the originand of radius ρ with ρ (cid:29)
1. Then, by using the divergence theorem, we calculate(4.8) lim ρ →∞ (cid:90) B ρ (cid:18) ΨΨΨ
Tik L i WWW i − WWW i L (cid:63)i ΨΨΨ ik (cid:19) dy = lim ρ →∞ (cid:90) ∂B ρ (cid:18) ΨΨΨ
Tik ∂ ρ WWW i − WWW Ti ∂ ρ ΨΨΨ ik (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) ρ = ρ ρ d Θ , where Θ denotes the solid angle for the unit sphere.To proceed, we use the following simple identities given in terms of the Kronecker symbol δ kl :(4.9) (cid:90) B ρ y k f ( ρ ) dy = 0 , (cid:90) B ρ y k y l f ( ρ ) dy = δ kl π (cid:90) ρ ρ f ( ρ ) dρ , for l, k = 1 , , . Since L (cid:63)i ΨΨΨ ik = 0, we can use (4.3a) and (4.9) to calculate the left-hand side of (4.8) aslim ρ →∞ (cid:90) B ρ ΨΨΨ
Tik L i WWW i dy = lim ρ →∞ (cid:18) − (cid:88) l =1 x (cid:48) il ( σ ) (cid:90) B ρ y k y l P i ( ρ ) V (cid:48) iε ( ρ ) ρ dy + 1 D (cid:90) B ρ y k P i ( ρ ) U iε ( ρ ) ρ dy (cid:19) = − π x (cid:48) ik ( σ ) (cid:90) ∞ P i ( ρ ) V (cid:48) iε ( ρ ) ρ dρ . (4.10) .00 0.05 0.10 0.15 0.20 S ( S ) Figure 7: Plot of the numerically-computed multiplier γ ( S ) as defined in the slow gradient flow dynamics (4.14). Next, in calculating the right-hand side of (4.8) by using the far-field behavior (4.5) and (4.6), we observethat only b iε and b iε terms play a role in the limit. In particular, in the D = O (1) regime we calculate interms of the components of b iεl of the vector b iε , as given in (4.5), thatlim ρ →∞ (cid:90) ∂B ρ ΨΨΨ
Tik ∂ ρ WWW i (cid:12)(cid:12) ρ = ρ ρ d Θ = lim ρ →∞ (cid:88) l =1 b iεl (cid:90) ∂B ρ y k y l ρ d Θ = 4 π b iεk , lim ρ →∞ (cid:90) ∂B ρ WWW Ti ∂ ρ ΨΨΨ ik (cid:12)(cid:12) ρ = ρ ρ d Θ = − ρ →∞ (cid:88) l =1 b iεl (cid:90) ∂B ρ y k y l ρ d Θ = − π b iεk . (4.11)¿From (4.8), (4.10), and (4.11), we conclude for the D = O (1) regime that(4.12) x (cid:48) ik ( σ ) = − γ ( S iε ) b iεk , where γ ( S iε ) ≡ (cid:90) ∞ P i ( ρ ) V (cid:48) i ( ρ, S iε ) ρ dρ , which holds for each component k = 1 , , i = 1 , . . . , N . From symmetry considerations wesee that the constant contribution to the far-field behaviour, as given by the first term in (4.5), is eliminatedwhen integrated over the boundary. In an identical way, we can determine x (cid:48) ik for the D = D /(cid:15) regime.In summary, in terms of the gradients of the Green’s functions and γ iε ≡ γ ( S iε ), as defined in (4.12), weobtain the following vector-valued ODE systems for the two distinguished ranges of D : dx i dσ = − πγ iε (cid:18) S iε ∇ R ( x i ; x i ) + (cid:80) j (cid:54) = i S jε ∇ G ( x i , x j ) (cid:19) , for D = O (1) , (cid:18) S iε ∇ R ( x i ; x i ) + (cid:80) j (cid:54) = i S jε ∇ G ( x i , x j ) (cid:19) , for D = D /ε . (4.13)Since only the symmetric N -spot configurations can be stable on an O (1) time scale (see Proposition3.1), it suffices to consider the ODE systems in (4.13) when S iε = S (cid:63) + O ( ε ) in the D = O (1) regime andwhen S iε = S c + O ( ε ), where S c solves (2.17), in the D = ε − D regime. In particular, we find that toleading order, where the O ( ε ) corrections to the source strengths are neglected, the ODE systems in (4.13)can be reduced to the gradient flow dynamics(4.14a) dx i dσ = − πSγ ( S ) ∇ x i H ( x , . . . , x N ) , with γ ( S ) = (cid:90) ∞ P ( ρ ) V ( ρ, S ) ρ dρ , .5 1.0 1.5 2.0 D h ( D ) Hopf Bifurcation Threshold for N = 1 (a) v ( , t ) D =0.60 =0.69=1.15 0 10 20 300.02.55.07.5 D =1.00 =0.49=0.810 10 20 30 t v ( , t ) D =1.40 =0.42=0.70 0 10 20 30 t D =1.80 =0.38=0.64 (b)Figure 8: (a) Leading order Hopf bifurcation threshold for a one-spot pattern. (b) Plots of the spot height v (0 , t )from numerically solving (1.1) using FlexPDE6 [6] in the unit ball with ε = 0 .
05 at the indicated τ and D values. where S = S (cid:63) or S = S c depending on whether D = O (1) or D = ε − D , respectively. In (4.14) the discreteenergy H , which depends on the instantaneous spot locations, is defined by(4.14b) H ( x , . . . , x N ) ≡ (cid:40)(cid:80) Ni =1 R ( x i ) + 2 (cid:80) Ni =1 (cid:80) j>i G ( x i , x j ) , for D = O (1) , (cid:80) Ni =1 R ( x i ) + 2 (cid:80) Ni =1 (cid:80) j>i G ( x i , x j ) , for D = ε − D . In accounting for the factor of two between (4.14) and (4.13), we used the reciprocity relations for theGreen’s functions. In this leading order ODE system, the integral γ ( S ) is the same for each spot, since P ( ρ ) is computed numerically from the homogeneous adjoint problem (4.7) using the core solution V ( ρ, S )and U ( ρ, S ) to calculate the matrix Q Tiε in (4.7). In Figure 7 we plot the numerically-computed γ ( S ), wherewe note that γ ( S ) >
0. Since γ ( S ) >
0, local minima of H are linearly stable equilibria for (4.14).We remark that this gradient flow system (4.14) differs from that derived in [16] for the 3-D Schnaken-berg model only through the constant γ ( S ). Since this parameter affects only the time-scale of the slowdynamics we deduce that the equilibrium configurations and stability properties for the ODE dynamics willbe identical to those of the Schnakenberg model. As such, we do not analyze (4.14) further and instead referto [16] for more detailed numerical investigations. Finally we note that the methods employed here andin [16] should be applicable to other 3-D RD systems yielding similar limiting ODE systems for slow spotdynamics. The similarity between slow dynamics for a variety of RD systems in 2-D has been previouslyobserved and a general asymptotic framework has been pursued in [13] for the dynamics on the sphere.
5. Numerical Examples.
In this section we use FlexPDE6 [6] to numerically solve (1.1) when Ω is theunit ball. In particular, we illustrate the emergence of Hopf and competition instabilities, as predicted in § D = D /ε regimes.We begin by considering a single spot centered at the origin in the unit ball, for the D = ε − D regime.Since no competition instabilities occur for a single spot solution, we focus exclusively on the onset of Hopf
100 200 300 400 500 t v ( x , t ) , v ( x , t ) Activator Spike Heights D (a) (b)Figure 9: (a) Plots of the spot heights (solid and dashed lines) in a two-spot symmetric pattern at the indicatedvalues of D . Results were obtained by using FlexPDE6 [6] to solve (1.1) in the unit ball with ε = 0 .
05 and τ = 0 . v ( x, t ) for D = 0 . v = 0 . , . , . instabilities as τ is increased. In Figure 8a we plot the Hopf bifurcation threshold obtained from our linearstability theory, and indicate several sample points below and above the threshold. Using FlexPDE6 [6],we performed full numerical simulations of (1.1) in the unit ball with ε = 0 .
05 and parameters D and τ corresponding to the labeled points in Figure 8a. The resulting activator height at the origin, v (0 , t ),computed from FlexPDE6 is shown in Figure 8b for these indicated parameter values. We observe thatthere is good agreement with the onset of Hopf bifurcations as predicted by our linear stability theory.Next, we illustrate the onset of a competition instability by considering a symmetric two-spot configu-rations with spots centered at ( ± . , ,
0) in the unit ball and with τ = 0 . ε = 0 .
05. The critical value of κ c ≈ . | Ω | = 4 π/ D ≈ . / (3 N ) = 0 . D = 0 .
09 and D = 0 . D = 0 .
6. The Weak Interaction Limit D = O ( ε ) . In § D = O (1) and D = O ( ε − ) regimes that N -spot quasi-equilibria are not susceptible to locally non-radially symmetricinstabilities. Here we consider the weak-interaction regime D = D ε , where we numerically determinethat locally non-radially symmetric instabilities of a localized spot are possible. First, we let ξ ∈ Ω satisfydist( ξ, ∂ Ω) (cid:29) O ( ε ) and we introduce the local coordinates x = ξ + εy and the inner variables v ∼ ε V ( ρ )and u ∼ ε U ( ρ ). With this scaling, and with D = D ε , the steady-state problem for (1.1) becomes(6.1) ∆ ρ V − V + U − V = 0 , D ∆ ρ U − U + V = 0 , ρ = | y | > . D V (0) versus D (a) V Dominant Eigenvalue ( l ) l 0234 (b)Figure 10: (a) Bifurcation diagram for solutions to the core problem (6.1) in the D = ε D regime. (b) Dominanteigenvalue of the linearization of the core problem for each mode l = 0 , , ,
4, as computed numerically from (6.5).
For this core problem, we impose the boundary conditions V ρ (0) = U ρ (0) = 0 and ( V, U ) → ρ → ∞ . Unlike the D = O (1) and D = O ( ε − ) regimes, u and v are both exponentially small in the outerregion. Therefore, for any well-separated configuration x , . . . , x N , the inner problems near each spot centreare essentially identical and independent. In Figure 10a we plot V (0) versus D obtained by numericallysolving (6.1). From this figure, we observe that for all D (cid:39) . V and inhibitor U decay exponentially there are only exponentially weakinteractions between individual spots. As a result, it suffices to consider only the linear stability of the coreproblem (6.1). Upon linearizing (1.1) about the core solution we obtain the eigenvalue problem(6.2) ∆ ρ Φ − l ( l + 1) ρ Φ − Φ + 2 VU Φ − V U Ψ = λ Φ , D ∆ ρ Ψ − l ( l + 1) ρ Ψ − Ψ + 2 V Φ = 0 , for each l ≥ (cid:48) (0) = Ψ (cid:48) (0) = 0 and (Φ , Ψ) → ρ → ∞ .We reduce (6.2) to a single nonlocal equation by noting that the Green’s function G l ( ρ, ρ ) satisfying(6.3) D ∆ ρ G l − l ( l + 1) ρ G l − G l = − δ ( ρ − ρ ) ρ , is given explicitly by(6.4) G l ( ρ, ρ ) = 1 D √ ρ ρ (cid:40) I l +1 / ( ρ/ √ D ) K l +1 / ( ρ / √ D ) , ρ < ρ ,I l +1 / ( ρ / √ D ) K l +1 / ( ρ/ √ D ) , ρ > ρ , where I n ( · ) and K n ( · ) are the n th order modified Bessel Functions of the first and second kind, respectively.As a result, by proceeding as in § M l Φ = λ Φ where(6.5) M l Φ ≡ ∆ ρ Φ − l ( l + 1) ρ Φ − Φ + 2 VU Φ − V U (cid:90) ∞ G l ( ρ, ˜ ρ ) V ( ˜ ρ ) Φ( ˜ ρ ) ˜ ρ d ˜ ρ . In Figure 10b we plot the real part of the largest numerically-computed eigenvalue of M l as a function of V (0) for l = 0 , , ,
4. From this figure, we observe that the entire lower solution branch in the V (0) versus igure 11: Snapshots of FlexPDE6 [6] simulation of (1.1) in the unit ball with ε = 0 . D = 16 ε , and τ = 1 andwith initial condition given by a single spot solution in the weak interaction limit calculated from (6.1) with V (0) = 5.The snapshots show contour plots of the activator v ( x, t ) at different times where for each spot the outermost, middle,and innermost contours correspond to values of 0 . . , and 0 .
012 respectively. Note that the asymptotic theorypredicts a maximum peak height of v ∼ ε V (0) ≈ . D bifurcation diagram in Figure 10a is unstable. However, in contrast to the D = O (1) and D = O ( ε − )regimes, we observe from the orange curve in Figure 10b for the l = 2 mode that when D = ε D there isa range of D values for which a peanut-splitting instability is the only unstable mode.In previous studies of singularly perturbed RD systems supporting peanut-splitting instabilities it hastypically been observed that such linear instabilities trigger nonlinear spot self-replication events (cf. [16],[9], [13], and [3]). Recently, in [23] it has been shown using a hybrid analytical-numerical approach thatpeanut-splitting instabilities are subcritical for the 2-D Schnakenberg, Gray-Scott, and Brusselator models,although the corresponding issue in a 3-D setting is still an open problem. Our numerical computationsbelow suggest that peanut-splitting instabilities for the 3-D GM model in the D = ε D regime are alsosubcritical. Moreover, due to the exponentially small interaction between spots, we also hypothesize thata peanut-splitting instability triggers a cascade of spot self-replication events that will eventually pack thedomain with identical spots. To explore this proposed behaviour we use FlexPDE6 [6] to numerically solve(1.1) in the unit ball with parameters τ = 1, ε = 0 .
05 and D = 16 ε , where the initial condition is a singlespot pattern given asymptotically by the solution to (6.1) with V (0) = 5. From the bifurcation and stabilityplots of Figure 10 our parameter values and initial conditions are in the range where a peanut-splittinginstability occurs. In Figure 11 we plot contours of the solution v ( x, t ) at various times. We observe thatthe peanut-splitting instability triggered between t = 20 and t = 60 leads to a self-replication processresulting in two identical spots at t = 110. The peanut-splitting instability is triggered for each of thesetwo spots and this process repeats, leading to a packing of the domain with N = 8 identical spots.
7. General Gierer-Meinhardt Exponents.
Next, we briefly consider the generalized GM model(7.1) v t = ε ∆ v − v + u − q v p , τ u t = D ∆ u − u + ε − u − s v m , x ∈ Ω ; ∂ n v = ∂ n u = 0 , x ∈ ∂ Ω , where the GM exponents ( p, q, m, s ) satisfy the usual conditions p > q > m > s ≥
0, and ζ ≡ mq/ ( p − − ( s + 1) > ifferences as compared to the prototypical set ( p, q, m, s ) = (2 , , ,
0) considered in this paper, many ofthe qualitative properties resulting from the properties of µ ( S ) in Conjecture 2.1, such as the existence ofsymmetric quasi-equilibrium spot patterns in the D = O (1) regime, remain unchanged.Suppose that (7.1) has an N -spot quasi-equilibrium solution with well-separated spots. Near the i th spot we introduce the inner expansion v ∼ D α V i ( y ), u ∼ D β U i ( y ), and y = ε − ( x − x i ), where∆ V i − V i + D ( p − α − qβ U − qi V pi = 0 , ∆ U i − ε D − U i = − D mα − ( s +1) β − U − si V mi , y ∈ R . Choosing α and β such that ( p − α − qβ = 0 and mα − ( s + 1) β = 1 we obtain α = ν/ζ , β = 1 /ζ , ν = q/ ( p − , with which the inner expansion takes the form v ∼ D ν/ζ V ( ρ ; S iε ) and u ∼ D /ζ U ( ρ ; S iε ), where V ( ρ ; S )and U ( ρ ; S ) are radially-symmetric solutions to the D -independent core problem∆ ρ V − V + U − q V p = 0 , ∆ ρ U = − U − s V m , ρ > , (7.2a) ∂ ρ V (0) = ∂ ρ U (0) = 0 , V −→ U ∼ µ ( S ) + S/ρ , ρ → ∞ . (7.2b)By using the divergence theorem, we obtain the identity S = (cid:82) ∞ U − s V m ρ dρ > µ ( S ) retains severalof the key qualitative properties of the exponent set ( p, q, m, s ) = (2 , , ,
0) observed in § § §
3. To path-follow solutions, we proceed as in § S (cid:28)
1. For S (cid:28)
1, we use the identity S = (cid:82) ∞ U − s V m ρ dρ > S scaling law, and from this we readily calculate that(7.3) V ( ρ ; S ) ∼ (cid:18) Sb (cid:19) νζ +1 w c ( ρ ) , U ( ρ ; S ) ∼ (cid:18) Sb (cid:19) ζ +1 , µ ( S ) ∼ (cid:18) Sb (cid:19) ζ +1 , b ≡ (cid:90) ∞ w mc ρ dρ , where w c > ρ w c − w c + w pc = 0 , ρ > ∂ ρ w c (0) = 0 , w c → ρ → ∞ . With this approximate solution for S (cid:28)
1, we proceed as in § µ ( S ) in (7.2) for differentGM exponent sets by path-following in S . In Figure 12b we plot µ ( S ) when ( p, q, m, s ) = ( p, , p, p = 2 , ,
4, while a similar plot is shown in Figure 12a for other typical exponent sets in [19]. Foreach set considered, we find that µ ( S ) satisfies the properties in Conjecture 2.1. Finally, to obtain theNAS for the source strengths we proceed as in § D with D /ζ in (2.8). Then, by using the matching condition u ∼ D /ζ ( µ ( S jε ) + S jε ε/ | x − x j | ) as x → x j , for each j = 1 , . . . , N , we conclude that the NAS (2.14) still holdsfor a general GM exponent set provided that µ ( S ) is now defined by the generalized core problem (7.2).
8. Discussion.
We have used the method of matched asymptotic expansions to construct and studythe linear stability of N -spot quasi-equilibrium solutions to the 3-D GM model (1.1) in the limit of anasymptotically small activator diffusivity ε (cid:28)
1. Our key contribution has been the identification of twodistinguished regimes for the inhibitor diffusivity, the D = O (1) and D = O ( ε − ) regimes, for which weconstructed N -spot quasi-equilibrium patterns, analyzed their linear stability, and derived an ODE system .0 0.1 0.2 0.3 0.4 0.5 S ( S ) for GM Exponents ( p , q , m , s ) (p,q,m,s)(2,1,3,0)(3,2,3,1)(3,2,2,0)(2,1,2,0) (a) S ( S ) for GM Exponents ( p , 1, p , 0) p 234 (b)Figure 12: Left panel: Plot of µ ( S ), computed from the generalized GM core problem (7.2), for the indicated exponentsets ( p, q, m, s ). Right panel: µ ( S ) for exponent sets ( p, , p,
0) with p = 2 , ,
4. For each set, there is a unique S = S (cid:63) for which µ ( S (cid:63) ) = 0. The properties of µ ( S ) in Conjecture 2.1 for the protypical set (2 , , ,
0) still hold. governing their slow spot dynamics. We determined that in the D = O (1) regime all N -spot patterns are,to leading order in ε , symmetric and linearly stable on an O (1) time scale. On the other hand, in the D = O ( ε − ) regime we found the existence of both symmetric and asymmetric N -spot patterns. However,we demonstrated that all asymmetric patterns are unstable on an O (1) time scale, while for the symmetricpatterns we calculated Hopf and competition instability thresholds. These GM results are related to thosein [16] for the 3-D singularly perturbed Schnakenberg model, with one of the key new features being theemergence of two distinguished limits, and in particular the existence of localized solutions in the D = O (1)regime for the GM model. For D = O (1), concentration behavior for the Schnakenberg model as ε → D = O ( ε ), where weused a hybrid analytical-numerical approach to calculate steady-state solutions and determine their linearstability properties. In this small D regime we found that spot patterns are susceptible to peanut-splittinginstabilities. Finally, using FlexPDE6 we illustrated how the weak-interaction between spots together withthe peanut-splitting instability leads to a cascade of spot self-replication events.We conclude by highlighting directions for future work and open problems. First, although we haveprovided numerical evidence for the properties of µ ( S ) highlighted in Conjecture 2.1, a rigorous proofremains to be found. In particular, we believe that it would be significant contribution to rigorously provethe existence and uniqueness of the ground state solution to the core problem (2.1), which we numericallycalculated when S = S (cid:63) . A broader and more ambitious future direction is to characterize the reactionkinetics F ( V, U ) and G ( V, U ) for which the core problem(8.1) ∆ ρ V + F ( V, U ) = 0 , ∆ ρ U + G ( V, U ) = 0 , in ρ > , admits a radially-symmetric ground state solution for which V → U = O ( ρ − ) as ρ → ∞ . The existence of such a ground state plays a key role in determining the regimes of D forwhich localized solutions can be constructed. For example, in the study of the 3-D singularly perturbed .0 0.1 0.2 0.3S0.0000.0050.0100.0150.0200.0250.030 ( S ) for the GMS Model (a) ( S ) for the S/GS Model (b) ( S ) for the Brusselator Model f 0.20.40.60.8 (c)Figure 13: Plots of the far-field constant behaviour for the (a) Gierer-Meinhardt with saturation, (b) Schnakenbergor Gray-Scott, and (c) Brusselator models. See Table 1 for the explicit form of the kinetics F ( v, u ) and G ( v, u ) foreach model. A zero-crossing of µ ( S ) at some S >
Schnakenberg model it was found that the core problem does not admit such a solution and as a resultlocalized spot solutions could not be constructed in the D = O (1) regime (cf. [16]). To further motivatesuch an investigation of (8.1) we extend our numerical method from § µ ( S ) for the core problems associated with the GM model with saturation (GMS),the Schnakenberg/Gray-Scott (S/GS) model, and the Brusselator (B) model (see Table 1 for more details).Note that for the GMS model we can find values of S (cid:63) such that µ ( S (cid:63) ) = 0, but such a zero-crossing doesnot appear to occur for the (S/GS) and (B) models. As a consequence, for these three specific RD systems,localized spot patterns in the D = O (1) regime should only occur for the GMS model. Additionally,understanding how properties of µ ( S ), such as convexity and positiveness, are inherited from the reactionkinetics would be a significant contribution. In this direction, it would be interesting to try extend therigorous numerics methodology of [1] to try to establish Conjecture 2.1. Acknowledgments.
Daniel Gomez was supported by an NSERC Doctoral Fellowship. Michael Wardand Juncheng Wei gratefully acknowledge the support of the NSERC Discovery Grant Program.
REFERENCES
RD Model F ( V, U ) G ( V, U ) Decay behaviorGierer-Meinhardt with Saturation (GMS) − V + V U (1+ κU ) V U ∼ µ ( S ) + S/ρ
Schnakenberg or Gray-Scott (S/GS) − V + V U − V U U ∼ µ ( S ) − S/ρ
Brusselator (B) − V + f V U V − V U U ∼ µ ( S ) − S/ρ
Table 1: Core problems and far-field inhibitor behaviour for some common RD systems.
1] I. Bal´azs, J. van den Berg, J. Courtois, J. Dud´as, J. P. Lessard, A. V¨or¨os-Kiss, J. F. Williams, and X. Y. Yin. Computer-assisted proofs for radially symmetric solutions of PDEs.
J. Comput. Dynamics , 5(1 & 2):61–80, 2018.[2] C. N. Chen, Y. S. Choi, Y. Hu, and X. Ren. Higher dimensional bubble profiles in a sharp interface limit of theFitzHugh-Nagumo system.
SIAM J. Math. Anal. , 50(5):5072–5095, 2018.[3] W. Chen and M. J. Ward. The stability and dynamics of localized spot patterns in the two-dimensional Gray-Scottmodel.
SIAM J. Appl. Dyn. Sys. , 10(2):582–666, 2011.[4] A. Doelman and H. Van der Ploeg. Homoclinic stripe patterns.
SIAM J. Appl. Dyn. Sys. , 1(1):65–104, 2002.[5] S. I. Ei and S. Y. Tzeng. Spike solutions for a mass conservation reaction-diffusion system.
DCDS Series A , 40(6):3357–3374, 2020.[6] P. FlexPDE. Solutions inc. , 2015.[7] A. Gierer and H. Meinhardt. A theory of biological pattern formation.
Kybernetik , 12(1):30–39, Dec 1972.[8] T. Kolokolnikov, M. J. Ward, and J. Wei. Pulse-splitting for some reaction-diffusion systems in one-space dimension.
Studies App. Math. , 114(2):115–165., 2005.[9] T. Kolokolnikov, M. J. Ward, and J. Wei. Spot self-replication and dynamics for the Schnakenberg model in a two-dimensional domain.
J. Nonlinear Sci. , 19(1):1–56, 2009.[10] M. Leda, V. K. Vanag, and I. R. Epstein. Instabilities of a three-dimensional localized spot.
Phys. Rev. E , 80:066204,2009.[11] Y. Nishiura.
Far-from Equilibrium dynamics: Translations of mathematical monographs , volume 209. AMS Publications,Providence, Rhode Island, 2002.[12] R. Straube and M. J. Ward. Intraceulluar signalling gradients arising from multiple compartments: A matched asymptoticexpansion approach.
SIAM J. Appl. Math. , 70(1):248–269, 2009.[13] P. H. Trinh and M. J. Ward. The dynamics of localized spot patterns for reaction-diffusion systems on the sphere.
Nonlinearity , 29(3):766–806, 2016.[14] J. Tzou. private communication.[15] J. C. Tzou, M. J. Ward, and J. C. Wei. Anomalous scaling of Hopf bifurcation thresholds for the stability of localizedspot patterns for reaction-diffusion systems in two dimensions.
SIAM J. Appl. Dyn. Syst. , 17(1):982–1022, 2018.[16] J. C. Tzou, S. Xie, T. Kolokolnikov, and M. J. Ward. The stability and slow dynamics of localized spot patterns for the3-D Schnakenberg reaction-diffusion model.
SIAM J. Appl. Dyn. Syst. , 16(1):294–336, 2017.[17] M. J. Ward. Spots, traps, and patches: Asymptotic analysis of localized solutions to some linear and nonlinear diffusiveprocesses.
Nonlinearity , 31(8):R189 (53), 2018.[18] M. J. Ward and J. Wei. Asymmetric spike patterns for the one-dimensional Gierer-Meinhardt model: Equilibria andstability.
European J. Appl. Math. , 13(3):283–320, 2002.[19] M. J. Ward and J. Wei. Hopf bifurcation of spike solutions for the shadow Gierer-Meinhardt model.
European J. Appl.Math. , 14(6):677–711, 2003.[20] J. Wei. Existence and stability of spikes for the Gierer-Meinhardt system. In M. Chipot, editor,
Handbook of DifferentialEquations, Stationary Partial Differential Equations , volume 5, pages 489–581. Elsevier, 2008.[21] J. Wei and M. Winter. Spikes for the two-dimensional Gierer-Meinhardt system: The weak coupling case.
J. NonlinearSci. , 11(6):415–458, 2001.[22] J. Wei and M. Winter.
Mathematial aspects of pattern formation in biological systems , volume 189. Applied MathematicalSciences Series, Springer, 2014.[23] T. Wong and M. J. Ward. Weakly nonlinear analysis of peanut-shaped deformations for localized spots of singularlyperturbed reaction-diffusion systems. to appear, SIAM J. Appl. Dyn. Sys. , 2020., 2020.