Stable Boundary Spike Clusters for the Two-Dimensional Gierer-Meinhardt System

We consider the Gierer-Meinhardt system with small inhibitor diffusivity and very small activator diffusivity in a bounded and smooth two-dimensional domain. For any given positive integer $k$ we construct a spike cluster consisting of $k$ boundary spikes which all approach the same nondegenerate local maximum point of the boundary curvature. We show that this spike cluster is linearly stable. The main idea underpinning these stable spike clusters is the following: due to the small inhibitor diffusivity the interaction between spikes is repulsive and the spikes are attracted towards a nondegenerate local maximum point of the boundary curvature. Combining these two effects can lead to an equilibrium of spike positions within the cluster such that the cluster is linearly stable.


Introduction
Turing in his pioneering work in 1952 [14] proposed that a patterned distribution of two chemical substances, called the morphogens, could trigger the emergence of cell structures. He also gives the following explanation for the formation of the morphogenetic pattern: It is assumed that one of the morphogens, the activator, diffuses slowly and the other, the inhibitor, diffuses much faster. In the mathematical framework of a coupled system of reaction-diffusion equations with hugely different diffusion coefficients it is shown by linear stability analysis that the homogeneous state may be unstable. In particular, a small perturbation of spatially homogeneous initial data may evolve to a stable spatially complex pattern of the morphogens.
Since the work of Turing, many different reaction-diffusion system in biological modelling have been proposed and the occurrence of pattern formation has been investigated by studying what is now called Turing instability. One of the most popular models in biological pattern formation is the Gierer-Meinhardt system [5], see also [10]. In two dimensions in a special case after rescaling it can be stated as follows: (1.1) The unknown functions u = u(x, t) and v = v(x, t) represent the concentrations of the activator and inhibitor, respectively, at the point x ∈ Ω ⊂ R 2 and at a time t > 0. Here ∆ is the Laplace operator in R 2 , Ω is a smooth bounded domain in R 2 , ν = ν(x) is the outer unit normal at x ∈ ∂Ω. Throughout this paper, we assume that 0 < ε ≪ 1, 0 < D ≪ 1, (1.2) τ ≥ 0 is a fixed constant independent of ε, D and x. Further, the diffusivities ε and D do not depend on x but they are both small constants. In this paper, we further assume that This means that ǫ is much smaller than D. On the other hand, ǫ cannot be exponentially small compared to √ D. In this paper, we study the Gierer-Meinhardt system in a bounded and smooth two-dimensional domain. We prove the existence and stability of a cluster consisting of k boundary spikes near a nondegenerate local maximum point P 0 of the boundary curvature h(P ).
The main idea underpinning these stable spike clusters is the following: due to the small inhibitor diffusivity the interaction between spikes is repulsive and the spikes are attracted towards a nondegenerate local maximum point of the boundary curvature. Combining these two effects can lead to an equilibrium of spike positions within the cluster such that the cluster is linearly stable.
The repulsive nature of spikes has been shown in [4]. Existence and stability of a spike cluster made up of two boundary spikes has been established in [3].
Before we state our main results, let us mention some previous ones concerning various regimes for the asymptotic behavior of D.
For the strong coupling case, i.e. D ∼ 1, the second and third authors constructed singleinterior spike solutions [20]. In [22], they continued the study, and proved the existence of solutions with k interior spikes.
Moreover, it is shown that this solution is linearly stable for τ = 0. For the weak coupling case D → ∞, in [21] the second and third authors proved the existence of multiple interior spike solutions. Further, they showed that there are stability thresholds D 1 (ε) > D 2 (ε) > · · · > D k (ε) > · · · such that if lim ε→0 D k (ε) D > 1, the k-peak solution is stable and if lim ε→0 D k (ε) D < 1, the kpeak solution is unstable. Multiple spikes for the Gierer-Meinhardt system in a one-dimensional interval have been studied in [7,13,23] and on the real line in [2].
In [17] the existence, uniqueness and spectral properties of a boundary spike solution have been studied for the shadow Gierer-Meinhardt system (i.e. after formally taking the limit D → ∞.) In [24] the existence and stability of N −peaked steady states for the Gierer-Meinhardt system with precursor inhomogeneity has been explored. The spikes in the patterns can vary in amplitude. In particular, the results imply that a precursor inhomogeneity can induce instability. Single-spike solutions for the Gierer-Meinhardt system with precursor including spike dynamics have been studied in [15].
Previous results on stable spike clusters include a stable spike cluster for a consumer chain model [25] and a stable spike cluster for the one-dimensional Gierer-Meinhardt system with precursor inhomogeneity [27].
For more background, modelling, analysis and computation on the Gierer-Meinhardt system, we refer to [26] and the references therein.
Theorem 2.1. Let k be a positive integer, and P 0 be a nondegenerate local maximum point of the curvature h(P ) of the boundary ∂Ω.
Further, we have Remark 2.1. The spike cluster is established by a balance between repelling spikes and attracting boundary point of local maximum curvature.
Theorem 2.2. The k-boundary spike cluster solution given in Theorem 2.1 is linearly stable if τ is small enough.
Remark 2.2. There are eigenvalues of two different orders: n − 1 eigenvalue related to repelling of neighbouring spikes are of order ǫ 3 log ξσ ǫD , and one eigenvalue stemming from the curvature of the boundary (corresponding to synchronous motion of all spikes) is of order ǫ 3 .
This paper is organised as follows. In sections 3-5 we show existence of the spike cluster steady state by using Liapunov-Schmidt reduction. In section 3 we introduce an approximation to the spike cluster steady state. In section 4 we use the Liapunov-Schmidt method to reduce the problem to finite dimensions. In section 5 we solve this reduced problem. In sections 6-7 we study the stability of this spike cluster steady state. In section 6 we consider large eigenvalues. Finally, in section 7 we study small eigenvalues. In two appendices we show some technical results: in appendix A (section 8) we prove Proposition 4.1 and in appendix B (section 9) we compute the small eigenvalues.

