Spike Solutions to the Supercritical Fractional Gierer-Meinhardt System

TL;DR

This paper investigates localized spike solutions in a one-dimensional supercritical fractional Gierer-Meinhardt system, employing asymptotic methods and stability analysis, and establishes rigorous existence results related to fractional Green's functions.

Contribution

It introduces a novel analysis of multi-spike solutions in a supercritical fractional GM system, combining formal asymptotics with rigorous existence and stability proofs.

Findings

Construction of multi-spike solutions via nonlinear algebraic systems

Stability analysis through eigenvalue problems

Rigorous existence of ground-state solutions near critical fractional order

Abstract

Localized solutions are known to arise in a variety of singularly perturbed reaction-diffusion systems. The Gierer-Meinhardt (GM) system is one such example and has been the focus of numerous rigorous and formal studies. A more recent focus has been the study of localized solutions in systems exhibiting anomalous diffusion, particularly with L\'evy flights. In this paper we investigate localized solutions to a one-dimensional fractional GM system for which the inhibitor's fractional order is supercritical. Using the method of matched asymptotic expansions we reduce the construction of multi-spike solutions to solving a nonlinear algebraic system. The linear stability of the resulting multi-spike solutions is then addressed by studying a globally coupled eigenvalue problem. In addition to these formal results we also rigorously establish the existence and stability of ground-state solutions when the inhibitor's fractional order is nearly critical. The fractional Green's function, for which we present a rapidly converging series expansion, is prominently featured throughout both the formal and rigorous analysis in this paper. Moreover, we emphasize that the striking similarities between the one-dimensional supercritical GM system and the classical three-dimensional GM system can be attributed to the leading order singular behaviour of the fractional Green's function.

Equations

$$(-\Delta {)}^s\phi (x)\equiv C_s{\int }_{-\infty }^{\infty }\frac{\phi (x)-\phi (\bar{x})}{|x-\bar{x}{|}^{1+2s}}d\bar{x},C_s\equiv \frac{2^{2s}s\Gamma (s+2^{-1})}{\sqrt{\pi }\Gamma (1-s)}.$$

$$[left=\text{empheqlbrace}]$$

$$\tau v_t+D(-\Delta {)}^{s_2}v+v-u^2=0,$$

$$u(x+2,t)=u(x,t),v(x+2,t)=v(x,t),$$

$$\begin{split}&(-\Delta)^{s}\varphi(x)\equiv C_{s}\int_{-1}^{1}K_{s}(x-\bar{x})\bigl{(}\varphi(x)-\varphi(\bar{x})\bigr{)}d\bar{x},\\ &K_{s}(z)\equiv\frac{1}{|z|^{1+2s}}+\sum_{j=1}^{\infty}\biggl{(}\frac{1}{|z+2j|^{1+2s}}+\frac{1}{|z-2j|^{1+2s}}\biggr{)}.\end{split}$$

$$(-\Delta {)}^{s_2}G+D^{-1}{G}_D=\delta (x),0<x<1,G_D(x+2)=G_D(x),$$

$$1/4<s_1<1,0<s_2<1/2.$$

$$[left=\text{empheqlbrace}]$$

$$U,V>0,$$

$$U,V\to 0,$$

$$\underset{{\epsilon }_0\to 0}{\lim }\left|{\tau }_s^{-1}U(x)-w(x)\right|=0,\underset{{\epsilon }_0\to 0}{\lim }\left|{\tau }_s^{-1}V(x)-1\right|=0,$$

$$(-\Delta {)}^{\frac{1}{2}}w+w-w^2=0,$$

$${\tau }_s={\left(\frac{\Gamma (1-2s)\sin (s\pi )}{\pi }{\int }_{\mathbb{R}}{w}^2(x)dx\right)}^{-1}.$$

$$[left=\text{empheqlbrace}]$$

$$(-\Delta {)}^s\psi -2U\varphi +\tau {\lambda }_s\psi =0,$$

$$[left=\text{empheqlbrace}]$$

$$D(-\Delta {)}^{s_2}v+v-u^2=0,$$

$$u(x+2)=u(x),v(x+2)=v(x).$$

$$[left=\text{empheqlbrace}]$$

$$U_c\sim \nu (S)|y{|}^{-(1+2s_1)},V_c\sim \mu (S)+S|y{|}^{2s_2-1},$$

$$(-\Delta {)}^s\varphi +\varphi =f,\text{for }-\infty <y<\infty ;\varphi \to 0,\text{as }|y|\to \infty ,$$

$$(-\Delta {)}^s\varphi =f,\text{for }-\infty <y<\infty ;\varphi \to 0,\text{as }|y|\to \infty ,$$

$$\underset{|x|\to \infty }{\lim }G(x)|x{|}^{1+2s}=C,$$

$$G_0(x,y)=\frac{1}{\pi }\Gamma (1-2s)\sin (s\pi )|x{|}^{2s-1}.$$

$$V_c(y;S)=C+{\mathfrak{a}}_{s_2}{\int }_{-\infty }^{\infty }|y-\bar{y}{|}^{2s_2-1}{U}_c(\bar{y};S{)}^{2}d\bar{y},{\mathfrak{a}}_{s_2}\equiv -2{\pi }^{-1}s\Gamma (-2s_2)\sin (\pi s_2).$$

$$S={\mathfrak{a}}_{s_2}{\int }_{-\infty }^{\infty }{U}_c(\bar{y};S{)}^{2}d\bar{y}.$$

$$\nu(S)=\lim_{y\rightarrow\infty}|y|^{1+2s_{1}}U_{c}(y,S),\qquad\mu(S)=\lim_{y\rightarrow\infty}\bigl{(}V_{c}(y;S)-S|y|^{2s_{2}-1}\bigr{)}.$$

$$U_c(y)\sim \sqrt{\frac{S}{{\mathfrak{b}}_{s_1}{\mathfrak{a}}_{s_2}}}{w}_{s_1}(y),V_c(y)\sim \sqrt{\frac{S}{{\mathfrak{b}}_{s_1}{\mathfrak{a}}_{s_2}}},\mu (S)\sim \sqrt{\frac{S}{{\mathfrak{b}}_{s_1}{\mathfrak{a}}_{s_2}}},$$

$${\mathfrak{b}}_{s_1}\equiv {\int }_{-\infty }^{\infty }{w}_{s_1}(y{)}^{2}dy,$$

$$[left=\text{empheqlbrace}]$$

$$w_{s_1}(y)=O(|y{|}^{-(1+2s_1)}),$$

$$u(x_{i}+\varepsilon y)\sim D\varepsilon^{-2s_{2}}\bigl{(}U_{i}(y)+o(1)\bigr{)},\quad v(x_{i}+\varepsilon y)\sim D\varepsilon^{-2s_{2}}\bigl{(}V_{i}(y;S_{i})+o(1)\bigr{)}.$$

$$(-\Delta {)}^{s_1}{U}_i+U_i-V_i^{-1}{U}_i^2=0,(-\Delta {)}^{s_2}{V}_i+D^{-1}{\epsilon }^{2s_2}{V}_i-U_i^2=0-(1+x_i)<\epsilon y<1-x_i.$$

$$U_i(y)\sim U_c(y;S_i)+o(1),V_i(y)\sim V_c(y;S_i)+o(1),$$

$$u^2\to {\epsilon }^{1-4s_2}{D}^2\sum _{i=1}^N{\int }_{-\infty }^{\infty }{U}_c(y;S_i{)}^{2}dy\delta (x-x_i)={\epsilon }^{1-4s_2}{D}^2{\mathfrak{a}}_{s_2}^{-1}\sum _{i=1}^N{S}_i\delta (x-x_i),$$

Full text

Keywords: Gierer-Meinhardt system, fractional Laplacian, Lévy flights, localized solutions, singular perturbation.

1. Introduction

Reaction diffusion systems have consistently been at the forefront of pattern formation research since Alan Turing’s seminal paper in 1952 [23] in which he demonstrated that sufficiently large differences in the diffusivities of reacting agents can lead to the formation of spatial patterns. By specifying reaction-diffusion systems either phenomenologically or from first principles, studies have used linear stability analysis to explore pattern formation in complex systems with applications to a variety of biological phenomena [16]. While these studies have traditionally assumed that individual agents undergo Brownian motion in which the mean squared displacement (MSD) is a linear function of the elapsed time, a recent growing body of literature has considered pattern formation in the context of anomalous diffusion in which there is alternative nonlinear relationships between the MSD and the elapsed time [9, 5, 31, 13]. Such anomalous diffusion may be better suited for describing the spatial distribution of agents in complex biological environments such as those found within individual cells [1, 19].

Of particular importance, and relevance to the present paper, is the case of anomalous superdiffusion with Lévy flights for which a heavy-tailed step-length distribution leads to an unbounded MSD. In this case the resulting fractional reaction-diffusion system features the fractional Laplacian which for one-dimensional problems is given by

[TABLE]

where $0<s<1$ and $\Gamma(z)$ is the Gamma function. A growing number of studies have considered such fractional-reaction diffusion systems with different reaction kinetics and using linear stability analysis have demonstrated that the introduction of anomalous diffusion can have a pronounced effect on pattern formation[5, 13]. Furthermore by considering parameters near the Turing stability threshold the authors in [13] used a weakly nonlinear analysis to investigate the resulting near equilibrium patterns that emerge after the onset of linear instabilities. An increasing number of studies have also analysed far from equilibrium solutions arising in the singularly perturbed limit where one of the diffusivities is asymptotically small [17, 30, 7, 15] and it is this latter thread of inquiry which we continue in this paper.

First proposed by Gierer and Meinhardt in 1972 [4], the Gierer-Meinhardt (GM) system is a canonical reaction diffusion system that addresses the consequences of short range activation and long range inhibition on pattern formation. In this paper we will be interested in the fractional GM system which is given by

[TABLE]

where $u(x,t)$ and $v(x,t)$ correspond to the activator and inhibitor concentrations respectively, $\tau>0$, $D>0$, and $0<\varepsilon\ll 1$. In this paper we impose periodic boundary conditions,

[TABLE]

which allows us to avoid the technical difficulties in assigning Dirichlet or Neumann boundary conditions [14] and with which (1.1) simplifies to

[TABLE]

When $s_{1}=s_{2}=1$ in (1.2) we recover the classical one-dimensional GM system which has been the focus of numerous studies. Specifically, in the classical case both rigorous [29] and formal asymptotic methods [12, 25] have been used to study the existence and stability of spike solutions. The analogous system in two- and three-dimensions have also been studied in [27, 28, 6], while the system posed on a two-dimensional Riemannian manifold was considered in [22]. In both rigorous and formal approaches the activator’s asymptotically small diffusivity leads to a separation of spatial scales which yields a tractable reduction of the full system. More recently spike solutions have been analysed for $s_{1}\in[1/2,1)$ and $s_{2}=1$ [17], $s_{1}=s_{2}\in[1/2,1)$ [30], $s_{1}\in(1/4,1)$ and $s_{2}\in(1/2,1)$ [7], and for $s_{1}=s_{2}=1/2$ [15]. The case of anomalous subdiffusion, in which the time derivative is of fractional order, was also considered in [18]. The individual treatment of each range of $s_{1}$ and $s_{2}$ values stems from differences in the singular behaviour of the outer inhibitor solution in a leading order theory. This is closely related to the limiting behaviour of the Green’s function $G(x)$ satisfying

[TABLE]

as $x\rightarrow 0$. In particular as discussed further in §2.4 and Appendix A below, the Green’s function remains bounded for subcritical values of $s_{2}\in(1/2,1]$ but has a logarithmic singularity at the critical value of $s_{2}=1/2$. There is a suggestive analogy here with the singular behaviour of the classical free space Green’s function in one- and two-dimensions which we elaborate in §2.4 below.

In this paper we use the method of matched asymptotic expansions to study the existence and stability of spike solutions to (1.2) in the parameter regime

[TABLE]

As discussed in §2.4 and Appendix A, the Green’s function in this regime has an algebraic (and in some cases an additional logarithmic) singularity as $x\rightarrow 0$ and for this reason we refer to the resulting fractional GM system as being supercritical. Consequently the asymptotic analysis of spike solutions in this regime is analogous to that found in the classical three dimensional Schnakenberg [24] and Gierer-Meinhard systems [6]. Using the method of matched asymptotic expansions we thus construct multi-spike solutions by deriving an appropriate nonlinear algebraic system (NAS). From the NAS we identify two distinguished parameter regimes: the $D=O(1)$ regime and the $D=O(\varepsilon^{2s_{2}-1})$ regime. Whereas in the former the NAS admits only symmetric solutions, we find that in the latter it admits both symmetric and asymmetric solutions. The stability of the resulting multi-spike solutions can then be determined by studying a globally coupled eigenvalue problem (GCEP). From the GCEP we deduce that asymmetric solutions are always linearly unstable while symmetric solutions are susceptible to two types of bifurcations: competition instabilities (see Figure 1(a) for an example), and Hopf bifurcations. In addition to these bifurcation which occur over an $O(1)$ timescale, otherwise stable multi-spike solutions can also undergo drift motion over an $O(\varepsilon^{2s_{2}-3})$ timescale (see Figure 1(b) for an example). This paper thus fully characterizes the equilibrium solutions to (1.2) and their linear stability, while also identifying key parameter regimes for the diffusivity $D$.

