Spot Patterns in the 2-D Schnakenberg Model with Localized Heterogeneities
SSpot Patterns in the 2-D Schnakenberg Model with LocalizedHeterogeneities
Tony Wong ∗ Michael J. Ward ∗ †
September 18, 2020
Abstract
A hybrid asymptotic-numerical theory is developed to analyze the effect of different types of localized heterogeneitieson the existence, linear stability, and slow dynamics of localized spot patterns for the two-component Schnakenbergreaction-diffusion model in a 2-D domain. Two distinct types of localized heterogeneities are considered: a stronglocalized perturbation of a spatially uniform feed rate and the effect of removing a small hole in the domain, throughwhich the chemical species can leak out. Our hybrid theory reveals a wide range of novel phenomena such as, saddle-node bifurcations for quasi-equilibrium spot patterns that otherwise would not occur for a homogeneous medium, a newtype of spot solution pinned at the concentration point of the feed rate, spot self-replication behavior leading to thecreation of more than two new spots, and the existence of a creation-annihilation attractor with at most three spots.Depending on the type of localized heterogeneity introduced, localized spots are either repelled or attracted towards thelocalized defect on asymptotically long time scales. Results for slow spot dynamics and detailed predictions of variousinstabilities of quasi-equilibrium spot patterns, all based on our hybrid asymptotic-numerical theory, are illustratedand confirmed through extensive full PDE numerical simulations. keywords —
Pattern formation, reaction-diffusion systems, spots, localized heterogeneites, pinning.
Localized spot patterns, in which a solution component becomes spatially localized near certain time-varying discretepoints within a bounded multi-dimensional domain, is a well-known “far-from-equilibrium” spatial pattern that occursfor certain two-component reaction-diffusion (RD) systems in the singular limit of a large diffusivity ratio. This class oflocalized pattern is observed in many chemical and biological systems, such as the chlorine-dioxide-malonic acid reaction[9], the ferrocyanide-iodate-sulphite reaction [25, 26], and the initiation of plant root hair cells mediated by the planthormone auxin [2], among others (see [36] and [14] for surveys). In a spatially homogeneous 2-D medium, and for variousspecific RD systems, the slow dynamical behavior of quasi-equilibrium spot patterns, together with their various typesof bifurcations that trigger a range of different instabilities of the pattern such as spot-annihilation, spot-replication, andtemporal oscillations of the spot amplitude, have been well-studied [7, 17, 31, 32, 33, 40, 41, 42, 47]. The primary focusof this article is to investigate, for one prototypical RD system, how certain spatial heterogeneities in the model affect thedynamics and instabilities of quasi-equilibrium spot patterns, and lead to new dynamical phenomena that would otherwisenot occur in a medium free of defects. For tractability of our analysis, and as we describe below, we will focus only oncertain types of spatially localized heterogeneities.There is a growing literature, primarily in a 1-D setting, of analyzing the effect of a spatial heterogeneity in eitherthe diffusivity or reaction kinetics on pattern-formation behavior for two-component RD systems with regards to bothsmall amplitude patterns (see [29], [30], [23] and the references therein) and for localized far-from-equilibrium spike-typepatterns (cf. [2], [3], [4], [5], [10], [18], [19], [20], [21], [33], [39], [43], [44], [45]). In particular, the analysis in [21] and [43]has revealed that a precursor gradient in the reaction kinetics can lead to the existence of stable asymmetric spike patternsfor the Gierer-Meinhardt (GM) model, which would otherwise not occur in a homogeneous medium. A precursor field forthe GM model can also lead to stable steady-states consisting of spike clusters near critical points of the precursor. In [22] ∗ Dept. of Mathematics, Univ. of British Columbia, Vancouver, B.C., Canada. † corresponding author, [email protected] a r X i v : . [ n li n . PS ] S e p t was shown that a different type of smooth heterogeneity in the 1-D GM model can lead to the formation of a creation-annihilation attractor, which consists of periodically repeating cycles of spike formation, propagation, and annihilationagainst a domain boundary. In the limit of a large number of spikes that are confined by a spatial heterogeneity, a meanfield equation for the spike density was derived in [20] and [18] for the 1-D GM and Schnakenberg models, respectively,and in [19] for 2-D spot clusters for the GM model. For the 1-D Schnakenburg model, the mean field limiting equationin [18] revealed the existence of a creation-annihilation attractor in which spikes undergo self-replication in the interior ofthe spike cluster, while other spikes are annihilated at the edges of the cluster. For an extended Klausmeir RD model ofspatial ecology, coarsening and pinning behavior of 1-D spike patterns with various spatial and temporal heterogeneitieswere analyzed in [3] and [4]. Spike dynamics and pinning effects for a 1-D RD model where the nonlinearities have smallspatial support, as is typical for catalytic reactions, was studied in [10].For a different class of localized pattern consisting of either a propagating pulse-type or a transition-layer solution,there has been much effort at analyzing the effect of a small step function barrier on pulse propagation properties forthe three-component Fitzhugh-Nagumo RD system (cf. [6], [35], [34], [48]) and for two-component bistable RD systems(cf. [11], [27]). For this type of jump-type spatial heterogeneity, the focus has been to analytically determine parameterranges where a 1-D propagating pulse will either be reflected, transmitted, or pinned by the barrier. In a 2-D setting, [28]provides a numerical study of similar propagation and collision properties for a single localized spot in the presence of a1-D step-function line barrier.In contrast to the simpler 1-D case, there are relatively few analytical studies of the effect of spatial heterogeneities forRD systems in higher spatial dimensions. For a generalized Schnakenberg-type RD system modeling the initiation of roothair profusion in plant cells, a spatially inhomogeneous auxin gradient in a 2-D rectangular domain was shown to lead tothe alignment of localized spots in the direction of the gradient (cf. [2]). In a 2-D rectangular domain, it was shown in[16] for a generalized Klausmeir RD system, modeling vegetation patterns in semi-arid environments, that an anisotropicdiffusivity can stabilize a localized stripe pattern to transverse perturbations. With isotropic diffusion, the homoclinicstripe would be unstable to either breakup into spots or zigzag deformations. Spot-pinning behavior for RD systems onclosed manifolds of non-constant curvature, which can be viewed as intrinsic spatial heterogeneities, has been analyzed forthe Schnakenberg model in [13]. In [33], which is most closely related to our study, the role of Robin boundary conditionsand boundary fluxes on the slow dynamics and instabilities of quasi-equilibrium spot patterns for the Brusselator RDmodel was analyzed.The goal of this paper is to analyze the effect of various types of localized heterogeneities for the singularly perturbedSchnakenberg model in a bounded 2-D domain Ω, formulated as v t = ε ∆ v − v + uv , τ u t = D ∆ u + a − ε − uv , x ∈ Ω ; ∂ n v = ∂ n u = 0 , x ∈ ∂ Ω , (1.1)where 0 < ε (cid:28)
1, while
D > τ > O (1) constants. One heterogeneity will be introduced through strong, butlocal, perturbations in the feed rate a = a ( x ), which characterizes the amount of material that is introduced from thesubstrate. Another localized heterogeneity that we will consider is to analyze the effect of perturbing (1.1) by removinga small hole in the domain, which thereby allows leakage of the chemical species out of the domain.For these types of localized heterogeneities we will extend the hybrid asymptotic-numerical framework of [17] and [33]to analyze the existence, linear stability, and slow dynamics of quasi-equilibrium spot patterns. Depending on the type oflocalized heterogeneity introduced, spot patterns are either repelled or attracted towards the defect on a long time scaleof order O ( ε − ). By formulating and analyzing various spectral problems arising from the linear stability analysis forinstabilities of the quasi-equilibrium pattern on short O (1) time-scales, we will show how peanut-splitting and competitioninstabilities that trigger either spot self-replication or spot-annihilation events, respectively, are affected by the type oflocalized heterogeneity. For a localized heterogeneity where there is a slowly moving localized source of feed in the domain,we will combine our linear stability theory for quasi-equilibrium spot patterns with our derived ODE system for slow spotdynamics to construct a novel attractor consisting of spot-replication and spot-annihilation events that has a maximumof three spots in the domain at any time.To both illustrate and validate our asymptotic theory for various types of localized heterogeneities, throughout thispaper we will compare our predictions for slow spot dynamics and spot amplitude instabilities with full PDE simulationsof (1.1). The full simulations are done using the open source finite element software FEniCS [1], which automates themesh generation and finite element assembly from user inputs. Our choice of node sizes range, approximately, from 20000to 80000. For time-stepping we used either a Backward-Euler time stepping scheme or a BDF-2 (backward differentiationformula), the latter of which is preferable for computing spot amplitude temporal oscillations due to a Hopf bifurcation.2he outline of this paper is as follows. In § §
2. In § § O ( ε ) from the domain, while imposing a homogeneous Dirichlet condition onthe boundary of the small hole. With this type of localized heterogeneity, which allows for the possibility of both chemicalspecies to leak out of the 2-D domain, we show that a one-spot quasi-equilibrium solution exists only if the feed rate islarge enough or if the spot is sufficiently far enough away from the hole. More specifically, in contrast to the scenario fora homogeneous medium, spot quasi-equilibria are shown to exhibit a novel saddle-node bifurcation structure in the unitdisk in terms of either the feed rate or the distance from the hole. In addition, from the derivation of a modified systemfor slow spot dynamics, we show that localized spots are dynamically repelled from the small hole (see Fig. 5 and Fig. 6below). Moreover, we show that a significantly larger threshold value for the feed rate, as compared to the case for ahomogeneous medium, is needed to initiate spot self-replication events (see Fig. 7 below). Although perforated domainshave been well-studied in the context of narrow capture mean first passage time problems for Brownian particles (see [24]and the references therein), the effect of a perforated domain, resulting in an open reaction-diffusion system (cf. [33]), onlocalized pattern formation problems has to our knowledge not been analyzed previously. The analysis in § § § §
4, in § § N + 1 spot patterns that consist of N unpinned spots together with anadditional spot centered at the concentration point of the feed rate. By deriving a globally coupled eigenvalue problem,we formulate a criterion for which this pattern undergoes a competition instability, triggering a spot-annihilation event,that is due to a zero-eigenvalue crossing of the linearization. Finally, in § § In § N -spot patterns for (1.1) andto characterize heir slow dynamics. In the limit ε → N -spot quasi-equilibrium solution for (1.1) with spots centered at x , . . . , x N .We assume that the spots are well-separated in the sense that | x i − x j | = O (1) for i (cid:54) = j and dist( x j , ∂ Ω) = O (1) for j = 1 , . . . , N . We assume that the quasi-equilibrium pattern is linearly stable on O (1) time intervals.In the inner region near the j th spot, we let x j = x j ( σ ) where σ = ε t is the slow time scale (cf. [17]). We introducethe inner variables v = √ DV j ( y ) , u = 1 √ D U j ( y ) , where y ≡ ε − ( x − x j ( σ )) , and ρ = | y | , (2.1a)together with the inner expansion V j = V j ( ρ ) + V j + · · · , U j = U j ( ρ ) + U j + · · · . (2.1b)Upon substituting (2.1) into (1.1), we collect powers of ε to obtain, at leading order, the radially symmetric core problem ∆ ρ V j − V j + U j V j = 0 , ∆ ρ U j − U j V j = 0 , < ρ < ∞ , (2.2a) V (cid:48) j (0) = U (cid:48) j (0) = 0 ; V j → , U j ∼ S j log ρ + χ ( S j ) , as ρ → ∞ , (2.2b)where ∆ ρ ≡ ∂ ρρ + ρ − ∂ ρ and S j is called the spot source strength. At next order, we find that v ≡ ( V j , U j ) T satisfies∆ y v + M j v = f j , y ∈ R , (2.3a)where ∆ y denote derivatives in y , and where we have defined M j ≡ (cid:18) − U j V j V j − U j V j − V j (cid:19) , f j ≡ (cid:18) − V (cid:48) j ( e φ · ˙x j ) (cid:19) . (2.3b)Here e φ ≡ (cos θ, sin θ ) T and ˙x j ≡ d x j /dσ . For (2.3a) we can impose that V j → | y | → ∞ . However, the far-fieldcondition of U j is determined only after asymptotic matching to an appropriate outer solution.In Fig. 1 we plot the numerical solution V j ( ρ ) to the core problem (2.2) for several source strengths. We show thatthere is a unique spot height V j (0) and a unique χ for each source strength.To derive the outer problem for u , we integrate the equation for U j in (2.2a) over 0 < ρ < ∞ to obtain the identity S j = (cid:90) ∞ U j V j ρ dρ. (2.4)4 (a) V j with different S (b) Spot height versus S (c) χ ( S ) Figure 1: Left panel: The core solution component V j ( ρ ), computed from (2.2), for S = 0 . , . , . , . V j (0) given approximately by 0 . , . , . , . , . S , showing a monotone increase on S < S (cid:63) ≈ .
83 and a unique spot height for each
S > χ versus S .Then, in the limit ε →
0, we use (2.4) to obtain, in the sense of distributions, that ε − uv → N (cid:88) i =1 (cid:18)(cid:90) R ( D − / U i )( D / V i ) d y (cid:19) δ ( x − x i ) = 2 π √ D N (cid:88) i =1 S i δ ( x − x i ) . (2.5)By using this distributional limit in (1.1), we obtain that the outer problem, defined away from the spots, is∆ u + aD − π √ D N (cid:88) i =1 S i δ ( x − x i ) = 0 in Ω , ∂ n u = 0 on ∂ Ω . (2.6)By integrating (2.6) over Ω and using the divergence theorem, we get N (cid:88) i =1 S i = a | Ω | π √ D ≡ p a , (2.7)where | Ω | denotes the area of Ω. The solution to (2.6) is represented as u ( x ) = − π √ D N (cid:88) i =1 S i G ( x ; x i ) + ¯ u , (2.8)where ¯ u is an undetermined additive constant and G ( x ; z ) is the unique Neumann Green’s function satisfying∆ G = 1 | Ω | − δ ( x − z ) in Ω , ∂ n G = 0 on ∂ Ω ; (cid:90) Ω G ( x ; z ) d x = 0 ,G ∼ − π log | x − z | + R ( z ; z ) + ∇ x R ( x , z ) | x = z · ( x − z ) + O ( | x − z | ) , as x → z , (2.9)where R ( z ; z ) is called the regular part of G .To determine a nonlinear algebraic system for the source strengths, and a DAE system for slow spot dynamics, wemust match the near-field behavior as x → x j of the outer solution (2.8) to the far-field behavior of the two-term innersolution, which is given from (2.2b) and (2.1) by u ∼ √ D ( S j log ρ + χ ( S j ) + εU j + · · · ) = 1 √ D (cid:18) S j log | x − x j | + S j ν + χ ( S j ) (cid:19) + εU j √ D + · · · , (2.10)5s ρ = ε − | x − x j | → ∞ , where we have defined ν ≡ − / log ε . Then, by Taylor-expanding (2.8) as x → x j , and replacing x − x j = ε y , we obtain after some algebra that u ∼ S j √ D log | x − x j | − π √ D S j R j,j + N (cid:88) i (cid:54) = j S i G j,i + ¯ u − ε √ D ( βββ j · y ) + · · · , (2.11)where βββ j is defined by βββ j ≡ π S j ∇ x R j,j + N (cid:88) i (cid:54) = j S i ∇ x G j,i . (2.12)Here we have labeled R j,j ≡ R ( x j ; x j ), G j,i ≡ G ( x j ; x i ), ∇ x R j,j ≡ ∇ x R ( x ; x j ) | x = x j , and ∇ x G j,i ≡ ∇ x G ( x ; x i ) | x = x j .By comparing the leading terms in (2.10) and (2.11), and recalling (2.7), we obtain in matrix form that s + 2 πν G s + ν χ = ν ¯ u √ D e , e T s = p a ≡ a | Ω | π √ D , (2.13a)where we have defined s , χ , e , and the Neumann Green’s matrix G ∈ R N × N for x , . . . , x N by s = ( S , . . . , S N ) T , χ ≡ ( χ ( S ) , . . . , χ ( S N )) T , e ≡ (1 , . . . , T ∈ R N , ( G ) i j = (cid:40) R j,j if i = j ,G i,j if i (cid:54) = j . (2.13b)By left-multiplying (2.13a) by e T , and by using e T s = p a , we can isolate ¯ u . Then, by substituting ¯ u back into (2.13a),we can decouple (2.13a) to obtain that s satisfies the nonlinear algebraic system (NAS) s + 2 πν ( I − E ) G s + ν ( I − E ) χ = p a N e , with ¯ u = p a + 2 πν e T G s + ν e T χ ν √ DN . (2.14)Here
E ≡ N − ee T ∈ R N × N and I ∈ R N × N is the identity matrix.To determine the slow dynamics, we proceed to next order and match the O ( ε ) terms in (2.10) and (2.11). This yieldsthat the far-field behavior for the solution U j to (2.3) is U j ∼ − βββ j · y , as ρ = | y | → ∞ , (2.15)where βββ j is defined in (2.12). The ODE system for the spot locations is obtained by imposing a solvability condition onthe solution to (2.3) with far-field behavior (2.15). By differentiating the core problem (2.2) with respect to y and y , itfollows that the homogeneous problem ∆ y Φ + M j Φ = has two non-trivial solutions. As such, there are two solutionsto the corresponding homogeneous adjoint problem ∆ y Ψ + M Tj Ψ = . These two solutions have the form Ψ c = P ( ρ ) cos φ , Ψ s = P ( ρ ) sin φ , (2.16)where P ( ρ ) ≡ ( P ( ρ ) , P ( ρ )) T is the normalized nontrivial solution to∆ ρ P − ρ P + M Tj P = , with P ∼ (cid:18) /ρ (cid:19) as ρ → ∞ . (2.17)In [17] the solvability condition is obtained by multiplying (2.3a) by Ψ c and Ψ s and applying Green’s second identityon a sufficiently large circle where the far-field conditions (2.15) and (2.17) are imposed. This yields the following ODEsystem for x j ( σ ), for j = 1 , . . . , N , with σ = ε t , that characterize the slow spot dynamics: d x j dσ = − γ ( S j ) βββ j , γ ( S j ) ≡ − (cid:82) ∞ P V (cid:48) j ρ dρ . (2.18)Here βββ j is defined in (2.12), while S , . . . , S N satisfies the NAS (2.14). The plot in Fig. 3 of [17] of the numericallycomputed γ ( S j ) shows that γ ( S j ) >
0. 6 .2 Linear stability analysis
The slow spot dynamics (2.18) is valid only when the quasi-equilibrium solution is linearly stable on O (1) time-scales. Inthis subsection we analyze the linear stability of the quasi-equilibrium solution, denoted by v = v e and u = u e . To do so,we introduce the perturbation v = v e + e λt φ and u = u e + e λt η into (1.1), and upon linearizing we obtain ε ∆ φ − φ + 2 u e v e φ + v e η = λφ , D ∆ η + a − ε − (cid:0) u e v e φ + v e η (cid:1) = τ λη , in Ω , (2.19)with ∂ n φ = ∂ n η = 0 on ∂ Ω.In the inner region near the j th spot we have to leading order that v e ∼ √ D V j ( ρ ) and u e ∼ U j ( ρ ) / √ D , where ρ = ε − | x − x j | , with V j and U j satisfying the core problem (2.2). By letting φ ∼ e imθ Φ j ( ρ ) and η ∼ e imθ N j ( ρ ) /D , forthe integer angular mode m ≥
0, we obtain the leading order inner eigenvalue problem∆ ρ Φ j − m ρ Φ j − Φ j + 2 U j V j Φ j + V j N j = λ Φ j , ∆ ρ N j − m ρ N j − U j V j Φ j − V j N j = 0 . (2.20)We first consider non-radially symmetric perturbations for which m >
0. The case m = 1 corresponds, trivially, to thetranslation mode (Φ j , N j ) = ( U (cid:48) j , V (cid:48) j ) with λ = 0. For angular modes with m ≥
2, we impose Φ j → ρ →
0. In addition, owing to the m N j /ρ term in (2.20), we impose the far-field decay condition N j ∼ O ( ρ − m ) as ρ → ∞ . The eigenvalue λ max in (2.20) with the largest real part has been numerically calculated in [17] for a range of S j . For each m ≥
2, it was found that λ max is real and negative (positive) when S j < Σ m ( S j > Σ m ) (see Fig. 4 of[17]). Moreover, as shown numerically in [17], the ordering principle Σ < Σ < . . . holds for the stability thresholdsfor non-radially symmetric perturbations. As such, the mode m = 2, referred as to the peanut-splitting mode, is thefirst to lose stability when S j is increased. The critical threshold for this mode is Σ ≈ . m = 0 is derived by globally coupling local problems near each spot.To derive this globally coupled eigenvalue problem (GCEP), we set m = 0 in (2.20) and impose that N j ∼ c j log ρ as ρ → ∞ , where c j is an unknown constant. We then writeΦ j = c j ˜Φ j , N j = c j ˜ N j , (2.21)so as to obtain from (2.20) that∆ ρ ˜Φ j − ˜Φ j + 2 U j V j ˜Φ j + V j ˜ N j = λ ˜Φ j , ∆ ρ ˜ N j − U j V j ˜Φ j − V j ˜ N j = 0 , ρ > , (2.22a)˜Φ (cid:48) j (0) = ˜ N (cid:48) j (0) = 0 ; ˜Φ j → , ˜ N j ∼ log ρ + ˜ B ( S j ; λ ) + o (1) , as ρ → ∞ , (2.22b)where ˜ B ( S j ; λ ) must be calculated numerically from (2.22). However, by differentiating the core problem (2.2) with respectto S j , we observe that ∂ S V j and ∂ S U j satisfy (2.22) when λ = 0. As a result, we have the identity that ˜ B ( S j ; 0) = χ (cid:48) ( S j ).By integrating the ˜ N j equation in (2.22), and using (2.21), we obtain the identity c j = (cid:90) ∞ (cid:0) U j V j Φ j + V j N j (cid:1) ρ dρ . (2.23)Then, in the limit ε →
0, we use (2.23) to derive, in the sense of distributions, that ε − (cid:0) u e v e φ + v e η (cid:1) → π N (cid:88) i =1 c i δ ( x − x i ) . (2.24)We use (2.24), together with the asymptotic matching condition η ∼ c j ˜ N j /D , where ˜ N j has the far-field behavior as ρ → ∞ in (2.22), to obtain the following outer problem for η , defined away from the spots:∆ η − τ λD η = 2 πD N (cid:88) i =1 c i δ ( x − x i ) in Ω , ∂ n η = 0 in ∂ Ω (2.25a) η ∼ c j D (cid:20) log | x − x j | + 1 ν + ˜ B ( S j ; λ ) (cid:21) , as x → x j , for j = 1 , . . . , N . (2.25b)7or λ (cid:54) = 0, we represent the solution to (2.25a) as η = − πD N (cid:88) i =1 c i G λ ( x ; x i ) , (2.26)where G λ ( x , z ) is the eigenvalue-dependent Green’s function satisfying∆ G λ − τ λD G λ = − δ ( x − z ) in Ω , ∂ n G λ = 0 on ∂ Ω , (2.27a) G λ ∼ − π log | x − z | + R λ ( z ; z ) + o (1) as x → z . (2.27b)By Taylor-expanding η in (2.26) as x → x j , and then equating the resulting expression with (2.25b), we conclude that c j + 2 πν c j R λ j,j + N (cid:88) i (cid:54) = j c i G λj,i + ν c j ˜ B ( S j ; λ ) = 0 , j = 1 , . . . , N , (2.28)where R λ j,j ≡ R λ ( x j ; x j ) and G λj,i ≡ G λ ( x j ; x i ). In matrix form, and with c ≡ ( c , . . . , c N ) T , (2.28) is equivalent to M ( λ ) c = , where M ( λ ) ≡ I + 2 πν G λ + ν ˜ B . (2.29a)Here I ∈ R N × N is the identity matrix, while the symmetric Green’s matrix G λ and diagonal matrix ˜ B are defined by (cid:0) G λ (cid:1) ij ≡ (cid:40) R λ j,j if i = j ,G λi,j if i (cid:54) = j , (cid:0) ˜ B (cid:1) ij = (cid:40) ˜ B ( S j ; λ ) if i = j , i (cid:54) = j . (2.29b)The homogeneous matrix system (2.29a) for c , referred to as the GCEP, has a nontrivial solution if and only ifdet M ( λ ) = 0 . (2.29c)A discrete root λ to (2.29c) for which Re( λ ) > c characterizes the small-scale perturbation of the spot amplitudes.In this way, the linear stability properties associated with locally radially symmetric perturbations near the spots isreduced to the problem of determining the number N of roots of det M ( λ ) = 0 in the right-half Re( λ ) > C ζ that consists of the semi-circle | λ | = ζ >
0, for − π/ ≤ arg λ ≤ π/
2, and the imaginary segment { λ = iλ I : − ζ ≤ λ I ≤ ζ } . However, since M ( λ ) is undefined at λ = 0, we need to first find the behavior of det M ( λ ) as λ → λ → G λ = D | Ω | τ λ ee T + Q , where Q ≡ G + O ( τ λ ) as λ → . (2.30)Here G is the Neumann Green’s matrix and e ≡ (1 , . . . , T . Since ee T is a rank one matrix, we substitute (2.30) into(2.29a) for M and, by using the well-known matrix determinant lemma, we obtaindet M ( λ ) = det (cid:16) I + 2 πν Q + ν ˜ B (cid:17) + 2 πνD | Ω | τ λ (cid:104) e T adj (cid:16) I + 2 πν Q + ν ˜ B (cid:17) e (cid:105) , (2.31)where adj( A ) denotes the adjugate of a matrix A . From (2.31) it follows that det M ( λ ) has a simple pole at λ = 0. As aresult, it is convenient to introduce the function T ( λ ) defined by T ( λ ) ≡ λ det M ( λ ), which has a removable singularityat λ = 0 and has the same number N of zeroes in Re( λ ) > M ( λ ). The argument principle for T yields that N = P + 12 π lim ζ →∞ [arg T ( λ )] C ζ , where T ( λ ) ≡ λ det M ( λ ) . (2.32)8ere P is the number of poles of T ( λ ) in Re( λ ) >
0. Since G λ is analytic in Re( λ ) >
0, any such pole can only arise fromthe diagonal matrix ˜ B as defined by (2.29b). However, from a numerical computation of the local problem (2.22), we findthat ˜ B is analytic in Re( λ ) > P = 0 in (2.32). To determine N in the examples below, the change [arg T ( λ )] C ζ in the argument of T over the contour C ζ is computed numerically.Next, we study zero-eigenvalue crossings. Since ˜ B ( S j ; 0) = χ (cid:48) ( S j ), the outer problem (2.25a) when λ = 0 becomes∆ η = 2 πD N (cid:88) i =1 c i δ ( x − x i ) in Ω , ∂ n η = 0 in ∂ Ω (2.33a) η ∼ c j D (cid:20) log | x − x j | + 1 ν + χ (cid:48) ( S j ) (cid:21) , as x → x j , j = 1 , . . . , N . (2.33b)From the divergence theorem we conclude that (cid:80) Ni =1 c i = 0. With this constraint, we represent the solution to (2.33a) interms of the Neumann Green’s function G of (2.9) as η = − πD N (cid:88) i =1 c i G ( x ; x i ) + ¯ η , (2.34)where ¯ η is an additive constant to be determined. Then, we Taylor-expand (2.34) as x → x j by recalling the local behaviorof G in (2.9). By equating the resulting expression with the required singularity condition (2.33b), we obtain a matrixsystem for c = ( c , . . . , c N ) T and ¯ η of the form( I + 2 πν G + ν ˜ B ) c = ν ¯ η e , e T c = 0 , (2.35)where G is the Neumann Green’s matrix and where the diagonal matrix ˜ B is defined by (cid:16) ˜ B (cid:17) i j = (cid:40) χ (cid:48) ( S j ) if i = j , i (cid:54) = j , −→ ˜ B ≡ diag ( χ (cid:48) ( S ) , . . . , χ (cid:48) ( S N )) . (2.36)By left-multiplying (2.35) by e T , and using e T c = 0, we find that ¯ η = N − (cid:16) π e T G c + e T ˜ B c (cid:17) . By substituting thisexpression back into the first equation in (2.35) we derive that M c = , where M ≡ I + 2 πν ( I − E ) G + ν ( I − E ) ˜ B , (2.37)where E ≡ N − ee T . We conclude that a zero-eigenvalue crossing associated with locally radially symmetric perturbationsnear the spots occurs if and only if det M = 0. Since the corresponding nontrivial eigenmode c satisfies e T c = 0, it isreferred to as a competition mode as it locally preserves the sum of all the spot amplitudes.Finally, we relate the zero-eigenvalue crossing condition det M = 0 to the local solvability of the NAS (2.14). Suppose,for a particular fixed parameter set, that s = s e is a non-degenerate solution to the NAS (2.14) in the sense that theJacobian matrix of the NAS is invertible at s = s e . Upon introducing the perturbation s = s e + c into (2.14) where | c | (cid:28)
1, we linearize the NAS to readily determine that this Jacobian matrix is in fact the GCEP matrix M of (2.37),in which ˜ B ≡ diag ( χ (cid:48) ( S e ) , . . . , χ (cid:48) ( S Ne )). As a result, if s e is a non-degenerate solution to the NAS (2.14) we must havedet M (cid:54) = 0, and so λ = 0 is not an eigenvalue of the GCEP. Therefore, it is only at a bifurcation point of the NAS (2.14)where a zero-eigenvalue crossing of the GCEP can occur. This correspondence is summarized asdet M = 0 ⇐⇒ s e is at a bifurcation point of the NAS (2.14) . (2.38) N -spot ring pattern An N -spot ring pattern is a pattern of N equally-spaced spots located on a ring of radius r , with 0 < r <
1, that isconcentric within the unit disk Ω. For j = 1 , . . . , N , the locations of the spots on the ring can be taken as x j = r e θ j , e θ j ≡ (cos θ j , sin θ j ) T , θ j ≡ π ( j − N , for j = 1 , . . . , N . (2.39)9or a ring pattern, the symmetric Neumann Green’s matrix G is also circulant, and so it has the eigenvector e = (1 , . . . , T .As a result, the NAS (2.14) admits a symmetric solution where the spots have the common source strength S j = S c ≡ p a /N , for j = 1 , . . . , N , where p a is given in (2.7).As shown in Appendix A, with x j = r ( σ ) e θ j the spot dynamics given in (2.18) can be reduced to the scalar ODE dr dσ = γ ( S c ) S c (cid:32) N − r − N r N − − r N − N r (cid:33) , with S c = a | Ω | πN √ D , (2.40)for the ring radius, where σ = ε t . On 0 < r <
1, this ODE (2.40) has a globally stable equilibrium point r e , given bythe unique root to N − N − r = r N − r N . (2.41)From § N -spot ring pattern is linear stable to locally non-radially symmetric perturbations near the spots onlywhen S c < Σ ≈ . S c = a | Ω | / [2 πN √ D ] with | Ω | = π . In terms of the feed rate a , this stability conditionholds when a < a f ≡ √ DN ≈ . DN .Next, we study the linear stability associated with radially-symmetric perturbations near the spots. For a ring pattern,the GCEP (2.29a) becomes M c = , where M = (1 + ν ˜ B c ) I + 2 πν G λ . (2.42)Here ˜ B c ≡ ˜ B ( S c ; λ ) is to be calculated from the inner problem (2.22) with S j = S c . Owing to the cyclic structure of thering pattern, the symmetric Green’s matrix G λ is also a circulant matrix and, as a result, it has the matrix spectrum (seeAppendix B) G λ e = ω e , G λ q j = ω j q j , j = 2 , . . . , N ; e T q j = 0 , q Ti q j = 0 , i (cid:54) = j , (2.43)where q j for j = 2 , . . . , N are given in (B.1b). The matrix eigenvalues ω j are given in terms of the first row of G λ by(B.1b), while the entries in G λ can be evaluated numerically from the infinite series result in (A.5) of Appendix A for theeigenvalue-dependent Green’s function of (2.27).Since M in (2.42) represents an update to G λ by a multiple of the identity matrix, the eigenspace of M is the same as G λ . As a result, we simply substitute c = e and c j = q j into (2.42) to obtain the root finding problems F j = 0, whichare defined in terms of ω j in (2.43) by F j ≡ ν ˜ B ( S c ; λ ) + 2 πν ω j , j = 1 , . . . , N . (2.44)We refer to c = e and c j = q j , for j = 2 , . . . , N , as the synchronous mode and asynchronous modes, respectively. We begin by analyzing the zero-eigenvalue crossing in the GCEP for an N -spot ring pattern. The criterion (2.37) becomes M c = , where M = (1 + νχ (cid:48) ( S c )) I + 2 πν ( I − E ) G − νχ (cid:48) ( S c ) E . (2.45)The matrix M shares the same eigenspace as the symmetric and circulant matrix G , and so has eigenvectors e , q , . . . , q N as in (2.43). Since G e = σ e , we use E e = e to calculate M e = e . Therefore, the synchronous mode c = e can never bea nullvector for M . In contrast, with c = q j for j = 2 , . . . , N , we use E q j = 0 to obtain that M c = if and only if1 + νχ (cid:48) ( S c ) + 2 πνσ j = 0 , where G q j = σ j q j , j = 2 , . . . , N . (2.46)From (2.38), det M = 0 can only occur at a bifurcation point for the NAS (2.14).As an example, we investigate competition instabilities for a two-spot equilibrium ring pattern in the unit disk with ε = 0 . D = 1, and ring radius r = 0 . q = (1 , − T is the only competition mode,we use S c = a/ [2 N √ D ] = a/ j = N = 2 to write (2.46) as a nonlinear algebraic equation in the feed rate a . Thisequation is solved numerically to obtain the competition threshold a comp ≈ .
45. To interpret this bifurcation point, wedetermine asymmetric branches of two-spot ring patterns from the NAS (2.14) with N = 2. Labeling S and S as the10ource strengths of the two spots, we set S = p a − S in (2.14) with N = 2 to obtain a scalar nonlinear algebraic equationfor S given by(1 + 2 πν ( R − G )) S + ν χ ( S ) − χ ( p a − S )] = p a πν ( R − G )) , where p a ≡ a . (2.47)By solving (2.47) numerically, in Fig. 2a we show the bifurcation structure of S versus a . The symmetric branchcorresponds to the common source strength S = S ≡ S c = a/
4. It undergoes a pitchfork bifurcation at a = a comp ≈ . S (cid:54) = S exist for a > a comp . In the same figure, we superimpose PDE simulation data computed from(1.1) with a slowly decreasing feed rate a = max(4 . − . t , . a comp , only one spotsurvives and there is a fast transition to the one-spot branch where S = a/ (a) Pitchfork bifurcation: two-spot ring pattern (b) Imperfect bifurcation: perturbed two-spot pattern Figure 2: Left panel: Source strength S versus a for a two-spot equilibrium ring pattern with ε = 0 . D = 1, andring radius r = 0 . S data interpolated from the PDE simulation with slowly decreasing feed rate a = max(4 . − . t , . S jumps to the one-spot branch where S = a/ x = (0 . ,
0) and x = ( − . , a is swept with a = max(4 , − ( ε/ t ) below the saddle-node point, only onespot survives. The sum of squares of the source strength jumps to the one-spot branch S = a / x = (0 . ,
0) and x = ( − . ,
0) in the unit disk with ε = 0 .
02 and D = 1. Throughnumerical continuation of the NAS (2.14) with bifurcation parameter a using COCO [8], in Fig. 2b we observe twoisolated branches of S + S , with one branch having a saddle-node bifurcation at a ≈ . a = max(4 , − ( ε/ t ) with ε = 0 .
02, in Fig. 2b we show that as a sweepsbelow the saddle-node point for two-spot quasi-equilibria, one spot gets annihilated while the remaining spot jumps tothe stable one-spot branch.Next, we illustrate how a pair of unstable eigenvalues emerge from a Hopf bifurcation as τ is increased. We considera two-spot equilibrium ring pattern in the unit disk with ε = 0 .
02 and D = 1. The two spots are centered at ( ± r , r ≈ . N = 2. By varying the feedrate a , on the range 4 . ≈ a comp < a < a f ≈ . τ for the synchronous mode ( j = 1) and the asynchronous mode ( j = 2). This is doneby using Newton’s method to solve for ( λ ( j ) I , τ ( j ) H ) inRe (cid:104) F j ( iλ ( j ) I , τ ( j ) H ) (cid:105) = 0 , Im (cid:104) F j ( iλ ( j ) I , τ ( j ) H ) (cid:105) = 0 , for j = 1 , . (2.48)The results for τ ( j ) H and λ ( j ) I for j = 1 , a are shown in the left and right panels of Fig. 3, respectively.From this figure, we observe that the mode that synchronizes the temporal oscillations in the spot amplitudes is the first11o go unstable as τ is increased. A numerical implementation of the winding number criterion in (2.32) yields that thetwo-spot ring pattern is linearly stable when τ < min( τ (1) H , τ (2) H ). (a) Hopf bifurcation threshold for τ (b) Imaginary eigenvalue at the Hopf bifurcation Figure 3: Left panel: The Hopf bifurcation value of τ for the synchronous ( j = 1) and asynchronous mode ( j = 2),as computed from (2.48), for the linearization of a two-spot ring steady-state solution with ε = 0 . D = 1, and ringradius r = 0 . a is increased. Right panel: Thecorresponding imaginary eigenvalue for the two modes.To confirm the Hopf bifurcation threshold, as calculated from (2.48), we compute full numerical solution to the PDE(1.1) for ε = 0 . D = 1, using as an initial condition a two-spot ring pattern with ring radius r = 0 . a = 6,we have τ (1) H ≈ .
56 and τ (2) H ≈ .
28 from (2.48). With the choice τ = 54, for which τ (1) H < τ < τ (2) H , we predict fromthe GCEP that the amplitudes of the two spots will oscillate in phase. In the PDE simulation results of Fig. 4a we showthat there are synchronous oscillations of the spot amplitudes, which eventually leads to the disappearance of both spots.By increasing the feed rate to a = 7 .
2, we have τ (1) H ≈ .
56 and τ (2) H ≈ .
11 from (2.48). With the choice τ = 137, wepredict that the two spots will be unstable to both synchronous and asynchronous perturbations in the spot amplitudes.In the PDE simulation results of Fig. 4b we show that, although initially the spot amplitudes oscillate synchronously. astime increases these oscillations become asynchronous, and eventually one of the two spots is annihilated. Sp o t h e i g h t (a) a = 6 and τ (1) H < τ < τ (2) H . Sp o t h e i g h t (b) a = 7 . τ (2) H < τ (1) H < τ . Figure 4: PDE simulation results of (1.1) for the spot amplitudes versus time starting from a two-spot steady-state ringpattern with ε = 0 . D = 1, and ring radius r = 0 . a = 6 and τ = 54. Synchronous oscillations occur,leading to the annihilation of both spots. Right panel: a = 7 . τ = 137. Eventually asynchronous spot amplitudeoscillations occur, leading to the annihilation of only one spot.12 A perforation of the domain as a localized defect
In this section we analyze how the existence, linear stability, and slow dynamics of quasi-equilibrium spot patterns areaffected by removing a small circular hole of radius O ( ε ) from Ω, given byΩ ε = (cid:8) x ∈ Ω : | x − x | ≤ Cε (cid:9) , where C > O (1) parameter controlling the size of the hole. In the perforated domain, the Schnakenberg model is v t = ε ∆ v − v + uv , τ u t = D ∆ u + a − uv ε , in ¯Ω ≡ Ω \ Ω ε , (3.1a) ∂ n v = ∂ n u = 0 on ∂ Ω ; v = u = 0 on ∂ Ω ε . (3.1b)The homogeneous Dirichlet boundary conditions on ∂ Ω ε models the leakage of the activator v and the substrate u throughthe boundary of the small hole. N -spot pattern and slow dynamics We begin by constructing a quasi-equilibrium N -spot pattern with spots located at x , . . . , x N in the perforated domain.We assume, initially, that this pattern is linearly stable on O (1) time intervals. In our analysis below, we assume that | x i − x j | = O (1) for i (cid:54) = j , and that dist( x i , ∂ Ω) = O (1) and dist( x i , ∂ Ω ε ) = O (1) for i = 1 , . . . , N .Following the derivation in § u + aD − π √ D N (cid:88) i =1 S i δ ( x − x i ) = 0 in ¯Ω , ∂ n u = 0 on ∂ Ω ; u = 0 on ∂ Ω ε , (3.2)where S , . . . , S N denote the spot source strengths. However, this outer problem is of singular perturbation type since u must satisfy the extra conditon u = 0 on ∂ Ω ε . To proceed, we will use strong localized perturbation theory to replacethe effect of the hole with a Dirac singularity. To do so, near the hole centered at x we introduce local coordinates y = ε − ( x − x ) and u ∼ U ( y ) / √ D . From (3.2), we obtain to leading order that∆ y U = 0 , | y | ≥ C ; U = 0 , on | y | = C , (3.3)which has the solution U = S log ( | y | /C ), where S is to be determined. This yields the matching condition u ∼ U √ D ∼ S √ D (cid:18) log | x − x | + 1 ν − log C (cid:19) , as x → x , (3.4)where ν = − / log ε . Owing to the identity (cid:90) ∂ Ω ε − D∂ n u | ∂ Ω ε ds ∼ πS √ D , (3.5)where ∂ n denotes the outward normal derivative to ¯Ω, the constant S is proportional to the diffusive flux of inhibitorthrough the hole. The strength of this leakage term, mediated by S , is calculated below in a self-consistent way.By superimposing the Dirac singularity 2 πS √ D δ ( x − x ) on the outer problem to account for the logarithmic singularityin (3.4), we replace (3.2) with the modified outer problem∆ u + aD − π √ D N (cid:88) i =0 S i δ ( x − x i ) = 0 in Ω ; ∂ n u = 0 on ∂ Ω , (3.6)which is defined at O (1) distances from the spot locations and from the center of the hole.13he solution to (3.6) is represented in terms of the Neumann Green’s function of (2.9) as u ( x ) = − π √ D N (cid:88) i =0 S i G ( x ; x i ) + ¯ u , (3.7)where ¯ u is a constant to be determined. By applying the divergence theorem to (3.6) we get N (cid:88) i =0 S i = a | Ω | π √ D ≡ p a . (3.8)We let x → x in (3.7) in order to asymptotically match the local behavior of u with the far-field behavior (3.4) for thesolution near the hole. This matching yields the algebraic equation S + 2 πν (cid:32) S R , + N (cid:88) i =1 S i G ,i (cid:33) − νS log C = ν √ D ¯ u , (3.9)where R , ≡ R ( x ; x ) and G ,i ≡ G ( x ; x i ).Next, we match the local behavior of the outer solution in (3.7) near each spot with the far-field behavior (2.10) of thecorresponding inner solution. Letting x → x j in (3.7) we obtain that u ∼ S j √ D log | x − x j | − π √ D S j R j,j + N (cid:88) i (cid:54) = j S i G j,i + S G j, + ¯ u − π S j ∇ x R j,j + N (cid:88) i (cid:54) = j S i ∇ x G j,i + S ∇ x G j, · ( x − x j ) + O ( | x − x j | ) , j = 1 , . . . , N . (3.10)By matching the O (1) terms in (2.10) and (3.10), we obtain that S j + 2 πν S j R j,j + N (cid:88) i (cid:54) = j S i G j,i + S G j, + νχ ( S j ) = ν √ D ¯ u , j = 1 , . . . , N . (3.11)We write the nonlinear algebraic system (3.8), (3.9), and (3.11) for S , . . . , S N and ¯ u in matrix form as S = p a − e T s , s + 2 πν ( G s + S g ) + ν χ = (cid:16) ν √ D ¯ u (cid:17) e , θS = ν √ D ¯ u − πν g T s , (3.12a)where we have defined s ≡ ( S , . . . , S N ) T , g ≡ ( G , , . . . , G ,N ) T , e ≡ (1 , . . . , T ∈ R N , χ ≡ ( χ ( S ) , . . . , χ ( S N )) T , θ ≡ πνR , − ν log C . (3.12b)Here
G ∈ R N × N is the Neumann Green’s matrix characterizing inter-spot interactions for spots centered at x , . . . , x N .By eliminating S between the first and third equations in (3.12a), we can solve for ¯ u as¯ u = θp a + s T (2 πν g − θ e ) ν √ D . (3.13)By substituting (3.13) together with S = p a − e T s into the middle equation of (3.12a) we obtain the following nonlinearalgebraic system for the vector s of spot strengths: s + 2 πν G s + ( e T s ) ( θ e − πν g ) + ν χ = p a ( θ e − πν g ) . (3.14)14ext, to derive the DAE system for slow spot dynamics, we match (2.10) with (3.10) for the O ( ε ) gradient terms.Denoting y = ε − ( x − x j ), and using S = p a − e T s , this yields the following far-field behavior for the correction U j tothe leading order core solution, as defined in (2.1): U j ∼ − ( βββ j + 2 πS ∇ x G j, ) · y = − (cid:34) βββ j + 2 π (cid:32) p a − N (cid:88) i =1 S i (cid:33) ∇ x G j, (cid:35) · y , as | y | → ∞ . (3.15)Here βββ j is defined in (2.12). Following the derivation in § d x j dσ = − γ ( S j ) (cid:34) βββ j + 2 π (cid:32) p a − N (cid:88) i =1 S i (cid:33) ∇ x G j, (cid:35) , j = 1 , . . . , N , (3.16)where σ = ε t and s ≡ ( S , . . . , S N ) T satisfies the nonlinear algebraic system (3.14). Here γ ( S j ) is defined in (2.18). (a) t = 0 (b) t = 999 x (c) x -coordinate of spots y (d) y -coordinate of spots Figure 5: For ε = 0 . , D = τ = 1 , a = 16 and a hole at x = (0 . , .
5) with radius ε ( C = 1), two spots initially locatedat (0 . ,
0) and (0 , . S ≈ . v at t = 0 and t = 999, respectively. In (c) and (d), we show the very close agreement of spot trajectoriesobtained by the PDE simulation (black dots) and the DAE (3.16) and (3.14) (red solid line). (a) t = 0 (b) t = 1998 x (c) x -coordinate of spots y (d) y -coordinate of spots Figure 6: For ε = 0 . , D = τ = 1 , a = 20 and a hole at x = (0 . ,
0) with radius ε ( C = 1), three spots are initially locatedat (0 . ,
0) ( S ≈ . − . , ± .
3) ( S ≈ . v at t = 0 and t = 1998, respectively. In (c) and (d), we show the very close agreement of spot trajectories obtained bythe PDE simulation (black dots) and the DAE (3.16) and (3.14) (red solid line). We note that the x -coordinates of twospots on the left of the hole almost coincide and their trajectories in the x -direction are indistinguishable.In the unit disk, in Fig. 5 and Fig. 6 we show a very favorable comparison between the spot trajectories as computedfrom the DAE system (3.16) and (3.14) and from the full PDE system (3.1) for the case of two or three spots, respectively.The hole location and radius, and the other parameter values, are given in the figure captions. From Fig. 5 and Fig. 6, weobserve that there is a repulsive interaction between the spots and the small hole. By increasing the feed-rate parameter a , in Fig. 7 we show that a one-spot solution will exhibit spot self-replication when the spot source strength exceeds thepeanut-splitting threshold Σ ≈ . a c = 2Σ ≈ .
6, and is independent of the spot location,we observe from Fig. 7 that a much larger feed rate is needed to trigger a peanut-splitting instability when the domaincontains a hole. Moreover, the required threshold of the feed rate depends on the relative locations of the spot and thecenter of the hole.
In this subsection we analyze the linear stability on an O (1) time-scale of the quasi-equilibria, denoted by v e and u e , asconstructed in § v = v e + e λt φ and u = u e + e λt η into (3.1a) and (3.1b), and linearize to obtain ε ∆ φ − φ + 2 u e v e φ + v e η = λφ , D ∆ η + a − ε − (cid:0) u e v e φ + v e η (cid:1) = τ λη , in ¯Ω ,∂ n φ = ∂ n η = 0 on ∂ Ω ; φ = η = 0 on ∂ Ω ε . (3.17)Following the analysis in § § j th spot is linearly unstable to the peanut-splitting mode when S j > Σ ≈ . S j is obtained from the nonlinear algebraic system (3.14) that depends on the location of the hole.We focus on deriving a GCEP associated with radially symmetric perturbation near a spot, in which m = 0 in the localproblem (2.20). Using the distributional limit (2.24), we obtain for λ (cid:54) = 0 that the outer problem for η away from thespots is ∆ η − τ λD η − πD N (cid:88) i =1 c j δ ( x − x i ) = 0 in ¯Ω , ∂ n η = 0 on ∂ Ω ; η = 0 on ∂ Ω ε . (3.18) (a) t = 0 (b) t = 28 (c) t = 80 Figure 7: For ε = 0 . , τ = D = 1 , a = 18 and a hole at the center with radius ε ( C = 1), a spot located at x = (0 . , S ≈ . ≈ .
302 in § η on thehole boundary by a Dirac Delta forcing of undetermined strength 2 πc D − δ ( x − x ). In this way, the modified outerproblem for η defined at O (1) distances from the spots and the hole is∆ η − τ λD η − πD N (cid:88) i =0 c i δ ( x − x i ) = 0 in Ω , ∂ n η = 0 on ∂ Ω , (3.19)which is subject to the matching condition η ∼ c D (cid:18) log | x − x | + 1 ν − log C (cid:19) , as x → x . (3.20)16he solution to (3.19) is represented in terms of the eigenvalue-dependent Green’s function G λ of (2.27) by η = − πD N (cid:88) i =0 c i G λ ( x ; x i ) . (3.21)We let x → x in (3.21) and equate the resulting O (1) limiting behavior with (3.20). This matching condition yields that c + 2 πν (cid:32) c R λ , + N (cid:88) i =1 c i G λ ,i (cid:33) − νc log C = 0 , (3.22)where R λ , ≡ R λ ( x ; x ) and G λ ,i ≡ G λ ( x ; x i ).Next, by expanding (3.21) as x → x j , we have for each j = 1 , . . . , N that η ∼ c j D log | x − x j | − πD c j R λ ( x j ; x j ) + N (cid:88) i (cid:54) = j c i G λ ( x j ; x i ) + c G λ ( x j ; x ) . (3.23)Upon matching (3.23) with the far-field behavior (2.25b) of the inner problem we obtain c j + 2 πν c j R λ j,j + N (cid:88) i =1 ,i (cid:54) = j c i G λj,i + c G λ j, + νc j ˜ B ( S j ; λ ) = 0 , (3.24)where R λ j,j ≡ R λ ( x j ; x j ) and G λj,i ≡ G λ ( x j ; x i ).We write (3.22) and (3.24) in matrix form as θ λ c + 2 πν g Tλ c = , c + 2 πν ( G λ c + c g λ ) + ν ˜ B c = , (3.25a)where the matrices G λ and ˜ B are defined in (2.29b). In (3.25a) we have defined c ≡ ( c , . . . , c N ) T , g λ ≡ (cid:0) G λ , , . . . , G λ ,N (cid:1) T , e ≡ (1 , . . . , T ∈ R N , θ λ ≡ πνR λ , − ν log C . (3.25b)The GCEP is obtained by eliminating c in (3.25a). In this way, we conclude that a discrete eigenvalue λ of the linearizationmust be such that M c = 0 , where M ( λ ) ≡ θ λ (cid:16) I + 2 πν G λ + ν ˜ B (cid:17) − π ν g λ g Tλ , (3.26a)has a nontrivial solution c (cid:54) = . Here I is the N × N identity matrix. Any such λ (cid:54) = 0 satisfyingdet M ( λ ) = 0 , (3.26b)for which Re( λ ) >
0, corresponds to an instability associated with locally radially symmetric perturbations near the spots.As similar to the analysis in § λ = 0.When λ = 0, the solution to the modified outer problem (3.19) is η = − πD N (cid:88) i =0 c i G ( x ; x i ) + ¯ η , where N (cid:88) i =0 c i = 0 , (3.27)and where ¯ η is an additive constant to be found. Here, G is the Neumann Green’s function satisfying (2.9). By matchingthe local behavior of η to the far-field behavior (3.20) near the hole as well as to the far field behavior (2.25b) near thespots, we obtain in matrix form that θc + 2 πν g T c = νD ¯ η , c + 2 πν ( G c + c g ) + ν ˜ B c = νD ¯ η e , (3.28a)where G ∈ R N × N is the Neumann Green’s matrix and ˜ B = diag ( χ (cid:48) ( S ) , . . . , χ (cid:48) ( S N )), as is given in (2.36). In (3.28a) wehave defined c ≡ ( c , . . . , c N ) T , g ≡ ( G , , . . . , G ,N ) T , e ≡ (1 , . . . , T ∈ R N , θ ≡ πνR ( x ; x ) − ν log C , (3.28b)17here G ,i ≡ G ( x ; x i ). Since (cid:80) Ni =0 c i = 0, we can write c = − e T c . Upon eliminating c in (3.28a), we conclude that λ = 0 is an eigenvalue of the linearization if and only if M c = 0 , where M ≡ I + θN E + 2 πν G + ν ˜ B − πν (cid:0) ge T + eg T (cid:1) , (3.29)has a nontrivial solution c (cid:54) = . Parameter values corresponding to zero-eigenvalue crossings are where det M = 0. N -spots with leakage at the center We consider a ring pattern of N -spots, with spots centered at (2.39), in the perforated unit disk ¯Ω that has a hole ofradius Cε at the origin. Since the N spots have a common source strength S c , we let s = S c e in (3.14). Upon using G e = p ( r ) e /N from (A.3), where r is the ring radius, together with θ = 1 + 2 πνR ( ; ) − ν log C = 1 − ν (cid:18) log C + 34 (cid:19) , g = G ( x j ; ) e = 12 π (cid:18) − log r + r − (cid:19) e , (3.30)as calculated from (A.1), we obtain from (3.14) and (3.8) that S c satisfies the scalar nonlinear equation S c + νS c N + 1 log (cid:34) r N +10 N C N (1 − r N ) (cid:35) + νχ ( S c ) N + 1 = p a N + 1 (cid:20) ν (cid:18) log (cid:16) r C (cid:17) − r (cid:19)(cid:21) , p a = a √ D . (3.31)Next, by using (A.2), we calculate for a ring pattern that2 π (cid:32) p a − N (cid:88) i =1 S i (cid:33) ∇ x G j, = ( p a − N S c ) (cid:18) r − r (cid:19) e θ j , where e θ j is defined in (2.39). Upon using this result, together with the expression (A.4) for β j for a ring pattern, theODE system (3.16) for slow spot dynamics reduces to the following scalar ODE for the ring radius r : dr dσ = γ ( S c ) (cid:34) p a (cid:18) r − r (cid:19) − S c (cid:32) N + 12 r + N r N − − r N (cid:33)(cid:35) , (3.32)where σ = ε t . Here S c = S c ( r ) is determined from the nonlinear constraint (3.31). It follows that the equilibrium ringradius r = r e of (3.32) with common source strength S c is a root of S c (cid:32) N + 12 r e + N r N − e − r N e (cid:33) = p a (cid:18) r e − r e (cid:19) , (3.33)where S c = S c ( r e ) satisfies (3.31).Next, the GCEP (3.26a) for a ring pattern reduces to finding values of λ for which there are nontrivial solutions to M c = , with M ≡ θ λ (cid:16) ν ˜ B c I + 2 πν G λ (cid:17) − N π ν β λ E , E ≡ N ee T , (3.34)where ˜ B c ≡ ˜ B ( S c ; λ ) is calculated from (2.22) and where β λ ≡ G λ ( x ; ) = . . . = G λ ( x N ; ). Since G λ is a cyclicsymmetric matrix, it has the eigenspace c = e and c = q j , where e T q j = 0 and q Tj q i = 0 for i (cid:54) = j and i, j = 2 , . . . , N .In this way, from (3.34), the discrete eigenvalues λ for the synchronous ( c = e ) mode and competition modes ( c = q j , j = 2 , . . . , N ) are the roots of F ≡ θ λ (1 + ν ˜ B c + 2 πν ω ) − N π ν β λ = 0 , (3.35a) F j ≡ θ λ (1 + ν ˜ B c + 2 πν ω j ) = 0 , j = 2 , . . . , N , (3.35b)where the matrix eigenvalues ω i = ω i ( λ ) of G λ are defined by G λ e = ω e and G λ q i = ω i q i for i = 2 , . . . , N .18ext, we derive the threshold condition on the parameters for which there is a zero-eigenvalue crossing in the GCEP.We use ˜ B = ˜ B ( S c ; 0) I ≡ χ (cid:48) ( S c ) I , together with (3.30), to obtain that (3.29) reduces to M c = 0 , where M ≡ [1 + νχ (cid:48) ( S c )] I + N (cid:20) − ν (cid:18) log Cr + r − (cid:19)(cid:21) E + 2 πν G . (3.36)By using (A.3), we conclude from (3.36) that a zero-eigenvalue crossing for the mode c = e occurs if and only if S c satisfies N + 1 + νχ (cid:48) ( S c ) + ν log (cid:32) r N +10 N C N (1 − r N ) (cid:33) = 0 . (3.37)We now show that this zero-eigenvalue threshold condition (3.37) occurs precisely at the value of S c for which the root S c = S c ( r ) to (3.31) has a saddle-node bifurcation. To see this, we differentiate (3.31) with respects to S c to obtain νr (cid:18) p a ( C − r ) r − S c ( N + 1 + ( N − r N ) r (1 − r N ) (cid:19) dr dS c = N + 1 + νχ (cid:48) ( S c ) + ν log (cid:32) r N +10 N C N (1 − r N ) (cid:33) . (3.38)At a saddle-point point ( r f , S cf ) we have dr /dS c = 0, and so the right-hand side of (3.38) must vanish at that point,which yields (3.37). We conclude that a zero-eigenvalue crossing of the GCEP can only occur at the location of asaddle-node bifurcation point for a quasi-equilibrium ring pattern.Next, to determine the threshold condition on the parameters for a zero-eigenvalue crossing for the competition modes,we substitute c = q i for i = 2 , . . . , N into (3.36), and use E q i = 0 to obtain1 + νχ (cid:48) ( S c ) + 2 πνσ i = 0 , i = 2 , . . . , N . (3.39)Here σ i are eigenvalues of the Neumann Green’s matrix for which G q i = σ i q i for i = 2 , . . . , N . Roots of the coupledproblem (3.39) and (3.30) correspond to the threshold values ( S ( i ) c , r ( i )0 ), for i = 2 , . . . , N , where a zero-eigenvalue crossingof the GCEP occurs. We first consider a one-spot quasi-equilibrium solution in the perforated unit disk ¯Ω. In this subsection, we fix ε = 0 . D = τ = 1. By taking the ring radius r as a bifurcation parameter, in Fig. 8a we show that (3.31) has a foldbifurcation structure for the source strength of the spot. From this figure, we observe that a one-spot quasi-equilibriumsolution does not exist when the spot is too close to the center of the hole located at the origin. In contrast, when thereis no hole, a one-spot quasi-equilibrium solution exists for all r ≥ r = 0 . S c versus r bifurcation diagram, computed from (3.31), changes with respectto the feed-rate parameter a and the parameter C > εC of the hole. We observe that as either a increases or C decreases (smaller hole radius), a one-spot quasi-equilibrium solution can exist closer to the hole.Next, we use numerical continuation on (3.31) and the saddle-node condition (3.37) to determine how the saddle-nodepoint r f for the ring radius depends on the feed-rate parameter a when C = 1. A similar numerical continuation of (3.31)and the steady-state ring radius condition (3.33), also reveals a saddle-node bifurcation structure of r e . These results,presented in Fig. 10a, show that a one-spot quasi-equilibrium solution exists only when a is greater than the saddle-nodevalue a f ≈ . a > a f , there are two fold-point values of r f for quasi-equilibria: one near the boundary ofthe unit disk (not shown in Fig. 8a) while the other is closer to the hole. For each a > . (a) Bifurcation diagram Sp o t h e i g h t (b) Spot height Figure 8: We fix C = 1, D = τ = 1, ε = 0 .
02, and a = 10. Left panel: S c versus r for a one-spot quasi-equilibriumsolution, as computed from (3.31). The saddle-node bifurcation is at ( r f , S cf ) ≈ (0 . , . r = 0 . S ≈ . S ≈ . v , with these twoinitial source strengths, as computed from the full PDE (3.1). The bottom (solid) curve shows that the spot on the lowerbranch is rapidly annihilated. (a) Outside to inside : a = 12 , , (b) Outside to inside : C = 0 . , , . Figure 9: In both panels, the middle curve is the same as in Fig. 8a, which corresponds to ε = 0 . , C = 1 and a = 10. Alllower branches have an unstable eigenvalue for the GCEP (3.34). Left panel: We fix C = 1. The saddle-node bifurcationfor a = 8 and a = 12 occurs at r f ≈ . r f ≈ . a = 10. The saddle-nodebifurcation for C = 0 . C = 1 . r f ≈ . r f ≈ . (a) Fix C = 1 (b) Fix a = 10 Figure 10: The saddle-node structures of r f (bigger U-shape) and the equilibrium ring radius r e (smaller U-shape) withrespect to the feed rate a (left panel) and the hole radius parameter C (right panel) for a one-spot solution. Along thedashed portion of the r e branch, the GCEP (3.34) has an unstable eigenvalue. For each feed-rate a exceeding a threshold,there is only one stable equilibrium location for the one-spot solution. a = 10, we show a similar saddle-node bifurcation structure for r f and r e versus the parameter C , which controls theradius of the hole. We observe that there is no quasi-equilibrium one-spot solution if the hole radius exceeds a certainthreshold.In Fig. 11a, we show full PDE results computed from (3.1) for a one-spot quasi-equilibrium solution, initially locatedat r = 0 .
57, in which the feed-rate parameter is slowly decreased in time according to a = max(7 . − . t, . t ≈
20. Thisrapid decay of the spot amplitude is due to the non-existence of one-spot quasi-equilibria for r = 0 .
57 when a decreasesbelow the saddle-node value a f . Alternatively, in Fig. 11a, the full PDE simulation results shows that the one-spotquasi-equilibrium persists when the feed rate is fixed at a = 7 . > a f . To motivate a further, but more delicate, PDEsimulation result, we observe from Fig. 10a that the saddle-node value for r e occurs at a e ≈ . a f ≈ . a between a f and a e , a quasi-equilibrium one-spot solution exists for some range of r , but there is no steady-state equilibrium value r e . In Fig. 11b we show results from a full PDE simulation of (3.1)for a one-spot quasi-equilibrium initially located at r = 0 .
57 and with feed-rate a = 7 .
48, which satisfies a f < a < a e .We observe that the one-spot quasi-equilibrium survives only until t ≈ a = 7 . > a e shows that the one-spot quasi-equilibrium solution persists, and slowly drifts away from thehole towards its stable equilibrium location at around t ≈ Next, we demonstrate the occurrence of a Hopf bifurcation in the spot amplitude for a one-spot quasi-equilibrium solutionin the perforated unit disk. By fixing ε = 0 . , a = 10, and C = D = 1, in Fig. 12 we plot the Hopf bifurcation thresholdvalue τ = τ H on the range r ∈ [0 . , . τ H , λ I ) fromRe [ F ( τ H , iλ I )] = 0 , Im [ F ( τ H , iλ I )] = 0 , (3.40)where F is defined in (3.35a). In particular, when r = 0 .
6, we compute that τ H ≈ .
6. To confirm this thresholdvalue, in Fig. 13 we plot the spot amplitude for a one-spot quasi-equilibrium solution with r = 0 . τ = 162 < τ H , τ = 168 > τ H , and for τ = 170, as computed from a full PDE simulation of (3.1). For τ = 168 we observe a small-scaleperiodic oscillation of the spot amplitude, suggesting that the Hopf bifurcation is supercritical. However, for the largervalue τ = 170, we observe that the temporal oscillation in the spot amplitude can grow and lead to spot annihilation.21
20 40 60 80 100t0.00.20.4 Sp o t h e i g h t (a) a = 7 . a = max(7 . − . t, . Sp o t h e i g h t (b) a = 7 . a = 7 . Figure 11: We fix C = 1. Left panel: Short-time evolution of the amplitude of a one-spot quasi-equilibrium solution for aconstant feed rate a = 7 . a = max(7 . − . t, .
4) (dashed line). Rightpanel: Longer time evolution of the spot amplitude for a ≡ . a = 7 .
48 (dashed line). When a = 7 .
6, theone-spot solution has become close to its equilibrium value when t ≈ r = 0 .
57, and the numerical results were computed from the full PDE (3.1).
Figure 12: Hopf bifurcation threshold τ H versus r for a one-spot quasi-equilibrium solution, as computed from (3.40),for a = 10, D = C = 1, and ε = 0 . Sp o t h e i g h t (a) τ = 162 Sp o t h e i g h t (b) τ = 168 Sp o t h e i g h t (c) τ = 170 Figure 13: For ε = 0 . , a = 10 , C = D = 1, we choose three values of τ near the Hopf bifurcation threshold τ H ≈ . r = 0 .
6. (a) τ = 162 < τ H : the spot amplitude has decayingoscillations. (b) τ = 168 > τ H : small amplitude oscillations indicating a supercritical Hopf bifurcation. (c) τ = 170: spotamplitude oscillations grow and trigger an oscillatory collapse of the spot.22 .3.3 Competition instability of a two-spot pattern Here we consider a two-spot quasi-equilibrium pattern in the perforated unit disk, with parameters ε = 0 . , C = D = τ = 1, and a = 10. In Fig. 14, we plot the bifurcation diagram of S c versus r for N = 2 spots, as computed from (3.31),showing a saddle-node bifurcation behavior. We calculate that the saddle-node point occurs at r = r (1)0 ≈ . r = r (2)0 ≈ . r (1)0 < r < r (2)0 , and the upper branch on r > r (2)0 . On the lower branch, we compute that there isa root to F = 0 to (3.35a) with Re( λ ) >
0, and so the GCEP (3.34) has an unstable eigenvalue. This indicates that,on the lower branch, the two-spot pattern is unstable to synchronous locally radially-symmetric perturbations near thespots. Along the upper branch with r < r (2)0 there is a root to F = 0 in (3.35b) with Re( λ ) >
0, and so this segmentof the bifurcation diagram is unstable to asynchronous locally radially-symmetric perturbations. Finally, on the upperbranch with r > r (2)0 , there is no root to (3.35b) in Re( λ ) >
0, and so this segment is linearly stable. These linearstability predictions are validated in Fig. 15 from full PDE simulations of (3.1) with initial conditions chosen in thesethree segments of the bifurcation diagram in Fig. 14.
Figure 14: The saddle-node bifurcation point and the competition threshold for a two-spot ring solution are shown asblack circle and square markers, respectively. The dashed (solid) segment of upper branch corresponds where r < r (2)0 ( r > r (2)0 ). Here r (2)0 is where there is a zero-eigenvalue crossing of the GCEP (3.36) for the competition mode. Theparameters are C = D = τ = 1, a = 10, and ε = 0 . Sp o t h e i g h t (a) Two spots survive. Sp o t h e i g h t (b) One spot survives. Sp o t h e i g h t (c) Both spots die. Figure 15: The evolution of the spot amplitudes for a two-spot quasi-equilibrium solution, as computed from the fullPDE (3.1). In (a) and (b), the two spots have initial condition on the upper branch of Fig. 14. Their initial locationsare at ( ± r , r = 0 . > r (2)0 and r = 0 . < r (2)0 in (a) and (b), respectively. In (c), the two-spots have initialcondition on the lower branch, with r = 0 .
29. The parameters are C = D = τ = 1, a = 10, and ε = 0 .
02. The PDEresults are in agreement with linear stability predictions. 23
Pinning effects from a spatially localized feed-rate
In this section, we analyze slow spot dynamics for the case where the localized heterogeneity consists of a localized sourceof feed from the substrate of the form a ( x ) = a + ε − a Φ (cid:0) | x − ξξξ | /ε (cid:1) , Φ( r ) ≡ exp( − r / / (2 π ) , (4.1)where a > a > ξξξ ∈ Ω is the location of the concentration of the feed.
We first modify our asymptotic construction of N -spot quasi-equilibria given in § § § u + a ( x ) D − π √ D N (cid:88) i =1 S i δ ( x − x i ) = 0 in Ω , ∂ n u = 0 on ∂ Ω , (4.2)where S , . . . , S N are the source strengths of the N spots. By applying the divergence theorem to (4.2) we get N (cid:88) i =1 S i = (cid:82) Ω a d x π √ D ≡ p a . (4.3)We decompose the solution to (4.2) as u ( x ) = u ( x ) D − π √ D N (cid:88) i =1 S i G ( x ; x i ) + ¯ u , (4.4)where ¯ u is a constant and G is the Neumann Green’s function of (2.9). Here u ( x ) is the unique solution to∆ u = − a ( x ) + (cid:82) Ω a ( x ) d x | Ω | in Ω , ∂ n u = 0 on ∂ Ω ; (cid:90) Ω u dx = 0 , (4.5)which is given in terms of G by u ( x ) = (cid:90) Ω a ( z ) G ( z ; x ) d z . (4.6)As in § x → x j for j = 1 , . . . , N between the outer solution and innersolutions to derive a nonlinear algebraic system for ¯ u and the source strengths. Letting x → x j in (4.4), we obtain that u ∼ u ( x j ) D + S j √ D log | x − x j | − π √ D S j R j,j + N (cid:88) i (cid:54) = j S i G j,i + ¯ u + D ∇ x u ( x j ) − π √ D S j ∇ x R j,j + N (cid:88) i (cid:54) = j S i ∇ x G j,i · ( x − x j ) + O ( | x − x j | ) , j = 1 , . . . , N , (4.7)where R j,j ≡ R ( x j ; x j ) and G j,i ≡ G ( x j ; x i ).Upon matching (4.7) with (2.10) for the O (1) terms, we write the resulting equations in matrix form as s + 2 πν G s + ν χ = ν (cid:18) √ D u + ¯ u √ D e (cid:19) , e T s = p a , (4.8a)24here G is the Neumann Green’s matrix, and where we have defined s ≡ ( S , . . . , S N ) T , χ ≡ ( χ ( S ) , . . . , χ ( S N )) T , e ≡ (1 , . . . , T ∈ R N , u ≡ ( u ( x ) , . . . , u ( x N )) T . (4.8b)Upon left-multiplying (4.8a) by e T , we can isolate ¯ u as¯ u = p a + 2 πν e T G s + ν e T χ νN √ D − e T u N D . (4.9)By using (4.9) to eliminate ¯ u in (4.8a), we obtain a nonlinear algebraic system for the vector of source strengths s , s + 2 πν ( I − E ) G s + ν ( I − E ) χ = ν √ D ( I − E ) u + p a N e , where E = 1 N ee T , (4.10)and p a is defined in (4.3).To derive the DAE system for slow spot dynamics we must match (2.10) with (4.7) for the O ( ε ) gradient terms. Thismatching yields the far-field behavior for the inner correction term U j , as defined in (2.1), given by U j ∼ (cid:18) √ D ∇ x u ( x j ) − βββ j (cid:19) · y as | y | → ∞ , (4.11)where y = ε − ( x − x j ) and βββ j is defined in (2.12). Following the derivation in § d x j dσ = γ ( S j ) (cid:18) √ D ∇ x u ( x j ) − βββ j (cid:19) , j = 1 , . . . , N , (4.12)where σ = ε t and s ≡ ( S , . . . , S N ) T satisfies the nonlinear algebraic system (4.10). Here γ ( S j ) is defined in (2.18).As ε →
0, we can approximate, in the sense of distributions, the heterogeneous feed rate in (4.1) as a ( x ) → a + a δ ( x − ξξξ ) . (4.13)In this way, u in (4.6) can be calculated explicitly, by using Green’s reciprocity and (cid:82) Ω G ( z ; x ) d x = 0, as u ( x ) = (cid:90) Ω a ( z ) G ( z ; x ) d x = a G ( x ; ξ ) . (4.14) For a one-spot solution, we use (4.13) in (4.3) to calculate S . Then, by using (4.14) in (4.12), together with the explicitexpressions (A.2) for the gradients of the Neumann Green’s function for the unit disk, we obtain from (4.12) that theslow dynamics of a one-spot quasi-equilibrium solution is d x dσ = − a γ ( S )2 π √ D H ( x ) , with S = a π + a π √ D , (4.15a)where σ = ε t and H is defined by H ( x ) ≡ a a (cid:20) x − ξ | x − ξ | + x | ξ | − ξ | x | | ξ | − x · ξ + 1 (cid:21) + x − | x | (cid:20) a a + π (2 − | x | ) (cid:21) . (4.15b)Without loss of generality we let ξξξ = ( ξ,
0) with 0 < ξ <
1. By symmetry, any equilibrium to (4.15) lies on the line thatconnects the origin and ξ . As such, we let x = ( r ,
0) and obtain from (4.15) that r = r ( σ ) satisfies the scalar ODE dr dσ = − a γ ( S )2 π √ D K ( r ) , where K ( r ) ≡ a a (cid:18) r − ξ + r − ξ (1 − r )(1 − ξr ) (cid:19) + πr (2 − r )1 − r . (4.16)25ince γ ( S ) > K ( r ) > ξ < r <
1, it follows that dr /dσ < ξ < r < r e for (4.16), satisfying K ( r e ) = 0, must be on the range 0 < r < ξ . The effect of therelative magnitude of the localized feed to the background feed appears in (4.16) in the form of their ratio a /a . Takingthis ratio as a bifurcation parameter, in Fig. 16a we plot the bifurcation diagram of the roots to K ( r ) = 0 for ξ = 0 . r (1)0 e < r (2)0 e provided that a /a < . a /a > . K (cid:48) ( r (1)0 e ) >
0, we conclude that r e is a stable equilibrium point of (4.16), while r e is an unstable equilibrium. To furtherdemonstrate the saddle-node bifurcation value of a /a , in Fig. 16b we plot K ( r ) on 0 < r < ξ for the four values a /a = 0 . , . , .
72 and 0 .
8. For a /a < . dr /dσ > r (2)0 e < r < ξ and dr /dσ < ξ < r <
1. Moreover, since dr /dσ = O [1 / ( r − ξ )] as r → ξ , this implies that a spot initially located at some r (0)with r (0) > r (2)0 e will get pinned at the concentration point ξ of the feed rate at a finite time. Moreover, if a /a > . any initial point r (0) in 0 < r (0) < r = r (1)0 e for any r (0) < r (2)0 e when a /a < . r = ξ if r (0) > r (2)0 e and a /a < . r = ξ for any r (0) in 0 < r (0) < a /a > . a /a is independent of the inhibitor diffusivity D . Although our asymptoticanalysis, leading to the ODE (4.15), is only valid when the spot is well-separated from the concentration point the feedrate, i.e. when | x − ξ | (cid:29) O ( ε ), the prediction of finite-time pinning phenomena provides a motivation for the analysis in § (a) Bifurcation diagram (b) K Figure 16: The concentration point for the feed rate is ξ = (0 . , r e of K ( r ) = 0, as defined in (4.16), versus a /a . A saddle-node bifurcation occurs at a /a ≈ . K ( r ) for a /a = 0 . , . , .
72 and 0 .
8, respectively.To illustrate these results we compare predictions based on the scalar ODE (4.16) with full PDE simulations of (1.1)with the feed rate (4.1) in the unit disk with D = τ = 1, and ε = 0 .
03. We set ξ = 0 . a = 6 and a = 4, for which a /a < . r (1)0 e ≈ . r (2)0 e ≈ . r (0) = 0 . < r (1)0 e slowly drifts to r (1)0 e . In contrast, for the same a and a , but with initial value r (0) = 0 . > r (2)0 e , we observe fromFig. 17b that the spot approaches ξ = 0 .
7. The full PDE and ODE results are found to agree well until the spot is near ξ = 0 .
7. We remark that the velocity field in the ODE becomes singular as r → ξ owing to the Dirac delta functionapproximation of the localized feed rate. Finally, if we increase the relative strength of the concentration of the feed rateso that a = 6 and a = 5, for which a /a > . r (0) = 0 . ξ owing to the absence of any equilibrium for this ratio a /a . Next, we consider a ring pattern of N -spots in the unit disk with localized feed rate concentrated at the origin, so that ξξξ = . By using ∇ x u = a ∇ x G ( x ; ), together with (A.2) and (A.4) for ∇ x G ( x ; ) and β j , respectively, we obtain from26
500 1000t0.200.250.30 x (a) Drift to r (1)0 e x (b) Pinned at ξ x (c) Pinned at ξ Figure 17: The x -coordinates of the spot trajectory computed from the full PDE (1.1) with (4.1) (black dots) and thescalar ODE (4.16) with ξ = 0 .
7. Left panel: a = 6, a = 4, and r (0) = 0 .
2. Middle panel: a = 6, a = 4, and r (0) = 0 .
53. Right panel: a = 6, a = 5, and r (0) = 0 . r satisfies the scalar ODE dr dσ = − a γ ( S c )2 πr √ D D ( r ) , where D ( r ) ≡ N + 12 N (cid:20) a a − π (cid:18) N − N + 1 (cid:19)(cid:21) + πr + (cid:18) π + a a (cid:19) r N − r N . (4.17)From (4.3) and (4.13), the common spot source strength is S c = ( a π + a ) / [2 πN √ D ].The equilibrium ring radius r e is a root to D ( r ) = 0. Since D (cid:48) ( r ) > D → + ∞ as r → − , the ODE (4.17)must have an equilibrium point in 0 < r < D (0) = N + 12 N (cid:20) a a − π (cid:18) N − N + 1 (cid:19)(cid:21) < , which implies a a < π (cid:18) N − N + 1 (cid:19) . (4.18)For N = 2, in Fig. 18 we plot the bifurcation diagram of the equilibrium ring radius r e versus the ratio a /a . Onthe range 0 ≤ a /a < π/ ≈ . r e → a /a → π/ ≈ . a /a in (4.18) for N = 2.Next, we fix a = 4 .
3, and D = τ = 1. The analysis of competition instabilities and the derivation of the GCEPfor two-spot equilibria with feed concentration at the origin is exactly the same as in § S c =( a π + a ) / [4 π √ D ] with D = 1 for the common source spot strength. This leads to the root finding criterion (2.44)with j = N = 2 for the GCEP (2.42) and the zero-eigenvalue crossing condition (2.46) with j = N = 2. When a = 0(no feed concentration), Fig. 2a showed that there is a competition instability for a steady-state two-spot ring patternif a < .
45. From a numerical computation of the winding number (2.32) and the zero-eigenvalue crossing condition(2.46) with j = N = 2, we obtain that the dashed portions in the bifurcation diagram in Fig. 18 for the equilibrium ringradius correspond to where the two-spot equilibrium solution is unstable to a competition instability. As expected, since a = 4 . < .
45, we observe that the two-spot equilibrium is unstable if a is sufficiently small. Moreover, the two-spotequilibrium is unstable near a /a ≈ π/ r e is too small). However, thekey new qualitative feature of Fig. 18 is that there is a range of a /a where a concentration of feed at the origin stabilizes a two-spot equilibrium solution, which without the concentration of feed would be unstable to a competition stability.To illustrate this linear stability prediction for a = 4 . D = τ = 1, we take ε = 0 .
02 and perform full PDEsimulations of (1.1) with (4.1) for a two-spot equilibrium ring pattern with spots located at ( ± r e , a /a = 0 . a /a = 0 . a /a ,we confirm from these figures that a competition instability occurs, which triggers the annihilation of a spot. In contrast,for a /a = 0 . ± . , v from the full PDE numerical solution for the parameter set in Fig. 19b. Thisfigure shows that after the competition instability triggers a spot-annihilation event, the surviving spot ultimately get27 Figure 18: Bifurcation diagram of equilibrium ring radius versus the ratio a /a , as computed from setting D ( r ) = 0 in(4.17), for a two-spot pattern in the unit disk with feed rate concentration at the origin. Fixing a = 4 . D = τ = 1,on the range 0 . < r e < . . .
8, as shown in Fig. 1b, for aconventional spot solution that is not near a concentration point of the feed. This observation motivates the analysis in § Sp o t h e i g h t (a) r e = 0 . , a /a = 0 . Sp o t h e i g h t (b) r e = 0 . , a /a = 0 .
25 50 75 100t0.460.480.500.52 Sp o t h e i g h t (c) r e = 0 . , a /a = 0 . Figure 19: Full PDE simulations of (1.1) with (4.1) of the spot amplitudes for three ratios of a /a . The initial conditionfor the PDE is a two-spot equilibrium ring pattern with spots located at ( ± r e , ε = 0 . D = τ = 1,and a = 4 .
3. The competition instability occurring in (a) and (b), leads to spot annihilation. In (c), the two-spotequilibrium is linearly stable.
In this section we consider the Schnakenberg model (1.1) with D = τ = 1 and with localized feed rate (4.1), given by v t = ε ∆ v − v + uv , u t = ∆ u + a + ε − (cid:0) a Φ (cid:0) ε − | x − ξξξ | (cid:1) − uv (cid:1) in Ω , (5.1)with ∂ n v = ∂ n v = 0 on ∂ Ω. For the choice Φ( r ) ≡ exp( − r / / (2 π ), we construct a new type of spot solution that ispinned at the site ξξξ ∈ Ω of the localization of the feed rate. Novel dynamical behaviors associated with including thisnew type of spot solution in a quasi-equilibrium spot pattern are analyzed.28 a) t = 0 (b) t = 25 (c) t = 29 (d) t = 99 Figure 20: Contour plots of v , from full PDE solutions of (1.1) with (4.1), corresponding to the parameter values shownin Fig. 19b. A competition instability triggers spot annihilation, and the surviving spot drifts to the origin where it ispinned by the localized feed rate. We construct the asymptotic profile of a pinned spot solution and we study its linear stability properties with respect tonon-radially symmetric perturbations near the spot. We then consider the effect of a time-varying localized concentrationof the feed rate.
We begin by constructing an asymptotic quasi-equilibrium solution for (5.1) corresponding to a single spot pinned at ξξξ .The quasi-equilibrium problem is ε ∆ v e − v e + u e v e = 0 , ∆ u e + a + ε − (cid:2) a Φ (cid:0) ε − | x − ξξξ | / (cid:1) − u e v e (cid:3) = 0 , (5.2)with ∂ n v e = ∂ n u e = 0 on ∂ Ω and Φ( r ) ≡ exp( − r / / (2 π ). In the inner region near the pinned spot, we look for a locallyradially symmetric solution of the form v e ∼ V ( ρ ) and u e ∼ U ( ρ ) where ρ = ε − | x − ξξξ | . From (5.2) we get that U and V satisfy a new core problem∆ ρ V − V + U V = 0 , ∆ ρ U + a Φ( ρ ) − U V = 0 , < ρ < ∞ , (5.3a) V (cid:48) (0) = U (cid:48) (0) = 0 ; V → , U ∼ S log ρ + χ ( S ; a ) as ρ → ∞ . (5.3b)The quantity χ ( S ; a ) is an O (1) nonlinear function of S and concentration intensity a of the feed rate. In Fig. 21 weplot the numerically computed spot profile V ( ρ ) for various S and a .By integrating the U equation in (5.3) on ρ >
0, we use (cid:82) ∞ Φ( ρ ) ρ dρ = 1 / (2 π ) to obtain the integral identity S + a π = (cid:90) ∞ U V ρ dρ . (5.4)With the identity (5.4), we derive in the sense of distributions that, for ε → ε − (cid:2) a Φ − u e v e (cid:3) → (cid:20) a − π (cid:18)(cid:90) ∞ U V ρ dρ (cid:19)(cid:21) δ ( x − ξξξ ) = − πS δ ( x − ξξξ ) . (5.5)Upon using (5.5) in (5.2), we obtain that the outer problem for u e , defined away from ξξξ , is∆ u e = − a + 2 πS δ ( x − ξξξ ) in Ω , ∂ n u = 0 on ∂ Ω , (5.6)which has the solution u = − πS G ( x ; ξξξ ) + ¯ u . (5.7)29ere G is the Neumann Green’s function satisfying (2.9) and ¯ u is an undetermined constant. By applying the divergencetheorem to (5.6), we obtain that the source strength for the pinned spot is S = a | Ω | π . (5.8)To determine ¯ u , we let x → ξξξ in (5.7) to obtain u ∼ S log | x − ξξξ | − πS R , + ¯ u , where R , = R ( ξξξ ; ξξξ ). Upon matchingthis expression with (5.3b) we obtain that ¯ u = ν − [ S + 2 πνS R , + νχ ( S ; a )]. a =0a =1a =5a =10 (a) V ( ρ ) with S = 3 for several a a =0a =1a =5a =10 (b) V ( ρ ) with S = 5 for several a Figure 21: Solution profiles V ( ρ ) with different a for two values of S , as computed numerically from (5.3). The spotheight increases as the strength a of the feed concentration increases.We now use this construction to account for the spot height of the pinned spot observed in the PDE simulations shownin Fig. 19b, in which a = 4 . a = 4 . a , (5.8) yields that S = 2 .
15. Then, by computingthe solution to the new core problem (5.3) with S = 2 .
15 and a = 4 . V (0) ≈ . . Next, we analyze the linear stability of a pinned spot. We let v e and u e denote the quasi-equilibrium solution and weintroduce the perturbation v = v e + e λt φ , u = u e + e λt η , into (5.1) and linearize. This yields the eigenvalue problem ε ∆ φ − φ + 2 u e v e φ + v e η = λφ , ∆ η − ε (cid:0) u e v e φ + v e η (cid:1) = λη . (5.9)To examine the possibility of locally non-radially symmetric instabilities near the spot, we let φ ∼ e imθ Φ ( ρ ) and η ∼ e imθ N ( ρ ) in (5.9) for integer modes m ≥
2, where ρ = ε − | x − ξξξ | . Then, upon using v e ∼ V ( ρ ) and u e ∼ U ( ρ ), toleading order we obtain an eigenvalue problem in the inner region∆ ρ Φ − m ρ Φ − Φ + 2 U V Φ + V N = λ Φ , ∆ ρ N − m ρ N − U V Φ − V N = 0 , (5.10)where ∆ ρ = ∂ ρρ + ρ − ∂ ρ . For the non-radially symmetric modes with m ≥
2, we can impose that Φ → ρ → ∞ and impose the algebraic decay condition N ∼ O ( ρ − m ) as ρ → ∞ . We remark that the eigenvalue problem(5.10) depends on S and a through the solution V and U to the new core problem (5.3).By discretizing (5.10), we obtain a generalized matrix eigenvalue problem. For each mode m ≥
2, we numericallycompute the eigenvalue λ of the discretization of (5.10) with the largest real part as a function of a and the sourcestrength S . The instability threshold occurs when Re( λ ) = 0. In Fig. 22, we plot Re( λ ) versus S for modes m = 2 , , a . We define Σ m ( a ) to be the spot source strength corresponding to the stability threshold30e( λ ) = 0 for angular mode m and concentrated feed intensity a . When a = 0, where there is no concentration of thefeed rate, we have from [17] (see the summary in § (0) < Σ (0) < Σ (0) < . . . for the mode instability thresholds. Therefore, when a = 0, the peanut-splitting mode m = 2 is the first mode tolose stability as S is increased. However, a qualitatively new result for our pinned spot solution is that this orderingprinciple can be violated if the feed intensity a is large enough. In particular, if a = 20, we observe from Fig. 22d thatΣ (20) < Σ (20), which implies that the m = 3 mode is the first to lose stability as S is increased.To illustrate this instability we compute full numerical solutions to the PDE (5.1) in the unit disk with ε = 0 . a = 20,and concentrated feed rate at the origin ξξξ = (0 , T . We choose a = 17, and so from (5.8) with | Ω | = π we get S = 8 . m = 2 and m = 3 modes are unstable since S = 8 . > Σ (20) > Σ (20), withthe m = 3 mode having the larger positive eigenvalue. In the numerical PDE results shown in Fig. 23 at times t = 212and t = 231, we observe a mode m = 3 instability for the pinned spot that triggers a nonlinear spot-splitting process,but with ultimately only one new spot surviving by time t = 300. However, by increasing the value of a to a = 18 and a = 19 for which S = 9 and S = 9 .
5, we observe from Fig. 24 and Fig. 25, respectively, that the most unstable mode m = 3 mode can trigger the creation of two or even three new spots by a nonlinear spot-splitting event. k=2k=3k=4 (a) a = 1 k=2k=3k=4 (b) a = 10 k=2k=3k=4 (c) a = 16 k=2k=3k=4 (d) a = 20 Figure 22: Numerically computed eigenvalue λ of (5.10) with the largest real part versus S for modes m = 2 , , a . The critical thresholds Σ m ( a ) are the values of S where Re( λ ) = 0. Top leftpanel: a = 1, Σ (1) ≈ . , Σ (1) ≈ . , Σ (1) ≈ . a = 10, Σ (10) ≈ . , Σ (10) ≈ . , Σ (10) ≈ . a = 16, Σ (16) ≈ . , Σ (16) ≈ . , Σ (16) ≈ . a = 20, Σ (20) ≈ . , Σ (20) ≈ . , Σ (20) ≈ . We have shown in § ξξξ for the feed rate, it will get pinned to ξξξ in finite time. This suggests that if the concentration point ξξξ is movingwith time, the spot will pursue ξξξ and remain pinned, provided that the dynamics of ξξξ is slow enough. To examine thisconjecture, we perform a full PDE simulation of (5.1) for a = 8, a = 5, and ξξξ = ξξξ ( ε t ) with ε = 0 .
03, where we choose ξξξ = ( ξ , ξ ) T , ξ = 0 . πε t ) , ξ = 0 . πε t ) . (5.11)31 a) t = 0 (b) t = 212 (c) t = 231 (d) t = 300 Figure 23: PDE simulation results of (5.1) for v in the unit disk with ε = 0 . a = 20, and concentrated feed rate atthe origin ξξξ = (0 , a = 17 the pinned spot exhibits a mode m = 3 instability by time t = 212, but ultimatelyonly one spot persists by t = 300. (a) t = 0 (b) t = 50 (c) t = 64 (d) t = 100 Figure 24: Same caption as in Fig. 23 except that a is increased to a = 18. The mode m = 3 instability of the pinnedspot leads to two new spots. (a) t = 0 (b) t = 31 (c) t = 39 (d) t = 99 Figure 25: Same caption as in Fig. 23 except that a is increased further to a = 19. The mode m = 3 instability of thepinned spot now leads to three new spots. 32n Fig. 26 we show that the trajectory of the pinned spot aligns closely with the motion of the rotating concentrationpoint ξξξ ( ε t ). This supports the conjecture that a spot will follow the trajectory of the concentration point of the feedrate. x , y Figure 26: The concentration point for the feed is rotating on the ring ξξξ = ( ξ ( t ) , ξ ( t )) T = 0 . πε t ) , sin(2 πε t )) T with ε = 0 .
03, and we choose a = 5 and a = 8. The x and y coordinates of the pinned spot, as computed numericallyfrom the full PDE (5.1), are shown by the black dots and square, respectively. We observe a close agreement betweenthe spot coordinates and the coordinates ξ ( t ) (solid line) and ξ ( t ) (dashed line) of the concentration point ξξξ of the feedrate. In this subsection we analyze the slow dynamics and linear stability of quasi-equilibrium spot patterns that have a pinnedspot, such as shown in Fig. 23–25.
We construct a quasi-equilibrium spot pattern, with N spots centered at x j ∈ Ω for j = 1 , . . . , N , and with an additionalpinned spot at the concentration point ξξξ ∈ Ω of the feed rate. We assume that the spots and the pinned-spot arewell-separated in the sense that | x i − x j | = O (1) , i (cid:54) = j , | x i − ξξξ | = O (1) , i = 1 , . . . , N . (5.12)Near the j th spot centered at x j , for j = 1 , . . . , N , we substitute the expansion (2.1) with D = 1 into (5.1). Due to theassumption (5.12), the term Φ ( | x j − ξξξ | /ε ) is exponentially small as ε →
0, and therefore absent to all algebraic ordersin ε . We retrieve the core problem (2.2) and the integration identity (2.4). Likewise, near the the pinned spot at ξξξ , wesubstitute v ∼ V ( ρ ) and u ∼ U ( ρ ) into (5.1) to obtain the new core problem (5.3). Upon using the distributional limits(2.5) (with D = 1) and (5.5), the outer problem for u , defined away from all the spots, is∆ u + a − π N (cid:88) i =1 S i δ ( x − x j ) − πS δ ( x − ξξξ ) = 0 , in Ω , ∂ n u = 0 , on ∂ Ω . (5.13)In terms of the Neumann Green’s function of (2.9), the solution to (5.13) is u = − πS G ( x ; ξξξ ) − π N (cid:88) i =1 S i G ( x ; x i ) + ¯ u , (5.14)where ¯ u is an undetermined constant. By using the divergence theorem on (5.13), we conclude that N (cid:88) i =0 S i = a | Ω | / (2 π ) . (5.15)33ext, we let x → x j , for j = 1 , . . . , N , in (5.14) to obtain that u ∼ S j log | x − x j | − π S j R j,j + S G ( x j ; ξξξ ) + N (cid:88) i (cid:54) = j S i G j,i + ¯ u − π S j ∇ x R j,j + S ∇ x G ( x j ; ξξξ ) + N (cid:88) i (cid:54) = j S i ∇ x G j,i · ( x − x j ) + O ( | x − x j | ) , j = 1 , . . . , N , (5.16)where R j,j ≡ R ( x j ; x j ) and G j,i ≡ G ( x j ; x i ). Upon matching the O (1) terms in (5.16) with the far-field behavior (2.2b)of the leading order core solution, we find that S j + 2 πν S j R j,j + S G ( x j ; ξξξ ) + N (cid:88) i (cid:54) = j S i G j,i + νχ ( S j ) = ν ¯ u , j = 1 , . . . , N . (5.17)Then, we expand (5.14) as x → ξξξ to get u ∼ S log | x − ξξξ | − πS R , − π N (cid:88) i =1 S i G j, + ¯ u + O ( | x − x j | ) , (5.18)where R , ≡ R ( ξξξ ; ξξξ ) and G j, ≡ G ( x j ; ξξξ ). Upon matching (5.18) with the far-field behavior (5.3b) of the new coreproblem, we conclude that S + 2 πν (cid:32) S R , + N (cid:88) i =1 S i G j, (cid:33) + νχ ( S ; a ) = ν ¯ u , (5.19)Next, we write (5.17), (5.19), and (5.15) in matrix form as s + 2 πν G s + ν χχχ = ν ¯ u e , e T s = p a ≡ a | Ω | π , (5.20a)where we have defined s ≡ ( S , S , . . . , S N ) T , χ ≡ ( χ ( S ; a ) , χ ( S ) , . . . , χ ( S N )) T , e ≡ (1 , . . . , T ∈ R N +1 . (5.20b)Here χ ( S ; a ) is defined by the new core problem (5.3) for the pinned spot, while G ∈ R ( N +1) × ( N +1) is the NeumannGreen’s matrix of ξξξ, x , . . . , x N . Upon eliminating ¯ u in (5.20b), we obtain that the nonlinear algebraic system for thevector s of source strengths is s + 2 πν ( I − E ) G s + ν ( I − E ) χχχ = p a N + 1 e , with ¯ u = p a + 2 πν e T G s + ν e T χχχν ( N + 1) . (5.21)Here E = N − ee T ∈ R ( N +1) × ( N +1) and I ∈ R ( N +1) × ( N +1) is the identity matrix.To derive the DAE system for slow spot dynamics we must match (2.10) (setting D = 1) with (5.16) for the O ( ε )gradient terms. This matching condition yields the far-field behavior for the inner correction term U j in (2.1): U j ∼ − π S j ∇ x R j,j + N (cid:88) i (cid:54) = j S i ∇ x G j,i + S ∇ x G ( x j ; ξξξ ) · y as | y | → ∞ , (5.22)where y = ε − ( x − x j ). Following the derivation in § d x j dσ = − γ ( S j ) ( βββ j + 2 πS ∇ x G ( x j ; ξξξ )) , j = 1 , . . . , N , (5.23)where σ = ε t . Here, βββ j and γ ( S j ) are defined in (2.12) and (2.18), respectively, while s ≡ ( S , S , . . . , S N ) T satisfies thenonlinear algebraic system (5.21). 34e now compare the DAE dynamics (5.23) and (5.21) with full numerical results computed from the PDE (5.1) in theunit disk. We set ε = 0 .
03, and for the localized feed rate we choose a = 15, a = 5, and ξξξ = . The initial quasi-equilibrium pattern has a pinned spot at the origin ξξξ = , and two additional spots centered at (0 . , T and (0 , . T . Asshown in Fig. 27c, the pinned spot remains at the origin while the other two spots move apart to form an almost colinearpattern. The spot trajectories computed from the full PDE simulation agree well with those from the DAE system. x (a) x-coordinate y (b) y-coordinate −1 0 1−101 (c) t = 1000 Figure 27: Left and middle panels: The spot trajectories for an initial quasi-equilibrium pattern with spots centered at(0 . , . T and (0 . , . T , and with a pinned spot centered at the origin ξξξ = . The full PDE results from (5.1) and DAEdynamics (5.23) and (5.21) are represented by the black dots and red solid line, respectively. Right panel: spot locationsat t = 1000 form a colinear pattern (near the steady-state). Parameters are ε = 0 . a = 15, and a = 5.For our second experiment in the unit disk, we set ε = 0 .
03 and consider a localized feed rate with a = 10 and a = 5,where the concentration point ξξξ moves slowly in time according to (5.11). At time t = 0 the quasi-equilibrium patternconsists of the pinned spot centered at ξξξ (0) = (0 . , T with an additional spot centered at ( − . , T . In Fig. 28 we showa favorable comparison between the spot trajectories obtained from the DAE dynamics (5.23) and (5.21) and from thefull PDE computations of (5.1) on 0 < t < ξξξ ( t ) along a circular trajectory. The other spot moves along a nearly circular trajectory in the unit disk. x (a) x-coordinate y (b) y-coordinate Figure 28: A two-spot quasi-equilibrium pattern with the moving concentration point ξξξ ( ε t ) of the feed rate given in(5.11), and with ε = 0 . a = 10, and a = 5. At time t = 0, there is a pinned spot centered at ξξξ (0) = (0 . , . T and anadditional spot at ( − . , . T . The solid red curve is the spot trajectory from the DAE dynamics (5.23) and (5.21). Thedashed red curve is the moving concentration point ξξξ ( ε t ). The black dots are the locations of the spot and pinned-spot,as computed numerically from the PDE (5.1). We now analyze the linear stability of a quasi-equilibrium pattern v e and u e that consists of spots centered at x , . . . , x N with an additional pinned spot at ξξξ , for which the source strengths are S , . . . , S N and S , respectively. For instabilities35ssociated with non-radially symmetric perturbations near the spots, our previous results in § § S j < Σ (0) ≈ . , for j = 1 , . . . , N and S < min m ≥ Σ m ( a ) . (5.24)Here the symmetry-breaking stability threshold Σ m ( a ) for the local angular mode m was defined in § v = v e + e λt φ and u = u e + e λt η into (5.1), we linearize to get ε ∆ φ − φ + 2 u e v e φ + v e η = λφ , ∆ η − ε − (cid:0) u e v e φ + v e η (cid:1) = λη , in Ω , (5.25)with ∂ n φ = ∂ n η = 0 on ∂ Ω. From the leading-order construction of the quasi-equilibrium pattern in § v e ∼ (cid:40) V (cid:0) ε − | x − ξξξ | (cid:1) , near ξξξ ,V j (cid:0) ε − | x − x j | (cid:1) , near x j , and u e ∼ (cid:40) U (cid:0) ε − | x − ξξξ | (cid:1) , near ξξξ ,U j (cid:0) ε − | x − x j | (cid:1) , near x j , (5.26)Here ( V j , U j ) is the solution to the core problem (2.2) while ( V , U ) is the solution to the new core problem (5.3) nearthe pinned spot, which depends on the feed intensity parameter a .In the inner region near a spot at x j , for j = 1 , . . . , N , we let φ ∼ c j ˜Φ j ( ρ ) and η ∼ c j ˜ N j ( ρ ) in (5.25), where ρ = ε − | x − x j | . Upon using (5.26), we retrieve the inner problem (2.22) for ˜Φ j and ˜ N j for each j = 1 , . . . , N . Similarly,by setting φ ∼ c ˜Φ ( ρ ) and η ∼ c ˜ N ( ρ ) in (5.25), where ρ = ε − | x − ξξξ | , we obtain the following inner problem for thepinned spot: ∆ ρ ˆΦ − ˆΦ + 2 U V ˆΦ + V ˆ N = λ ˆΦ , ∆ ρ ˆ N − U V ˆΦ − V ˆ N = 0 , ρ > , (5.27a)ˆΦ (cid:48) (0) = ˆ N (cid:48) (0) = 0 ; ˆΦ → , ˆ N ∼ log ρ + ˜ B ( S j ; a , λ ) + o (1) , as ρ → ∞ . (5.27b)Here ˜ B ( S j ; a , λ ) depends on the feed intensity a through the pinned core solution ( V , U ). By differentiating (5.3) withrespect to S , and then comparing the resulting system with (5.27) when λ = 0, we identify ˜ B ( S ; a ,
0) = ∂ S χ ( S ; a ).As in § η we first derive the distributional limit ε − (cid:0) u e v e φ + v e η (cid:1) → πc δ ( x − ξξξ ) + 2 π N (cid:88) i =1 c i δ ( x − x i ) , as ε →
0. By using this limit in (5.25), and by enforcing the asymptotic matching condition to the inner solutions nearthe spots, we obtain that the outer problem for η , defined away from the spots, is∆ η − λη − π (cid:34) N (cid:88) i =1 c i δ ( x − x i ) + c δ ( x − ξξξ ) (cid:35) = 0 in Ω , ∂ n η = 0 on ∂ Ω . (5.28a) η ∼ c (cid:16) log | x − ξξξ | + 1 /ν + ˜ B ( S ; a , λ ) (cid:17) , as x → ξξξ , (5.28b) η ∼ c j (cid:16) log | x − x j | + 1 /ν + ˜ B ( S j ; λ ) (cid:17) , as x → x j , j = 1 , . . . , N . (5.28c)For λ (cid:54) = 0, the solution to (5.28a) is represented as η = − πc G λ ( x ; ξξξ ) − π N (cid:88) i =1 c i G λ ( x ; x i ) , (5.29)where G λ is the eigenvalue-dependent Green’s function defined by (2.27). By matching the near-field behavior of (5.29)as x → ξξξ and as x → x j , for j = 1 , . . . , N , to the required singularity behavior in (5.28b) and (5.28c), respectively, wederive a new GCEP for c ≡ ( c , c , . . . , c N ) T , which we write in matrix form as M c ≡ where M ≡ I + 2 πν G λ + ν ˜ B . (5.30a)36ere the entries of the Green’s matrix G λ ∈ R ( N +1) × ( N +1) and the diagonal matrix ˜ B ∈ R ( N +1) × ( N +1) are given by( G λ ) i +1 j +1 = (cid:40) G λ ( x i ; x j ) i (cid:54) = j ,R λ ( x j ; x j ) i = j , for i, j = 0 , . . . , N , ( ˜ B ) = ˜ B ( S ; a , λ ) , ( ˜ B ) j +1 j +1 = ˜ B ( S j ; λ ) , for j = 1 , . . . , N , (5.30b)where, for convenience of notation, we have defined x ≡ ξξξ . We conclude that the N -spot quasi-equilibrium solution withan additional pinned spot at ξξξ is linearly stable on O (1) time-scales to locally radially symmetric perturbations near thespots when there is no root in Re( λ ) > M ( λ ) = 0 . (5.31)Next, we formulate the GCEP for zero-eigenvalue crossings where λ = 0 in (5.28). For λ = 0, the solution to (5.28a) is η = − πc G ( x ; ξξξ ) − π N (cid:88) i =1 c i G ( x ; x j ) + ¯ η , (5.32)where G is the Neumann Green’s function of (2.9), and ¯ η is an undetermined additive constant. By applying the divergencetheorem to (5.28a) we obtain that e T c = 0. Then, by matching the near-field behavior of (5.32) as x → ξξξ and as x → x j ,for j = 1 , . . . , N , to the required singularity behavior in (5.28b) and (5.28c), respectively, and by recalling the identities˜ B ( S ; a ,
0) = ∂ S χ ( S ; a ) and ˜ B ( S j ; 0) = χ (cid:48) ( S j ), we obtain in matrix form that (cid:16) I + 2 πν G + ν ˜ B (cid:17) c = ν ¯ η e , e T c = 0 , (5.33)where c ≡ ( c , c , . . . , c N ) T and e = (1 , . . . , T ∈ R N +1 . Here G is the Neumann Green’s matrix of ξξξ , x , . . . , x N , and˜ B ∈ R ( N +1) × ( N +1) is a diagonal matrix with diagonal entries (cid:16) ˜ B (cid:17) = ∂ S χ ( S ; a ) , (cid:16) ˜ B (cid:17) j +1 j +1 = χ (cid:48) ( S j ) , j = 1 , . . . , N . (5.34)By left-multiplying (5.33) by e T , we use e T c = 0 to calculate that ¯ η = N − (cid:16) π e T G c + e T ˜ B c (cid:17) . By substituting ¯ η backinto the first equation in (5.33) we obtain the following GCEP for detecting zero-eigenvalue crossings: M c = , where M ≡ I + 2 πν ( I − E ) G + ν ( I − E ) ˜ B , (5.35)where E ≡ N − ee T . In summary, a zero-eigenvalue crossing associated with locally radially symmetric perturbations nearthe spots occurs if and only if det M = 0. Since e T c = 0, this criterion detects the initiation of an inter-spot competitioninstability. In this subsection we show a PDE simulation of (5.1) that involves a repeating loop of spot replication and annihilation.We choose ε = 0 .
03 and consider a feed rate with a = 19 and a = 5, where the concentration point ξξξ of the feed evolvesslowly in time according to (5.11). We consider an initial two-spot quasi-equilibrium pattern where the pinned spot isinitially at ξξξ (0) = (0 . , T and with an additional unpinned spot initially centered at ( − . , T . The PDE simulationresults are shown in the right panels of Figs. 29–34 at the indicated times. In the PDE results, we observe that when t ≈
41 the unpinned spot splits into two spots. The resulting two spots evolve dynamically and remain separated fromthe approaching pinned spot. However, as the pinned spot becomes close enough to one of the two unpinned spots, at t ≈
685 a competition instability is triggered and one of these spots is annihilated, leaving a pattern with only one spotand the pinned spot. Later at t ≈ § O (1)time-scale instability of the quasi-equilibrium pattern, we need to augment the DAE solver with a numerical detection37 efore predictedsplitting ( t = 9) After predictedsplitting ( t = 10) v at t = 10Before actualsplitting ( t = 40) After actualsplitting ( t = 41) v at t = 41 Figure 29: Left and middle panels: The spot locations obtained from PDE simulation of (5.1) are represented by theblack circles. The star markers represent the spot locations from the DAE simulation of (5.23) and (5.21). Right panel:contour plot of v at the indicated time, as computed numerically from the PDE. Parameters are ε = 0 . a = 19, a = 5,with the concentration point for the feed rate evolving dynamically by (5.11).38trategy for the initiation of spot-replication or spot-annihilation events, and the subsequent addition or removal of newlycreated or annihilated spots. At the end of each time step in the DAE solver, we first use (5.35) and the conditiondet M = 0 to detect zero-eigenvalue crossings in GCEP. In practice, in our algorithm we identify a zero-eigenvaluecrossing if | det M | ≤ TOL , where TOL (cid:28) . (5.36)This zero-eigenvalue crossing corresponds to an inter-spot competition instability, and triggers the annihilation of the spotwith the smallest source strength. Once the criterion (5.36) is met, we eliminate that particular spot with S j = min ≤ i ≤ N S i from the DAE system. Next, to detect a peanut-splitting instability, we choose a number Σ eff2 that is slightly largerthan the peanut splitting threshold Σ ≈ .
302 for the unpinned spots. If there is an unpinned spot with S j > Σ eff2 , thepeanut-splitting instability triggers a nonlinear spot creation process that divides the spot into two separate spots. Tomodel this process in our algorithm, we replace this spot at x j , with two new spots located at x new = x j + δ ( v , − v ) , x new = x j + δ ( − v , v ) , (5.37)where δ (cid:28) v = ( v , v ) is the normalized velocity field ( | v | = 1 ) as computed from the DAE system (5.23) and(5.21). The choice in (5.37) for this two newly created spot locations is motivated from the result in [17], which showedthat the direction of spot-splitting is perpendicular to the direction of motion of the spot.We use this algorithm for augmenting the DAE solver with parameters TOL = 10 − , Σ eff2 = 4 . δ = 0 .
1, withresults shown in the left and middle panels in Figs. 29– 34. At t = 9, the algorithm detects a peanut-splitting instabilityof the spot. The spot is replaced with two spots given in (5.37) (see Fig. 29). Then at t = 583, the DAE solver detectsa zero-eigenvalue crossing based on (5.36). The spot with the minimum source strength is removed from the DAEsystem. Right after the removal, the DAE solver detects a peanut-splitting instability, and two new spots are created. Inconclusion, the DAE solver predicts a spot creation-annihilation event at t = 583. The PDE simulation confirms a spotannihilation event at t ≈ t ≈ t = 1248 and t = 1952, respectively. These events are all confirmed by the full PDEsimulation (see Fig. 32 and Fig. 33). We have developed a hybrid asymptotic-numerical theory to analyze the effect of several types of localized heterogeneitieson the existence, linear stability, and slow dynamics of spot patterns for the prototypical two-component Schnakenbergmodel (1.1) in a bounded 2-D domain. Our analysis has focused on distinct types of localized heterogeneities: a stronglocalized perturbation of a spatially uniform feed rate and the effect of removing a small hole in the domain, through whichthe chemical species can leak out. Although our overall approach relies on the theoretical framework first introduced in[17], and later extended in [33] for analyzing the effect of heterogeneities in the Brusselator model, our analysis of localizedheterogeneities for the Schnakenberg model has revealed a wide range of novel phenomena such as, saddle-node bifurcationsfor quasi-equilibrium spot patterns that otherwise would not occur for a homogeneous medium, a new type of spot solutionpinned at the concentration point of the feed rate, spot self-replication behavior that generates more than two new spots,and the existence of a creation-annihilation attractor with at most three spots. The hybrid approach presented herein canbe readily extended to other well-known RD models such as the Gray-Scott, Gierer-Meinhardt, and Brusselator models.We conclude by briefly discussing a few problems that warrant further investigation. One interesting direction would beto extend the algorithm, introduced in § N -spot quasi-equilibrium patterns over much longer time scales.This would involve coupling the ODE-DAE system for slow spot dynamics occuring on long O ( ε − ) time scales withsudden “surgeries”, resulting in either spot-creation or spot-annihilation events, that are informed through monitoringlinear stability thresholds at each time step as the quasi-equilibrium spot pattern evolves. In particular, in this generalsetting, it would be interesting to classify whether spot-annihilation events, triggered by a competition instability dueto a zero-eigenvalue crossing of the globally coupled eigenvalue problem, can be interpreted more geometrically in termsof crossing through a saddle-node bifurcation point of manifolds of spot quasi-equilibria. Such manifolds depend on theinstantaneous spatial configuration of spots, and they evolve slowly in time. For a two-spot pattern in the unit disk,and with a feed rate that is slowly ramped in time, such a fold-point crossing was observed in Fig. 2. In a 1-D setting,39 efore predicted killing andsplitting ( t = 583) After predicted killing andsplitting ( t = 584) v at t = 584Before actual killing ( t = 684) After actual killing ( t = 685) v at t = 684Before actual splitting ( t = 719) After actual splitting ( t = 720) v at t = 720 Figure 30: Our approximation algorithm detects spot-annihilation and spot-splitting at the same time. On the otherhand, the spot-tracking algorithm based on the PDE simulation detects killing and splitting, separately at later time.
After one cycle ( t = 1112) v at t = 1112 Figure 31: When the pinned-spot finishes a full cycle, the predicted spot location from the DAE and the actual spotlocation (from the PDE spot-tracking algorithm) have a good agreement. This shows that the spots can catch up withthe prediction from the augmented DAE algorithm. 40 efore predicted killing andsplitting ( t = 1248) After predictedkilling and splitting( t = 1249) v at t = 1249Before actual killing ( t = 1350) After actual killing ( t = 1351) v at t = 1350Before actual splitting( t = 1384) After actual splitting ( t = 1385) v at t = 1385 Figure 32: The second spot creation-annihilation event.41 efore predicted killing andsplitting ( t = 1952) After predicted killing andsplitting ( t = 1953) v at t = 1953Before actualkilling ( t = 2016) After actualkilling ( t = 2017) v at t = 2016Before actualsplitting ( t = 2049) After actualsplitting ( t = 2050) v at t = 2050 Figure 33: The third spot creation-annihilation event t = 2500 v at t = 2500 Figure 34: The approximation algorithm provides a reasonably close prediction for the unpinned spot trajectory at t = 2500 after three spot creation-annihilation events. 42pike-annihilation events have been recently interpreted in [4] for the extended Klausmeir RD model as arising from rapidtransitions between manifolds of spike quasi-equilibria as the pattern evolves.Another open problem is to determine whether a creation-annihilation attractor for spot quasi-equilibria, which involvesonly a few spots, can occur for a time-independent feed rate that has a smooth (not localized) spatial variation. For the1-D Schnakenberg model, but for a very large number of spikes, such a creation-annihilation attractor has been predictedand observed in [18] through the analysis of a limiting mean-field equation for the spike density.Finally, it would be interesting to study how variations in the domain geometry or the domain boundary conditioninfluence the slow spot dynamics and the linear stability properties of quasi-equilbrium spot patterns, leading to to newtypes of spot-pinning behavior. Some work in this direction for the Brusselator model with a Robin boundary condition isgiven in [33] for the disk. For general domains, the determination of the Neumann Green’s function and the reduced-waveGreen’s function would be central to this study. For an elliptical domain of arbitrary eccentricity, spot-pinning behaviorcan readily be studied using the new analytical result in [12] for the Neumann Green’s function for the ellipse. Acknowledgements
Tony Wong was supported by a UBC Four-Year Graduate Fellowship. Michael Ward gratefully acknowledges the financialsupport from the NSERC Discovery Grant program. We thank Prof. Colin B. MacDonald for helpful suggestions regardingthe numerical PDE computations.
References [1] M. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, and G. N.Wells. The fenics project version 1.5.
Archive of Numerical Software , 3(100), 2015.[2] D. Avitabile, V. F. Brena, and M. J. Ward. Spot dynamics in a reaction-diffusion model of plant root hair initiation.
SIAM J. Appl. Math. , 78(1):291–319, 2018.[3] R. Bastiaansen, M. Chirilus-Bruckner, and A. Doelman. Pulse solutions for an extended Klausmeier model withspatially varying coefficients.
SIAM J. Appl. Dyn. Sys. , 19(1):1–57, 2020.[4] R. Bastiaansen and A. Doelman. The dynamics of disappearing pulses in a singularly perturbed reaction-diffusionsystem with parameters that vary in space and time.
Physica D , 388:45–72, 2019.[5] V. F. Brena, A. Champneys, C. Grierson, and M. J. Ward. Mathematical modeling of plant root hair initiation:Dynamics of localized patches.
SIAM J. Appl. Dyn. Sys. , 13(1):210–248, 2014.[6] C. N. Chen, S. I. Ei, and S. Tzeng. Heterogeneity-induced effects for pulse dynamics in Fitzhugh-Nagumo-typesystems.
Physica D , 382-383(1):22–32, 2018.[7] 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.[8] H. Dankowicz and F. Schilder.
Recipes for continuation , volume 11 of
Computational Science & Engineering . Societyfor Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2013.[9] P. Davies, P. Blanchedeau, E. Dulos, and P. De Kepper. Dividing blobs, chemical flowers, and patterned islands ina reaction- diffusion system.
J. Phys. Chem. A , 102(43):8236–8244, 1998.[10] A. Doelman, P. van Heijster, and J. Shen. Pulse dynamics in reaction-diffusion systems with strong spatially localizedimpurities.
Phil. Trans. Roy. Soc. A. , 376:20170183, 2018.[11] H. Ikeda and S. I. Ei. Front dynamics in heterogeneous diffusive media.
Physica D. , 239:1637–1649, 2010.[12] S. Iyaniwura, T. Wong, C. B. MacDonald, and M. J. Ward. Optimization of the mean first passage time in near-diskand elliptical domains in 2-D with small absorbing traps.
SIAM Review , 2020, submitted.[13] L. T. J.C. Tzou. Spot patterns of the Schnakenberg reaction-diffusion system on a curved torus.
Preprint , 2019.4314] E. Knobloch. Spatial localization in dissipative systems.
Annu. Rev. Cond. Mat. Phys. , 6:325–359, 2015.[15] T. Kolokolnikov, M. S. Titcombe, and M. J. Ward. Optimizing the fundamental Neumann eigenvalue for the Laplacianin a domain with small traps.
European J. of Appl. Math. , 16(2):161–200, 2005.[16] T. Kolokolnikov, M. Ward, J. Tzou, and J. Wei. Stabilizing a homoclinic stripe.
Phil. Trans. Roy. Soc. A. ,376(2135):20180110, 2018.[17] T. Kolokolnikov, M. J. Ward, and J. Wei. Spot self-replication and dynamics for the Schnakenburg model in atwo-dimensional domain.
J. Nonlinear Science , 19(1):1–56, 2009.[18] T. Kolokolnikov and J. Wei. Pattern formation in a reaction-diffusion system with space-dependent feed rate.
SIAMReview , 60(3):626–645, 2018.[19] T. Kolokolnikov and J. Wei. Hexagonal spike clusters for some pde’s in 2-d.
DCDS-B , 25(10):4057–4070, 2020.[20] T. Kolokolnikov and S. Xie. Spike density distribution for the Gierer-Meinhardt model with precursor.
Physica D ,31:132247, 2019.[21] T. Kolokolniov, F. Paquin-Lefebvre, and M. J. Ward. Stable asymmetric spike equilibria for the Gierer-Meinhardtmodel with a precursor field.
IMA J. of Appl. Math , 2020, (to appear).[22] A. L. Krause, V. Klika, T. E. Woolley, and E. A. Gaffney. Heterogeneity induces spatiotemporal oscillations inreaction-diffusion systems.
Phys. Rev. E , 97(5):052206, 2018.[23] A. L. Krause, V. Klika, T. E. Woolley, and E. A. Gaffney. From one pattern into another: Analysis of Turing patternsin heterogeneous domains via WKBJ.
J. Roy. Soc. Interface , 17:20190621, 2020.[24] V. Kurella, J. C. Tzou, D. Coombs, and M. J. Ward. Asymptotic analysis of first passage time problems inspired byecology.
Bull. Math. Biol. , 77(1), 2015.[25] K.-J. Lee, W. D. McCormick, J. E. Pearson, and H. L. Swinney. Experimental observation of self-replicating spotsin a reaction–diffusion system.
Nature , 369(6477):215–218, 1994.[26] K. J. Lee and H. L. Swinney. Lamellar structures and self-replicating spots in a reaction-diffusion system.
Phys. Rev.E , 51(3):1899, 1995.[27] K. Nishi, Y. Nishiura, and T. Teramoto. Dynamics of two interfaces in a hybrid system with jump-type heterogeneity.
Japan J. Indust. Appl. Math. , 30:351–395, 2013.[28] Y. Nishiura, T. Teramoto, and X. Yuan. Heterogeneity-induced spot dynamics for a three-component reaction-diffusion system.
Comm. Pure and Appl. Anal. , 11(1):307–338, 2012.[29] K. Page, P. K. Maini, and N. A. M. Monk. Pattern formation in spatially heterogeneous turing reaction-diffusionmodels.
Physica D , 181(1-2):80–101, 2003.[30] K. Page, P. K. Maini, and N. A. M. Monk. Complex pattern formation in reaction-diffusion systems with spatiallyvarying parameters.
Physica D , 202(1-2):95–115, 2005.[31] I. Rozada, S. J. Ruuth, and M. Ward. The stability of localized spot patterns for the Brusselator on the sphere.
SIAM J. Appl. Dyn. Sys. , 13(1):564–627, 2014.[32] P. H. Trinh and M. J. Ward. The dynamics of localized spot patterns for reaction-diffusion systems on the sphere.
Nonlinearity , 29(3):766, 2016.[33] J. C. Tzou and M. J. Ward. The stability and slow dynamics of spot patterns in the 2D Brusselator model: Theeffect of open systems and heterogeneities.
Physica D , 373:13–37, 2018.[34] P. van Heijster, C. N. Chen, Y. Nishiura, and T. Teramoto. Pinned solutions in a heterogeneous three-componentFitzHugh-Nagumo model.
Journal of Dynamics and Differential Equations , 31(1):153–203, 2019.[35] P. van Heijster, A. Doelman, T. J. Kaper, Y. Nishiura, and K. I. Ueda. Pinned fronts in heterogeneous media ofjump type.
Nonlinearity , 24(1):127–157, 2010.[36] V. K. Vanag and I. R. Epstein. Localized patterns in reaction-diffusion systems.
Chaos , 17(3):037110, 2007.[37] M. J. Ward, W. D. Henshaw, and J. B. Keller. Summing logarithmic expansions for singularly perturbed eigenvalueproblems.
SIAM J. Appl. Math. , 53(3):799–828, 1993. 4438] M. J. Ward and J. B. Keller. Strong localized perturbations of eigenvalue problems.
SIAM J. Appl. Math. , 53(3):770–798, 1993.[39] M. J. Ward, D. McInerney, H. P., G. D., and P. Maini. The dynamics and pinning of a spike for a reaction-diffusionmodel.
SIAM J. Appl. Math. , 62(4):1297–1328, 2002.[40] J. Wei and M. Winter. Asymmetric spotty patterns for the Gray–Scott model in R . Studies in Appl. Math. ,110(1):63–102, 2003.[41] J. Wei and M. Winter. Existence and stability of multiple-spot solutions for the Gray–Scott model in R . PhysicaD: Nonlinear Phenomena , 176(3-4):147–180, 2003.[42] J. Wei and M. Winter. Stationary multiple spots for reaction–diffusion systems.
Journal of mathematical biology ,57(1):53–89, 2008.[43] J. Wei and M. Winter. On the Gierer-Meinhardt system with precursors.
DCDS-A , 25(1):363–398, 2009.[44] J. Wei and M. Winter. Spikes for the Gierer-Meinhardt system with discontinuous diffusion coefficients.
J. NonlinearSci. , 12(3):301–339, 2009.[45] J. Wei and M. Winter. Stable spike clusters for the one-dimensional Gierer-Meinhardt system.
European J. of Appl.Math. , 28(4):576–635, 2017.[46] T. Wong and M. J. Ward. Weakly nonlinear analysis of peanut-shaped deformations for localized spots of singularlyperturbed reaction-diffusion systems.
SIAM J. Appl. Dyn. Sys. , 2020, to appear.[47] S. Xie and T. Kolokolnikov. Moving and jumping spot in a two-dimensional reaction–diffusion model.
Nonlinearity ,30(4):1536, 2017.[48] X. Yuan, T. Teramoto, and Y. Nishiura. Heterogeneity-induced defect bifurcation and pulse dynamics for a three-component reaction-diffusion system.
Phys. Rev. E. , 75(036220), 2007.
A The Green’s functions for the unit disk
The Neumann Green’s function and its regular part, satisfying (2.9), have explicit formulae for the unit disk (cf. [17]): G ( x ; z ) = 12 π (cid:18) − log | x − z | − log (cid:12)(cid:12)(cid:12)(cid:12) | z | x − | z | z (cid:12)(cid:12)(cid:12)(cid:12) + 12 ( | x | + | z | ) − (cid:19) , (A.1a) R ( z ; z ) = 12 π (cid:18) − log (cid:0) − | z | (cid:1) + | z | − (cid:19) . (A.1b)Their gradients are given by ∇ x G = − π (cid:32) ( x − z ) | x − z | + | z | (cid:0) | z | x − z (cid:1)(cid:12)(cid:12) | z | x − z (cid:12)(cid:12) − x (cid:33) , ∇ x R = 12 π (cid:18) − | z | − | z | (cid:19) z . (A.2)For a ring pattern, where x , . . . , x N are equally-spaced on a ring of radius r concentric within the unit disk as givenin (2.39), we have from Proposition 4.3 of [15] that G e = p ( r ) N e , p ( r ) ≡ π (cid:18) − N log( N r N − ) − N log(1 − r N ) + r N − N (cid:19) . (A.3)As such, for a ring pattern, there is a symmetric solution to the NAS (2.14) given by S j = S c = p a /N for j = 1 , . . . , N ,where p a is defined in (2.7). Then, upon defining ∇ x R j,j ≡ ∇ x R ( x ; x j ) | x = x j , and ∇ x G j,i ≡ ∇ x G ( x ; x i ) | x = x j , we thenuse the reciprocity property of the Green’s function to calculate βββ j in (2.12) as βββ j = 2 πS c ∇ x R j,j + N (cid:88) i (cid:54) = j ∇ x G j,i = 2 πS c (cid:18) p (cid:48) ( r )2 N (cid:19) e θ j = S c (cid:32) − N − r + N r N − − r N + N r (cid:33) e θ j , (A.4)45here e θ j ≡ (cos θ j , sin θ j ) T and θ j = 2 π ( j − /N . By substituting (A.4) in (2.18), and using x j = r ( σ ) e θ j , we obtainthe scalar ODE (2.40) for the ring radius r .For the unit disk, and for λ (cid:54) = 0, the eigenvalue-dependent Green’s function G λ ( x ; x ), as defined by (2.27), can beexpressed as an infinite series as (cf. Appendix A.1 of [7]) G λ ( x ; x ) = 12 π (cid:20) K ( θ λ | x − x | ) − K (cid:48) ( θ λ ) I (cid:48) ( θ λ ) I ( θ λ r ) I ( θ λ r ) (cid:21) − π ∞ (cid:88) n =1 cos [ n ( ψ − ψ )] K (cid:48) n ( θ λ ) I (cid:48) n ( θ λ ) I n ( θ λ r ) I n ( θ λ r ) . (A.5a)Here x = r (cos( ψ ) , sin( ψ )) , x = r (cos( ψ ) , sin( ψ )), I n and K n are the n th order modified Bessel functions of the firstand second kind, respectively, and θ λ is the principal branch of θ λ ≡ (cid:112) τ λ/D . The regular part of G λ is R λ ( x ; x ) = 12 π (cid:20) log 2 − γ e − log( D/τ )2 − log λ − K (cid:48) ( θ λ ) I (cid:48) ( θ λ ) I ( θ λ r ) (cid:21) − π ∞ (cid:88) n =1 K (cid:48) n ( θ λ ) I (cid:48) n ( θ λ ) I n ( θ λ r ) , (A.5b)where γ e ≈ . B Spectrum of circulant matrices