Introduction of the approximate solutions
In this paper, we consider the following problem: where Ω ⊂ R 2 is a bounded and smooth two-dimensional domain. Let P 0 be a nondegenerate maximum point of the boundary curvature h(P ) on the boundary of Ω. For P ∈ ∂Ω, where ∂ ∂τ (P ) denotes the tangential derivative with respect to P at P ∈ ∂Ω. We will sometimes drop the variable P if this can be done without causing confusion.
In this section, we construct an approximation to a spike cluster solution to (3.1) which concentrates at P 0 .
The approximate cluster consists of spikes σ −2 ξ σ,i w( x−P i ε ) which are centred at the points P i for i = 1, · · · , k, where σ = ε √ D and the amplitude ξ σ,i satisfies (see (3.38)). Let P 1 , · · · , P k be k points distributed along the boundary ∂Ω such that we have for i = 2, · · · , k where ν 2 is given in (5.10) below. Further, η > 0 is a small constant independent of ε and D. The reason for assuming (3.2) and (3.3) will become clear in Section 5 when we solve the reduced problem.
Remark 3.1. By (3.2), the distance of neighbouring spikes satisfies Since we want to construct multiple boundary spikes which collapse at one point, we require that assumptions (1.2) and (1.3) hold.
After re-scaling,û if we drop the hat and still denote solutions by (u, v), equation (3.1) is equivalent to where Ω ε = ε −1 Ω.
From now on we will deal with (3.4). Before introducing an approximation to the spike solutions, we first define some notation.
3.1. The analysis of the projection P w q (z − q). Before calculating the heights of the spikes, we need some preliminaries of the projection P w p i (z − p i ) defined in (3.7) which are rather standard by now. Some of these results have been derived in [18,19]. Let P ∈ ∂Ω. We define a diffeomorphism straightening the boundary. We may assume that the inward normal to ∂Ω at P is pointing in the direction of the positive x 2 axis. Denote B ′ (R) = {x ∈ R 2 ||x 1 | ≤ R}. Then since ∂Ω is smooth, we can find a constant R such that ∂Ω can be represented by the graph of a smooth function ρ P : B ′ (R) → R, where ρ P (0) = 0, and ρ ′ P (0) = 0. From now on, we omit the use of P in ρ P and write ρ if this can be done without causing confusion. So near P , ∂Ω can be represented by (x 1 , ρ(x 1 )). The curvature of ∂Ω at p is h(P ) = ρ ′′ (0). Let After rescaling, it follows that near p = P ε , the boundary ∂Ω ε can be represented by (z 1 − p 1 , ε −1 ρ(ε(z 1 − p 1 ))), where (z 1 , z 2 ) = ε −1 (x 1 , x 2 ) and p = (p 1 , p 2 ). By Taylor expansion, we have For z ∈ Ω 1ε = 1 ε Ω 1 we set Under this transformation, the Laplace operator and the boundary derivative operator become Let v (1) be the unique solution of where R 2 + is the upper half plane, namely R 2 (3.12) Note that v (1) , v (2) are even functions in y 1 and v (3) is an odd function in y 1 . Moreover, it is easy to see that |v i (y)| ≤ Ce −µ|y| for any 0 < µ < 1. Let χ be a smooth cut-off function, such that χ(a) = 1 for a ∈ B(0, R 0 √ D log ξσ εD ), and χ(a) = 0 for x ∈ B(0, 2R 0 √ D log ξσ εD ) c for some suitable R 0 such that |p i − p 0 | < R 0 σ log ξσ ǫD , and (3.14) Then we have the following estimate: Proposition 3.1 was proved in [19] by Taylor expansion including a rigorous bound for the remainder using estimates for elliptic partial differential equations. Moreover, it has been shown that that |ξ p (z)| ≤ Ce −µ|z−p| for any 0 < µ < 1.
Similarly we know from [19] that [

16)
where η is the unique solution of the following equation: (3.17) It follows that η(y) is an odd function in y 1 . It can be seen that |η 1 (y)| ≤ Ce −µ|y| for some 0 < µ < 1.
Remark 3.2. In the following sections, we will denote by y i = (y i 1 , y i 2 ) the transformation defined by (3.10) centred at the point p i and let v For D > 0, let G D (x, y) be the Green's function given by where y ∈ ∂Ω √ D and let G 0 be the Green's function of the upper half plane: (3.23) and let η 2 be the solution of It can be seen easily that η 1 is even in y 1 and η 2 is odd in y 1 . Then one can get the following result. Proof. First we compute on ∂Ω √ D , for any 0 < l < 1.
Since we have From the expansion above, we can get the asymptotic behaviour of G √ D .
Next we have the following expansion of G 0 : The following expansion of G 0 holds: for 0 < r < 1, where ψ is a smooth function with ψ(0) = ψ ′ (0) = 0 and c 1 , c 2 are universal constants.
Proof. By an even extension in y 2 , one can get the Green's function in the whole space R 2 . For the expansion of fundamental solution, see Lemma 4.1 in [1]. Then the above expansion follows.
We set From the estimates above, and for points p ε ∈ Λ k , we have Generally, we have For the derivatives, we estimate Calculating the heights of the spikes. In this section, we are going to determine the heights of spikes ξ σ,i to leading order. In the sequel, by T [h] we denote the unique solution of the equation (3.32) As mentioned before, we will choose the approximate solution to be First we calculate the heights of the peaks: Thus We assume that the heights of the spikes are asymptotically equal as ε, D → 0, i.e.
Then we get that The analysis in this subsection calculates the heights of the spikes under the assumption that their shape is given. In the next two sections, we provide the rigorous proof for the existence.