The bulk of this paper uses formal asymptotic methods to characterize localized solutions as discussed in the preceding paragraph. While a rigorous justification of these results remains open for the full range of values $0<s_{2}<1/2$ there are some results that we can rigorously prove when $s_{2}<1/2$ is close to $s_{2}=1/2$. Specifically, for such values of $s_{2}$ we can rigorously prove the existence and stability of ground state solutions to the core problem considered in §2. Specifically we have the following theorem.

Theorem 1.1**.**

There exists an $\varepsilon_{0}>0$ such that for each $s\in\left(\frac{1}{2}(1-\varepsilon_{0}),\frac{1}{2}\right)$ the core problem

[TABLE]

admits a solution $(U,V)$ such that

[TABLE]

uniformly in compact sets in $x$. Here $U$ is the ground state solution of

[TABLE]

and

[TABLE]

Interestingly, we can also study the stability and instability of the ground state solution constructed in Theorem 1.1. Writing the associated eigenvalue problem for the system as

[TABLE]

where $\lambda_{s}\in\mathbb{C}.$ Here we say $(U,V)$ is linearly stable if the real part of each eigenvalue is negative, while $(U,V)$ is called linear unstable if there exists a $\lambda_{s}$ such that its real part $\Re(\lambda_{s})>0.$

Theorem 1.2**.**

Let $(U,V)$ be the solution constructed in Theorem 1.1. There exists $\tau_{1}$ such that the solution is linearly stable for any $\tau<\tau_{1}.$

The remainder of this paper is organized as follows. In §2 we construct multi-spike quasi-equilibrium solutions by first considering the relevant core problem in §2.1 and then deriving the NAS in §2.2. This is followed by §2.3 where we specifically consider symmetric and asymmetric solutions in the $D=O(\varepsilon^{2s_{2}-1})$ regime and by §2.4 where we discuss in more detail the singular behaviour of the Green’s function and its connection with higher dimensional problems. In §3 we study the linear stability of multi-spike solutions by deriving the GCEP and focusing in particular on the $D\ll O(\varepsilon^{2s_{2}-1})$ and $D=O(\varepsilon^{2s_{2}-1})$ regimes in §3.1 and §3.2 respectively. This is followed by §4 where we derive an ordinary differential equation (ODE) system governing the slow dynamics of multi-spike solutions. In §5 we then perform full numerical simulations of (1.2) to validate our asymptotic theory. In §6 we prove Theorems 1.1 and 1.2. Finally in §7 we summarize our results and make some concluding remarks.

2. Asymptotic Approximation of $N$-Spike Quasi-Equilibrium Solutions

In this section we will use the method of matched asymptotic expansions to calculate asymptotic approximations of $N$-spike solutions to

[TABLE]

The successful use of the method of matched asymptotic expansions relies on the asymptotically small activator diffusivity $\varepsilon^{2s_{1}}\ll 1$ which leads to the emergence of two distinct length scales. Specifically the activator concentrates at $N$ points $-1<x_{1}<...<x_{N}<1$ that are well separated in the sense that $|x_{i}-x_{j}|\gg\varepsilon$ for all $i\neq j$ as well as $x_{1}+1\gg\varepsilon$ and $1-x_{N}\gg\varepsilon$. Over an $O(\varepsilon)$ length scale centred at each $x_{1},...,x_{N}$ the system (2.1) is approximated by a core problem in $\mathbb{R}$ whose solutions yields the local profile of the activator and inhibitor. This core problem depends on an undetermined spike strength parameter whose value determines the far-field behaviour of the core solution. On the other hand over an $O(1)$ length scale away from each spike location $x_{1},...,x_{N}$ the nonlinear term appearing in (2.1b) can be approximated, in the sense of , by a sum of appropriately weighted Dirac delta functions centred at each $x_{1},..,x_{N}$. As a consequence the inhibitor can be approximated as a weighted sum of Green’s functions over an $O(1)$ length scale. By matching the behaviour of this sum of Green’s functions as $x$ approaches each spike location with the far-field behaviour of each core solution we can then derive a NAS of $N$ equations in the $N$ undetermined spike strength parameters. The method of matched asymptotic expansions therefore reduces the original PDE system (2.1) to a finite number of nonlinear algebraic equations whose solutions yields an asymptotic approximation of an $N$-spike solution.

Guided by the preceding discussion, in the remainder of this section we will first discuss the core problem and highlight some of its key properties. We will then use the method of matched asymptotic expansions as outlined above to derive the NAS. The remainder of the section will then be dedicated to a discussion on the existence of symmetric and asymmetric $N$-spike solutions as well as to some of the peculiarities of the fractional Gierer-Meinhardt system which distinguish it from the classical one- and three-dimensional Gierer-Meinhardt systems.

2.1. The Core Problem

The core problem is one of the key ingredients in deriving an asymptotic approximation of an $N$-spike quasi-equilibrium solutions as outlined above. It is given by

[TABLE]

where $S>0$ is a parameter which we refer to as the spike strength while $\nu(S)$ and $\mu(S)$ are two $S$-dependent constants. Solutions to (2.2) will be denoted by $U_{c}(y;S)$ and $V_{c}(y;S)$ to make explicit the dependence on the parameter $S$. The core problem (2.2) is a leading order approximation of (2.1) after the rescaling $y=\varepsilon^{-1}(x-x_{i})$ and its solutions yield the local profile of each spike in an $N$-spike quasi-equilibrium solution of (2.1). The far-field behaviour of $U_{c}(y;S)$ and $V_{c}(y;S)$ is a consequence of the following lemma

Lemma 2.1**.**

Let $0<s<1/2$ and suppose that $f(y)=O(|y|^{-\sigma})$ as $|y|\rightarrow\infty$ for $\sigma>0$.

1. (1)

If $\sigma>1+2s$ then the solution to

[TABLE]

decays like $\phi\sim C|y|^{-1-2s}$ as $|y|\rightarrow\infty$. 2. (2)

If $\sigma>1$ then the solution to

[TABLE]

decays like $\phi\sim C|y|^{2s-1}$ as $|y|\rightarrow\infty$.

Proof.

The conclusion follows easily by using classical potential analysis and the decay properties of the Green’s functions associated with the operators $(-\Delta)^{s}+I$ and $(-\Delta)^{s}$. Specifically the Green’s function $G(x,y)$ of $(-\Delta)^{s}+I$ has the asymptotic behaviour

[TABLE]

for some constant $C>0$, while the Green’s function $G_{0}(x,y)$ of $(-\Delta)^{s}$ has the form

[TABLE]

We refer the readers to [30, Section 2] and [20, Section 1.12] for the proof of (2.3) and (2.4) respectively. ∎

We can in fact be more explicit about the solution $V_{c}$ of (2.2) by taking the Fourier transform of the second equation in (2.2a) to get

[TABLE]

Taking the limit as $|y|\rightarrow\infty$ then yields as a special case of Lemma 2.1 the limiting behaviour $V_{c}(y;S)\sim C+\mathfrak{a}_{s_{2}}|y|^{2s_{2}-1}\int_{-\infty}^{\infty}U_{c}(\bar{y};S)^{2}d\bar{y}$. Comparing this with the far-field behaviour of the core solution given in (2.2b) we deduce the useful identity

[TABLE]

which in particular reinforces our assumption that $S>0$ since $\mathfrak{a}_{s_{2}}>0$ for $s_{2}<1/2$.

In light of the above discussion the specification of the parameter $S$ is equivalent to fixing the $L^{2}(\mathbb{R})$ norm of $U_{c}$. By solving (2.2) for a fixed value of $S$ we can then extract the values of the far-field constants $\nu(S)$ and $\mu(S)$ by taking the limits

[TABLE]

The nonlinearity in the first equation of (2.2a) implies that we must have $V_{c}(y;S)>0$ for all $y\in\mathbb{R}$ and this leads us to the constraint $\mu(S)\geq 0$. We next have to determine whether there any values of $S>0$ for which this constraint holds. To address this we first consider the small $S$-asymptotics. Specifically if $S\ll 1$ then (2.6) implies that $U_{c}(y;S)=O(\sqrt{S})$ and by balancing terms in (2.2) we also deduce that $V_{c}(y;S)=O(\sqrt{S})$ and $\mu(S)=O(\sqrt{S})$. It is then straightforward to see that to leading order in $S\ll 1$ we have the asymptotic expansions

[TABLE]

and $w_{s_{1}}(y)$ is the fractional homoclinic solution satisfying

[TABLE]

We refer the reader to Section 4 in [30] for further properties of the nonlinear problem (2.9). The small-$S$ asymptotics (2.8a) imply that $\mu(S)>0$ for $0<S\ll 1$. A numerical continuation in $S$ then further extends the range of $S$ values for which $\mu(S)>0$ holds (see Appendix B.2 for details).

Plots of the numerically calculated far-field constants $\mu(S)$ and $\nu(S)$ are shown in Figure 2. These plots indicate that there exists a value of $S=S_{\star}>0$ beyond which $\mu(S)>0$ no longer holds. In Figure 3(a) we plot $S_{\star}$ as a function of $s_{2}$ for select values of $s_{1}$. In addition the plots in Figure 2 indicate that $\mu(S)$ attains a unique global maximum in $0<S<S_{\star}$ at some value $S=S_{\text{crit}}$ which we plot for distinct values of $s_{1}$ in Figure 3(b). This critical value of $S=S_{\text{crit}}$ plays a crucial role in the leading order stability theory of multiple spike solutions as will be further discussed in §3 below. Finally, in Figure 3(c) we plot the profiles of the core solutions for select values of $S>0$ when $s_{1}=0.5$ and $s_{2}=0.4$.

We conclude by remarking that our preceding discussion has thus far been limited to numerical calculations of solutions to the core problem (2.2). In §6 we rigorously prove the existence and stability of ground state solutions, i.e. those for which $\mu(S)=0$, when $s_{2}\approx 1/2$. For more general values of $s_{2}<1/2$ the rigorous justification of such solutions remains an open problem. We remark also that the existence of a ground state is not guaranteed as can be seen, for example, in the case of the core problems associated with the three-dimensional Gray-Scott, Schnakenberg, and Brusselator systems [6].

2.2. Asymptotic Matching and the Nonlinear Algebraic System

We now consider the asymptotic construction of an $N$ spike solution to (2.1). Assuming that the $N$-spikes concentrate at $N$ well separated (in the sense made precise above) points $-1<x_{1}<...<x_{N}<1$ we begin by making the ansatz that

[TABLE]

A simple change of variables then yields that $U_{i}$ and $V_{i}$ must satisfy

[TABLE]

Approximating the domain $-\varepsilon^{-1}(1+x_{i})<y<\varepsilon^{-1}(1-x_{i})$ with $-\infty<y<\infty$ and dropping the $D^{-1}\varepsilon^{2s_{2}}$ term in the $V_{i}$ equation we deduce that

[TABLE]

where $U_{c}(y;S)$ and $V_{c}(y;S)$ are the solutions of the core problem (2.2) discussed and $S_{i}>0$ is an as-of-yet undetermined constant. Implicit in the asymptotic approximation (2.12) is the assumption that the inner profiles interact only through the far-field behaviour constants $S_{i}$, the nature of which is revealed by formulating the outer problem and deriving an appropriate matching condition.

Next we derive an outer problem valid for values of $-1<x<1$ such that $|x-x_{i}|\gg\varepsilon$ for all $i=1,...,N$. We first make note of the limit

[TABLE]

as $\varepsilon\rightarrow 0^{+}$ which is to be understood in the sense of distributions and for which we have used (2.6) in the equality. The outer inhibitor solution must then be a $2$-periodic function satisfying

[TABLE]

We let $G_{D}(x)$ be the the $2$-periodic fractional Green’s function satisfying

[TABLE]

which can be written as (see Appendix A)

[TABLE]

where $\mathfrak{a}_{ks_{2}}\equiv-2ks_{2}\pi^{-1}\Gamma(-2ks_{2})\sin(\pi ks_{2})$ and $R_{D}(x)$ is given explicitly by (A.7b). In terms of this Green’s function the solution to (2.13a) can be explicitly written as

[TABLE]

Comparing the limiting behaviour of (2.16) as $x\rightarrow x_{i}$ with the limiting behaviour (2.13b) from the inner solution yields the algebraic equation

[TABLE]

The $S_{i}|x-x_{i}|^{2s_{2}-1}$ term on the left-hand-side cancels the $k=1$ term in the sum on the right-hand-side while the remaining singular terms corresponding to $k=2,...,k_{\max}$ are cancelled out by higher order corrections to the inner solution. On the other hand the constant term $\mu(S)$ on the left-hand-side must be balanced with the constant terms appearing on the right-hand-side. Since this must hold for each value of $i=1,...,N$ we are thus led to the NAS

