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

TL;DR

This paper constructs and analyzes stable boundary spike clusters in a two-dimensional Gierer-Meinhardt system with small diffusivities, demonstrating their stability and the underlying mechanisms of spike interactions.

Contribution

It introduces a method to construct stable boundary spike clusters in the Gierer-Meinhardt system with small diffusivities, highlighting the role of boundary curvature and spike interactions.

Findings

Existence of boundary spike clusters approaching boundary curvature maxima

Proof of linear stability of these spike clusters

Identification of repulsive spike interactions due to small inhibitor diffusivity

Abstract

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.

Equations

$$\left\{\begin{array}[]{l}u_{t}=\varepsilon^{2}\Delta u-u+\frac{u^{2}}{v},u>0\mbox{ in }\Omega,\\[5.69054pt] \tau v_{t}=D\Delta v-v+u^{2},v>0\mbox{ in }\Omega,\\[5.69054pt] \frac{\partial u}{\partial\nu}=\frac{\partial v}{\partial\nu}=0\mbox{ on }\partial\Omega.\end{array}\right.$$

$$0<\epsilon \ll 1,0<D\ll 1,$$

$$e^{-\frac{1}{\sqrt{D}}}\ll \epsilon \ll \sqrt{D}.$$

$$D_1(\epsilon )>D_2(\epsilon )>\cdots >D_k(\epsilon )>\cdots$$

$$\left\{\begin{array}[]{l}\Delta w-w+w^{2}=0,\\[5.69054pt] w(0)=\max_{y\in{\mathbb{R}}^{2}}w(y),\ w(y)\to 0\mbox{ as }|y|\to\infty.\end{array}\right.$$

$$w(y)\sim |y{|}^{-\frac{1}{2}}{e}^{-|y|}\text{mbox}as|y|\to \infty$$

$$w^{\prime }(y)=-(1+o(1))w(y)\text{mbox}as|y|\to \infty ,$$

$$u_{ϵ}\sim \frac{D{\xi }_{\sigma }}{{ϵ}^2}\sum _{i=1}^{k}w(\frac{x-P_{i,ϵ}}{ϵ}),$$

$${\xi }_{\sigma }\sim \frac{1}{\log \frac{\sqrt{D}}{ϵ}}$$

$$|P_i-P_{i-1}|\sim \sqrt{D}\log \frac{{\xi }_{\sigma }}{ϵD},i=2,\dots ,k.$$

$$\left\{\begin{array}[]{l}\varepsilon^{2}\Delta u-u+\frac{u^{2}}{v}=0\mbox{ in }\Omega,\\[5.69054pt] D\Delta v-v+u^{2}=0\mbox{ in }\Omega,\\[5.69054pt] \frac{\partial u}{\partial\nu}=\frac{\partial v}{\partial\nu}=0\mbox{ on }\partial\Omega,\end{array}\right.$$

$${\nabla }_{\tau (P)}:=\frac{\partial }{\partial \tau (P)},$$

$${\xi }_{\sigma ,i}\sim \frac{1}{\frac{1}{\pi }\log \frac{1}{\sigma }{\int }_{{\mathbb{R}}_{+}^2}{w}^{2}dx}$$

$$\left|\frac{P_i-P_{i-1}}{\sqrt{D}}-\log \frac{{\xi }_{\sigma }}{\epsilon D}+\frac{3}{2}\log \log \frac{{\xi }_{\sigma }}{\epsilon D}+\log (-\frac{\frac{{\partial }^2}{\partial {\tau }^2}h(P^0){\nu }_1}{2{\nu }_2})+\log [(i-1)(k+1-i)]\right|\le \eta$$

$$\left|\frac{1}{k}\sum _{i=1}^k{P}_i-P^0\right|\le \eta \sqrt{D}\log \frac{{\xi }_{\sigma }}{\epsilon D},$$

$$|P_i-P_{i-1}|\sim C\sqrt{D}\log \frac{{\xi }_{\sigma }}{\epsilon D}.$$

$$\hat{u}(z)={\sigma }^{2}u(\epsilon z),\hat{v}(z)={\sigma }^{2}v(\epsilon z),$$

$$\left\{\begin{array}[]{l}\Delta u-u+\frac{u^{2}}{v}=0\mbox{ in }\Omega_{\varepsilon},\\[5.69054pt] \Delta v-\sigma^{2}v+u^{2}=0\mbox{ in }\Omega_{\varepsilon},\\[5.69054pt] \frac{\partial u}{\partial\nu}=\frac{\partial v}{\partial\nu}=0\mbox{ on }\partial\Omega_{\varepsilon},\end{array}\right.$$

$${\Lambda }_k=\{{\mathbf{p}}^{\epsilon }={\epsilon }^{-1}{\mathbf{P}}^{\epsilon }:P_1,\cdots ,P_k\text{mbox}suchthat(\text{ref}e:2)\text{mbox}and(\text{ref}e:3)\text{mbox}hold\}.$$

$$\left\{\begin{array}[]{l}u(z)\sim\sum_{i=1}^{k}\xi_{\sigma,i}Pw_{p_{i}}(z-p_{i}),\\[5.69054pt] v(p_{i})\sim\xi_{\sigma,i},\end{array}\right.$$

$$\Delta u-u+w(\cdot -p_i{)}^2=0\text{mbox}in{\Omega }_{\epsilon ,p_i},\frac{\partial u}{\partial \nu }=0\text{mbox}on\partial {\Omega }_{\epsilon ,p_i}.$$

$${\epsilon }^{-1}\rho (\epsilon (z_1-{\mathbf{p}}_1))=\frac{1}{2}{\rho }^{\prime \prime }(0)\epsilon (z_1-{\mathbf{p}}_1{)}^2+\frac{1}{6}{\rho }^{(3)}(0){\epsilon }^2(z_1-{\mathbf{p}}_1{)}^3+O({\epsilon }^3(z_1-{\mathbf{p}}_1{)}^4).$$

$$\left\{\begin{array}[]{l}\Delta h_{p}(z)-h_{p}(z)=0,\mbox{ in }\Omega_{\varepsilon},\\[5.69054pt] \frac{\partial h_{p}}{\partial\nu}=\frac{\partial}{\partial\nu}w(z-p)\mbox{ on }\partial\Omega_{\varepsilon}.\end{array}\right.$$

$$\left\{\begin{array}[]{l}y_{1}=z_{1}-{\bf p}_{1},\\[5.69054pt] y_{2}=z_{2}-{\bf p}_{2}-\varepsilon^{-1}\rho(\varepsilon(z_{1}-{\bf p}_{1})).\end{array}\right.$$

$${\Delta }_z={\Delta }_y+({\rho }^{\prime }(\epsilon y_1){)}^2{\partial }_{y_2{y}_2}-2{\rho }^{\prime }(\epsilon y_1){\partial }_{y_1{y}_2}-\epsilon {\rho }^{\prime \prime }(\epsilon y_1){\partial }_{y_2},$$

$$(1+{\rho }^{\prime }(\epsilon y_1{)}^2{)}^{\frac{1}{2}}\frac{\partial }{\partial \nu }={\rho }^{\prime }(\epsilon y_1){\partial }_{y_1}-(1+{\rho }^{\prime }(\epsilon y_1{)}^2){\partial }_{y_2}.$$

$$\left\{\begin{array}[]{l}\Delta v-v=0\mbox{ in }{\mathbb{R}}^{2}_{+},\\[5.69054pt] \frac{\partial v}{\partial y_{2}}=\frac{w^{\prime}}{|y|}\frac{\rho^{\prime\prime}(0)}{2}y_{1}^{2}\mbox{ on }\partial{\mathbb{R}}^{2}_{+},\end{array}\right.$$

$$\left\{\begin{array}[]{l}\Delta v-v-2\rho^{\prime\prime}(0)y_{1}\frac{\partial^{2}v^{(1)}}{\partial_{y_{1}}\partial_{y_{2}}}=0\mbox{ in }{\mathbb{R}}^{2}_{+},\\[5.69054pt] \frac{\partial v}{\partial y_{2}}=-\rho^{\prime\prime}(0)y_{1}\frac{\partial v^{(1)}}{\partial y_{1}}\mbox{ on }\partial{\mathbb{R}}^{2}_{+}.\end{array}\right.$$

$$\left\{\begin{array}[]{l}\Delta v-v=0\mbox{ in }{\mathbb{R}}^{2}_{+},\\[5.69054pt] \frac{\partial v}{\partial y_{2}}=\frac{w^{\prime}}{|y|}\frac{1}{3}\rho^{(3)}(0)y_{1}^{3}\mbox{ on }\partial{\mathbb{R}}^{2}_{+}.\end{array}\right.$$

$${\chi }_{\epsilon }(z-p)=\chi (\epsilon (z-p))\text{mbox}forz\in {\Omega }_{\epsilon }.$$

$$h_p(z)=-(\epsilon v^{(1)}(y)+{\epsilon }^2(v^{(2)}(y)+v^{(3)}(y))){\chi }_{\epsilon }(z-p)+{\epsilon }^3{\xi }_p(z),z\in {\Omega }_{\epsilon }.$$

Full text

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.

Keywords: Pattern formation, reaction-diffusion system, spike, cluster, stability.

1. 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:

[TABLE]

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\in\Omega\subset{\mathbb{R}}^{2}$ and at a time $t>0$. Here $\Delta$ is the Laplace operator in ${\mathbb{R}}^{2}$, $\Omega$ is a smooth bounded domain in ${\mathbb{R}}^{2}$, $\nu=\nu(x)$ is the outer unit normal at $x\in\partial\Omega$.

Throughout this paper, we assume that

[TABLE]

$\tau\geq 0$ is a fixed constant independent of $\varepsilon$, $D$ and $x$. Further, the diffusivities $\varepsilon$ and $D$ do not depend on $x$ but they are both small constants. In this paper, we further assume that

[TABLE]

This means that $\epsilon$ is much smaller than ${D}$. On the other hand, $\epsilon$ cannot be exponentially small compared to $\sqrt{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\sim 1$, the second and third authors constructed single-interior 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 $\tau=0$.

For the weak coupling case $D\to\infty$, in [21] the second and third authors proved the existence of multiple interior spike solutions. Further, they showed that there are stability thresholds

[TABLE]

such that if $\lim_{\varepsilon\to 0}\frac{D_{k}(\varepsilon)}{D}>1$, the $k$-peak solution is stable and if $\lim_{\varepsilon\to 0}\frac{D_{k}(\varepsilon)}{D}<1$, the $k$-peak solution is unstable. Multiple spikes for the Gierer-Meinhardt system in a one-dimensional interval have been studied in [13, 7, 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\to\infty$.)

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.

2. Main Results: Existence and Stability

Let $\Omega\subset{\mathbb{R}}^{2}$ be bounded and smooth two-dimensional domain. Let $w$ be the unique solution in $H^{1}({\mathbb{R}}^{2})$ of the problem

[TABLE]

For the existence and uniqueness of (2.1), we refer to [8] and [11]. We also recall that $w$ is radially symmetric and

[TABLE]

and

[TABLE]

where $w^{\prime}$ is the radial derivative of $w$, i.e. $w^{\prime}=w_{r}(r)$.

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 $\partial\Omega$.

Then for $0<e^{-\frac{1}{\sqrt{D}}}\ll\varepsilon\ll\sqrt{D}\ll 1$, the Gierer-Meinhardt system (1.1) has a $k$-boundary spike cluster steady-state solution $(u_{\epsilon},v_{\epsilon})$ which concentrates near $P^{0}$. In particular, it satisfies

[TABLE]

where $P_{i,\epsilon}\to P^{0}$ as $\epsilon\to 0$ for $i=1,2,\ldots,k$.

Further, we have

[TABLE]

and

[TABLE]

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 $\tau$ 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 $\epsilon^{3}\log\frac{\xi_{\sigma}}{\epsilon D}$, and one eigenvalue stemming from the curvature of the boundary (corresponding to synchronous motion of all spikes) is of order $\epsilon^{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.

3. Introduction of the approximate solutions

In this paper, we consider the following problem:

[TABLE]

where $\Omega\subset{\mathbb{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 $\Omega$. For $P\in\partial\Omega$,

[TABLE]

where $\frac{\partial}{\partial\tau(P)}$ denotes the tangential derivative with respect to $P$ at $P\in\partial\Omega$. 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 $\sigma^{-2}\xi_{\sigma,i}w(\frac{x-P_{i}}{\varepsilon})$ which are centred at the points $P_{i}$ for $i=1,\cdots,k$, where $\sigma=\frac{\varepsilon}{\sqrt{D}}$ and the amplitude $\xi_{\sigma,i}$ satisfies

[TABLE]

(see (3.38)).

Let $P_{1},\cdots,P_{k}$ be $k$ points distributed along the boundary $\partial\Omega$ such that we have for $i=2,\cdots,k$

[TABLE]

and

[TABLE]

where $\nu_{2}$ is given in (5.10) below. Further, $\eta>0$ is a small constant independent of $\varepsilon$ and $D$. The reason for assuming (3) and (3.3) will become clear in Section 5 when we solve the reduced problem.

Remark 3.1**.**

By (3), the distance of neighbouring spikes satisfies

[TABLE]

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,

[TABLE]

if we drop the hat and still denote solutions by $(u,v)$, equation (3.1) is equivalent to

[TABLE]

where $\Omega_{\varepsilon}=\varepsilon^{-1}\Omega$.

From now on we will deal with (3.4). Before introducing an approximation to the spike solutions, we first define some notation.

Fixing ${\bf P}^{\varepsilon}=(P_{1},\cdots,P_{k})$ such that (3) and (3.3) hold, we set

[TABLE]

We are looking for multiple spike solutions to (3.4) of the form

[TABLE]

where $Pw_{p_{i}}(z-p_{i})$ is defined to be the unique solution of

[TABLE]

Here $\Omega_{\varepsilon,p_{i}}=\{z\,:\,z+p_{i}\in\Omega_{\varepsilon}\}$, the function $w$ has been defined in (2.1) and $\xi_{\sigma,i},\ i=1,\cdots,k$ are the heights of the spikes, which will be determined in (3.38).

3.1. The analysis of the projection $Pw_{q}(z-q)$

Before calculating the heights of the spikes, we need some preliminaries of the projection $Pw_{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\in\partial\Omega$. We define a diffeomorphism straightening the boundary. We may assume that the inward normal to $\partial\Omega$ at $P$ is pointing in the direction of the positive $x_{2}$ axis. Denote $B^{\prime}({R})=\{x\in{\mathbb{R}}^{2}||x_{1}|\leq{R}\}$. Then since $\partial\Omega$ is smooth, we can find a constant $R$ such that $\partial\Omega$ can be represented by the graph of a smooth function $\rho_{P}:B^{\prime}(R)\to{\mathbb{R}}$, where $\rho_{P}(0)=0$, and $\rho_{P}^{\prime}(0)=0$. From now on, we omit the use of $P$ in $\rho_{P}$ and write $\rho$ if this can be done without causing confusion. So near $P$, $\partial\Omega$ can be represented by $(x_{1},\rho(x_{1}))$. The curvature of $\partial\Omega$ at $p$ is $h(P)=\rho^{\prime\prime}(0)$. Let $\Omega_{1}=\Omega\cap B(P,{R})=\{(x_{1},x_{2})\in B(P,{R})|x_{2}-P_{2}>\rho(x_{1}-P_{1})\}$, where $B(P,R)=\{x\in{\mathbb{R}}^{2}||x-P|<R\}$.

After rescaling, it follows that near $p=\frac{P}{\varepsilon}$, the boundary $\partial\Omega_{\varepsilon}$ can be represented by $(z_{1}-{\bf p}_{1},\varepsilon^{-1}\rho(\varepsilon(z_{1}-{\bf p}_{1})))$, where $(z_{1},z_{2})=\varepsilon^{-1}(x_{1},x_{2})$ and $p=({\bf p}_{1},{\bf p}_{2})$. By Taylor expansion, we have

[TABLE]

Let $h_{p}(z)=w(z-p)-Pw_{p}(z-p)$. Then $h_{p}$ satisfies

[TABLE]

For $z\in\Omega_{1\varepsilon}=\frac{1}{\varepsilon}\Omega_{1}$ we set

[TABLE]

Under this transformation, the Laplace operator and the boundary derivative operator become

[TABLE]

Let $v^{(1)}$ be the unique solution of

[TABLE]

where ${\mathbb{R}}^{2}_{+}$ is the upper half plane, namely ${\mathbb{R}}^{2}_{+}=\{y=(y_{1},y_{2})\in{\mathbb{R}}^{2}|y_{2}>0\}.$

Let $v^{(2)}$ be the unique solution of

[TABLE]

Let $v^{(3)}$ be the unique solution of

[TABLE]

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)|\leq Ce^{-\mu|y|}$ for any $0<\mu<1$. Let $\chi$ be a smooth cut-off function, such that $\chi(a)=1$ for $a\in B(0,R_{0}\sqrt{D}\log\frac{\xi_{\sigma}}{\varepsilon D})$, and $\chi(a)=0$ for $x\in B(0,2R_{0}\sqrt{D}\log\frac{\xi_{\sigma}}{\varepsilon D})^{c}$ for some suitable $R_{0}$ such that $|p_{i}-p^{0}|<\frac{R_{0}}{\sigma}\log\frac{\xi_{\sigma}}{\epsilon D}$, and

[TABLE]

Set

[TABLE]

Then we have the following estimate:

Proposition 3.1**.**

[TABLE]

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 $|\xi_{p}(z)|\leq Ce^{-\mu|z-p|}$ for any $0<\mu<1$.

Similarly we know from [19] that

Proposition 3.2**.**

[TABLE]

where $\eta$ is the unique solution of the following equation:

[TABLE]

Moreover,

[TABLE]

It follows that $\eta(y)$ is an odd function in $y_{1}$. It can be seen that $|\eta_{1}(y)|\leq Ce^{-\mu|y|}$ for some $0<\mu<1$.

Finally, let

[TABLE]

We have

Lemma 3.1**.**

[TABLE]

where $H^{2}_{N}({\mathbb{R}}^{2}_{+})=\{u\in H^{2}({\mathbb{R}}^{2}_{+})\,:\,\frac{\partial u}{\partial y_{2}}=0\mbox{ on }\partial{\mathbb{R}}^{2}_{+}\}$.

Proof.

See Lemma 4.2 in [12]. ∎

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_{i}^{(j)}$ be the corresponding solutions in the expansion of $h_{p_{i}}$.

3.2. The analysis of the Green’s function

Next we introduce a Green’s function $G_{\sqrt{D}}$ which is needed to derive our main results.

For $D>0$, let $G_{D}(x,y)$ be the Green’s function given by

[TABLE]

where $y\in\partial\Omega_{\sqrt{D}}$ and let $G_{0}$ be the Green’s function of the upper half plane:

[TABLE]

Then $H(x)=G_{\sqrt{D}}(x)-G_{0}(x)$ will satisfy

[TABLE]

Let $\eta_{1}$ be the solution of

[TABLE]

and let $\eta_{2}$ be the solution of

[TABLE]

It can be seen easily that $\eta_{1}$ is even in $y_{1}$ and $\eta_{2}$ is odd in $y_{1}$. Then one can get the following result.

Lemma 3.2**.**

[TABLE]

Proof.

First we compute on $\partial\Omega_{\sqrt{D}}$,

[TABLE]

for any $0<l<1$.

Since we have

[TABLE]

we get

[TABLE]

From the expansion above, we can get the asymptotic behaviour of $G_{\sqrt{D}}$.

∎

Next we have the following expansion of $G_{0}$:

Lemma 3.3**.**

The following expansion of $G_{0}$ holds:

[TABLE]

for $0<r<1$, where $\psi$ is a smooth function with $\psi(0)=\psi^{\prime}(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 ${\mathbb{R}}^{2}$. For the expansion of fundamental solution, see Lemma 4.1 in [1]. Then the above expansion follows.

∎

We set

[TABLE]

From the estimates above, and for points ${\bf p}^{\varepsilon}\in\Lambda_{k}$, we have

[TABLE]

[TABLE]

Generally, we have

[TABLE]

For the derivatives, we estimate

[TABLE]

3.3. Calculating the heights of the spikes

In this section, we are going to determine the heights of spikes $\xi_{\sigma,i}$ to leading order. In the sequel, by $T[h]$ we denote the unique solution of the equation

[TABLE]

Then we know that

[TABLE]

As mentioned before, we will choose the approximate solution to be

[TABLE]

and

[TABLE]

First we calculate the heights of the peaks:

[TABLE]

Thus

[TABLE]

We assume that the heights of the spikes are asymptotically equal as $\varepsilon,D\to 0$, i.e.

[TABLE]

Then we get that

[TABLE]

where

[TABLE]

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.

4. 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:

[TABLE]

where we have used the notation

[TABLE]

On the other hand, we calculate for $x=p_{i}+z$

[TABLE]

where $R_{1,i}(z)=2w(z)[v_{i}^{(1)}(z)+\varepsilon v_{i}^{(2)}(z)]+\varepsilon(v_{1}^{(i)}(z))^{2}$, $R_{2,i}=2w(z)v_{i}^{(3)}(z)$. This implies

[TABLE]

where $R_{1}(z)$ is even in $z_{1}$ and $R_{2}(z)$ is odd in $z_{1}$ and

[TABLE]

and

[TABLE]

Thus we can get that

[TABLE]

By the above estimates, we have the following key estimate:

Lemma 4.1**.**

For $x=p_{i}+z$, $|z|<\frac{R_{0}}{\sigma}\log\frac{\varepsilon D}{\xi_{\sigma}}$, we have

[TABLE]

where

[TABLE]

[TABLE]

where $\tilde{R}_{1}(z)$ is even in $z_{1}$ with the property that $\tilde{R}_{1}(z)=O(\log(1+|z|))$, and $R_{2}(z),R_{2,i}(z)$ are defined above.

Further, $S_{1}(U,V)=\varepsilon^{\frac{R_{0}}{\sigma}}$ for $|x-p_{i}|\geq\frac{R_{0}}{\sigma}\log\frac{\varepsilon D}{\xi_{\sigma}}$ for all $i$.

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

[TABLE]

[TABLE]

We first define

[TABLE]

and define the approximate kernels by

[TABLE]

and choose the approximate cokernels as follows:

[TABLE]

We then define

[TABLE]

where $C_{\varepsilon,{\bf p}}^{\perp},K_{\varepsilon,{\bf p}}^{\perp}$ denote the orthogonal complements with the scalar product of $L^{2}(\Omega_{\varepsilon})$.

Let $\pi_{\varepsilon,{\bf p}}$ denote the projection in $L^{2}(\Omega_{\varepsilon})\times L^{2}(\Omega_{\varepsilon})$ onto $\mathcal{C}_{\varepsilon,{\bf p}}^{\perp}$. We are going to show that the equation

[TABLE]

has a unique solution \Sigma_{\varepsilon,{\bf p}}=\left(\begin{array}[]{c}\phi_{\varepsilon,{\bf p}}\\ \psi_{\varepsilon,{\bf p}}\end{array}\right)\in\mathcal{K}_{\varepsilon,{\bf p}}^{\perp} if $\max\{\sigma,D\}$ is small enough.

Set

[TABLE]

Written in components, we have

[TABLE]

and

[TABLE]

As a preparation we state a result on the invertibility of the corresponding linearised operator $\mathcal{L}_{\varepsilon,{\bf p}}$ whose proof is postponed to Appendix A.

Proposition 4.1**.**

There exist positive constants $\bar{\delta},C$ such that for $\max\{\sigma,D\}<\bar{\delta}$, the map $\mathcal{L}_{\varepsilon,{\bf p}}$ is surjective for arbitrary ${\bf p}\in\Lambda_{k}$. Moreover the following estimate holds:

[TABLE]

Now we are in the position to solve the equation

[TABLE]

Since $\mathcal{L}_{\varepsilon,{\bf p}}|_{\mathcal{K}_{\varepsilon,{\bf p}}^{\perp}}$ is invertible, we can write the above equation as

[TABLE]

where

[TABLE]

and

[TABLE]

and the operator $M_{\varepsilon,{\bf p}}$ is defined above for $\Sigma\in H_{N}^{2}(\Omega_{\varepsilon})\times H^{2}_{N}(\Omega_{\varepsilon})$. We are going to show that the operator $M_{\varepsilon,{\bf p}}$ is a contraction mapping on

[TABLE]

if $C_{0}$ is large enough. We have that

[TABLE]

where $C>0$ is independent of $\varepsilon>0$ and $c(\xi_{\sigma})\to 0$ as $\xi_{\sigma}\to 0$. Similarly we can show that

[TABLE]

where $c(\xi_{\sigma})\to 0$ as $\xi_{\sigma}\to 0$. If we choose $C_{0}$ large enough, then $M_{\varepsilon,{\bf p}}$ is a contraction mapping on $B_{\varepsilon}$. The existence of a fixed point $\Sigma_{\varepsilon,{\bf p}}$ together with an error estimate now follows from the contraction mapping principle. Moreover, $\Sigma_{\varepsilon,{\bf p}}$ is a solution. Thus we have proved

Lemma 4.2**.**

There exists $\bar{\delta}>0$ such that for every triple $(\varepsilon,D,{\bf p})$ with $\max\{\sigma,D\}<\bar{\delta}$, and ${\bf p}\in\Lambda_{k}$, there exists a unique $(\phi_{\varepsilon,{\bf p}},\psi_{\varepsilon,{\bf p}})\in\mathcal{K}_{\varepsilon,{\bf p}}^{\perp}$ satisfying

[TABLE]

and

[TABLE]

More refined estimates for $\phi_{\varepsilon,{\bf 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 $\phi_{\varepsilon,{\bf p}}$.

Lemma 4.3**.**

Let $\phi_{\varepsilon,{\bf p}}$ be defined by (4.17). Then for $x=p_{i}+z$, we have

[TABLE]

where $\phi_{\varepsilon,{\bf p},1}$ is even in $z_{1}$ which can be estimated by

[TABLE]

and $\phi_{\varepsilon,{\bf p},2}$ can be estimated by

[TABLE]

Proof.

Let $S(u)=S_{1}(u,T(u))$. We first solve

[TABLE]

for $\phi_{\varepsilon,{\bf p},1}\in K_{\varepsilon,{\bf p}}^{\perp}$. Then we solve

[TABLE]

for $\phi_{\varepsilon,{\bf p},2}\in K_{\varepsilon,{\bf p}}^{\perp}$.

Using the same proof as in Lemma 4.2, the above two equations are uniquely solvable for $\max\{\sigma,D\}\ll 1$. By uniqueness, $\phi_{\varepsilon,{\bf p}}=\phi_{\varepsilon,{\bf p},1}+\phi_{\varepsilon,{\bf p},2}$. Since $S_{1,2}=S_{1,2}^{0}+S_{1,2}^{\perp}$, where $S_{1,2}^{0}=O(\xi_{\sigma}\varepsilon^{2}\sqrt{D}\log\frac{\xi_{\sigma}}{\varepsilon D})$ and $S_{1,2}^{\perp}\in C_{\varepsilon,{\bf p}}^{\perp}$, it is easy to see that $\phi_{\varepsilon,{\bf p},1}$ and $\phi_{\varepsilon,{\bf p},2}$ have the required properties.

∎

5. 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 ${\bf p}\in\Lambda_{k}$, there exists a unique solution $(\phi,\psi)\in\mathcal{K}_{\varepsilon,{\bf p}}^{\perp}$ such that

[TABLE]

We need to determine ${\bf p}=(p_{1},\cdots,p_{k})\in\Lambda_{k}$ such that

[TABLE]

and therefore S_{\varepsilon}\left(\begin{array}[]{c}U+\phi\\ V+\psi\end{array}\right)=0.

To this end, we calculate the projection:

[TABLE]

where $I_{1},I_{2}$ are defined by the last equality.

For $I_{1}$, we have

[TABLE]

Note that by the estimates satisfied by $\phi$ in Lemma 4.3, we have

[TABLE]

where we have used the estimates (3.28)-(3.31). Further, we have

[TABLE]

We compute

[TABLE]

By (3.12), we have

[TABLE]

where the constant $\nu_{1}>0$ is defined by

[TABLE]

Now by (5.2)-(5.6),

[TABLE]

Next we estimate $I_{12}$:

[TABLE]

where we have used $\sqrt{D}\log\frac{\xi_{\sigma}}{\epsilon D}\ll 1$ and

[TABLE]

Thus by (5.8) and (5.9), we have

[TABLE]

Next we estimate $I_{2}$:

[TABLE]

By the equation for $\psi$, we have

[TABLE]

and therefore

[TABLE]

where $R(z)$ is an even function in $z_{1}$.

Thus we have

[TABLE]

Combining the estimates for $I_{1}$ and $I_{2}$, (5.11) and (5.12), we have

[TABLE]

Thus $W_{\varepsilon,i}=0$ is reduced to the following system:

[TABLE]

We first solve the limiting case:

[TABLE]

This system is uniquely solvable with

[TABLE]

Moreover, we have

[TABLE]

where we have used the notation $h^{\prime\prime}(p^{0})=\frac{\partial^{2}}{\partial\tau^{2}}h(P^{0})<0$.

From (5), we know

[TABLE]

To find $p_{i}$ such that $W_{\varepsilon,i}=0$, we expand $p_{i}=p_{i}^{0}+\tilde{p_{i}}$. Then adding the first $i$ equations, we have

[TABLE]

One can get that

[TABLE]

By the last equation

[TABLE]

Thus

[TABLE]

From the equation above, we can estimate $\tilde{p}_{1}$ by

[TABLE]

In conclusion, we solve $W_{\varepsilon,i}=0$ with

[TABLE]

Thus we have proved the following proposition:

Proposition 5.1**.**

For $\max\{\sigma,D\}$ small enough, there exists ${\bf p}^{\varepsilon}\in\Lambda_{k}$ with $P_{i}\to P^{0}$ such that $W_{\varepsilon,i}=0$.

Finally, we complete the proof of Theorem 2.1.

Proof. By Proposition 5.1, there exists ${\bf P}^{\varepsilon}\to{\bf P}^{0}$, such that $W_{\varepsilon}({\bf p}^{\varepsilon})=0$. In other words, we have

[TABLE]

Moreover, by the maximum principle, $(U,V)>0$ and the solution satisfies all the properties of Theorem 2.1.

∎

6. Study of the large eigenvalues

We consider the stability of the steady-state $(u,v)$ constructed in Theorem 2.1.

In this section, we first study the large eigenvalues which satisfy $\lambda_{\varepsilon}\to\lambda_{0}\neq 0$ as $\max\{\sigma,D\}\to 0$.

Linearizing the system around the equilibrium states $(u,v)$ obtained in Theorem 2.1, we obtain the following eigenvalue problem:

[TABLE]

for $(\phi,\psi)\in H^{2}_{N}(\Omega_{\varepsilon})\times H^{2}_{N}(\Omega_{\varepsilon})$.

In this section, since we study the large eigenvalues, we may assume that $|\lambda_{\varepsilon}|\geq c>0$ for $\max\{\sigma,D\}$ small enough. If $\mbox{Re}(\lambda_{\varepsilon})\leq-c<0$, then $\lambda_{\varepsilon}$ is a stable large eigenvalue, we are done. Therefore, we may assume that $\mbox{Re}(\lambda_{\varepsilon})\geq-c$ and for a subsequence $\max\{\sigma,D\}\to 0$, $\lambda_{\varepsilon}\to\lambda_{0}\neq 0$. We shall derive the limiting eigenvalue problem which is given by a coupled system of NLEPs.

The second equation of (6.1) is equivalent to

[TABLE]

We introduce the following notation:

[TABLE]

where in $\sqrt{1+\tau\lambda_{\varepsilon}}$, we take the principal part of the square root.

Let us assume that

[TABLE]

We cut off $\phi$ as follows:

[TABLE]

where the cutoff function $\chi_{\epsilon}$ has been defined in (3.13).

From (6.1) and the exponential decay of $w$, it follows that

[TABLE]

Then by a standard procedure (see [6], Section 7.12), we extend $\phi_{\varepsilon,j}$ to a function defined on ${\mathbb{R}}^{2}$ such that

[TABLE]

Since $\|\phi_{\varepsilon}\|_{H^{2}(\Omega_{\varepsilon})}=1$, $\|\phi_{\varepsilon,j}\|_{H^{2}({\mathbb{R}}^{2})}\leq C$ . By taking a subsequence, we may assume that $\phi_{\varepsilon,j}\to\phi_{j}$ as $\max\{\sigma,D\}\to 0$ in $H^{1}({\mathbb{R}}^{2})$ for some $\phi_{j}\in H^{1}({\mathbb{R}}^{2})$ for $j=1,\cdots,k$.

By (6.1), we have

[TABLE]

Substituting the above equation into the first equation of (6.1), letting $\max\{\sigma,D\}\to 0$, and using the expansion of $\xi_{\sigma,j}$, we arrive at the following nonlocal eigenvalue problem (NLEP):

[TABLE]

By Theorem 3.5 in [26], (6.5) has only stable eigenvalues if $\tau$ 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 $\tau$ is small enough.

7. Study of the small eigenvalues

Now we study the eigenvalue problem (6.1) with respect to small eigenvalues. Namely, we assume that $\lambda_{\varepsilon}\to 0$ as $\max\{\sigma,D\}\to 0$. We will show that the small eigenvalues in leading order are related to the matrix $\mathcal{M}({\bf p}^{0})$ given in (7.3) which is computed from the Green’s function. Our main result in this section says that if $\lambda_{\varepsilon}\to 0$, then in leading order

[TABLE]

where $\sigma_{0}$ is an eigenvalue of $\mathcal{M}({\bf p}^{0})$. We will show that all the eigenvalues of $\mathcal{M}({\bf 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 $\mathcal{M}({\bf 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 $\epsilon^{3}\log\frac{\xi_{\sigma}}{\epsilon D}$, the part from the boundary curvature is of order $\epsilon^{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

[TABLE]

where $\nu_{2}$ has been defined in (5.10). Further, we have

[TABLE]

where

[TABLE]

Using the estimates for $p^{0}_{i}$ in (5.15) and (5), we have

[TABLE]

This shows that if all the eigenvalues of $\mathcal{M}({\bf p}^{0})$ have positive real part, then the small eigenvalues are stable. On the other hand, if $\mathcal{M}({\bf 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\times k$ matrix $A$ defined by

[TABLE]

We have the following result from Lemma 16 in [27]:

Lemma 7.1**.**

The eigenvalues of the matrix $A$ are given by

[TABLE]

By Lemma 7.1, the eigenvalues of $\mathcal{M}({\bf p}^{0})$ are all positive except for a single eigenvalue zero with eigenvector $(1,1,\ldots,1)^{T}$.

Equation (7.2) shows that the small eigenvalues $\lambda_{\varepsilon}$ are

[TABLE]

We remark that the scaling of these small eigenvalues is

[TABLE]

for some $c_{5}<0$.

However, one of the eigenvalues of $\mathcal{M}({\bf p}^{0})$ is exactly zero, with eigenvector $(1,1,\ldots,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 $\partial\Omega$ we have computed in Appendix B that for this eigenvalue

[TABLE]

for some $c_{6}<0$.

To summarise, there are small eigenvalues of two orders which differ by the logarithmic factor $\log\frac{\xi_{\sigma}}{\epsilon D}$.

8. 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 $\varepsilon_{n},\sigma_{n},{\bf p}_{n}$ and $\Sigma_{n}$ with $\varepsilon_{n},\sigma_{n}\to 0$ and \Sigma_{n}=\left(\begin{array}[]{c}\phi_{n}\\ \psi_{n}\end{array}\right)\in\mathcal{K}_{\varepsilon_{n},{\bf p}_{n}}^{\perp} such that

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

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\{\sigma,D\}\to 0$.

We define

[TABLE]

and

[TABLE]

It is easy to see that $\phi_{\varepsilon,k+1}=o(1)$ in $H^{2}(\Omega_{\varepsilon})$, since it satisfies

[TABLE]

We define $\psi_{\varepsilon,i}$ by

[TABLE]

Note that since $\xi_{\sigma}^{-1}\|g_{n}\|_{L^{2}(\Omega_{\varepsilon})}\to 0$, we also have $\|\psi_{\varepsilon,k+1}\|_{H^{2}(\Omega_{\varepsilon})}=o(1)$.

Next by the equation satisfied by $\psi_{\varepsilon}$, we have

[TABLE]

So at $x=p_{i}$, we calculate

[TABLE]

Substituting the above equation into the first equation of (6.1), letting $\max\{\sigma,D\}\to 0$, we can show that

[TABLE]

and

[TABLE]

where $\phi_{i}$ is solution of the following nonlocal problem:

[TABLE]

where $C_{0}^{\perp},K_{0}^{\perp}$ denote the orthogonal complements with respect to the scalar product of $L^{2}({\mathbb{R}}^{2}_{+})$ in the space of $H^{2}({\mathbb{R}}^{2}_{+})$ and $L^{2}({\mathbb{R}}^{2}_{+})$ respectively.

By Theorem 1.4 in [16], we know that $\phi_{i}=0,i=1,\cdots,k$.

By taking the limit in the equation satisfied by $\psi_{\varepsilon}$, we see that this implies that $\psi_{\varepsilon}\to 0$ in $H^{2}(\Omega_{\varepsilon})$. This contradicts the assumption

[TABLE]

This proves the boundedness of the linear operator $\mathcal{L}_{\epsilon,{\bf p}}$.

To complete the proof of Proposition 4.1, we just need to show that conjugate operator to $\mathcal{L}_{\epsilon,{\bf p}}$ (denoted by $\mathcal{L}_{\epsilon,{\bf p}}^{*}$) is injective from $\mathcal{K}_{\epsilon,{\bf}}^{\perp}$ to $\mathcal{C}_{\epsilon,{\bf p}}^{\perp}$. The proof for $\mathcal{L}_{\epsilon,{\bf p}}^{*}$ follows almost the same process as for $\mathcal{L}_{\epsilon,{\bf p}}$ and therefore it is omitted.

The proof is complete.

9. 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

[TABLE]

where $u_{\epsilon}$ is the exact boundary cluster solution derived in Section 3-5 and $\chi_{\varepsilon}$ is the cutoff function given in (3.13). It is easy to see that

[TABLE]

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 $\partial\Omega$ 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 $\Omega$ 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:

[TABLE]

where $T_{R}[u_{0}^{2}]=v_{0}$ is the solution of the inhibitor equation

[TABLE]

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 $\alpha=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 $\frac{\partial}{\partial\alpha}$ in (9.1) and get

[TABLE]

As a preparation for this perturbation analysis, we represent the boundary $\partial\Omega$ 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 $\epsilon^{3}$.

Near the point $p_{i}$ we have expanded the boundary $\partial\Omega$ in Section 3.1 using Cartesian coordinates. We have derived that $\rho(0)=\rho^{\prime}(0)=0$ and $\rho^{\prime\prime}(0)$ is the curvature of $\partial\Omega$ at $p_{i}$. Recall that $\rho$ was used to flatten the boundary $\partial\Omega$ near $p_{i}$ to a line.

Using polar coordinates $(\phi,r)$ such that $x_{1}-p_{1}=r\sin\phi$, $x_{2}-p_{2}=R-r\sin\phi$, the boundary $\partial\Omega$ can be represented as

[TABLE]

The radius $R$ will be chosen such that $\rho^{\prime\prime}(0)=\frac{1}{R}.$ We now derive the function $f(\phi)$ from the function $\rho(x_{1}-p_{1})$ introduced in Section 3.1.

Substituting the Taylor expansion of $\rho$ given in equation (3.8) into $f(\phi)^{2}$ gives

[TABLE]

[TABLE]

We expand

[TABLE]

and compute the coefficients $\alpha,\,\beta,\,\gamma,\,\delta$ by matching powers $\phi^{i}$ for $i=1,2,3,4$. First we get

[TABLE]

Second we have

[TABLE]

which implies that

[TABLE]

and from now on we choose $R$ such that this condition is satisfied. Third we compute

[TABLE]

which gives

[TABLE]

Finally, we get

[TABLE]

which gives

[TABLE]

To summarise, we have

[TABLE]

where $f(\phi)$ denotes the radius and $\phi$ the angle and

[TABLE]

[TABLE]

The point $\phi=0$ with $f(0)=R$ and $f^{\prime}(0)=f^{\prime\prime}(0)=0$ corresponds to $p_{i}$. Note that $f(0)=R$ and $f^{\prime}(0)=0$ follow from $\rho(0)=\rho^{\prime}(0)=0$.

Further, we have $f^{\prime\prime}(0)=0$ due to the choice of the leading term $f(0)=R$. Since we have $f(\phi)=r$, we get $r=R+O(\phi^{3})$ which enables us to replace $r$ by $R$ in the derivation of $f(\phi)$.

In polar coordinates, for $r=R$ we get a point on $\partial B_{R}$. We change variables such that for the variables $(\phi,r^{\prime})$ at $r^{\prime}=R$ we get a point on $\partial\Omega$. Thus $r^{\prime}=R$ has to map into $r=f(\phi)$. This is achieved by defining $r^{\prime}=r-f(\phi)+R$ for each $\phi$.

Then the Laplacian is transformed as follows:

[TABLE]

using partial derivatives $\partial_{r}=\frac{\partial}{\partial r}$ etc. and

[TABLE]

[TABLE]

for the transformed variables $(\phi,r^{\prime})$.

We introduce rescaled variables $(\alpha,b)$ inside the spike such that

[TABLE]

Then in rescaled variables $(\alpha,b)$ we have

[TABLE]

This implies that in rescaled variables we get

[TABLE]

Using

[TABLE]

and

[TABLE]

we have

[TABLE]

and

[TABLE]

Thus second, third and fourth terms in the Laplacian are small and they can be estimated as follows: $O(\epsilon^{4}\alpha^{4})$, $O(\epsilon^{2}\alpha^{2})$ and $O(\epsilon^{2}\alpha)$, respectively.

Comparing with Cartesian coordinates $(z_{1},z_{2})$ for the $\epsilon$-scale inside the spike used in Section 3.1, by elementary trigonometry we get

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

Now we approximate the exact solution for the activator $u_{\epsilon}$ as follows:

[TABLE]

Near the centre $p_{i}$ of each spike we have

[TABLE]

where we have used the notation

[TABLE]

for $i=1,2,\ldots,k,\,j=1,2,3$ and $\chi_{\varepsilon}$ 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 $\rho^{\prime\prime}(0)=\frac{1}{R}$. Using polar coordinates, we can represent $u_{0}(\epsilon\alpha,R-r)$ and this function satisfies the Neumann boundary condition in a ball:

[TABLE]

However, the function $u_{0}$ does not satisfy the boundary condition at $\partial\Omega$ (and for $r=R$ we reach the boundary of the disk (circle) but not the domain boundary $\partial\Omega$). Recall that $r-f(\phi)=r^{\prime}-R$. Thus, to get a function with a better approximation to the boundary condition at $\partial\Omega$ (and such that for $r=R$ we reach $\partial\Omega$), we define

[TABLE]

Then for $r^{\prime}=R$, we have

[TABLE]

i.e. for $r^{\prime}=R$ the arguments of the function $\bar{w}$ are contained in $\partial\Omega$ and the arguments of $u_{0}$ are contained in $\partial B_{R}$. This means that the function $f$ deforms the boundary to a circle (in the same way as in Section 3.1 $\rho$ 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 of $\bar{w}$ as follows:

[TABLE]

The outward unit normal vector in polar coordinates is given by

[TABLE]

Using rescaled variables, this implies

[TABLE]

[TABLE]

[TABLE]

Subtracting the radial unit vector $e_{r}=(0,1)$, this implies

[TABLE]

Using that $\|\bar{w}_{i}-w\|_{H^{2}(\Omega_{\epsilon})}=O(\epsilon\alpha)$ and that $\bar{w}$ satisfies (9.7), the outward normal derivative of $\bar{w}$ is computed as

[TABLE]

Now we compute the terms $\epsilon^{3}\bar{v}^{(1)}$ (even around $\alpha=0$) and $\epsilon^{2}\bar{v}^{(2)}$ (odd around $\alpha=0$) in the expansion (9.5) near $p_{i}$ such that the solution satisfies the Neumann boundary condition to higher order.

Let $\bar{v}_{i}^{(1)}$ satisfy

[TABLE]

Let $\bar{v}_{i}^{(2)}$ be given by

[TABLE]

Substituting the expansion (9.4) of the solution in the activator equation, we get

[TABLE]

[TABLE]

We calculate for $x=p_{i}+z$

[TABLE]

where $\bar{R}_{1,i}(z)=2\bar{w}_{i}(z)v_{i}^{(1)}(z)$, $\bar{R}_{2,i}=2\bar{w}_{i}(z)v_{i}^{(2)}(z)$. Further, we recall that by (4.4) we have

[TABLE]

By the reduced problem in Section 5, we have

[TABLE]

Therefore we can add another contribution $\bar{v}_{i}^{(3)}$ to the solution such that $\varepsilon^{2}\bar{v}_{i}^{(3)}$ satisfies

[TABLE]

where

[TABLE]

[TABLE]

Setting

[TABLE]

and adding this part to the solution will cancel out the odd terms (with respect to $\alpha=0$) of order $\epsilon^{2}$ in the activator equation and we get

[TABLE]

[TABLE]

Taking the derivative $\frac{\partial}{\partial\alpha}$ in this relation near $x=p_{i}$, we compute

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

On the other hand,

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

Finally, by an expansion similar to (9.8) we have

[TABLE]

and

[TABLE]

9.2. Expansion of the eigenfunction

We define the approximate kernels to be

[TABLE]

Then we expand the eigenfunction as follows:

[TABLE]

[TABLE]

where $\phi_{\epsilon}^{\perp}\in\mathcal{K}_{\varepsilon,{\bf p}}^{\perp}$.

Suppose that $\|\phi_{\varepsilon}\|_{H^{2}(\Omega_{\varepsilon})}=1$. Then $|a_{j,\varepsilon}|\leq C$. Let us put

[TABLE]

Then for a subsequence and

[TABLE]

The decomposition of $\phi_{\varepsilon}$ in (9.11) implies that

[TABLE]

where $\psi_{i,\varepsilon}$ satisfies

[TABLE]

which we also write as $\psi_{i,\varepsilon}=T_{\tau\lambda_{\varepsilon}}[\phi_{i,\varepsilon}]$, and $\psi_{\varepsilon}^{\perp}$ is given by

[TABLE]

which we also represent as $\psi_{\varepsilon}^{\perp}=T_{\tau\lambda_{\varepsilon}}[\phi_{\varepsilon}^{\perp}]$.

Let us first consider the leading term of $\phi_{\varepsilon}$. For $\frac{\partial\bar{w}_{i}}{\partial\alpha}$ we get, using (9.8),

[TABLE]

Therefore, expanding the boundary condition as we did above for $\bar{w}_{i}$, we define $\bar{v}_{eig,i}^{(1)}$ as the unique solution of

[TABLE]

Similarly, let $\bar{v}_{eig,i}^{(2)}$ be the unique solution of

[TABLE]

Let us compare this with the derivative $\frac{\partial\bar{v}^{(1)}}{\partial\alpha}$ which satisfies

[TABLE]

Further, $\frac{\partial\bar{v}^{(2)}}{\partial\alpha}$ solves

[TABLE]

Using (9.10), we get

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

and the remainder $\bar{v}_{R,i}^{(1)}$ is given by the difference of the previous two contributions as follows:

$\bar{v}_{R,i}^{(1)}=\bar{v}_{eig,i}^{(1)}-\frac{\partial\bar{v}^{(1)}}{\partial\alpha}$ which satisfies

[TABLE]

We note that $\bar{v}_{R,i}^{(1)}$ is an even function around $\alpha=0$.

Substituting the decompositions of $\phi_{\varepsilon}$ and $\psi_{\varepsilon}$ into (6.1), we have

[TABLE]

Set

[TABLE]

[TABLE]

and

[TABLE]

We first compute

[TABLE]

Since

[TABLE]

we can estimate $I_{4}$ as follows:

[TABLE]

where we use the equation satisfied by $\psi_{i,\varepsilon}$.

9.3. Expansion of the small eigenvalues

Multiplying both sides of (9.2) by $\frac{\partial u_{l}}{\partial\tau(p_{l})}$ and integrating over $\Omega_{\varepsilon}$, we have

[TABLE]

Using the estimate $I_{4}$, we have

[TABLE]

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

[TABLE]

where $a_{i,\varepsilon}$ has been defined in (9.11).

Proof.

We first study the asymptotic behaviour of $\psi_{j,\varepsilon}$.

Note that for $l\neq k$, we have

[TABLE]

Next we compute $\psi_{l,\varepsilon}-\frac{\partial v}{\partial\tau(p_{l})}$ near $p_{l}$:

[TABLE]

so

[TABLE]

where we have used the relation $e^{-\frac{1}{\sqrt{D}}}\ll\varepsilon$.

Thus

[TABLE]

Combining (9.3) and (9.25), we have

[TABLE]

Similarly, we have the following:

[TABLE]

for $i\neq l$ and

[TABLE]

Thus we have for $J_{4,l}$,

[TABLE]

For $J_{6,l}$, we have from (9.26)

[TABLE]

Similarly, using (9.27), (9.28), (5.15) and (5), we get

[TABLE]

where we have used the relation

[TABLE]

Combining the above two estimates, we have

[TABLE]

and this proves (9.23).

Using

[TABLE]

and

[TABLE]

we can evaluate the integral $J_{5,l}$. A computation analogous to (5.6) gives

[TABLE]

where the constant $\nu_{1}>0$ has been defined in (5.7). Thus

[TABLE]

From the estimates on $I_{3}$ and $I_{4}$, we know that

[TABLE]

Next by the definition of $J_{1,l}$, using (9.30)

[TABLE]

We have obtained (9.19). Next we estimate $J_{2,l}$:

[TABLE]

By the equation satisfied by $\psi_{\varepsilon}^{\perp}$, we have

[TABLE]

Further,

[TABLE]

We have by the estimate for $\frac{\partial v}{\partial\tau(p_{l})}$,

[TABLE]

and

[TABLE]

Thus we have

[TABLE]

and (9.20) follows. Last we consider $J_{3,l}$,

[TABLE]

Thus (9.21) follows.

∎

10. 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.