Existence I: Reduction to finite dimensions
Let us start to prove Theorem 2.1. The first step is choosing a good approximate solution which was done in (3.6). The second step is using the Liapunov-Schmidt method to reduce the problem to a finite dimensions which we do in this section. The last step is solving the reduced problem which will be done in Section 5.
First we need to calculate the error terms caused by the approximate solution given in (3.6) to show that this is a good choice: where we have used the notation On the other hand, we calculate for Thus we can get that By the above estimates, we have the following key estimate: The above estimates will be very useful in the existence proof using the Liapunov-Schmidt reduction. In particular, they will imply an explicit formula for the positions of the spikes in Section 5. Now we study the linearised operator defined by and define the approximate kernels by , and choose the approximate cokernels as follows: We then define p denote the orthogonal complements with the scalar product of L 2 (Ω ε ). Let π ε,p denote the projection in L 2 (Ω ε ) × L 2 (Ω ε ) onto C ⊥ ε,p . We are going to show that the equation As a preparation we state a result on the invertibility of the corresponding linearised operator L ε,p whose proof is postponed to Appendix A.
Proposition 4.1. There exist positive constantsδ, C such that for max{σ, D} <δ, the map L ε,p is surjective for arbitrary p ∈ Λ k . Moreover the following estimate holds: Now we are in the position to solve the equation Since L ε,p | K ⊥ ε,p is invertible, we can write the above equation as where Σ = φ ψ and We are going to show that the operator M ε,p is a contraction mapping on If we choose C 0 large enough, then M ε,p is a contraction mapping on B ε . The existence of a fixed point Σ ε,p together with an error estimate now follows from the contraction mapping principle. Moreover, Σ ε,p is a solution. Thus we have proved Lemma 4.2. There existsδ > 0 such that for every triple (ε, D, p) with max{σ, D} <δ, and p ∈ Λ k , there exists a unique (φ ε,p , ψ ε,p ) ∈ K ⊥ ε,p satisfying and More refined estimates for φ ε,p are needed. We recall that from Lemma 4.1 that S 1 (U, V ) can be decomposed into two parts S 1,1 and S 1,2 , where S 1,1 is in leading order an even function in z 1 and S 1,2 is in leading order an odd function in z 1 . Similarly we can decompose φ ε,p .
where φ ε,p,1 is even in z 1 which can be estimated by 18) and φ ε,p,2 can be estimated by Proof. Let S(u) = S 1 (u, T (u)). We first solve for φ ε,p,1 ∈ K ⊥ ε,p . Then we solve for φ ε,p,2 ∈ K ⊥ ε,p . Using the same proof as in Lemma 4.2, the above two equations are uniquely solvable for max{σ, D} ≪ 1. By uniqueness, it is easy to see that φ ε,p,1 and φ ε,p,2 have the required properties.