[TABLE]

Note that this NAS must in general be solved numerically since $\mu(S)$ can only be computed numerically (see §2.1). It nevertheless provides a substantial reduction in the construction of multi-spike solutions to the equilibrium equation (2.1).

We remark that the NAS (2.18) is $\varepsilon$-dependent and yields distinct leading order approximations depending on whether $D=O(1)$ or $D=D_{0}\varepsilon^{2s_{2}-1}$ where $D_{0}=O(1)$. Indeed in the former case (2.18) implies that $\mu(S_{i})=0$ to leading order and hence $S_{i}\sim S_{\star}+O(D\varepsilon^{1-2s_{2}})$. On the other hand if $D=D_{0}\varepsilon^{2s_{2}-1}\gg 1$ then the asymptotics

[TABLE]

imply that $S_{1},..,S_{N}>0$ must solve the leading order system

[TABLE]

for each $i=1,...,N$ with the next order correction being $O(D_{0}^{-1}\varepsilon^{1-2s_{2}})$. The shape of $\mu(S)$ illustrated in Figure 2 suggests the possibility that $S_{1},...,S_{N}\in\{S_{l},S_{r}\}$ for some $0<S_{l}<S_{r}<S_{\star}$. Thus whereas the $D=O(1)$ regime supports solutions in which the profiles of each spike are identical, i.e. the $N$-spike solution is symmetric, the $D=D_{0}\varepsilon^{2s_{2}-1}$ regime may admit both symmetric and asymmetric $N$-spike solutions which we discuss further in §2.3 below.

While the leading order approximations discussed above are suggestive of the solutions we may encounter it is important to highlight that their associated errors are $O(D\varepsilon^{1-2s_{2}})$ when $D=O(1)$, and $O(D_{0}^{-1}\varepsilon^{1-2s_{2}})$ when $D=D_{0}\varepsilon^{2s_{2}-1}$. Although these errors are small in the limit $\varepsilon\rightarrow 0$ they may in practice be unacceptably large. For example if $\varepsilon=0.01$ and $s_{2}=0.4$ then $\varepsilon^{1-2s_{2}}\approx 0.4$. In contrast if we solve the $\varepsilon$-dependent NAS (2.18) directly then the next order correction to the inner problem can be deduced from the matching condition (2.17) and is either $O(D^{-1}\varepsilon^{2s_{2}})$ if $s_{2}<1/4$ or $O(D\varepsilon^{2-2s_{2}})$ if $1/4<s_{2}<1/2$. In particular this yields an $O(\varepsilon)$ error when $D=D_{0}\varepsilon^{2s_{2}-1}$ and For this reason we will be using the $D=D_{0}\varepsilon^{2s_{2}-1}$ regime when we perform numerical simulations of (1.2) in §5 below.

2.3. Symmetric and Asymmetric Solutions in the $D=D_{0}\varepsilon^{2s_{2}-1}$ Regime

As discussed above, the leading order equation of the NAS (2.18) in the $D=D_{0}\varepsilon^{2s_{2}-1}$ regime given by (2.20) admits both symmetric and asymmetric $N$-spike solutions. In the following section we will explore these two types of solutions in more detail while also drawing parallels to the analogous solutions encountered in the case of the three-dimensional Gierer-Meinhardt system [6].

Symmetric $N$-spike solutions are perhaps the easiest to analyze since in this case the spike strengths are all equal, $S_{1}=...=S_{N}=S_{c}$ and the leading order NAS (2.20) reduces to the scalar equation

[TABLE]

From the plots of $\mu(S)$ in Figure 2 it is clear that $S_{c}\rightarrow S_{\star}$ as $D_{0}\rightarrow 0$ thereby providing a connection between the $D=O(1)$ and $D=O(\varepsilon^{2s_{2}-1})$ regimes. On the other hand as $D_{0}\rightarrow\infty$ we obtain $S\rightarrow 0$ and in particular using the small-$S$ asymptotics (2.8a) we find that $S\sim(\mathfrak{b}_{s_{1}}\mathfrak{a}_{s_{2}}\kappa^{2})^{-1}$. In Figure 4(a) we plot $S_{c}$ versus $\kappa$ for a selection of $s_{1}$ and $s_{2}$ values.