Existence proof II: The reduced problem
In this section, we solve the reduced problem. This completes the proof for our main existence result given in Theorem 2.1.
By Lemma 4.2, for every p ∈ Λ k , there exists a unique solution (φ, ψ) ∈ K ⊥ ε,p such that To this end, we calculate the projection: where I 1 , I 2 are defined by the last equality.
For I 1 , we have Note that by the estimates satisfied by φ in Lemma 4.3, we have where we have used the estimates (3.28)-(3.31). Further, we have We compute By (3.12), we have where the constant ν 1 > 0 is defined by Now by (5.2)-(5.6), Next we estimate I 12 : where we have used √ D log ξσ ǫD ≪ 1 and Thus by (5.8) and (5.9), we have Next we estimate I 2 : By the equation for ψ, we have ∆ψ − σ 2 ψ + 2U φ + φ 2 = 0 and therefore where R(z) is an even function in z 1 .
Thus we have Combining the estimates for I 1 and I 2 , (5.11) and (5.12), we have Thus W ε,i = 0 is reduced to the following system: (5.13) We first solve the limiting case: (5.14) This system is uniquely solvable with where we have used the notation h ′′ (p 0 ) = ∂ 2 ∂τ 2 h(P 0 ) < 0. From (5.16), we know To find p i such that W ε,i = 0, we expand p i = p 0 i +p i . Then adding the first i equations, we have One can get thatp By the last equation From the equation above, we can estimatep 1 bỹ In conclusion, we solve W ε,i = 0 withp Thus we have proved the following proposition: Finally, we complete the proof of Theorem 2.1. Proof. By Proposition 5.1, there exists P ε → P 0 , such that W ε (p ε ) = 0. In other words, we have Moreover, by the maximum principle, (U, V ) > 0 and the solution satisfies all the properties of Theorem 2.1.

Study of the large eigenvalues
We consider the stability of the steady-state (u, v) constructed in Theorem 2.1.
Linearizing the system around the equilibrium states (u, v) obtained in Theorem 2.1, we obtain the following eigenvalue problem: . In this section, since we study the large eigenvalues, we may assume that |λ ε | ≥ c > 0 for max{σ, D} small enough. If Re(λ ε ) ≤ −c < 0, then λ ε is a stable large eigenvalue, we are done. Therefore, we may assume that Re(λ ε ) ≥ −c and for a subsequence max{σ, D} → 0, λ ε → λ 0 = 0. We shall derive the limiting eigenvalue problem which is given by a coupled system of NLEPs.
Since φ ε H 2 (Ωε) = 1, φ ε,j H 2 (R 2 ) ≤ C . By taking a subsequence, we may assume that φ ε,j → φ j as max{σ, D} → 0 in H 1 (R 2 ) for some φ j ∈ H 1 (R 2 ) for j = 1, · · · , k. By (6.1), we have Substituting the above equation into the first equation of (6.1), letting max{σ, D} → 0, and using the expansion of ξ σ,j , we arrive at the following nonlocal eigenvalue problem (NLEP): By Theorem 3.5 in [26], (6.5) has only stable eigenvalues if τ is small enough. In conclusion, we have shown that the large eigenvalues of the k-peaked solutions given in Theorem 2.1 are all stable if τ is small enough.

Study of the small eigenvalues
Now we study the eigenvalue problem (6.1) with respect to small eigenvalues. Namely, we assume that λ ε → 0 as max{σ, D} → 0. We will show that the small eigenvalues in leading order are related to the matrix M(p 0 ) given in (7.3) which is computed from the Green's function. Our main result in this section says that if λ ε → 0, then in leading order where σ 0 is an eigenvalue of M(p 0 ). We will show that all the eigenvalues of M(p 0 ) have negative real part provided that the eigenvector is orthogonal to (1, 1, ..., 1) T . However, for the eigenvector (1, 1, ..., 1) T the eigenvalue of M(p 0 ) is zero, the leading order term in the eigenvalue expansion vanishes and the next order term is needed to prove stability. To establish it we have to compute the contribution from the boundary curvature. It follows that for a local maximum point of the boundary curvature this eigenvalue has negative real part. Whereas the Green's function part is of order ǫ 3 log ξσ ǫD , the part from the boundary curvature is of order ǫ 3 . Thus the small eigenvalues of (6.1) are all stable.
To compute the small eigenvalues, we need to expand the spike cluster solution to higher order. Then we expand the eigenfunction and compute the small eigenvalues. This will be done in Appendix B. The key estimates are given in Lemma 9.1.
We compute the small eigenvalues using Lemma 9.1. Comparing l.h.s and r.h.s, we obtain where ν 2 has been defined in (5.10). Further, we have Using the estimates for p 0 i in (5.15) and (5.16), we have This shows that if all the eigenvalues of M(p 0 ) have positive real part, then the small eigenvalues are stable. On the other hand, if M(p 0 ) has eigenvalues with negative real part, then there are eigenfunctions and eigenvalues to make the system unstable. Next we study the spectrum of the k × k matrix A defined by a s,s = (s − 1)(k − s + 1) + s(k − s), s = 1, · · · , k, a s,s+1 = a s+1,s = −s(k − s), s = 1, · · · , k, a s,l = 0, |s − l| > 1.
We have the following result from Lemma 16 in [27]: Lemma 7.1. The eigenvalues of the matrix A are given by λ n = n(n + 1), n = 0, · · · , k − 1. Equation (7.2) shows that the small eigenvalues λ ε are We remark that the scaling of these small eigenvalues is for some c 5 < 0. However, one of the eigenvalues of M(p 0 ) is exactly zero, with eigenvector (1, 1, . . . , 1) T . To determine the sign of the real part of this eigenvalue, we have to expand to the next order. By considering contributions for the curvature of the boundary ∂Ω we have computed in Appendix B that for this eigenvalue λ ǫ ∼ c 6 ǫ 3 for some c 6 < 0.
To summarise, there are small eigenvalues of two orders which differ by the logarithmic factor log ξσ ǫD .