In addition to symmetric $N$-spike solutions the leading order NAS (2.20) and plots of $\mu(S)$ in Figure 2 further suggest the possibility of asymmetric $N$-spike solution. Specifically, recalling that $0<S_{\text{crit}}<S_{\star}$ is the value where $\mu(S)$ attains its unique maximum we deduce that for any $S_{r}\in[S_{\text{crit}},S_{\star})$ there is a unique $S_{l}(S_{r})\in(0,S_{\text{crit}}]$ which we plot for‘ $s_{1}=0.45$ and a selection of $s_{2}$ values in Figure 4(b). Notice from its definition that $S_{l}(S_{\text{crit}})=S_{\text{crit}}$ whereas $S_{l}(S_{\star})=0$. Moreover, by differentiating $\mu(S_{l}(S(r))=\mu(S_{r})$ we obtain $S_{l}^{\prime}(S_{r})=[\mu^{\prime}(S_{l}(S_{r}))]^{-1}\mu^{\prime}(S_{r})$ so that in particular $S_{l}^{\prime}(S_{r})\rightarrow 0$ as $S_{r}\rightarrow S_{\star}$ due to the small $S$ asymptotics (2.8a). Plots of $S_{l}^{\prime}(S_{r})$ in Figure 4(c) further indicate that $-1\leq S_{l}^{\prime}(S_{r})\leq 0$.

We next consider the construction of asymmetric $N$-spike solutions consisting of $1\leq n\leq N-1$ large and $N-n$ small spikes by letting

[TABLE]

where $\sigma$ is a permutation of $\{1,...,N\}$. With this assumption the leading order system (2.20) reduces to the scalar equation

[TABLE]

This scalar equation was previously encountered in the classical 3D Gierer-Meinhardt model [6]. For that model two key properties of $\mu(S)$ and $S_{l}(S_{r})$ allowed for a complete characterization of the bifurcation structure of (2.22), the first being that $\mu^{\prime}(S)<0$ for $S_{\text{crit}}<S<S_{\star}$, and the second that $-1<S_{l}^{\prime}(S_{r})<0$ for all $S_{\text{crit}}<S_{r}<S_{\star}$. Since these properties likewise hold for the $\mu(S)$ and $S_{l}(S_{r})$ in our present case we will simply state the results from [6], referring the interested reader to Section 2.3 of [6] for more details. The first result states that if

[TABLE]

then (2.22) has a unique solution for any $1\leq n\leq N-1$. In addition if $n\geq N-n$ then (2.22) does not have a solution for any $\kappa\geq\kappa_{c1}$. If on the other hand $n<N-n$ then (2.22) has exactly two distinct solutions for

[TABLE]

where $S_{\text{crit}}<S_{r}^{\star}<S_{\star}$ is the unique solution to

[TABLE]

and no solutions if $\kappa\geq\kappa_{c2}$.

2.4. On the Fractional Green’s Function

The preceding sections have highlighted the importance of the fractional Green’s function satisfying (2.14) in the asymptotic construction of quasi equilibrium solutions. We conclude this section by highlighting some of the key properties of the fractional Green’s function and relating them to the behaviour of the classical Green’s function in one-, two-, and three-dimensions.

The limiting behaviour of $G_{D}(x)$ as $x\rightarrow 0$ plays a crucial role in the existence and stability of multi-spike solutions. Interestingly this behaviour is markedly different when $s_{2}\in(1/2,1]$, $s_{2}\in(0,1/2)\setminus\{\tfrac{1}{2r}\,|\,r\in\mathbb{Z},\,r\geq 1\}$, and $s_{2}\in\{\tfrac{1}{2r}\,|\,r\in\mathbb{Z},\,r\geq 1\}$. In particular when $s_{2}\in(1/2,1)$ the Green’s function is not singular with $G_{D}(x)\sim\mathfrak{a}_{s_{2}}|x|^{2s_{2}-1}+O(1)$ as $x\rightarrow 0$ [7]. On the other hand, referring to Propositions A.1 and A.2 in Appendix A, we have

[TABLE]

where $k_{\max}=\lceil\tfrac{1}{2s_{2}}-1\rceil$. The singular behaviour in each of these cases has direct analogies with the singular behaviour of the non-fractional Green’s function in one-, two, and three-dimensions. Specifically, we may view the fractional Green’s function as analogous to the one-, two-, and three-dimensional non-fractional Green’s function when $s_{2}\in(1/2,1)$, $s_{2}=1/2$, and $s_{2}\in(0,1/2)\setminus\{\tfrac{1}{2r}\,|\,r\in\mathbb{Z},\,r\geq 1\}$ respectively. This analogy further extends to the methods used in the analysis of spike solutions as is evident by the similarities between the analysis in [7] for $s_{2}\in(1/2,1)$ and the classical Gierer-Meinhardt system (e.g. in [12]), that in [15] for $s_{2}=1/2$ and the two-dimensional Gierer-Meinhardt system [27, 28], and that in the present paper with the analysis of spike solutions in the three-dimensional Schnakenberg [24] and Gierer-Meinhardt [6] systems. For the remaining values of $s_{2}\in\{\tfrac{1}{2r}\,|\,r\in\mathbb{Z},\,r\geq 2\}$ the mixing between logarithmic and algebraic singularities leads to problems which don’t appear to have a clear classical analog. The analysis of the fractional Gierer-Meinhardt system for these remaining parameter values is not addressed in this paper but is an interesting direction for future research.

We now consider the regular part of the Green’s function $R_{D}(x)$ which can be computed using the series expansion (A.7b). Numerical calculations indicate that $R_{D}(0)>0$ for all $D>0$ when $s_{2}\in(\tfrac{1}{2(r+1)},\tfrac{1}{2r})$ for even values of $r\geq 1$, whereas there is a threshold $D_{R}(s_{2})>0$ for which $R_{D}(0)<0$ for all $D<D_{R}(s_{2})$ when $s_{2}\in(\tfrac{1}{2(r+1)},\tfrac{1}{2r})$ for odd values of $r\geq 1$. The threshold $D_{R}(s_{2})$ can be numerically computed using (A.7b) and is plotted in Figure 5. Note that special care must be taken when using the series (A.7b) as $s_{2}\rightarrow\tfrac{1}{2r}^{+}$ for any integer $r\geq 2$. Specifically, in this case $k_{\max}=\lceil\tfrac{1}{2s_{2}}-1\rceil<r$ so that in the limit the second term in (A.7b) does not converge. This is easily fixed by letting $k_{\max}=r$ when $s_{2}$ is sufficiently close to $\tfrac{1}{2r}$ and we use this in our numerical computations. In addition our numerical calculations indicate that $D_{R}(s_{2})\rightarrow+\infty$ as s_{2}\rightarrow\bigl{(}\tfrac{1}{2(r+1)}\bigr{)}^{+} and s_{2}\rightarrow\bigl{(}\tfrac{1}{2r}\bigr{)}^{-} for odd values of $k\geq 1$. This diverging behaviour can be explicitly characterized by balancing dominant terms in the series (A.7b). Specifically by noting that $\Gamma(-z)\sim-(1-z)^{-1}$ as $z\rightarrow 1$ we deduce that $\mathfrak{a}_{rs_{2}}\sim\pi^{-1}(1-2rs_{2})^{-1}$ as $s_{r}\rightarrow\tfrac{1}{2r}$. Assuming that $D=D_{R}(s_{2})\gg 1$ in (A.7b) and balancing dominant terms then implies that

[TABLE]

which we plot using dashed lines in Figure 5. This diverging behaviour of the threshold $D_{R}(s_{2})$ is suggestive of an alternative scaling that arises in this limit. Indeed the appearance of a logarithmic singularity at values of $s_{2}=\tfrac{1}{2r}$ for integer $r\geq 1$ suggests that an additional small parameter $\nu=-\tfrac{1}{\log\varepsilon}$ must be incorporated into the asymptotic theory, likely leading to alternative distinguished asymptotic regimes for the diffusivity.

We conclude by remarking that the negativity of $R_{D}(0)$ for $D<D_{R}(0)$ and $s_{2}\in(\tfrac{1}{2(r+1)},\tfrac{1}{2r})$ for even values of $r\geq 1$ poses a challenge to the application of our asymptotic theory. Indeed, recalling the NAS (2.18) we observe that $R_{D}(0)<0$ in this regime contradicts the positivity of $\mu(S)$. By restricting our attention to the case $D=D_{0}\varepsilon^{2s_{2}-1}$ this difficulty can be circumvented at least in theory when $\varepsilon\ll 1$ and for which $D>D_{R}(s_{2})$. However to validate our asymptotic theory with numerical simulations we have to use a finite value of $\varepsilon>0$ which may lead to $D<D_{R}(s_{2})$ especially as $s_{2}$ approaches any of the values for which $D_{R}(s_{2})$ diverges. Such behaviour is not in the range of validity of our asymptotic theory and we will henceforth ignore it though we would be remiss to not at least suggest approaches for handling this issue. One possibility is to develop a higher order asymptotic theory though this falls out of the scope of this paper. An alternative approach is to consider a $\varepsilon$-dependent core problem (2.2) posed on the truncated domain $|y|<L/\varepsilon$ for some $L>0$ in which case negative values of $\mu(S)$ are permissible provided that the solution $V_{c}$ remains positive. This approach however has two major shortcomings: it requires an appropriate assignment for $V_{c}$ in $|y|\geq L/\varepsilon$ in order to have a well-posed problem, and the core problem will need to be recomputed anew for different values of $\varepsilon$.

3. Linear Stability: The Large, $O(1)$, Eigenvalues

In this section we consider the linear stability on an $O(1)$ timescale of the $N$-spike equilibrium solution $u_{e}$ and $v_{e}$ constructed in Section 2 above. We proceed by substituting into (1.2) the perturbed solutions $u=u_{e}+e^{\lambda t}\phi$ and $v=v_{e}+e^{\lambda t}\psi$ where $|\phi|,|\psi|\ll 1$. To linear order in $\phi$ and $\psi$ we then have the spectral problem

[TABLE]

for which we seek $\lambda=O(1)$ eigenvalues. If $\Re(\lambda)>0$ (resp. $\Re(\lambda)<0$) then the $N$-spike equilibrium solution is linearly unstable (resp. stable) and we will commonly refer to such an eigenvalue as being unstable (resp. stable). Noting that the diffusivity $\varepsilon^{2s_{1}}\ll 1$ appearing in (3.1a) is asymptotically small we will use the method of matched asymptotic expansions to derive a globally coupled eigenvalue problem (GCEP) from which distinct modes of instabilities and their respective thresholds can be determined.

For each $i=1,..,N$ and $y=O(1)$ we begin by substituting

[TABLE]

into (3.1) to get

[TABLE]

Assuming that $D\gg O(\varepsilon^{2s_{2}})$ and exploiting the homogeneity of this system we obtain the leading order asymptotic expansion

[TABLE]

where $\Phi_{c}^{\lambda}(y;S)$ and $\Psi_{c}^{\lambda}(y;S)$ satisfy

[TABLE]

where we assume the general far-field behaviour

[TABLE]

The undetermined constants $c_{1},...,c_{N}$ correspond to distinct instability modes and moreover yield additional degrees of freedom with which we can normalize the behaviour of solutions to (3.3). We can solve (3.3b) explicitly as

[TABLE]

Substituting this back into (3.3) then results in the inhomogeneous equation

[TABLE]

Observe that if $B(\lambda,S)=0$ then $\lambda$ is an eigenvalue of $\mathscr{M}(S)$ and $\Phi_{c}^{\lambda}(y;S)$ the corresponding eigenfunction. If $\lambda$ is not an eigenvalue of $\mathscr{M}$ then we can uniquely solve (3.6a) for $\Phi_{c}$ which gives

[TABLE]

In addition we make note of the far-field behaviour

[TABLE]

in which the second term vanishes if $\Phi_{c}^{\lambda}(\cdot;S)$ is odd. On the other hand if $\Phi_{c}^{\lambda}(y;S)$ is not odd then using the additional degrees of freedom granted by $c_{1},...,c_{N}$ we can normalize $\Phi_{c}^{\lambda}(y;S)$ such that

[TABLE]

with which we get the far-field behaviour

[TABLE]

Note in addition that such a normalization fixes $B(\lambda,S)$which we obtain by multiplying (3.7) by $U_{c}(y;S)$ and integrating to get

[TABLE]

We next consider the distributional limit

[TABLE]

from which we observe that any $i\in\{1,...,N\}$ corresponding to odd-valued $\Phi_{c}^{\lambda}(y;S_{i})$ will not contribute to the outer problem. A modification of the proceeding calculations in which we keep track of such odd-valued $\Phi_{c}(\cdot,S_{i})$ reveals that such terms do not contribute to the linear stability over an $O(1)$ timescale, though they do contribute to drift instabilities considered in §4 below. Without loss of generality we therefore assume that none of the $\Phi_{c}^{\lambda}(y;S_{i})$ ($i=1,...,N$) are odd-valued. Using the normalization (3.8) we thus obtain the outer problem

[TABLE]

for each $i=1,...,N$. The solution to (3.11a) can then be expressed in terms of the Green’s function satisfying (2.14) as

[TABLE]

Using (2.15) the matching condition (3.11b) then becomes

[TABLE]

as $x\rightarrow x_{i}$ for each $i=1,...,N$. The leading order singular behaviour immediately balances whereas balancing the leading order constants for each $i=1,...,N$ yields the GCEP

[TABLE]

In the following subsections we consider the leading order behaviour of the GCEP (3.13) when $D\ll O(\varepsilon^{2s_{2}-1})$ and when $D=D_{0}\varepsilon^{2s_{2}-1}$. This leading order behaviour will provide insights into the modes of instabilities arising in each of these asymptotic regimes. However, as in the case of the NAS (2.18) analyzed in §2 we remind the reader that the errors in such leading order approximations will typically be unacceptably large for moderately small value of $\varepsilon>0$. Therefore when we perform full numerical simulations of (1.2) to support our asymptotic predictions in §5 we will be numerically computing the relevant stability thresholds from the $\varepsilon$-dependent GCEP (3.13) directly.

3.1. Linear Stability in the $D\ll O(\varepsilon^{2s_{2}-1})$ Regime

We consider first perhaps the simplest case which is when $D\ll O(\varepsilon^{2s_{2}-1})$ or $D=O(1)$ in particular. From our discussion in §2 we know that in this case all $N$-spike solutions are symmetric to leading order in $\varepsilon\ll 1$ with $S_{1}=...=S_{N}=S_{\star}$. Moreover in this regime the GCEP (3.13) reduces to the single scalar equation $B(\lambda,S_{\star})=0$. Therefore $\lambda$ must be an eigenvalue of the operator $\mathscr{M}(S_{\star})$ defined in (3.6b) above with the far-field asymptotics (3.4).

By numerically calculating the spectrum of $\mathscr{M}$ as outlined in Appendix B we have observed that the dominant eigenvalue is always stable when $S=S_{\star}$. In Figures 6(a)-6(c) we plot the three largest eigenvalues of $\mathscr{M}$. Note that $\lambda=0$ is always an eigenvalue of $\mathscr{M}$ but that this corresponds to the translational mode $\Phi=\partial U_{c}/\partial y$ and $\Psi=\partial V_{c}/\partial y$ whose analysis is deferred to §4 below. Therefore $\lambda_{1}$ is the appropriate eigenvalue for which (3.4) is satisfied when $S=S_{\star}$ and the plot of $\Re\lambda_{1}$ at $S=S_{\star}$ versus $0.2<s_{2}<0.5$ for select values of $s_{1}$ in Figure 6(d) indicates that this eigenvalue is always stable. In summary, when $D\ll O(\varepsilon^{2s_{2}-1})$ all $N$-spike solutions are linearly stable to leading order in $\varepsilon\ll 1$.

3.2. Linear Stability in the $D=O(\varepsilon^{2s_{2}-1})$ Regime

In this section we consider the case when $D=D_{0}\varepsilon^{2s_{2}-1}$ and for which we will consider the case $D_{0}\rightarrow\infty$ as a special case. Using the large $D$ asymptotics of the Green’s function (2.19) the GCEP (3.13) becomes to leading order in $\varepsilon\ll 1$

[TABLE]

where $\boldsymbol{e}=(1,\cdots,1)^{T}$ and where we remind the reader that $\kappa=ND_{0}/(2\mathfrak{a}_{s_{2}})$. In this section we will consider the linear stability of both the symmetric and asymmetric solutions described in §2.3. We demonstrate that the symmetric $N$-spike solutions are susceptible to two types of instabilities: oscillatory instabilities arising through a Hopf bifurcation, and competition instabilities arising through a zero eigenvalue crossing. On the other hand we will show that asymmetric solutions are always linearly unstable with respect to competition instabilities. The proceeding analysis closely follows previous work done on the three-dimensional Gierer-Meinhardt model [6] with its successful adaptation to the present one-dimensional fractional case being due to the properties of $\mu(S)$ described in §2.1.

3.2.1. The Shadow Limit $D_{0}\rightarrow\infty$

Before analyzing (3.14) in general we first consider the shadow limit obtained by letting $D_{0}\rightarrow\infty$. As discussed in §2.3 all $N$-spike solutions are then symmetric with $S_{c}\sim(\mathfrak{b}_{s_{1}}\mathfrak{a}_{s_{2}}\kappa^{2})^{-1}\ll 1$. Moreover by using the small $S$ asymptotics (2.8a) and the definition of $\mathscr{M}$ given in (3.6b) we readily deduce that

[TABLE]

From (3.10) we obtain

[TABLE]

with which (3.14) becomes

[TABLE]

Note that the shadow limit case is independent of $s_{2}$. If $N\geq 2$ then this equation is satisfied if $\boldsymbol{c}$ is any competition mode satisfying $c_{1}+...+c_{N}=0$ and $\lambda$ is the dominant eigenvalue of $\mathscr{L}$. Since the dominant eigenvalue, $\Lambda_{0}$, of $\mathscr{L}$ has a positive real part (see Section 4 of [3]) we therefore deduce that multi-spike solutions in the $D_{0}\rightarrow\infty$ are always linearly unstable. We refer to the resulting instabilities as competition instabilities since the condition $c_{1}+...+c_{N}=0$ leads to the growth of some spikes at the expense of the decay of others. If on the other hand $N=1$ then (3.16) becomes the scalar nonlocal eigenvalue problem (NLEP)

[TABLE]

Following the arguments used in the classical one-dimensional Gierer-Meinhardt system in [26] it can be shown that there is a Hopf bifurcation threshold $\tau_{h}^{\infty}$ such that all eigenvalues are stable if $\tau<\tau_{h}^{\infty}$ whereas there is exactly one complex conjugate pair of unstable eigenvalues when $\tau>\tau_{h}^{\infty}$. To calculate this Hopf bifurcation threshold we substitute the purely imaginary eigenvalue $\lambda=i\lambda_{I}$ into (3.17) and isolate real and imaginary parts to get the system

[TABLE]

We can then find the Hopf bifurcation threshold by numerically solving the first equation for $\lambda_{I}=\lambda_{h}^{\infty}(s_{1})$ and then substituting into the second equation to get a value for the Hopf bifurcation threshold $\tau=\tau_{h}^{\infty}(s_{1})$ (for plots of $\tau_{h}^{\infty}$ and $\lambda_{h}^{\infty}$ see Figure 1A of [7]).

In summary, when $D_{0}\rightarrow\infty$ multi-spike solutions are always linearly unstable due to competition instabilities whereas single spike solutions are linearly stable provided $\tau$ does not exceed the numerically calculated Hopf bifurcation threshold $\tau=\tau_{h}^{\infty}(s_{1})$. We now address the question of what happens to these competition instability and Hopf bifurcation thresholds for symmetric $N$-spike solutions when $D_{0}$ is finite.

3.2.2. Stability Threshold for Symmetric Solutions

We now consider the linear stability of symmetric $N$-spike solutions for which we remind the reader that $S_{1}=...=S_{N}=S_{c}$ where $S_{c}$ satisfies (2.21). To determine the linear stability of these solutions with respect to competition modes we first let $\boldsymbol{c}$ satisfy $c_{1}+...+c_{N}=0$. It follows that (3.14) reduces to $B(\lambda,S_{c})=0$ so that $\lambda$ is an eigenvalue of $\mathscr{M}$ whose eigenfunction satisfies the far-field behaviour (3.4). Numerical calculations of the spectrum of $\mathscr{M}$ indicate that this eigenvalue is positive if $S<S_{\text{crit}}$ whereas it is negative if $S_{\text{crit}}<S<S_{\star}$ (see Figure 6(a)-6(c)). From (2.21) and the plots of $\mu(S)$ in Figure 2 we therefore conclude that symmetric $N$-spike solutions are linearly stable with respect to competition instabilities if $\kappa<\kappa_{c_{1}}$ and linearly unstable otherwise. Recall that $\kappa_{c1}=\mu(S_{\text{crit}})/S_{\text{crit}}$ was previously encountered in (2.23) when considering the existence of asymmetric solutions. From the definition of $\kappa$ we can alternatively express this threshold for $\kappa$ as a threshold for the diffusivity

[TABLE]

As in the classical Gierer-Meinhardt model (and other singularly perturbed reaction diffusion systems) the stability of multi-spike solutions decreases as the number of spikes increases. In Figure 7(a) we plot the leading order competition instability threshold $\kappa_{c_{1}}$ versus $s_{2}$ for several values of $s_{1}$. From which we observe that the competition instability threshold is monotone decreasing in $s_{1}$. Moreover the threshold decreases monotonically with $s_{2}$ for $s_{1}>0.5$ whereas we see that for $s_{1}<0.5$ it is non-monotone, increasing for smaller values of $s_{2}$ and then decreasing.

Since $c_{1}+...+c_{N}=0$ spans an $(N-1)$-dimensional subspace of $\mathbb{R}^{N}$ it remains only to consider the synchronous modes $\boldsymbol{c}$ for which $c_{1}=...=c_{N}$. By substituting such a synchronous mode $\boldsymbol{c}$ into (3.14) we get

[TABLE]

First we show that $\lambda=0$ is not a solution of (3.20). Differentiating the core problem (2.2) with respect to $S$ we first make the observation that $B(0,S_{c})=\mu^{\prime}(S_{c})$ so that after solving (2.21) for $\kappa$ (3.20) becomes

[TABLE]

for which we claim the left-hand-side is strictly negative. This is clearly true for $S_{c}\geq S_{\text{crit}}$ since $\mu^{\prime}(S_{c})<0$ (see Figure 2). On the other hand for $S_{c}<S_{\text{crit}}$ the claim follows by observing that the derivative of the left-hand-side is $\mu^{\prime\prime}(S_{c})<0$ whereas the small $S$ asymptotics (2.8a) imply $S_{c}\mu^{\prime}(S_{c})-\mu(S_{c})\sim-\tfrac{1}{2}\sqrt{S_{c}/(\mathfrak{b}_{s_{1}}\mathfrak{a}_{s_{2}})}<0$ as $S_{c}\rightarrow 0^{+}$. Therefore instabilities with respect to the synchronous mode must arise through a Hopf bifurcation. Seeking purely imaginary eigenvalues $\lambda=i\lambda_{I}$ and separating the real and imaginary parts of (3.20) we obtain the system

[TABLE]

which we can numerically solve for $\lambda_{I}=\lambda_{h}(S_{c},s_{1},s_{2})$ from the first equation and then calculate the Hopf bifurcation threshold $\tau=\tau_{h}(S,s_{1},s_{2})$ from the second equation. The first equation is numerically solved using Newton’s method by slowly increasing $S_{c}$ starting from a small value for which the shadow-limit value $\lambda_{h}^{\infty}(s_{1},s_{2})$ provides an accurate initial guess. The resulting (leading order) Hopf bifurcation thresholds $\tau_{h}(S_{c},s_{1},s_{2})$ and associated eigenvalue $\lambda_{h}(S_{c},s_{1},s_{2})$ are plotted in Figures 7(b) and 7(c). In all cases we observe that the Hopf bifurcation threshold diverges toward $+\infty$ as $S_{c}\rightarrow S_{\text{crit}}^{-}$ and this is a consequence of the nonlocal operator $\mathscr{M}$ having a zero eigenvalue for this value of $S_{c}$. As discussed in §5 the Hopf bifurcation threshold can be extended beyond this critical value of $S_{c}=S_{\text{crit}}$ but this requires calculating the Hopf bifurcation threshold from the $\varepsilon$-dependent GCEP (3.13) directly.

3.2.3. Asymmetric $N$-Spike Solutions are Always Unstable

We conclude this section on the leading order stability of multi-spike solutions by adapting the analysis for the three-dimensional Gierer-Meinhardt model [6] to show that the asymmetric solutions of §2.3 are linearly unstable. The analysis follows closely that previously done in [6] so we provide only an outline, highlighting the key properties of $\mu(S)$ which allow the adaptation of the analysis in [6].

The key idea in showing that the asymmetric solutions are always linearly unstable is to construct specific modes $\boldsymbol{c}$ for which an instability is guaranteed. Assuming without loss of generality that $S_{1}=...=S_{n}=S_{r}>S_{\text{crit}}$ and $S_{n+1}=...=S_{N}=S_{l}(S_{r})$ the leading order GCEP (3.14) becomes

[TABLE]

where $\mathcal{I}_{n}$ is the $n\times n$ identity matrix and $\mathcal{O}_{n,m}$ is the $n\times m$ zero matrix. If $1\leq n\leq N-2$ then the mode $\boldsymbol{c}$ with $c_{1}=...=c_{n}=0$ and $c_{n+1}+...+c_{N}=0$ is immediately seen to be unstable since (3.22) reduces to $B(\lambda,S_{l}(S_{r}))=0$ and $S_{l}(S_{r})<S_{\text{crit}}$ implies the dominant eigenvalue of $\mathscr{M}$ with the far-field behaviour (3.4) is unstable (see Figures 6(a)-6(c)). Thus, competition between the $N-n$ small spikes is always destabilizing.

Since the modes considered above are trivial when $n=N-1$ a different argument must be used to show the instability of asymmetric solutions in this case. In particular if $n\geq N-n$ then it it was shown in [6] that unstable modes of the form $c_{1}=...=c_{n}=c_{r}$ and $c_{n+1}=...=c_{N}=c_{l}$ can always be found. The argument used in [6] relies on the some key properties of $\mu(S)$ and the spectrum of $\mathscr{M}$. First it requires that $\mu^{\prime}(S_{l})>0$, $\mu^{\prime}(S_{r})<0$, and $S_{l}^{\prime}(S_{r})>-1$ all of which our numerical calculations indicate are satisfied in the present case (see Figures 2 and 4(c)). Second it requires that the eigenvalues of $\mathscr{M}$ satisfying the appropriate far-field behaviour (3.4) are stable for $S_{r}>S_{\text{crit}}$. Since this condition is likewise satisfied (see for example Figures 6(a)-6(c)) we are able to adapt the argument from [6] and therefore conclude that all asymmetric $N$-spike solutions are linearly unstable.

4. Slow Spike Dynamics and the Equilibrium Configurations

It is well known that localization solutions to a variety of singularly perturbed reaction diffusion systems exhibit slow dynamics [11, 24, 6]. Similar behaviour has likewise been observed for the fractional Gierer-Meinhardt system in one-dimension when $s_{2}>1/2$ [7]. In this section we establish that these slow dynamics persist in the case $1/4<s_{2}<1/2$ albeit at a different time scale. The dynamics in this parameter regime share qualitative similarities with their classical counterparts in one-, two-, and three-dimensions. Specifically the dynamics are determined by the gradient of the Green’s function which leads to a mutual repulsion between spikes. The derivation of the slow dynamics however more closely resembles that for the three-dimensional Gierer-Meinhardt and Schnakenberg systems [6, 24] owing to, as in previous sections, the strong coupling between the activator and inhibitor in the inner region. In this section we formally derive the equations governing the slow dynamics of a multi-spike solution and in §5.3 we validate our theory with numerical examples of two-spike solutions (see also Figure 1(b) in §1).

Slow spike dynamics are the result of higher order corrections so we begin by first substituting $x=x_{i}+\varepsilon y$ into (2.16) to obtain the higher order expansion

[TABLE]

where we have defined

[TABLE]

Note that due to the periodic boundary conditions we have $R_{D}^{\prime}(0)=0$. This implies, as subsequent calculations will show, that the dynamics of individual spikes are independent of their absolute position in the interval $-1<x<1$ but are due solely to interactions between spikes. Next we refine the inner expansion (2.10) by letting

[TABLE]

where $U_{i}^{\varepsilon}\equiv U_{c}(y;S_{i}^{\varepsilon})$, $V_{i}^{\varepsilon}\equiv V_{c}(y;S_{i}^{\varepsilon})$, $|\Phi_{i}^{\varepsilon}|\ll U_{i}^{\varepsilon}$ and $|\Psi_{i}^{\varepsilon}|\ll V_{i}^{\varepsilon}$, and where h.o.t. refers to higher order terms whose order will become evident after the asymptotic expansions are carried out. Substituting (4.2) into (1.2) we find that $\boldsymbol{\Phi}_{i}^{\varepsilon}\equiv(\Phi_{i}^{\varepsilon},\Psi_{i}^{\varepsilon})^{T}$ satisfies $\mathcal{L}_{i}^{\varepsilon}\boldsymbol{\Phi_{i}}^{\varepsilon}=\boldsymbol{f}_{i}^{\varepsilon}$, where

[TABLE]

We observe that $(\tfrac{d}{dy}U_{i}^{\varepsilon},\tfrac{d}{dy}V_{i}^{\varepsilon})^{T}$ is in the kernel of $\mathcal{L}_{i}^{\varepsilon}$ and assume that the kernel of $(\mathcal{L}_{i}^{\varepsilon})^{T}$ is likewise one-dimensional and spanned by $\boldsymbol{P}_{i}^{\varepsilon}\equiv(P_{i}^{\varepsilon},Q_{i}^{\varepsilon})^{T}$. We can then impose the solvability condition

[TABLE]

Numerical calculations indicate that $P_{i}^{\varepsilon}$ and $Q_{i}^{\varepsilon}$ are odd so that the final term vanishes and therefore

[TABLE]

Together with the NAS (2.18) this constitutes a differential algebraic system for the $N$ spike locations $x_{1},...,x_{N}$ and their strengths $S_{1}^{\varepsilon},...,S_{N}^{\varepsilon}$.

We immediately observe that (4.4) implies that the slow dynamics occur over a slow $O(\varepsilon^{2s_{2}-3})$ timescale. Furthermore since $P_{i}^{\varepsilon}$ is odd in $y$ we also have

[TABLE]

If $x_{i}\gtrless x_{j}$ then $G_{D}^{\prime}(x_{i}-x_{j})\lessgtr 0$ and we thus conclude that spikes are mutually repulsing. In particular it is easy to see from (4.4) that two-spike solutions are stationary if and only if $|x_{1}-x_{2}|=1$. In §5.3 we compare the slow-dynamics predicted by the differential algebraic system (4.4) and (2.18) with numerical simulations of (1.2) for two spike solutions that are initially separated by a distance $|x_{1}(0)-x_{2}(0)|<1$.

We conclude this section by outlining how to calculate the function $P_{i}^{\varepsilon}$ needed to evaluate the coefficient appearing in (4.4). Following the analysis of §3 we write

[TABLE]

from which we deduce that $P_{i}^{\varepsilon}$ solves $\mathscr{M}^{\star}(S_{i}^{\varepsilon})P_{i}^{\varepsilon}=-2C_{i}^{\varepsilon}U_{i}^{\varepsilon}$ where we define the adjoint operator $\mathscr{M}^{\star}=\mathscr{M}^{\star}(S)$ by

[TABLE]

Numerical calculations (not included) indicate that the adjoint operator $\mathscr{M}^{\star}$, like $\mathscr{M}$ in §3, has exactly one zero eigenvalue for all $0<S<S_{\star}$ except at $S=S_{\text{crit}}$ for which it has exactly two zero eigenvalues. In particular assuming $S_{i}^{\varepsilon}\neq S_{\text{crit}}$ we may set $C_{i}^{\varepsilon}=0$ and thus deduce that $P_{i}^{\varepsilon}$ is in the kernel of $\mathscr{M}^{\star}(S_{i}^{\varepsilon})$.

5. Numerical Simulations

In this section we numerically simulate the fractional Gierer-Meinhardt system (1.2) to support our asymptotic calculations in the preceding section. Using the asymptotically constructed solutions from §2 as initial conditions we choose parameter values to support the stability thresholds found by numerically solving (2.18). We proceed in three parts. In the first we consider Hopf bifurcations of single spike solutions, in the second we consider competition instabilities of two-spike solutions, and in the third and final part we consider the slow dynamics of two-spike solutions. In the first two parts we will first numerically compute the corresponding $\varepsilon$-dependent stability thresholds and compare them with their leading order counterparts. As emphasized in §3, due to the fractional powers of $\varepsilon$ in the asymptotic expansions of the stability thresholds we anticipate that the leading order thresholds deviate substantially from those obtained by solving (2.18) directly. Finally, when considering the slow-dynamics of two-spike solutions in the third part we will choose parameter values for which the two-spike solutions are linearly stable with respect to Hopf and competition instabilities.

5.1. Hopf Bifurcation of One-Spike Solutions

We first verify the Hopf bifurcation threshold for a single spike solution centred, without loss of generality, at $x=0$. With $N=1$ the NAS (2.18) and GCEP (3.13) become

[TABLE]

where we remind the reader that $R_{D}(x)$ is given by (A.7b). For a given value of $D_{0}$ we first solve (5.1a) for $S_{c}=S_{c}^{\varepsilon}$. Separating real and imaginary parts in (5.1b) we can then numerically solve for the Hopf bifurcation threshold $\tau=\tau_{h}^{\varepsilon}$ and accompanying eigenvalue $\lambda_{I}=\lambda_{h}^{\varepsilon}$. Specifically we solve the resulting system with Newton’s method starting with a large value of $D_{0}$ for which the shadow limit solutions $\tau_{h}^{\infty}$ and $\lambda_{I}^{\infty}$ are good initial guesses. Using $\varepsilon=0.01$ the resulting $\varepsilon$-dependent Hopf bifurcation thresholds are shown in Figures 8(a)-8(d) which illustrate the persistence of the Hopf bifurcation threshold for $S>S_{\text{crit}}$ not captured by the leading order theory. To support our asymptotically predicted threshold we performed several numerical simulations of the full system (1.2) with $\varepsilon=0.01$ and using a single spike solution centred at the origin as an initial condition. In Figures 8(e)-8(h) we plot $u(0,t)$ when $s_{1}=0.5$ and $s_{2}=0.34$ for select values of $D_{0}$ and values of $\tau$ slightly below and slightly above the Hopf bifurcation threshold, all of which validate the Hopf bifurcations thresholds from the asymptotic theory.

5.2. Competition Instabilities of Two-Spike Solutions

Turning our attention now to the case of a symmetric $N=2$-spike solution we perform numerical simulations to verify the onset of competition instabilities as predicted by our stability theory. We assume that $|x_{1}-x_{2}|=1$ so that there are no small eigenvalues or, equivalently, there are no slow dynamics as discussed in §4. In this case the NAS (2.18) and GCEP (3.13) with $\lambda=0$ become

[TABLE]

We can numerically solve this system for $D_{0}$ as a function of $s_{2}$ at select values of $s_{1}$. Doing so with $\varepsilon=0.01$ we obtain the higher order competition instability threshold shown in Figures 9(a)-9(d). In contrast to the leading order competition threshold which can be calculated as in §3.2.2 there is an upper limit to the value of $s_{2}$ for which we can compute the higher order $\varepsilon$-dependent threshold from (5.2). This is a consequence of the change in sign of $R_{D}(0)$ for smaller values of $D$ as described in §2.4. For sufficiently small values of $\varepsilon$ the value of $D=D_{0}\varepsilon^{2s_{2}-1}$ will always exceed this threshold and a competition instability threshold $D_{0,\text{comp}}^{\varepsilon}$ can be calculated for values of $s_{2}$ closer to $1/2$. Otherwise higher order correction terms need to be calculated or the inhibitor in the numerical discretization of the core problem (2.2) needs to be allowed to become negative as described in §2.4. We will not address these additional technical difficulties further.

To support our asymptotically calculated higher order competition instability threshold we performed several numerical experiments. In each experiment we use the methods of §2 to asymptotically construct a symmetric two-spike solution with spikes centred at $x_{1}=-0.5$ and $x_{2}=0.5$. Using this solution as the initial condition we then solve (1.2) numerically for select values of $s_{1}$ and $s_{2}$ and with small value of $\tau=0.05$ (so that there are no Hopf bifurcations) as well as values of $D=D_{0}\varepsilon^{2s_{2}-1}$ such that $D_{0}$ is either slightly below or slightly above the numerically calculated competition instability threshold $D_{0,\text{comp}}^{\varepsilon}$. In each case we found good agreement with the higher order calculated threshold $D_{0,\text{comp}}^{\varepsilon}$ and in Figure 9(e)-9(h) we show a sampling of numerically calculated values of the spike heights $u(x_{1},t)$ and $u(x_{2},t)$ for values of $D_{0}=0.95D_{0,\text{comp}}^{\varepsilon}$ (top) and $D_{0}=1.05D_{0,\text{comp}}^{\varepsilon}$ (bottom).

5.3. Slow Dynamics of Two-Spike Solutions

We conclude the numerical validation of our asymptotic theory by considering the slow dynamics of symmetric two-spike solutions. Using the translational invariance granted by the periodic boundary conditions we reduce the differential algebraic system (4.4) and (2.18) to the pair of scalar equations

[TABLE]

We remind the reader that the NAS (second equation) determines the common spike strength $S_{c}^{\varepsilon}$ for a given spike separation distance $|x_{2}-x_{1}|$. The common spike strength is then used to solve (2.2) for $U_{c}^{\varepsilon}$ and $V_{c}^{\varepsilon}$ as well as to solve for $P_{c}^{\varepsilon}$ in the adjoint problem of §4. We implement this system numerically by pre-computing $S_{c}^{\varepsilon}$ as a function of $0\leq|x_{2}-x_{1}|\leq 2$ and then computing each of $U_{c}^{\varepsilon}$, $V_{c}^{\varepsilon}$, and $P_{c}^{\varepsilon}$ as functions of $0<S_{c}^{\varepsilon}<S_{\star}$. The differential algebraic system can then be easily solved with any standard ordinary differential equation library (we used solve_ivp from the SciPy integrate library).

To validate our asymptotic theory we performed multiple numerical simulations of (1.2) with an initial condition consisting of a symmetric two-spike solution constructed using the methods in §2 where the spikes are concentrated at $x_{1}=-0.2$ and $x_{2}=0.2$. For each of our simulations we set $\varepsilon=0.01$ and used values of $\tau=0.1$ and $D=0.8D_{0,\text{comp}}^{\varepsilon}\varepsilon^{2s_{2}-1}$ with which we can avoid Hopf bifurcations and competition instabilities (see Sections 5.1 and 5.2). These simulations were completed for the pairs $(s_{1},s_{2})=(0.4,0.35),(0.4,0.3),(0.46,0.3),$ and $(0.46,0.35)$ and the resulting spike trajectories are shown as solid blue lines in Figure 10. In each of these plots the trajectories predicted by solving (5.3) are shown as dashed orange lines. We see that in each case the asymptotics provide good qualitative agreement of the spike trajectories. Finally we direct the reader to Figure 1(b) where have plotted in more detail the time evolution of the activator to accompany Figure 10(a).

6. Rigorous Results for $s_{2}\approx 1/2$

In this section we shall rigorously study the existence and stability of the ground state solution to the core problem

[TABLE]

We proceed by first presenting in § 6.1 several known results which will be used throughout this section. Then in § 6.2 and § 6.3 we provide the rigorous study on the existence and stability of the ground state solution respectively.

6.1. Preliminaries

Lemma 6.1**.**

Let $s<\frac{1}{2}$ and $G(x)$ be the Green’s function of the equation

[TABLE]

Then

[TABLE]

Proof.

Using the Fourier transform, we can write the (6.2) as

[TABLE]

Therefore, we have

[TABLE]

where we used $\int_{0}^{\infty}\frac{\cos\xi}{\xi^{2s}}d\xi=\Gamma(1-2s)\sin(s\pi)$ if $2s<1.$ ∎

We introduce the transformation

[TABLE]

where $w$ is the unique ground state solution to

[TABLE]

In this case we can give the explicit form of $w$ and the integral of $w^{2}$ on the real line

[TABLE]

Based on (6.5) we can write (6.1) as

[TABLE]

We look for a solution to (6.8) in the form $U=w+\phi$ with $\phi$ being a lower order term. Denoting by $T(h)$ the unique solution of the equation

[TABLE]

for $h\in L^{\infty}({\mathbb{R}})$, then formally we have

[TABLE]

where $h.o.t.$ indicates the higher order terms. We denote $v_{w}=T(w^{2})$ so that using the Green’s function given in Lemma 6.1 we have

[TABLE]

Expanding the Green function in the following way

[TABLE]

and then using the fact $\Gamma(1-2s)\sim(1-2s)^{-1}$ as $s\to\frac{1}{2}$ we get that

[TABLE]

As a consequence when $x$ is bounded we have

[TABLE]

then the nonlinear term of the first equation in (6.8) can be written as

[TABLE]

Substituting it into the first equation of (6.8) we get

[TABLE]

where

[TABLE]

and

[TABLE]

and represents the higher order terms in $\phi.$

Concerning the ground state $w$ and the non-local linearized operator $L$, we have the following result

Proposition 6.1**.**

Let $w$ be the unique, positive, radially symmetric solution to (6.6).

• (a)

Let $L_{0}=(-\Delta)^{\frac{1}{2}}+(1-2w)id$. Then we have

[TABLE]

• (b)

*Let $L$ be the lineazied operator defined in (6.16) and *

[TABLE]

Then

[TABLE]

Proof.

The proof of part (a) is given by [3, Proposition 1.1 and Theorem 2.3]. To prove part (b) we first notice that $L_{0}w=-w^{2}$. If $\phi\in\mbox{Ker}(L)$ then

[TABLE]

Therefore, by conclusion (a) we get $\phi-c(\phi)w\in\mbox{Ker}(L_{0})$ and in particular

[TABLE]

for some constant $\beta$. As a consequence, we have

[TABLE]

which implies $c(\phi)=0$. Hence, $\phi\in\mbox{Ker}(L_{0})$ and we get that $\phi\in\mathrm{Span}\left\{\frac{dw}{dx}\right\}$. Similarly, if $\phi\in\mathrm{Ker}(L^{*})$ then

[TABLE]

Using the fact

[TABLE]

we have

[TABLE]

Then

[TABLE]

where we used

[TABLE]

Thus $c_{1}(\phi)=0$ and $\phi\in\mathrm{Span}\left\{\frac{dw}{dx}\right\}$ which proves the third conclusion. ∎

In the end of this subsection, we provide the analysis of the linear operator $L$ in a framework of weighted $L^{\infty}$ spaces. For this purpose we consider the following norm for a function defined on ${\mathbb{R}}$. We define

[TABLE]

Given a function $h$ with $\|h\|_{*}<\infty$, due to the fact that $\mathrm{Ker}(L)={{Span}}\left\{\frac{dw}{dx}\right\}$, we need to study the related linear problem in the following form

[TABLE]

Our aim is to find $(\phi,c)$ such that (6.21) holds. Concerning (6.21) we have the following existence result and a-priori estimate for which a proof can be found in [30, Theorem 4.2].

Theorem 6.1**.**

If $h$ satisfies $\|h\|_{*}<\infty$ then problem (6.21) has an unique solution $\phi=\mathcal{T}(h)$ and $c=c(h)$. Moreover there exists a constant $C>0$ such that for any such $h$

[TABLE]

6.2. The rigorous proof of the existence results

In this section we shall give rigorous proof of Theorem 1.1.

6.2.1. Error estimates

We begin by studying $v_{w}(x)$ for which we prove improved estimates. By the definition (6.11) we see that

[TABLE]

where $\tau_{s}$ is given in (6.5). We will consider $v_{w}(x)$ in the two disjoint regions $x\in I_{s}$ and $x\in\mathbb{R}\setminus I_{s}$ where we define the interval

[TABLE]

Starting with $x\in I_{s}$ we use the Green representation formula (6.11) and the asymptotics (6.12) to get

[TABLE]

Next we define

[TABLE]

and readily deduce that $H_{i}(x)$ is even since $w(y)$ is even. Furthermore, as $|x|$ is sufficiently large, by standard potential analysis we can write

[TABLE]

where $f$ is an even function and itself with its first derivative are uniformly bounded.

While for $|x|\geq 100(1-2s)^{-1}$, using the potential analysis, we get that

[TABLE]

Summarizing the above estimates, we have the following conclusion.

Lemma 6.2**.**

Letting $v_{w}$ be defined as in (6.11) we have the following estimates:

• (a).

If $x\in I_{s}$, then

[TABLE]

• (b).

If $x\in{\mathbb{R}}\setminus I_{s}$, then

[TABLE]

for some constant $c>0$.

We now focus on estimating the quantity $S(w)=-(-\Delta)^{\frac{1}{2}}w-w+v_{w}^{-1}w^{2}$ which, using (6.6), can be rewritten as

[TABLE]

Let us first analyze the term $S(w)$ in the interval $I_{s}$ introduced in (6.24). It is easy to see that in this region we have

[TABLE]

where $\delta$ is any small positive number, and therefore

[TABLE]

Writing

[TABLE]

we deduce that for $x\in I_{s}$

[TABLE]

On the other hand, by Lemma 6.2 we find that for $x\in\mathbb{R}\setminus I_{s}$

[TABLE]

In conclusion, we have

Lemma 6.3**.**

Let $\mu=2$ in the definition of $\|\cdot\|_{*}$. If $1-2s$ is sufficiently small then we have

[TABLE]

where $C$ is some constant independent of $\varepsilon$ and $\delta$ is any small positive number independent of $\varepsilon$.

6.2.2. The existence of solution

Recall that the original problem was cast in the form

[TABLE]

Rather than solving (6.34) directly we consider instead the problem of finding $A$ satisfying

[TABLE]

for a certain constant $c$, and such that $\langle A-w,Z\rangle=0$. Rewriting $A=w+\phi$ we get that this problem is equivalent to

[TABLE]

and

[TABLE]

Using the operator $\mathcal{T}$ introduced in Theorem 6.1, we see that the problem is equivalent to finding a $\phi\in\mathcal{H}$ so that

[TABLE]

We shall show that this fixed point problem has a unique solution in the region of the form

[TABLE]

for any small positive constant $\delta$, provided that $1-2s$ is sufficiently small. Here

[TABLE]

We have already proved that $\|S(w)\|_{*}\leq C(1-2s)^{1-\delta}$. In the following lemma we estimate the higher order error term $N(\phi)$.

Lemma 6.4**.**

Assume that $\phi\in\mathcal{D}$, then for $1-2s$ sufficiently small, we have

[TABLE]

where $\sigma(1-2s)\leq C(1-2s)^{1-\delta}$ as $1-2s\to 0$.

Proof.

Let us assume first $x\in{\mathbb{R}}\setminus I_{s}$. In this region we have $w(x)\leq C\rho(x)$. Combined with the standard potential analysis one can show that

[TABLE]

As a consequence,

[TABLE]

Therefore we have

[TABLE]

provided $s\to\frac{1}{2}$.

Considering next the case $x\in I_{s}$ we decompose $N(\phi)$ in the form

[TABLE]

where

[TABLE]

and

[TABLE]

It is known that

[TABLE]

and

[TABLE]

and in particular $|T(w\phi)|=O(\|\phi\|_{*})$. Likewise, $T(\phi^{2})=O(\|\phi\|_{*}^{2})$. Combining these facts we obtain

[TABLE]

A similar analysis yields

[TABLE]

and therefore

[TABLE]

for $x\in I_{s}$. Together with (6.41) this proves the lemma. ∎

With Lemma 6.4 we are able to give the proof of Theorem 1.1.

Proof of Theorem 1.1..

Using the definition of the corresponding norms, repeating almost the same arguments as Lemma 6.4 one can prove that if $\|\phi_{i}\|_{*}\leq C(1-2s)^{1-\delta}$ for $i=1,2$, then, given any small $\kappa\in(0,1)$, we have the following inequality

[TABLE]

provided $1-2s$ is sufficiently small. As a consequence, we get that the operator $Q$ is a contraction mapping in the set $\mathcal{D}$ defined in (6.38). On the other hand, we also get from Lemma 6.4 that $Q$ maps $\mathcal{D}$ into itself. Thus, by using the Banach fixed point theorem, we get the existence of a unique fixed point of $Q$ in $\mathcal{D},$ that is,

[TABLE]

Next, we notice that $S(w)$ is an even function and the linearized problem can be solved in the even symmetric function class. Without loss of generality, we can pose the further restriction on the set $\mathcal{H}$ that all the perturbations $\phi$ are even symmetric functions. As a consequence, we see that apart for the term $\frac{dw}{dx}$ all the remaining terms are even symmetric and this implies that $c=0.$ Hence $w+\phi$ is a solution to the original Gierer-Meinhardt system (6.1). ∎

6.3. Stability Analysis: large and small eigenvalues

In this section we characterize the linear stability of the ground state solution constructed in §6.2 above by considering both large and small eigenvalues.

6.3.1. Large eigenvalue

Linearizing (6.1) about the equilibrium solution $(u,v)$ we obtain the following eigenvalue problem

[TABLE]

where $\lambda_{s}\in\mathbb{C}$, $\phi\in H^{1}({\mathbb{R}})$, and $\psi\in H^{2s}({\mathbb{R}})$. Let

[TABLE]

Then (6.44) can be rewritten as

[TABLE]

Our aim is to study the large eigenvalues, i.e. those for which we may assume that there exists $c>0$ such that $|\lambda_{s}|\geq c>0$ for $1-2s$ is small. If $\Re(\lambda_{s})<-c$ then we are done and we therefore may assume that $\Re(\lambda_{s})\geq-c$. For a subsequence $1-2s\to 0$ and $\lambda_{s}\to\lambda_{0}$ we shall derive a limiting NLEP satisfied by $\lambda_{0}$.

To simplify our argument, we shall assume $\tau=0$ and the general case can be proved by a perturbation argument. When $x\in I_{s}$, we calculate

[TABLE]

Substituting this into (6.45a), and letting $2s-1\rightarrow 0$, we derive the following nonlocal eigenvalue problem

[TABLE]

By Theorem 3.2 in [7] we see that $\lambda_{0}<0$, which implies that the large eigenvalues are stable.

6.3.2. Small eigenvalue

We next consider the small eigenvalues of (6.45), i.e. those for which $\lambda_{s}\to 0$ as $s\to\frac{1}{2}$. In last section, we have already shown the existence of solutions $(\hat{U},\hat{V})$ to (6.8). We notice that this equation is translation invariant. By differentiating (6.8) we derive that

[TABLE]

This suggests that $(\phi,\psi)$ of (6.45) can be written as

[TABLE]

where $\phi^{\perp}\perp\frac{d\hat{U}}{dx}$ and $\psi^{\perp}$ satisfy

[TABLE]

As $s\to\frac{1}{2}$, we know that

[TABLE]

Multiplying (6.50a) by $\phi^{\perp}$ we have

[TABLE]

From Lemma A.2 in [7] we have that

[TABLE]

where $a_{1}>0$ and

[TABLE]

Since $\phi^{\perp}\perp\frac{d\hat{U}}{dx}$ and $\hat{U}$ is well approximated by $w$, we get from (6.51) that

[TABLE]

Hence, we have shown all the small eigenvalues are stable. Thus, Theorem 1.2 follows by combining the conclusions of the last two sections.

7. Discussion

In this paper we have used formal asymptotic methods to study the existence and linear stability of localized solutions for the fractional Gierer-Meinhardt system where the fractional order of the inhibitor is $s_{2}\in(0,1/2)$. These results extend those previously obtained in [7] and [15] for $s_{2}\in(1/2,1)$ and $s_{2}=1/2$ respectively. Using the method of matched asymptotic expansions the construction of localized solutions was reduced to solving a system of nonlinear algebraic equations while the study of their linear stability was reduced to analyzing a globally coupled eigenvalue problem. We found that when $D=O(\varepsilon^{2s_{1}-1})$ both symmetric and asymmetric multi-spike solutions can be constructed though the latter were found to always be linearly unstable. On the other hand symmetric spikes were found to have stability regions outside of which they may undergo either a competition instability or a Hopf bifurcation. Using a leading order theory we found that the competition instability threshold is monotone decreasing in $s_{1}$ and it is either monotone decreasing in $s_{2}$ when $s_{1}>0.5$ or non-monotonic (first increasing and then decreasing) when $s_{1}<0.5$. In addition we found that the Hopf bifurcation threshold increases with $1/4<s_{1}<1$ provided $s_{2}$ and $\kappa$ are large enough, whereas it decreases with $0<s_{2}<1/2$ for all values of $s_{1}$ and $\kappa$. We also computed higher-order stability thresholds for specific cases of one- and two-spike solutions and these were supported by full numerical simulations of the system (1.2). Finally, in addition to the linear stability over an $O(1)$ timescale we also determined that spike solutions may be susceptible to drift instabilities leading to mutual repulsion between spikes, though these arise over a much slower $O(\varepsilon^{3-2s_{2}})$ timescale.

A key component in the formal construction of multi-spike solutions is the core problem (2.2) which was considered in detail numerically in §2.1 for general $s_{2}\in(0,1/2)$ and rigorously in §6 for $s_{2}\approx 1/2$. We found that the behaviour of the far-field constant $\mu(S)$ shares some properties with its counterpart in the three-dimensional Gierer-Meinhardt system previously studied in [6]. In particular we used numerical continuation to deduce the existence of a value $S=S_{\star}$ for which the core problem admits a ground state solution (i.e. one for which $\mu(S_{\star})=0$). The existence and linear stability of such a ground state was then rigorously established in §6 for $s_{2}\approx 1/2$.

Finally, throughout our paper we have highlighted the similarities between both the analysis and structure of localized solutions for the one-dimensional fractional Gierer-Meinhardt system when $s_{2}\in(0,1/2)$ and the corresponding localized solutions in the three-dimensional Gierer-Meinhardt system [6]. This connection is a result of the leading order algebraic singularity of the Green’s function which in particular fixes the far-field behaviour of solutions to the core problem (2.2) and also plays a key role in the asymptotic matching. In Appendix A we provide an expression for the Green’s function which makes explicit its singular behaviour, showing in particular that the singular behaviour consists of multiple algebraic singularities when $s_{2}\in(0,1/2)\setminus\{\tfrac{1}{2r}\,|\,r\in\mathbb{Z},r\geq 1\}$ (see Proposition A.1) as well as logarithmic singularities for $s_{2}=\tfrac{1}{2r}$ for $r\in\mathbb{Z}$ with $r\geq 1$ (see Proposition A.2) We believe that these expressions for the Green’s function will be particularly useful for future studies of localized solutions in one-dimensional fractional reaction-diffusion systems.

We conclude by highlighting some outstanding problems and suggestions for future research. One of the first outstanding problems is to derive a higher-order asymptotic theory in the case when $s_{2}=\tfrac{1}{2r}$ for $r=1,2,...$. The key hurdle in this direction is the emergence of both logarithmic and algebraic singularities in the Green’s functions and we believe that a resolution of this would spark some interesting mathematics. Additionally, it would be interesting to provide a rigorous justification for the existence and linear stability results we have formally derived for general values of $0<s_{2}<1/2$. Extensions of the current model to incorporate non-periodic boundary conditions as well as different reaction-kinetics would also be an interesting direction for future research. Moreover the consideration of such fractional problems in two- and three-dimensional domains will also lead to interesting mathematical questions.

Acknowledgement

D. Gomez is supported by NSERC and the Simons Foundation, M. Medeiros is partially supported by NSERC, J. Wei is partially supported by NSERC, and W. Yang is partially supported by NSFC No.11801550 and 1187147.

Conflicts of Interest

The authors don’t have any financial or non-financial conflicts of interest to disclose in relation to the contents of this paper.

Data Availability

The data generated during and/or analysed during the current study is available from the corresponding author upon a reasonable request.

Appendix A A Rapidly Converging Series for the Green’s Function

In this appendix we derive a rapidly converging series expansion of the periodic fractional Green’s function satisfying (2.14). We begin by formally computing the Fourier series

[TABLE]

In contrast to the classical Green’s function, the fractional Green’s function contains multiple singular terms. By identifying and removing these singular terms from the series expansion (A.1) we can therefore obtain a rapidly converging series expansion.

We will first restrict $s_{2}$ to be such that $s_{2}\neq\tfrac{1}{2r}$ for any integers $r\geq 1$, returning to the remaining cases at the end of this appendix. We first note that for any $l>1$ and sufficiently large values of $n>0$ such that $D(\pi n)^{2s_{2}}>1$ we have

[TABLE]

Substituting into (A.1) then gives

[TABLE]

where we choose $k_{\max}$ to be the smallest positive integer such that $2(1+k_{\max})s_{2}>1$, i.e. $k_{\max}=\lceil\tfrac{1}{2s_{2}}-1\rceil$. This choice of $k_{\max}$ guarantees the second term in (A.2) converges. We remark that other choices are also possible. For example, if second order derivatives at $x=0$ are needed then it will be more convenient to choose $k_{\max}$ to be the smallest integer such that $k_{\max}>\frac{3}{2s_{2}}-1$.

To determine the singular terms from the remaining sum in (A.2) we first let $\beta>0$ and consider the series expansion

[TABLE]

Repeated integration by parts yields the identities

[TABLE]

Letting $q\geq 1$ be the largest positive integer such that $\beta-2(q-1)>0$ the above identities give

[TABLE]

Note that if $\beta>0$ is an integer then (A.4) yields an explicit expression for $c_{\beta,n}$ Indeed if $\beta$ is an odd integer then $\beta-2(q-1)=1$ so that the product in the second term of (A.4) vanishes whereas if $\beta$ is an even integer then $\beta-2(q-1)=2$ and the integral in the second term of (A.4) evaluates to $1-(-1)^{n}$. In particular when $\beta=3$ we obtain the useful series

[TABLE]

If instead $\beta>0$ is not an integer then we first write

[TABLE]

Since in this case $-1<\beta-2q+1<1$ we can use standard properties of the Gamma function (see equations (5.9.7) and (5.5.3) in [2]) to write

[TABLE]

where $\mathfrak{a}_{\beta/2}=-\pi^{-1}\beta\Gamma(-\beta)\sin(\pi\beta/2)$. Moreover since $\beta-2q<0$ we can integrate by parts to get

[TABLE]

where we remark that the last term is $O(n^{\beta-2q-2})$. In summary, for non-integer values of $\beta>0$ we have the expression

[TABLE]

The constant term appearing in (A.6a) allows us to relate $|x|^{\beta-1}$ with the series (A.2) whereas the decay $a_{\beta,n}=O(n^{\beta-2q-2})$ from (A.6b) yields a quickly converging series. Specifically, since our choice of $k_{\max}=\lceil\tfrac{1}{2s_{2}}-1\rceil$ implies that $\beta\leq 1-2s_{2}$ we deduce that $\beta=2ks_{2}\leq 1-2s_{2}$ is not an integer for all $1\leq k\leq k_{\max}$. Therefore for any $0<s_{2}<1/2$ we can rewrite the summands in the rightmost term of (A.2) as

[TABLE]

Substituting this back into (A.2) we reach our final result.

Proposition A.1**.**

Let $s_{2}\in(0,\tfrac{1}{2})\setminus\{\tfrac{1}{2r}\,|\,r\in\mathbb{Z},\,r\geq 1\}$. Then the Green’s function satisfying (2.14) is given by

[TABLE]

We conclude this appendix by revisiting our restriction that $s_{2}\neq\tfrac{1}{2r}$ for all integers $r\geq 1$. If instead $s_{2}=\tfrac{1}{2r}$ for some integer $r\geq 1$, then $k_{\max}=r$ and we see that $\mathfrak{a}_{k_{\max}s_{2}}=-2k_{\max}s_{2}\pi\Gamma(-1)\sin(\pi/2)$ is undefined. This extends beyond a technical difficulty and is closely tied to the singular behaviour of the Green’s function (2.15). Indeed when $s_{2}>1/2$ the Green’s function has no singular terms, but if $1/4<s_{2}<1/2$ then it has one singular term, when $1/6<s_{2}<1/2$ it has two, and so on. At the transition points, i.e. $s_{2}=2^{-1},4^{-1},6^{-1}$, etc., the Green’s function has an additional logarithmic singularity. One way to see this formally is to observe that if $s_{2}\rightarrow\tfrac{1}{2r}$ then $|x|^{2rs_{2}-1}\sim(2rs_{2}-1)\log|x|$ whereas $\mathfrak{a}_{rs_{2}}\sim\pi^{-1}(1-2rs_{2})^{-1}$ so that the $k=k_{\max}=r$ term in (A.7) behaves like $(-1)^{r}D^{1-r}\pi^{-1}\log|x|$. To make this more precise we can calculate the Fourier series of $\log|x|$ and integrate by parts to get

[TABLE]

where $\text{Si}(z)\equiv\int_{0}^{z}\tfrac{\sin t}{t}dt$ is the Sine integral. Note that the rightmost term converges for all $-1<x<1$ since $\text{Si}(z)\sim\pi/2+O(z^{-1})$ as $z\rightarrow+\infty$. Comparing (A.8) with the rightmost term in (A.2) we readily deduce our next result.

Proposition A.2**.**

Let $s_{2}=\tfrac{1}{2r}$ for some integer $r\geq 1$. Then the Green’s function satisfying (2.14) is given by

[TABLE]

where $\text{Si}(z)$ is the sine integral and where $a_{k/r,n}$ and $\mathfrak{a}_{k/(2r)}$ are defined as in (A.7).

Appendix B Numerical Implementation

In this appendix we outline the numerical methods used for the numerous computations in this paper. Specifically we first describe the numerical implementation of the fractional Laplacian in both $\mathbb{R}$ with a prescribed far-field behaviour as well as in the interval $-1<x<1$ with periodic boundary conditions. We then outline the key steps in the numerical continuation used to construct solutions to the core problem (2.2). This is followed by a discussion of the numerical computation of the spectrum of the nonlocal operator $\mathscr{M}$ defined in (3.6b) as well as its adjoint $\mathscr{M}^{\star}$ from §4. We then outline the implementation of IMEX methods for numerically simulating (1.2).

B.1. Numerical Computation of The Fractional Laplacian

To numerically solve the core problem (2.2) as well as to determine the spectrum of the nonlocal operator $\mathscr{M}$ defined in (3.6b) and its adjoint $\mathscr{M}^{*}$ considered in §4 we must first calculate an appropriate discretization of the fractional Laplacian in $\mathbb{R}$. On the other hand to numerically simulate the full system (1.2) we must discretize the fractional Laplacian on $-1<x<1$ with periodic boundary conditions.

B.1.1. The Fractional Laplacian in $\mathbb{R}$

We use the finite-difference method of Huang and Oberman [10] to calculate the discretized fractional Laplacian. For completeness we include here the most important details in the implementation. To simplify our presentation we consider the problem of discretizing $(-\Delta)^{s}\varphi$ for $0<s<1/2$ when $\varphi$ satisfies the far-field behaviour $\varphi\sim a_{\pm}+b_{\pm}|y|^{-\beta}$ for constants $a_{\pm}$, $b_{\pm}$, and $\beta>0$. We first introduce the truncated domain $-L<y<L$ and consider the approximate boundary conditions $\varphi(y)\approx a_{\pm}+b_{\pm}|y|^{-\beta}$ for $\pm y\geq L$. Next we introduce the discretization $y_{i}=ih$ for $i=-2N,...,2N$ where $h=L/N$. This leads to the computational domain $-2L<y<2L$ which has been expanded from the original truncated domain to account for nonlocal contributions. For each $i>N$ we impose a fixed value for $\varphi_{\pm i}\equiv\varphi(y_{\pm i})$ that depends on its value at $i=\pm N$ in one of two ways

[TABLE]

These two cases account, respectively, for the activator and inhibitor in the core problem (2.2). We next let $\nu_{s}(y)\equiv C_{s}|y|^{-(1+2s)}$ and for each $|i|\leq N$ we decompose (1.1) as

[TABLE]

The first integral in (B.2) is approximated by performing a piecewise quadratic interpolation of $\varphi$ which gives us

[TABLE]

where $w_{j}^{Q}$ ($j=-2N,...,2N$) are the quadratic interpolation weights calculated using Definition 3.2 of [10] with $\alpha=2s$ and explicitly given by

[TABLE]

Note that some of the entries in (B.3) will have $|i-j|>N$ and for these we use (B.1).

Next, by using the definition of $\nu_{s}(y)$ we readily calculate

[TABLE]

To calculate IIIi we first consider the portion of the integral for which $y>2L$ and hence $y_{i}-y<-L$ so that by using (B.1) we calculate

[TABLE]

The rightmost integral in both cases can be written in terms of the Gauss Hypergeometric function ${}_{2}F_{1}(a,b,c,z)$ which we recall has the integral representation (see (9.111) in [8])

[TABLE]

A simple change of variables then immediately yields that for any $\beta$ and $\alpha$ with $\beta+\alpha>0$

[TABLE]

with an analogous result for the integration over $-\infty<y<-2L$. Thus for Case 1 and Case 2 we respectively obtain

[TABLE]

Substituting equations (B.3), (B.5), and (B.7) into (B.2) thus yields a discretization of the fractional Laplacian incorporating the relevant far-field behaviour for the core problem and spectrum calculations. We conclude by remarking that in Case 1 this discretization leads to multiplication by a dense matrix whereas for Case 2 the discretization involves both multiplication by a dense matrix as well as the addition of an inhomogeneous term that captures the far-field behaviour.

B.1.2. The Periodic Fractional Laplacian in $-1<x<1$

We now consider the discretization of $(-\Delta)^{s}\varphi(x)$ in $-1<x<1$ with periodic boundary conditions. In contrast to the discretization of the fractional Laplacian in $\mathbb{R}$ considered above, periodic boundary conditions lead to a significantly simpler implementation. We discretize the domain $-1<x<1$ by letting $x_{i}=-1+2ih$ for $i=0,...,N-1$ where $h=1/N$. By approximating $\varphi$ with a piecewise quadratic interpolator and letting $\varphi_{i}\equiv\varphi(x_{i})$ we then get

[TABLE]

where $w_{i}^{Q}$ are again the quadratic weight functions found in [10]. This discretization is readily implemented and leads to a dense matrix. To compute $W_{n}^{Q}$ we truncate the sum by taking the first $10^{3}$ terms.

B.2. Solving the Core Problem

In this section we outline the key steps for calculating $\mu(S)$ and $\nu(S)$ for $S>0$. This computation has two main steps. In the first step we numerically compute the fractional homoclinic $w_{s_{1}}$ satisfying (2.9) for values of $1/4<s_{1}<1$. In the second step we use the $S\ll 1$ asymptotics (2.8a) to initiate a numerical continuation in $S$ to solve (2.2) from which $\mu(S)$ and $\nu(S)$ can then be computed from (2.7).

The fractional homoclinic solution to (2.9) was previously numerically computed for values of $1/4<s_{1}<1$ in Appendix B of [7] and therefore we provide only an outline here and refer the reader to that paper for more details. The first step is to discretize (2.9) using the method in Appendix B.1.1. This yields a nonlinear system for the discretized solution which we can solve using Newton’s method. Next we note that the exact solution to (2.9) when $s_{1}=1/2$ is $w_{1/2}(y)=2/(1+y^{2})$. Letting $(s_{1})_{0}=1/2<(s_{1})_{1}<...<(s_{1})_{M}<1$ for a large $M>0$ we can then numerically solve the discretized nonlinear system for $s_{1}=(s_{1})_{i}$ when $i=1,\cdots,M$ by using the solution for $s_{1}=(s_{1})_{i-1}$ as an initial guess where the exact solution is used for $i=1$. A similar continuation is likewise performed for values of $1/4<s_{1}<1/2$.

The numerical solution of the core problem (2.2) follows a similar procedure. Numerically discretizing both the fractional Laplacians appearing in (2.2) using the method in (B.1.1) yields a nonlinear system for the discretized solutions. Note that when performing the discretization the far-field condition for Case 1 and Case 2 (see Equation (B.1)) is used for $U_{c}$ and $V_{c}$ respectively. It is then straightforward to perform a numerical continuation in $S>0$ by slowly incrementing it and using the small $S$ asymptotics (2.8a) as the initializing guess.

B.3. Computing the Spectrum of $\mathscr{M}$

In this section we consider the numerical computation of the spectrum of the nonlocal operator $\mathscr{M}$ defined in (3.6b). This is done by discretizing $\mathscr{M}$ and then calculating the spectrum of the resulting matrix operator using standard eigenvalue libraries (we used the eig function from the SciPy linalg library). In this section it therefore suffices to describe the discretization of $\mathscr{M}$. The operator $\mathscr{M}$ consists of two nonlocal contributions: the fractional Laplacian $(-\Delta)^{s_{1}}$ and the convolution with the fractional Green’s function. The numerical discretization of the former was considered in detail in Appendix B.1.1 and so it remains only to discuss the discretization of the latter.

To simplify (and generalize) our presentation we henceforth focus on numerically approximating

[TABLE]

where $0<s<1/2$ and $\varphi$ is assumed to have the far-field behaviour

[TABLE]

where $\varphi_{\pm\infty}$ are unknown. Note that we use $s=s_{2}$ in the discretization of both $\mathscr{M}$ and $\mathscr{M}^{\star}$ whereas we use $\varphi(y)=U_{c}(y;S)\Phi_{c}^{\lambda}(y;S)$ so that $\beta=2+2s_{1}$ in the former and $\varphi(y)=(U_{c}(y;S)/V_{c}(y;S))^{2}P(y)$ so that $\beta=3+3s_{1}$ in the latter. We fix our notation by introducing the discretization $y_{i}=ih$ and letting $\Phi_{i}=\varphi(y_{i})$ for $i=-2N,...,2N$ where $h=L/N$. Then we decompose

[TABLE]

The first integral is easily approximated using a Taylor series as

[TABLE]

The approximation of the remaining two integrals proceeds as for the fractional Laplacian above. Specifically imposing that $\varphi_{j}=\varphi_{\pm N}|y_{j}/L|^{-\beta}$ for all $\pm j>N$ we first approximate $J_{2}$ with the finite sum

[TABLE]

where $\tilde{w}_{j}^{Q}$ ($j=-2N,...,2N$) are the quadratic interpolation weights given by

[TABLE]

where $G(t)$ is given by (B.4b). Observe that the multiplicative factor $C_{s}$ as well as the $(2-2s)^{-1}$ term for the $j=\pm 1$ case present in (B.4) are omitted in (B.14). The former is omitted because it is specific to the definition of the fractional Laplacian while the latter is omitted because it arises from a finite difference approximation for the singular part of the fractional Laplacian which in the present case is contained in the $J_{1}$ contribution.

The remaining integral contributions $J_{3}$ are then computed using (B.6) which yields

[TABLE]

B.4. Time Stepping of the Fractional Gierer-Meinhardt System

We performed full numerical simulations of (1.2) by first discretizing the system using the methods in §B.1.2 which yields a system of $2N$ ordinary differential equations of the form

[TABLE]

In this expression $\bm{\Phi}(t)=(u(x_{0},t),\cdots,u(x_{N-1},t),v(x_{0},t),\cdots,v(x_{N-1},t))^{T}$ approximates the solution at the discretization points $x_{i}=-1+2ih$ for $i=0,...,N-1$ as in §B.1.2. The $2N\times 2N$ block-diagonal matrix $\mathcal{A}$ then has entries corresponding to the discretization (B.8) while the $2N$ dimensional vector $\bm{\mathcal{N}}(\bm{\Phi})$ accounts for the nonlinearities in (1.2). We integrate the ODE system using a second-order semi-implicit backwards difference scheme (2-SBDF) [21] which leads to the linear system

[TABLE]

where $\Delta t>0$ is the time-step size and $\bm{\Phi}_{n}=\bm{\Phi}(n\Delta t)$. The initial condition is given by $\bm{\Phi}_{0}$ and since 2-SBDF is second-order we also need $\bm{\Phi}_{1}$. We obtain this second initial value by integrating (B.16) using the first-order semi-implicit backward difference scheme (1-SBDF) with a smaller time step $\Delta t/\tilde{N}$ for $\tilde{N}>0$ steps. Specifically, we do this by solving

[TABLE]

for $n=1,...,\tilde{N}$. For our numerical simulations we used $N=2000$, $\tilde{N}=5$, and $\Delta t=0.01$. Moreover in the spatial discretization of the fractional Laplacian we truncated the infinite sum appearing in (B.8) after $250$ terms.