Appendix A: Linear Theory
In this section we prove Proposition 4.1. We follow the Liapunov-Schmidt reduction method. Suppose that to the contrary, there exist sequences ε n , σ n , p n and Σ n with ε n , σ n → 0 and Σ n = φ n ψ n ∈ K ⊥ εn,pn such that and φ n H 2 (Ωε) + ψ n H 2 (Ωε) = 1.
We now show that this is impossible. To simplify notation, we omit the index n. In the first step we show that the linearised problem given above tends to a limit problem as max{σ, D} → 0. We define We define ψ ε,i by ∆ψ ε,i − σ 2 ψ ε,i + 2U φ ε,i = 0, in Ω ε , ∂ψ ε,i ∂ν = 0 on ∂Ω ε .
To complete the proof of Proposition 4.1, we just need to show that conjugate operator to L ǫ,p (denoted by L * ǫ,p ) is injective from K ⊥ ǫ, to C ⊥ ǫ,p . The proof for L * ǫ,p follows almost the same process as for L ǫ,p and therefore it is omitted.
The proof is complete.

Appendix B: Computation of the small eigenvalues
In this appendix we will compute the small eigenvalues. First we expand the solution to a higher degree of accuracy than in Section 3. Then we expand the eigenfunctions and finally we calculate the small eigenvalues. 9.1. Further expansion of the solution. In this subsection, we further improve our expansion to the solutions derived in Section 3.
First we define where u ǫ is the exact boundary cluster solution derived in Section 3-5 and χ ε is the cutoff function given in (3.13). It is easy to see that We will derive an approximation to u i which is more accurate than that given in Section 3.
Our main idea is to start with a single boundary spike solution of the Gierer-Meinhardt system in a disk B R of radius R such that the curvature at the centres of the boundary spikes on the disk and the domain ∂Ω agree. Further, the boundary spike solution in the disk is invariant under rotations of the solution which results in a zero eigenvalue. Then the small eigenvalue of the boundary spike solution in Ω can be computed by a perturbation analysis of the disk case.
The single boundary spike solution (u 0 , v 0 ) in the ball B R with max B R u 0 = u 0 (0, R) in polar coordinates solves the following system: It can be constructed from the ground state w following the approach in Section 3 by using a fixed-point argument in the space of even functions around α = 0 and no Liapunov-Schmidt reduction is needed. Note that the single-boundary spike solution u 0 is invariant under rotations of the solution. Therefore we can apply ∂ ∂α in (9.1) and get As a preparation for this perturbation analysis, we represent the boundary ∂Ω in a neighbourhood of the centre of the spike p i in polar coordinates and deform it to a circle with the same curvature. This will imply that the perturbation of the boundary will only be in order ǫ 3 .
Near the point p i we have expanded the boundary ∂Ω in Section 3.1 using Cartesian coordinates. We have derived that ρ(0) = ρ ′ (0) = 0 and ρ ′′ (0) is the curvature of ∂Ω at p i . Recall that ρ was used to flatten the boundary ∂Ω near p i to a line.
Second we have 2β = R(1 − ρ ′′ (0)R) which implies that β = 0 provided ρ ′′ (0) = 1 R and from now on we choose R such that this condition is satisfied. Third we compute Finally, we get To summarise, we have where f (φ) denotes the radius and φ the angle and Further, we have f ′′ (0) = 0 due to the choice of the leading term f (0) = R. Since we have f (φ) = r, we get r = R + O(φ 3 ) which enables us to replace r by R in the derivation of f (φ).
In polar coordinates, for r = R we get a point on ∂B R . We change variables such that for the variables (φ, r ′ ) at r ′ = R we get a point on ∂Ω. Thus r ′ = R has to map into r = f (φ). This is achieved by defining r ′ = r − f (φ) + R for each φ.
Comparing with Cartesian coordinates (z 1 , z 2 ) for the ǫ-scale inside the spike used in Section 3.1, by elementary trigonometry we get Now we approximate the exact solution for the activator u ǫ as follows: Near the centre p i of each spike we have where we have used the notation for i = 1, 2, . . . , k, j = 1, 2, 3 and χ ε is the cutoff function defined in (3.13).
Let us derive the terms in this expansion step by step. We start from the single boundary spike solution u 0 defined in (9.1) in a ball of radius R such that ρ ′′ (0) = 1 R . Using polar coordinates, we can represent u 0 (ǫα, R − r) and this function satisfies the Neumann boundary condition in a ball: ∂u 0 ∂r r=R = 0.
However, the function u 0 does not satisfy the boundary condition at ∂Ω (and for r = R we reach the boundary of the disk (circle) but not the domain boundary ∂Ω). Recall that r−f (φ) = r ′ −R.
Thus, to get a function with a better approximation to the boundary condition at ∂Ω (and such that for r = R we reach ∂Ω), we definē w(ǫα, r ′ ) = u 0 (ǫα, r).
Then for r ′ = R, we havew (ǫα, R) = u 0 (ǫα, f (ǫα)), i.e. for r ′ = R the arguments of the functionw are contained in ∂Ω and the arguments of u 0 are contained in ∂B R . This means that the function f deforms the boundary to a circle (in the same way as in Section 3.1 ρ deforms the boundary to a straight line). Since the circle also takes into account the curvature it gives a better approximation to the boundary than the straight line and the approximate spike solution will give a better approximation to the exact solution than the approximations in Section 3.1. Note that the boundary spike solution in the ball is invariant under rotation (in the same way as the boundary spike is translation invariant in half space). Now we calculate the radial derivative ofw as follows: The outward unit normal vector in polar coordinates is given by Using rescaled variables, this implies Subtracting the radial unit vector e r = (0, 1), this implies .
Using that w i − w H 2 (Ωǫ) = O(ǫα) and thatw satisfies (9.7), the outward normal derivative of w is computed as Now we compute the terms ǫ 3v(1) (even around α = 0) and ǫ 2v(2) (odd around α = 0) in the expansion (9.5) near p i such that the solution satisfies the Neumann boundary condition to higher order. Letv We calculate for x = p i + z i (z). Further, we recall that by (4.4) we have By the reduced problem in Section 5, we have Therefore we can add another contributionv to the solution such that ε 2v (3) and adding this part to the solution will cancel out the odd terms (with respect to α = 0) of order ǫ 2 in the activator equation and we get Taking the derivative ∂ ∂α in this relation near x = p i , we compute ∂ ∂α On the other hand, (9.10) where Finally, by an expansion similar to (9.8) we have 9.2. Expansion of the eigenfunction. We define the approximate kernels to be Then we expand the eigenfunction as follows: (9.11) where φ ⊥ ǫ ∈ K ⊥ ε,p . Suppose that φ ε H 2 (Ωε) = 1. Then |a j,ε | ≤ C. Let us put a ε := (a 1,ε , · · · , a k,ε ) T (9.12) Then for a subsequence and a i,0 = lim ε→0 a i,ε , a 0 := (a 1,0 , · · · , a k,0 ). (9.13) The decomposition of φ ε in (9.11) implies that where ψ i,ε satisfies ∆ψ i,ε − σ 2 (1 + τ λ ε )ψ i,ε + 2uφ i,ε = 0 in Ω ε ∂ψ i,ε ∂ν = 0 on ∂Ω ε (9.15) which we also write as ψ i,ε = T τ λε [φ i,ε ], and ψ ⊥ ε is given by Let us first consider the leading term of φ ε . For ∂w i ∂α we get, using (9.8), Therefore, expanding the boundary condition as we did above forw i , we definev (1) eig,i as the unique solution of ∂α g (4) (0)ǫ 3 α 3 on ∂B R . Similarly, letv (2) eig,i be the unique solution of Let us compare this with the derivative ∂v (1) ∂α which satisfies R,i is given by the difference of the previous two contributions as follows: v (1) R,i =v (1) eig,i − ∂v (1) ∂α which satisfies We note thatv R,i is an even function around α = 0.
Substituting the decompositions of φ ε and ψ ε into (6.1), we have R,i , and We first compute we can estimate I 4 as follows: where we use the equation satisfied by ψ i,ε . 9.3. Expansion of the small eigenvalues. Multiplying both sides of (9.18) by ∂u l ∂τ (p l ) and integrating over Ω ε , we have (1)).
Using the estimate I 4 , we have where J i,l are defined as the integrals in the last equality. We divide our proof into several steps. The following lemma contains the key estimates: Lemma 9.1. We have where a i,ε has been defined in (9.11).

Acknowledgements
The research of J. Wei is partially supported by NSERC of Canada. MW thanks the Department of Mathematics at UBC for their kind hospitality.