# Spike Solutions to the Supercritical Fractional Gierer-Meinhardt System

Daniel Gomez  Daniel Gomez, Department of Mathematics, University of
Pennsylvania, Philadelphia, PA USA 19104-6395<EMAIL_ADDRESS>,  Markus
De Medeiros  Markus Medeiros, Department of Mathematics, The University of
British Columbia, Vancouver, BC Canada V6T 1Z2<EMAIL_ADDRESS>,  Jun-cheng Wei  Jun-cheng Wei, Department of Mathematics, The University of
British Columbia, Vancouver, BC Canada V6T 1Z2<EMAIL_ADDRESS>and  Wen
Yang  Wen Yang, Wuhan Institute of Physics and Mathematics, Innovation Academy
for Precision Measurement Science and Technology, Chinese Academy of Sciences,
P.O. Box 71010, Wuhan 430071, P. R. China<EMAIL_ADDRESS>
###### 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évy 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.

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

$(-\Delta)^{s}\varphi(x)\equiv
C_{s}\int_{-\infty}^{\infty}\frac{\varphi(x)-\varphi(\bar{x})}{|x-\bar{x}|^{1+2s}}d\bar{x},\qquad
C_{s}\equiv\frac{2^{2s}s\Gamma(s+2^{-1})}{\sqrt{\pi}\Gamma(1-s)}.$ (1.1)

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

$\displaystyle[left=\empheqlbrace]$ $\displaystyle
u_{t}+\varepsilon^{2s_{1}}(-\Delta)^{s_{1}}u+u-v^{-1}u^{2}=0,$
$\displaystyle-1<x<1,$ (1.2a) $\displaystyle\tau
v_{t}+D(-\Delta)^{s_{2}}v+v-u^{2}=0,$ $\displaystyle-1<x<1,$ (1.2b)

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,

$u(x+2,t)=u(x,t),\qquad v(x+2,t)=v(x,t),$ (1.3)

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

$\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}$
(1.4)

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

$(-\Delta)^{s_{2}}G+D^{-1}G_{D}=\delta(x),\qquad 0<x<1,\qquad
G_{D}(x+2)=G_{D}(x),$

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.

(a)

(b)

Figure 1. (A) Numerically calculated profile of the activator in a symmetric
two-spike solution undergoing (A) a competition instability using parameters
$s_{1}=0.5$, $s_{2}=0.39$, $\varepsilon=0.01$, $\tau=0.1$, and
$D=1.095\varepsilon^{2s_{2}-1}$, and (B) slow dynamics over an
$O(\varepsilon^{2s_{2}-3})$ timescale using parameters $s_{1}=0.4$,
$s_{2}=0.35$, $\varepsilon=0.01$, $\tau=0.1$, and
$D=0.768\varepsilon^{2s_{2}-1}$. The solid blue and dashed orange curves in
the $(x,t)$ plane in (B) indicate the numerically and asymptotically
calculated spike locations.

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

$1/4<s_{1}<1,\qquad 0<s_{2}<1/2.$ (1.5)

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

$\displaystyle[left=\empheqlbrace]$
$\displaystyle(-\Delta)^{\frac{1}{2}}U+U-V^{-1}U^{2}=0,\quad(-\Delta)^{s}V-U^{2}=0,$
$\displaystyle-\infty<x<\infty,$ (1.6a) $\displaystyle U,V>0,$
$\displaystyle-\infty<x<\infty,$ (1.6b) $\displaystyle U,V\to 0,$
$\displaystyle\mbox{as}\quad|x|\to+\infty,$ (1.6c)

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

$\lim_{\varepsilon_{0}\to
0}\left|\tau_{s}^{-1}U(x)-w(x)\right|=0,\qquad\lim_{\varepsilon_{0}\rightarrow
0}\left|\tau_{s}^{-1}V(x)-1\right|=0,$ (1.7)

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

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

and

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

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

$\displaystyle[left=\empheqlbrace]$
$\displaystyle(-\Delta)^{\frac{1}{2}}\phi+\phi-2V^{-1}U\phi+V^{-2}U^{2}\psi+\lambda_{s}\phi=0,$
$\displaystyle-\infty<x<\infty,$ (1.8a)
$\displaystyle(-\Delta)^{s}\psi-2U\phi+\tau\lambda_{s}\psi=0,$
$\displaystyle-\infty<x<\infty,$ (1.8b)

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

$\displaystyle[left=\empheqlbrace]$
$\displaystyle\varepsilon^{2s_{1}}(-\Delta)^{s_{1}}u+u-v^{-1}u^{2}=0,$
$\displaystyle-1<x<1,$ (2.1a) $\displaystyle D(-\Delta)^{s_{2}}v+v-u^{2}=0,$
$\displaystyle-1<x<1,$ (2.1b) with periodic boundary conditions
$u(x+2)=u(x),\qquad v(x+2)=v(x).$ (2.1c)

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

$\displaystyle[left=\empheqlbrace]$
$\displaystyle(-\Delta)^{s_{1}}U_{c}+U_{c}-V_{c}^{-1}U_{c}^{2}=0,\quad(-\Delta)^{s_{2}}V_{c}-U_{c}^{2}=0,$
$\displaystyle-\infty<y<\infty,$ (2.2a) $\displaystyle
U_{c}\sim\nu(S)|y|^{-(1+2s_{1})},\quad V_{c}\sim\mu(S)+S|y|^{2s_{2}-1},$
$\displaystyle\text{as}\quad|y|\rightarrow\infty.$ (2.2b)

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

$(-\Delta)^{s}\phi+\phi=f,\quad\text{for
}-\infty<y<\infty;\qquad\phi\rightarrow 0,\quad\text{as
}|y|\rightarrow\infty,$

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

  2. (2)

If $\sigma>1$ then the solution to

$(-\Delta)^{s}\phi=f,\quad\text{for }-\infty<y<\infty;\qquad\phi\rightarrow
0,\quad\text{as }|y|\rightarrow\infty,$

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

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

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

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

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

$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},\qquad\mathfrak{a}_{s_{2}}\equiv-2\pi^{-1}s\Gamma(-2s_{2})\sin(\pi
s_{2}).$ (2.5)

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

$S=\mathfrak{a}_{s_{2}}\int_{-\infty}^{\infty}U_{c}(\bar{y};S)^{2}d\bar{y}.$
(2.6)

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

(a)

(b)

(c)

(d)

(e)

(f)

Figure 2. Plots of core problem far-field constants $\mu(S)$ and $\nu(S)$ for
distinct values of $1/4<s_{1}<1$. In each plot the darkest and lightest curves
corresponds to $s_{2}=0.2$ and $s_{2}=0.49$ respectively, with the
intermediate curves corresponding to $0.01$ increments in $s_{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

$\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{)}.$
(2.7)

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

$U_{c}(y)\sim\sqrt{\tfrac{S}{\mathfrak{b}_{s_{1}}\mathfrak{a}_{s_{2}}}}w_{s_{1}}(y),\quad
V_{c}(y)\sim\sqrt{\tfrac{S}{\mathfrak{b}_{s_{1}}\mathfrak{a}_{s_{2}}}},\qquad\mu(S)\sim\sqrt{\tfrac{S}{\mathfrak{b}_{s_{1}}\mathfrak{a}_{s_{2}}}},$
(2.8a) where
$\mathfrak{b}_{s_{1}}\equiv\int_{-\infty}^{\infty}w_{s_{1}}(y)^{2}dy,$ (2.8b)

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

$\displaystyle[left=\empheqlbrace]$
$\displaystyle(-\Delta)^{s_{1}}w_{s_{1}}+w_{s_{1}}-w_{s_{1}}^{2}=0,$
$\displaystyle-\infty<y<\infty,$ (2.9a) $\displaystyle
w_{s_{1}}(y)=O(|y|^{-(1+2s_{1})}),$
$\displaystyle\text{as}\quad|y|\rightarrow\infty.$ (2.9b)

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

(a)

(b)

(c)

Figure 3. Plots of the critical values (A) $S=S_{\star}$ and (B)
$S=S_{\text{crit}}$ at which the far-field constant $\mu(S)$ vanishes and
attains its global maximum respectively. In both plots the darkest and
lightest curves corresponds to values of $s_{1}=0.3$ and $s_{1}=0.7$
respectively with the intermediate curves corresponding intermediate values in
increments of $0.05$. (C) The core solution for $s_{1}=0.5$ and $s_{2}=0.4$ at
the indicated values of $S$.

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

$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{)}.$ (2.10)

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

$(-\Delta)^{s_{1}}U_{i}+U_{i}-V_{i}^{-1}U_{i}^{2}=0,\quad(-\Delta)^{s_{2}}V_{i}+D^{-1}\varepsilon^{2s_{2}}V_{i}-U_{i}^{2}=0\quad-(1+x_{i})<\varepsilon
y<1-x_{i}.$ (2.11)

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

$U_{i}(y)\sim U_{c}(y;S_{i})+o(1),\qquad V_{i}(y)\sim V_{c}(y;S_{i})+o(1),$
(2.12)

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

$u^{2}\rightarrow\varepsilon^{1-4s_{2}}D^{2}\sum_{i=1}^{N}\int_{-\infty}^{\infty}U_{c}(y;S_{i})^{2}dy\,\delta(x-x_{i})=\varepsilon^{1-4s_{2}}D^{2}\mathfrak{a}_{s_{2}}^{-1}\sum_{i=1}^{N}S_{i}\delta(x-x_{i}),$

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

$(-\Delta)^{s_{2}}v+D^{-1}v=\varepsilon^{1-4s_{2}}\mathfrak{a}_{s_{2}}^{-1}D\sum_{i=1}^{N}S_{i}\delta(x-x_{i}),\qquad
x\in(-1,1)\setminus\\{x_{1},...,x_{N}\\},$ (2.13a) and having the limiting
behaviour $v\sim
D\varepsilon^{-2s_{2}}\bigl{(}\varepsilon^{1-2s_{2}}S_{i}|x-x_{i}|^{2s_{2}-1}+\mu(S_{i})\bigr{)},\qquad\text{as}\quad
x\rightarrow x_{i},$ (2.13b) for each $i=1,...,N$ obtained from the far-field
behaviour of the inner solution (2.2b).

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

$(-\Delta)^{s_{2}}G_{D}+D^{-1}G_{D}=\delta(x),\quad-1<x<1,\qquad
G_{D}(x+2)=G_{D}(x),$ (2.14)

which can be written as (see Appendix A)

$G_{D}(x)=D\sum_{k=1}^{k_{\max}}\frac{(-1)^{k+1}\mathfrak{a}_{ks_{2}}}{D^{k}}|x|^{2ks_{2}-1}+R_{D}(x),\qquad
k_{\max}\equiv\lceil\tfrac{1}{2s_{2}}-1\rceil,$ (2.15)

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

$v(x)=\varepsilon^{1-4s_{2}}\mathfrak{a}_{s_{2}}^{-1}D\sum_{i=1}^{N}S_{i}G_{D}(x-x_{i}).$
(2.16)

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

$\begin{split}\mu(S_{i})+\varepsilon^{1-2s_{2}}S_{i}|x-x_{i}|^{2s_{2}-1}\sim\frac{\varepsilon^{1-2s_{2}}}{\mathfrak{a}_{s_{2}}}\biggl{(}DS_{i}\sum_{k=1}^{k_{\max}}\frac{(-1)^{k+1}\mathfrak{a}_{ks_{2}}}{D^{k}}|x-x_{i}|^{2ks_{2}-1}+S_{i}R_{D}(0)&\\\
+\sum_{j\neq i}S_{j}G_{D}(|x_{i}-x_{j}|)+O(|x-x_{i}|)&\biggr{)}.\end{split}$
(2.17)

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

$\mu(S_{i})=\frac{\varepsilon^{1-2s_{2}}}{\mathfrak{a}_{s_{2}}}\biggl{(}S_{i}R_{D}(0)+\sum_{j\neq
i}S_{j}G_{D}(|x_{i}-x_{j}|)\biggr{)},\qquad i=1,...,N.$ (2.18)

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

$R_{D}(0)\sim\frac{1}{2}D+O(1),\qquad
G_{D}(|x_{i}-x_{j}|)\sim\frac{1}{2}D+O(1)\quad\text{for }i\neq j\qquad(D\gg
1),$ (2.19)

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

$\mu(S_{i})=\frac{\kappa}{N}\sum_{j=1}^{N}S_{j},\qquad\kappa\equiv\frac{ND_{0}}{2\mathfrak{a}_{s_{2}}},$
(2.20)

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

(a)

(b)

(c)

Figure 4. Plots of (a) the common spike strength in a symmetric solution, (b)
the small spike strength value corresponding to the large spike strength value
in an asymmetric solution, and (c) its derivative. In each plot the darkest
and lightest curves corresponds to $s_{2}=0.2$ and $s_{2}=0.49$ respectively,
with the intermediate curves corresponding to $0.01$ increments in $s_{2}$

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

$\mu(S_{c})=\kappa S_{c},\qquad 0<S_{c}<S_{\star}.$ (2.21)

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

$S_{\sigma(1)}=...=S_{\sigma(n)}=S_{r},\qquad
S_{\sigma(n+1)}=...=S_{\sigma(N)}=S_{l}(S_{r}),\qquad
S_{\text{crit}}<S_{r}<S_{\star},$

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

$\mu(S_{r})=\kappa f(S_{r},\tfrac{n}{N}),\qquad f(S,\theta)\equiv\theta
S+(1-\theta)S_{l}(S).$ (2.22)

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

$0<\kappa<\kappa_{c1}\equiv\mu(S_{crit})/S_{\text{crit}},$ (2.23)

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

$\kappa_{c1}<\kappa<\kappa_{c2}\equiv\mu(S_{r}^{\star})/f(S_{r}^{\star},n/N),$
(2.24)

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

$f(S_{r}^{\star},n/N)\mu^{\prime}(S_{r}^{\star})=f^{\prime}(S_{r}^{\star},n/N)\mu(S_{r}^{\star}),$
(2.25)

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

$G_{D}(x)\sim\begin{cases}\sum_{k=1}^{k_{\max}}\tfrac{(-1)^{k-1}\mathfrak{a}_{ks_{2}}}{D^{k-1}}|x|^{2ks_{2}-1}+O(1),&s_{2}\in(0,1/2)\setminus\\{\tfrac{1}{2r}\,|\,r\in\mathbb{Z},\,r\geq
1\\},\\\
\sum_{k=1}^{r-1}\frac{(-1)^{k-1}\mathfrak{a}_{ks_{2}}}{D^{k-1}}|x|^{2ks_{2}-1}+\frac{(-1)^{r}}{\pi
D^{r-1}}\log|x|+O(1),&s_{2}=\tfrac{1}{2r},\,\text{for }r\in\mathbb{Z},\,r\geq
1,\end{cases}$ (2.26)

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.

Figure 5. Plots of the threshold $D_{R}(s_{2})$ (solid black) and its limiting
behaviour as $s_{2}\rightarrow\bigl{(}\tfrac{1}{2r}\bigr{)}^{-}$ for odd
values of $r\geq 1$ and as $s_{2}\rightarrow\bigl{(}\tfrac{1}{2r}\bigr{)}^{+}$
for even values of $r\geq 1$ (dashed red).

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

$D_{R}(s_{2})\sim\begin{cases}\bigl{(}\frac{2}{\pi(1-2rs_{2})}\bigr{)}^{1/r},&\text{as}\quad
s_{2}\rightarrow\bigl{(}\tfrac{1}{2r}\bigr{)}^{-}\quad\text{ for odd $r$},\\\
\bigl{(}\frac{2}{\pi(2rs_{2}-1)}\bigr{)}^{1/r},&\text{as}\quad
s_{2}\rightarrow\bigl{(}\tfrac{1}{2r}\bigr{)}^{+}\quad\text{for even
$r$}.\end{cases}$ (2.27)

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

$\displaystyle[left=\empheqlbrace]$
$\displaystyle\lambda\phi+\varepsilon^{2s_{1}}(-\Delta)^{s_{1}}\phi+\phi-2v_{e}^{-1}u_{e}\phi+v_{e}^{-2}u_{e}^{2}\psi=0,$
$\displaystyle-1<x<1,$ (3.1a)
$\displaystyle\tau\lambda\psi+D(-\Delta)^{s_{2}}\psi+\psi-2u_{e}\phi=0,$
$\displaystyle-1<x<1,$ (3.1b)
$\displaystyle\phi(x+2)=\phi(x),\qquad\psi(x+2)=\psi(x),$
$\displaystyle-1<x<1,$ (3.1c)

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

$\phi(x_{i}+\varepsilon
y)=\Phi_{i}^{\varepsilon}(y),\quad\psi(x_{i}+\varepsilon
y)=\Psi_{i}^{\varepsilon}(y),$

into (3.1) to get

$\displaystyle[left=\empheqlbrace]$
$\displaystyle\lambda\Phi_{i}+(-\Delta)^{s_{1}}\Phi_{i}+\Phi_{i}-2V_{i}^{-1}U_{i}\Phi_{i}+V_{i}^{-2}U_{i}^{2}\Psi_{i}=0,$
$\displaystyle-1+x_{i}<\varepsilon y<1-x_{i}$ (3.2a)
$\displaystyle\tau\lambda\varepsilon^{2s_{2}}D^{-1}\Psi_{i}+(-\Delta)^{s_{2}}\Psi+\varepsilon^{2s_{2}}D^{-1}\Psi_{i}-2U_{i}\Phi_{i}=0,$
$\displaystyle-1+x_{i}<\varepsilon y<1-x_{i}.$ (3.2b)

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

$\Phi_{i}^{\varepsilon}\sim
c_{i}\Phi_{c}^{\lambda}(y;S_{i})+o(1),\qquad\Psi_{i}^{\varepsilon}\sim
c_{i}\Psi_{c}^{\lambda}(y;S_{i})+o(1),$

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

$\displaystyle[left=\empheqlbrace]$
$\displaystyle(-\Delta)^{s_{1}}\Phi_{c}^{\lambda}+\Phi_{c}^{\lambda}-2V_{c}^{-1}U_{c}\Phi_{c}^{\lambda}+V_{c}^{-2}U_{c}^{2}\Psi_{c}^{\lambda}=-\lambda\Phi_{c}^{\lambda},$
$\displaystyle-\infty<y<\infty,$ (3.3a)
$\displaystyle(-\Delta)^{s_{2}}\Psi_{c}^{\lambda}-2U_{c}\Phi_{c}^{\lambda}=0,$
$\displaystyle-\infty<y<\infty,$ (3.3b)

where we assume the general far-field behaviour

$\Phi_{c}^{\lambda}\rightarrow 0,\quad\Psi_{c}^{\lambda}\sim
B(\lambda,S)+o(1),\qquad\text{as}\quad|y|\rightarrow\infty.$ (3.4)

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

$\Psi_{c}^{\lambda}(y;S)=B(\lambda,S)+2\mathfrak{a}_{s_{2}}\int_{-\infty}^{\infty}\frac{U_{c}(z;S)\Phi_{c}^{\lambda}(z;S)}{|y-z|^{1-2s_{2}}}dz.$
(3.5)

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

$\mathscr{M}\Phi_{c}^{\lambda}=\lambda\Phi_{c}^{\lambda}+B(\lambda,S)V_{c}^{-2}U_{c}^{2}$
(3.6a) where the nonlocal operator $\mathscr{M}=\mathscr{M}(S)$ is defined by
$\mathscr{M}\Phi\equiv-(-\Delta)^{s_{1}}\Phi-\Phi_{c}^{\lambda}+2V_{c}^{-1}U_{c}\Phi-2\mathfrak{a}_{s_{2}}V_{c}^{-2}U_{c}^{2}\int_{-\infty}^{\infty}\frac{U_{c}(z)\Phi(z)}{|y-z|^{1-2s_{2}}}dz.$
(3.6b)

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

$\Phi_{c}^{\lambda}(y,S)=B(\lambda,S)(\mathscr{M}-\lambda)^{-1}\bigl{(}V_{c}(y;S)^{-1}U_{c}(y;S)\bigr{)}^{2}.$
(3.7)

In addition we make note of the far-field behaviour

$\Psi_{c}^{\lambda}\sim
B(\lambda,S)+2\mathfrak{a}_{s_{2}}|y|^{2s_{2}-1}\int_{-\infty}^{\infty}U_{c}(z;S)\Phi_{c}^{\lambda}(z;S)dz\qquad\text{as}\quad|y|\rightarrow\infty,$

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

$\int_{-\infty}^{\infty}U_{c}(z;S)\Phi_{c}^{\lambda}(z;S)dz=\frac{1}{2\mathfrak{a}_{s_{2}}},$
(3.8)

with which we get the far-field behaviour

$\Psi_{c}^{\lambda}\sim
B(\lambda,S)+|y|^{2s_{2}-1}\qquad\text{as}\quad|y|\rightarrow\infty.$ (3.9)

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

$B(\lambda,S)=\biggl{(}2\mathfrak{a}_{s_{2}}\int_{-\infty}^{\infty}U_{c}(z;S)(\mathscr{M}-\lambda)^{-1}\bigl{(}V_{c}(z;S)^{-1}U_{c}(z;S)\bigr{)}^{2}dz\biggr{)}^{-1}.$
(3.10)

We next consider the distributional limit

$2u_{e}\phi\rightarrow
2\varepsilon^{1-2s_{2}}D\sum_{i=1}^{N}c_{i}\int_{-\infty}^{\infty}U_{c}(y;S_{i})\Phi_{c}^{\lambda}(y;S_{i})dy\delta(x-x_{i})$

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

$(-\Delta)^{s_{2}}\psi+\frac{1+\tau\lambda}{D}\psi=\varepsilon^{1-2s_{2}}\mathfrak{a}_{s_{2}}^{-1}\sum_{i=1}^{N}c_{i}\delta(x-x_{i}),\qquad
x\in(-1,1)\setminus\\{x_{1},...,x_{N}\\},$ (3.11a) together with the singular
behaviour $\psi(x)\sim
c_{i}\bigl{(}B(\lambda,S_{i})+\varepsilon^{1-2s_{2}}|x-x_{i}|^{2s_{2}-1}\bigr{)},\quad
x\rightarrow x_{i},$ (3.11b)

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

$\psi(x)=\frac{\varepsilon^{1-2s_{2}}}{\mathfrak{a}_{s_{2}}}\sum_{i=1}^{N}c_{i}G_{D_{\lambda}}(x-x_{i}),\qquad
D_{\lambda}\equiv\frac{D}{1+\tau\lambda}.$ (3.12)

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

$\begin{split}c_{i}\bigl{(}B(\lambda,S_{i})+\varepsilon^{1-2s_{2}}|x-x_{i}|^{2s_{2}-1}\bigr{)}\sim\frac{\varepsilon^{1-2s_{2}}}{\mathfrak{a}_{s_{2}}}\biggl{(}c_{i}\sum_{k=1}^{k_{\max}}\frac{(-1)^{k-1}\mathfrak{a}_{ks_{2}}}{D_{\lambda}^{k-1}}|x-x_{i}|^{2ks_{2}-1}+c_{i}R_{D_{\lambda}}&(0)\\\
+\sum_{j\neq
i}c_{j}G_{D_{\lambda}}(|x_{i}-x_{j}|)+O(|x-x_{i}|)&\biggr{)},\end{split}$

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

$\mathcal{B}(\lambda,\boldsymbol{S})\boldsymbol{c}=\varepsilon^{1-2s_{2}}\mathfrak{a}_{s_{2}}^{-1}\mathcal{G}_{D_{\lambda}}\boldsymbol{c},$
(3.13a) where $\boldsymbol{S}=(S_{1},...,S_{N})^{T}$,
$\boldsymbol{c}=(c_{1},...,c_{N})^{T}$, and
$\mathcal{B}(\lambda,\boldsymbol{S})$ and $\mathcal{G}_{D_{\lambda}}$ are
$N\times N$ matrices with entries
$(\mathcal{B}(\lambda,\boldsymbol{S}))_{ij}=\begin{cases}B(\lambda,S_{i}),&i=j,\\\
0,&i\neq
j,\end{cases}\qquad(\mathcal{G}_{D_{\lambda}})_{ij}=\begin{cases}R_{D_{\lambda}}(0),&i=j,\\\
G_{D_{\lambda}}(|x_{i}-x_{j}|),&i\neq j.\end{cases}$ (3.13b)

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.

(a)

(b)

(c)

(d)

Figure 6. (A)-(C) Plots of the real part of the eigenvalues of $\mathscr{M}$
versus $0<S<S_{\star}$ at the indicated values of $s_{1}$ and $s_{2}$. In each
plot the dashed vertical line corresponds to the value of $S=S_{\text{crit}}$
at which $\mu^{\prime}(S_{\text{crit}})=0$. (D) Plots of $\Re\lambda_{1}$ at
$S=S_{\star}$ versus $0.2<s_{2}<0.5$ for different values of $s_{1}$. The
darkest (uppermost) and lightest (lowermost) curves correspond to values of
$s_{1}=0.4$ and $s_{1}=0.49$ respectively, with the intermediate curves being
separated by intervals of $0.01$.

###  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$

$\mathcal{B}(\lambda,\boldsymbol{S})\boldsymbol{c}=\frac{\kappa}{1+\tau\lambda}\mathcal{E}_{N}\boldsymbol{c},\qquad\mathcal{E}_{N}=\frac{1}{N}\boldsymbol{e}\boldsymbol{e}^{T},$
(3.14)

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.

(a)

(b)

(c)

(d)

Figure 7. (A) Leading order competition instability threshold $\kappa_{c1}$
versus $s_{2}$. The darkest (uppermost) and lightest (lowermost) curves
correspond to values of $s_{1}=0.3$ and $s_{1}=0.7$ respectively, with the
intermediate curves corresponding to increments of $0.05$. The dashed line
corresponds to $s_{1}=0.5$. (B) The leading order Hopf bifurcation threshold
$\tau_{h}$ at $s_{1}=0.5$. The darkest and lightest curves correspond to
$s_{2}=0.3$ and $s_{2}=0.48$ respectively, with the intermediate curves
corresponding to $0.02$ increments in $s_{2}$. (C-D) The leading order Hopf
bifurcation threshold $\tau_{h}$ at indicated values of $s_{2}=0.45$ and
$s_{2}=0.40$. The darkest and lightest curves in both plots correspond to
$s_{1}=0.3$ and $s_{1}=0.7$ respectively, with the intermediate curves
corresponding to $0.05$ increments in $s_{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

$\mathscr{M}\Phi\sim\mathscr{L}\Phi+O(\kappa^{-1}),\qquad\mathscr{L}\Phi\equiv-(-\Delta)^{s_{1}}\Phi-\Phi+2w_{s_{1}}\Phi.$
(3.15)

From (3.10) we obtain

$B\bigl{(}\lambda,(\mathfrak{b}_{s_{1}}\mathfrak{a}_{s_{2}}\kappa^{2})^{-1}\bigr{)}\sim\frac{\kappa\int_{-\infty}^{\infty}w_{s_{1}}(y)^{2}dy}{2\int_{-\infty}^{\infty}w_{s_{1}}(y)\bigl{(}\mathscr{L}-\lambda\bigr{)}^{-1}w_{s_{1}}(y)^{2}dy},$

with which (3.14) becomes

$\frac{\int_{-\infty}^{\infty}w_{s_{1}}(y)^{2}dy}{2\int_{-\infty}^{\infty}w_{s_{1}}(y)\bigl{(}\mathscr{L}-\lambda\bigr{)}^{-1}w_{s_{1}}(y)^{2}dy}\boldsymbol{c}=\frac{1}{1+\tau\lambda}\mathcal{E}_{N}\boldsymbol{c}.$
(3.16)

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)

$1-\frac{2}{1+\tau\lambda}\frac{\int_{-\infty}^{\infty}w_{s_{1}}(y)\bigl{(}\mathscr{L}-\lambda\bigr{)}^{-1}w_{s_{1}}(y)^{2}dy}{\int_{-\infty}^{\infty}w_{s_{1}}(y)^{2}dy}=0.$
(3.17)

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

$\begin{cases}2\Re\bigl{(}\int_{-\infty}^{\infty}w_{s_{1}}(y)\bigl{(}\mathscr{L}-i\lambda_{I}\bigr{)}^{-1}w_{s_{1}}(y)^{2}dy\bigr{)}=\int_{-\infty}^{\infty}w_{s_{1}}(y)^{2}dy,&\\\
2\Im\bigl{(}\int_{-\infty}^{\infty}w_{s_{1}}(y)\bigl{(}\mathscr{L}-i\lambda_{I}\bigr{)}^{-1}w_{s_{1}}(y)^{2}dy\bigr{)}=\tau\lambda_{I}\int_{-\infty}^{\infty}w_{s_{1}}(y)^{2}dy.&\end{cases}$
(3.18)

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

$D_{0,\text{comp}}=\frac{2\mathfrak{a}_{s_{2}}}{N}\frac{\mu(S_{\text{crit}})}{S_{\text{crit}}}.$
(3.19)

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

$B(\lambda,S_{c})-\frac{\kappa}{1+\tau\lambda}=0.$ (3.20)

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

$S_{c}\mu^{\prime}(S_{c})-\mu(S_{c})=0,$

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

$\frac{|B(i\lambda_{I},S_{c})|^{2}}{\text{Re}[B(i\lambda_{I},S_{c})]}-\frac{\mu(S_{c})}{S_{c}}=0,\qquad\tau=-\frac{\text{Im}[B(i\lambda_{I},S_{c})]}{\lambda_{I}\text{Re}[B(i\lambda_{I},S_{c})]},$
(3.21)

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

$\begin{pmatrix}B(\lambda,S_{r})\mathcal{I}_{n}&\mathcal{O}_{n,N-n}\\\
\mathcal{O}_{N-n,n}&B(\lambda,S_{l}(S_{r}))\mathcal{I}_{N-n}\end{pmatrix}\boldsymbol{c}=\frac{\kappa}{1+\tau\lambda}\mathcal{E}_{N}\boldsymbol{c},$
(3.22)

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

$\displaystyle
v\sim\varepsilon^{1-4s_{2}}D\mathfrak{a}_{s_{2}}^{-1}\biggl{[}S_{i}^{\varepsilon}\biggl{(}\sum_{k=1}^{k_{\text{max}}}\frac{(-1)^{k+1}\mathfrak{a}_{ks_{2}}}{D^{k-1}}|y|^{2ks_{2}-1}\varepsilon^{2ks_{2}-1}+R_{D}(0)+\varepsilon
R_{D}^{\prime}(0)y+O(\varepsilon^{2})\biggr{)}$ $\displaystyle+\sum_{j\neq
i}S_{j}^{\varepsilon}G_{D}(x_{i}-x_{j})+\varepsilon\beta_{i,1}y+O(\varepsilon^{2})\biggr{]},$

where we have defined

$\beta_{i,1}=\sum_{j\neq i}S_{j}^{\varepsilon}G_{D}^{\prime}(x_{i}-x_{j}).$
(4.1)

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

$\displaystyle u(x_{i}+\varepsilon
y)\sim\varepsilon^{-2s_{2}}D\bigl{(}U_{i}^{\varepsilon}+\Phi_{i}^{\varepsilon}+\text{h.o.t.}\bigr{)},$
(4.2a) $\displaystyle v(x_{i}+\varepsilon
y)\sim\varepsilon^{-2s_{2}}D\bigl{(}V_{i}^{\varepsilon}+\Psi_{i}^{\varepsilon}+\varepsilon^{2-2s_{2}}\mathfrak{a}_{s_{2}}^{-1}\beta_{i,1}y+\text{h.o.t.}\bigr{)},$
(4.2b)

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

$\mathcal{L}_{i}^{\varepsilon}\equiv\begin{pmatrix}(-\Delta)^{s_{1}}+1-2\tfrac{U_{i}^{\varepsilon}}{V_{i}^{\varepsilon}}&\bigl{(}\tfrac{U_{i}^{\varepsilon}}{V_{i}^{\varepsilon}}\bigr{)}^{2}\\\
-2U_{i}^{\varepsilon}&(-\Delta)^{s_{2}}\end{pmatrix},\quad\boldsymbol{f}_{i}^{\varepsilon}\equiv\begin{pmatrix}\tfrac{1}{\varepsilon}\tfrac{dx_{i}}{dt}\tfrac{dU_{i}^{\varepsilon}}{dy}-\varepsilon^{2-2s_{2}}\mathfrak{a}_{s_{2}}^{-1}\bigl{(}\tfrac{U_{i}^{\varepsilon}}{V_{i}^{\varepsilon}}\bigr{)}^{2}\beta_{i,1}y\\\
-\tfrac{\varepsilon^{2s_{2}}}{D}V_{i}^{\varepsilon}\end{pmatrix}.$ (4.3)

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

$0=\int_{-\infty}^{\infty}(\boldsymbol{P}_{i}^{\varepsilon})^{T}\mathcal{L}_{i}^{\varepsilon}\boldsymbol{\Phi}_{i}^{\varepsilon}dy=\tfrac{1}{\varepsilon}\tfrac{dx_{i}}{dt}\int_{-\infty}^{\infty}P_{i}^{\varepsilon}\tfrac{dU_{i}^{\varepsilon}}{dy}dy-\varepsilon^{2-2s_{2}}\mathfrak{a}_{s_{2}}^{-1}\beta_{i,1}\int_{-\infty}^{\infty}yP_{i}^{\varepsilon}\bigl{(}\tfrac{U_{i}^{\varepsilon}}{V_{i}^{\varepsilon}}\bigr{)}^{2}dy-\tfrac{\varepsilon^{2s_{2}}}{D}\int_{-\infty}^{\infty}Q_{i}^{\varepsilon}V_{i}^{\varepsilon}dy.$

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

$\frac{dx_{i}}{dt}\sim\mathfrak{a}_{s_{2}}\varepsilon^{3-2s_{2}}\frac{\int_{-\infty}^{\infty}yP_{i}^{\varepsilon}\bigl{(}\tfrac{U_{i}^{\varepsilon}}{V_{i}^{\varepsilon}}\bigr{)}^{2}dy}{\int_{-\infty}^{\infty}P_{i}^{\varepsilon}\tfrac{dU_{i}^{\varepsilon}}{dy}dy}\sum_{j\neq
i}S_{j}^{\varepsilon}G_{D}^{\prime}(x_{i}-x_{j}),\qquad(i=1,...,N).$ (4.4)

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

$\frac{\int_{-\infty}^{\infty}yP_{i}^{\varepsilon}\bigl{(}\tfrac{U_{i}^{\varepsilon}}{V_{i}^{\varepsilon}}\bigr{)}^{2}dy}{\int_{-\infty}^{\infty}P_{i}^{\varepsilon}\tfrac{dU_{i}^{\varepsilon}}{dy}dy}\leq
0.$ (4.5)

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

$Q_{i}^{\varepsilon}=C_{i}^{\varepsilon}-\mathfrak{a}_{s_{2}}\int_{-\infty}^{\infty}\frac{(U_{i}^{\varepsilon}(z)/V_{i}^{\varepsilon}(z))^{2}P_{i}^{\varepsilon}(z)}{|y-z|^{1-2s_{2}}}dz,$

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

$\mathscr{M}^{\star}(S)P\equiv-(-\Delta)^{s_{1}}P-P+2\frac{U_{c}}{V_{c}}P-2\mathfrak{a}_{s_{2}}U\int_{-\infty}^{\infty}\frac{(U_{c}(z)/V_{c}(z))^{2}P(z)}{|y-z|^{1-2s_{2}}}dz.$
(4.6)

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

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 8. (A)-(D) Hopf bifurcation thresholds for a one-spike solution
obtained by numerically solving the $\varepsilon$-dependent system (5.1) with
$\varepsilon=0.01$. (E)-(F) Plots of $u(0,t)$ from numerically simulating the
fractional Gierer-Meinhardt system with $\varepsilon=0.01$, $s_{1}=0.5$,
$s_{2}=0.34$, and indicated values of $D_{0}$ with
$\tau=0.95\tau_{h}^{\varepsilon}(D_{0})$ (top) and
$\tau=1.05\tau_{h}^{\varepsilon}(D_{0})$ (bottom). In each case a single spike
solution (obtained using the asymptotics of §2) centred at $x=0$ with
multiplicative noise was used as the initial condition.

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

$\displaystyle[left=\empheqlbrace]\mu(S_{c})=\mathfrak{a}_{s_{2}}^{-1}\varepsilon^{1-2s_{2}}R_{D_{0}\varepsilon^{2s_{2}-1}}(0)S_{c},$
(5.1a) $\displaystyle
B(i\lambda_{I},S_{c})=\frac{\mathfrak{a}_{s_{2}}^{-1}\varepsilon^{1-2s_{2}}}{1+i\tau\lambda_{I}}R_{\tfrac{D_{0}\varepsilon^{2s_{2}-1}}{1+i\tau\lambda_{I}}}(0),$
(5.1b)

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

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 9. (A)-(D) Competition instability thresholds for a two-spike solution
obtained by numerically solving the $\varepsilon$-dependent system (5.2) with
$\varepsilon=0.01$. (E)-(F) Plots of $u(x_{1},t)$ (solid blue) and
$u(x_{2},t)$ (dashed orange) from numerically simulating the fractional
Gierer-Meinhardt system with $\varepsilon=0.01$, $s_{1}=0.5$, at the indicated
values of $s_{2}$ with $D_{0}=0.95D_{0,\text{comp}}^{\varepsilon}$ (top) and
$D_{0}=1.05D_{0,\text{comp}}^{\varepsilon}$ (bottom). In each case a two spike
solution (obtained using the asymptotics of §2) separated by a distance of
$|x_{1}-x_{2}|=1$ with multiplicative noise was used as the initial condition.

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

$\displaystyle[left=\empheqlbrace]\mu(S_{c})=\mathfrak{a}_{s_{2}}^{-1}\varepsilon^{1-2s_{2}}\bigl{(}R_{D_{0}\varepsilon^{2s_{2}-1}}(0)+G_{D_{0}\varepsilon^{2s_{2}-1}}(1)\bigr{)}S_{c},$
(5.2a)
$\displaystyle\mu^{\prime}(S_{c})=\mathfrak{a}_{s_{2}}^{-1}\varepsilon^{1-2s_{2}}\bigl{(}R_{D_{0}\varepsilon^{2s_{2}-1}}(0)-G_{D_{0}\varepsilon^{2s_{2}-1}}(1)\bigr{)}.$
(5.2b)

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

(a)

(b)

(c)

(d)

Figure 10. Time evolution of spike locations $x_{1}(t)$ and $x_{2}(t)$ in a
two-spike solution. In each example $\tau=0.1$, $\varepsilon=0.01$, and
$D=D_{0}\varepsilon^{2s_{2}-1}$ where
$D_{0}=0.8D_{0,\text{comp}}^{\varepsilon}$ and
$D_{0,\text{comp}}^{\varepsilon}$ is the corresponding two-spike competition
instability threshold.

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

$\displaystyle[left=\empheqlbrace]$
$\displaystyle\frac{d(x_{2}-x_{1})}{dt}=2\mathfrak{a}_{s_{2}}\varepsilon^{3-2s_{2}}\frac{\int_{-\infty}^{\infty}yP_{c}^{\varepsilon}\bigl{(}\tfrac{U_{c}^{\varepsilon}}{V_{c}^{\varepsilon}}\bigr{)}^{2}dy}{\int_{-\infty}^{\infty}P_{c}^{\varepsilon}\tfrac{dU_{c}^{\varepsilon}}{dy}dy}G_{D}^{\prime}(x_{2}-x_{1}),$
(5.3a)
$\displaystyle\mu(S_{c}^{\varepsilon})=\frac{\varepsilon^{1-2s_{2}}}{\mathfrak{a}_{s_{2}}}\bigl{(}R_{D}(0)+G_{D}(|x_{2}-x_{1}|)\bigr{)}S_{c}^{\varepsilon}.$
(5.3b)

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

$\displaystyle[left=\empheqlbrace]$
$\displaystyle(-\Delta)^{\frac{1}{2}}U+U-V^{-1}U^{2}=0,\quad(-\Delta)^{s}V-U^{2}=0,$
$\displaystyle-\infty<x<\infty,$ (6.1a) $\displaystyle U,V>0,$
$\displaystyle-\infty<x<\infty,$ (6.1b) $\displaystyle U,V\to 0,$
$\displaystyle\mbox{as}\quad|x|\to+\infty.$ (6.1c)

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

$(-\Delta)^{s}G(x)=\delta(x).$ (6.2)

Then

$G(x)=\dfrac{\Gamma(1-2s)\sin(s\pi)}{\pi}x^{2s-1}.$

###### Proof.

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

$|\xi|^{2s}\hat{G}(\xi)=1.$ (6.3)

Therefore, we have

$G(x)=\frac{1}{\pi}\int_{0}^{\infty}\frac{\cos(x\xi)}{\xi^{2s}}d\xi=x^{2s-1}\frac{1}{\pi}\int_{0}^{\infty}\frac{\cos\xi}{\xi^{2s}}d\xi=\frac{\Gamma(1-2s)\sin(s\pi)}{\pi}x^{2s-1},$
(6.4)

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

$U=\tau_{s}U,\quad
V=\tau_{s}V,\quad\tau_{s}=\left(\frac{\Gamma(1-2s)\sin(s\pi)}{\pi}\int_{{\mathbb{R}}}w^{2}(y)dy\right)^{-1}=\frac{1}{2\Gamma(1-2s)\sin(s\pi)},$
(6.5)

where $w$ is the unique ground state solution to

${(-\Delta)^{1/2}}w+w-w^{2}=0\quad\mbox{in}\quad{\mathbb{R}},\qquad w(x)\to
0\quad\mbox{as}\quad|x|\to\infty.$ (6.6)

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

$w(x)=\frac{2}{1+x^{2}}\quad\mathrm{and}\quad\int_{{\mathbb{R}}}w^{2}dx=2\pi.$
(6.7)

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

$\displaystyle[left=\empheqlbrace]$
$\displaystyle(-\Delta)^{\frac{1}{2}}U+U-V^{-1}U^{2}=0,\quad(-\Delta)^{s}V-\tau_{s}U^{2}=0,$
$\displaystyle-\infty<x<\infty,$ (6.8a) $\displaystyle U,V>0,$
$\displaystyle-\infty<x<\infty,$ (6.8b) $\displaystyle U,V\to 0,$
$\displaystyle\mbox{as}\quad|x|\to\infty.$ (6.8c)

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

$(-\Delta)^{s}V=\tau_{s}h\quad\mbox{in}\quad{\mathbb{R}},\qquad V(x)\to
0\quad\mbox{as}\quad|x|\to\infty,$ (6.9)

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

$T(U^{2})=T(w^{2})+2T(w\phi)+h.o.t.,$ (6.10)

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

$v_{w}=T(w^{2})=\tau_{s}\int_{{\mathbb{R}}}w^{2}(y)G(x-y)dy.$ (6.11)

Expanding the Green function in the following way

$\displaystyle G(x)={}$
$\displaystyle\frac{\Gamma(1-2s)\sin(s\pi)}{\pi}|x|^{2s-1}=e^{(2s-1)\log|x|+\log\Gamma(1-2s)\sin(s\pi)/\pi}$
(6.12) $\displaystyle={}$
$\displaystyle\frac{\Gamma(1-2s)\sin(s\pi)}{\pi}\left(1+(2s-1)\log|x|+(2s-1)^{2}(\log|x|)^{2}/2+\cdots\right),$

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

$T(w^{2})=\tau_{s}\frac{\Gamma(1-2s)\sin(s\pi)}{\pi}\int_{{\mathbb{R}}}w^{2}(y)dy+O(2s-1).$
(6.13)

As a consequence when $x$ is bounded we have

$v_{w}\equiv T(w^{2})=1+O(2s-1)\quad\mbox{and}\quad
T(w\phi)=\frac{1}{\int_{{\mathbb{R}}}w^{2}}\int_{{\mathbb{R}}}w\phi
dy+O(2s-1),$ (6.14)

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

$\frac{U^{2}}{V}=\frac{w^{2}+2w\phi+h.o.t.}{v_{w}+2T(w\phi)+h.o.t.}=\frac{w^{2}}{v_{w}}+2w\phi-2\frac{\int_{{\mathbb{R}}}w\phi
dy}{\int_{{\mathbb{R}}}w^{2}dy}w^{2}+h.o.t.+O(2s-1).$ (6.15)

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

$L(\phi)\equiv(-\Delta)^{\frac{1}{2}}\phi+(1-2w)\phi+2\frac{\int_{{\mathbb{R}}}w\phi
dy}{\int_{{\mathbb{R}}}w^{2}dy}w^{2}=S(w)+N(\phi),$ (6.16)

where

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

and

$N(\phi)=\frac{(w+\phi)^{2}}{T((w+\phi)^{2})}-\frac{w^{2}}{v_{w}}-2w\phi+2\frac{\int_{{\mathbb{R}}}w\phi
dy}{\int_{{\mathbb{R}}}w^{2}dy}w^{2},$

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

$\mathrm{Ker}(L_{0})=\mathrm{Span}\left\\{\frac{dw}{dx}\right\\}.$

  * (b)

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

$L^{*}\phi=(-\Delta)^{\frac{1}{2}}\phi+(1-2w)\phi+2\frac{\int_{{\mathbb{R}}}w^{2}\phi
dx}{\int_{{\mathbb{R}}}w^{2}dx}w.$

Then

$\mathrm{Ker}(L)=\mathrm{Ker}(L^{*})={\mathrm{Span}}\left\\{\frac{dw}{dx}\right\\}.$
(6.17)

###### 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

$L_{0}\phi=-c(\phi)w^{2},\quad\mbox{where}\quad
c(\phi)=2\frac{\int_{{\mathbb{R}}}w\phi dx}{\int_{{\mathbb{R}}}w^{2}dx}.$
(6.18)

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

$\phi=\beta\frac{dw}{dx}+c(\phi)w,$

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

$c(\phi)=2c(\phi)\frac{\int_{{\mathbb{R}}}w^{2}dx}{\int_{{\mathbb{R}}}w^{2}dx}=2c(\phi),$

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

$L^{*}\phi=-c_{1}(\phi)w,\quad\mbox{where}\quad
c_{1}(\phi)=2\frac{\int_{{\mathbb{R}}}w^{2}\phi
dx}{\int_{{\mathbb{R}}}w^{2}dx}.$

Using the fact

$L_{0}(w+x\cdot\partial_{x}w)=-w,$

we have

$\phi-2\frac{\int_{{\mathbb{R}}}w^{2}\phi
dx}{\int_{{\mathbb{R}}}w^{2}dx}(w+x\cdot\partial_{x}w)\in\mbox{Ker}(L_{0}).$
(6.19)

Then

$c_{1}(\phi)=2c_{1}(\phi)\frac{\int_{{\mathbb{R}}}(w+x\cdot\partial_{x}w)w^{2}dx}{\int_{{\mathbb{R}}}w^{2}dx}=2c_{1}(\phi)=2c_{1}(\phi)\frac{\frac{2}{3}\int_{{\mathbb{R}}}w^{3}dx}{\int_{{\mathbb{R}}}w^{2}dx}=2c_{1}(\phi),$

where we used

$\int_{{\mathbb{R}}}w^{3}dx=3\pi,\quad\int_{{\mathbb{R}}}w^{2}dx=2\pi.$

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

$\|\phi\|_{*}=\|\rho(x)^{-1}\phi\|_{L^{\infty}({\mathbb{R}})},\quad\mbox{where}\quad\rho(x)=\frac{1}{(1+|x|)^{\mu}},\qquad\frac{1}{2}<\mu\leq
2.$ (6.20)

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

$L\phi=h+c\frac{dw}{dx}\quad-\infty<x<\infty,\qquad\phi(x)\to 0\quad\mbox{as
}|x|\to\infty,\qquad\langle\phi,\tfrac{dw}{dx}\rangle=0.$ (6.21)

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$

$\|\mathcal{T}h\|_{*}\leq C\|h\|_{*}.$ (6.22)

###  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

$(-\Delta)^{s}v_{w}=\tau_{s}w^{2},\quad v_{w}(x)\to 0~{}\ \mbox{as}~{}\
|x|\to+\infty,$ (6.23)

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

$I_{s}\equiv\left[-100(1-2s)^{-1},~{}100(1-2s)^{-1}\right].$ (6.24)

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

$\begin{split}v_{w}(x)=&\tau_{s}\frac{\Gamma(1-2s)\sin(s\pi)}{\pi}\int_{{\mathbb{R}}}w^{2}(y)dy+(2s-1)\tau_{s}\frac{\Gamma(1-2s)\sin(s\pi)}{\pi}\int_{{\mathbb{R}}}\log|x-y|w^{2}(y)dy\\\
&+(2s-1)^{2}\tau_{s}\frac{\Gamma(1-2s)\sin(s\pi)}{2\pi}\int_{{\mathbb{R}}}w^{2}(y)(\log|x-y|)^{2}dy+o((2s-1)^{2})\\\
=&1+(2s-1)\frac{1}{\int_{{\mathbb{R}}}w^{2}(y)dy}\int_{{\mathbb{R}}}\log|x-y|w^{2}(y)dy\\\
&+\frac{(2s-1)^{2}}{2}\frac{1}{\int_{{\mathbb{R}}}w^{2}(y)dy}\int_{{\mathbb{R}}}(\log|x-y|)^{2}w^{2}(y)dy+o((2s-1)^{2}).\end{split}$
(6.25)

Next we define

$H_{i}(x)=\frac{\int_{{\mathbb{R}}}w^{2}(y)(\log|x-y|)^{i}dy}{\int_{{\mathbb{R}}}w^{2}(y)dy},$
(6.26)

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

$H_{i}(x)=(\log|x|)^{i}+f(x),$ (6.27)

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

$v_{w}(x)\geq c\tau_{s}|x|^{2s-1}\quad\mbox{for}\quad|x|\geq 100(1-2s)^{-1}.$
(6.28)

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

$v_{w}(x)=1+(2s-1)H_{1}(x)+\frac{(2s-1)^{2}}{2}H_{2}(x)+o((1-2s)^{3}).$ (6.29)

  * (b).

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

$v_{w}(x)\geq c\tau_{s}|x|^{2s-1}$ (6.30)

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

$S(w)=\frac{w^{2}}{v_{w}}-w^{2}.$

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

$v_{w}(x)=1+O\left((1-2s)^{1-\delta}\right),$

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

$\frac{1-v_{w}}{v_{w}}w^{2}=O\left((1-2s)^{1-\delta}w^{2}\right).$

Writing

$S(w)=\frac{1-v_{w}}{v_{w}}w^{2}.$ (6.31)

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

$|S(w)|=O\left((1-2s)^{1-\delta}{\rho(x)}\right)\quad\mbox{for}\quad|1-2s|\ll
1.$ (6.32)

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

$\displaystyle|S(w)|\leq
C(1-2s)^{-1}|x|^{1-2s}w^{2}=C(1-2s)^{-1}|x|^{-2}\rho(x)=O(1-2s)\rho(x).$
(6.33)

In conclusion, we have

######  Lemma 6.3.

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

$\|S(w)\|_{*}\leq C(1-2s)^{1-\delta},$

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

$(-\Delta)^{\frac{1}{2}}U+U-\frac{U^{2}}{T(V^{2})}=0.$ (6.34)

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

$(-\Delta)^{\frac{1}{2}}A+A-\frac{A^{2}}{T(A^{2})}=c\frac{dw}{dx},$ (6.35)

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

$\displaystyle(-\Delta)^{\frac{1}{2}}\phi+\phi-2w\phi+2w^{2}\frac{\int_{{\mathbb{R}}}w\phi
dx}{\int_{{\mathbb{R}}}w^{2}dx}$ (6.36)
$\displaystyle=-(-\Delta)^{\frac{1}{2}}w-w+\frac{w^{2}}{v_{w}}+\frac{(w+\phi)^{2}}{T((w+\phi)^{2})}-\frac{w^{2}}{v_{w}}-2w\phi+2w^{2}\frac{\int_{{\mathbb{R}}}w\phi
dx}{\int_{{\mathbb{R}}}w^{2}dx}+c\frac{dw}{dx}$
$\displaystyle=S(w)+N(\phi)+c\frac{dw}{dx}$

and

$N(\phi)=\frac{(w+\phi)^{2}}{T((w+\phi)^{2})}-\frac{w^{2}}{v_{w}}-2w\phi+2w^{2}\frac{\int_{\mathbb{R}}w\phi
dx}{\int_{\mathbb{R}}w^{2}dx}.$ (6.37)

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

${\phi=Q(\phi)\equiv\mathcal{T}(S(w)+N(\phi))}.$

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

$\mathcal{D}=\left\\{\phi\in\mathcal{H}\mid\|\phi\|_{*}\leq
C(1-s)^{1-\delta}\right\\},$ (6.38)

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

$\mathcal{H}=\left\\{\phi\in
L^{\infty}\bigr{|}\big{\langle}\phi,\tfrac{dw}{dx}\big{\rangle}=0\right\\}.$
(6.39)

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

$\|N(\phi)\|_{*}\leq C(\|\phi\|_{*}+\sigma(1-2s))\|\phi\|_{*},$ (6.40)

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

$\displaystyle T((w+\phi)^{2})\geq C(1-2s)|x|^{2s-1},$ $\displaystyle
T(w\phi)\leq C(1-2s)|x|^{2s-1}\|\phi\|_{*},$ $\displaystyle T(\phi^{2})\leq
C(1-2s)^{2-\delta}|x|^{2s-1}\|\phi\|_{*}.$

As a consequence,

$\displaystyle|N(\phi)|\leq$
$\displaystyle~{}\left(\frac{2wv_{w}\phi+v_{w}\phi^{2}-2w^{2}T(w\phi)-w^{2}T(\phi^{2})}{v_{w}T((w+\phi)^{2})}-2w\phi+2w^{2}\frac{\int_{{\mathbb{R}}}w\phi
dx}{\int_{{\mathbb{R}}}w^{2}dx}\right)$ $\displaystyle\leq$
$\displaystyle~{}C\left(\frac{\rho(x)}{(1-2s)(1+|x|)^{2s+1}}+\frac{\rho(x)}{(1-2s)(1+|x|)^{2s+1}}\|\phi\|_{*}\right)\|\phi\|_{*}+C\rho(x)^{2}\|\phi\|_{*}.$

Therefore we have

$|\rho^{-1}N(\phi)|\leq C(\|\phi\|_{*}+(1-2s)^{6s-2})\|\phi\|_{*},$ (6.41)

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

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

$N(\phi)=N_{1}(\phi)+N_{2}(\phi),$

where

$\displaystyle
N_{1}(\phi)=(w+\phi)^{2}\Big{[}\frac{1}{T((w+\phi)^{2})}-\frac{1}{v_{w}}+\frac{2T(w\phi)}{v_{w}^{2}}\Big{]}-(2w+\phi)\phi\frac{2T(w\phi)}{V^{2}}$

and

$\displaystyle N_{2}(\phi)=-2\phi
w\left(1-\frac{1}{v_{w}}\right)+2U^{2}\left(\frac{\int_{{\mathbb{R}}}w\phi
dx}{\int_{{\mathbb{R}}}w^{2}dx}-\frac{T(w\phi)}{v_{w}^{2}}\right)+\frac{\phi^{2}}{v_{w}}.$

It is known that

$\displaystyle v_{w}(x)=1+O((1-2s)^{1-\delta})$

and

$\displaystyle T(w\phi)=\frac{\int_{{\mathbb{R}}}w\phi
dx}{\int_{{\mathbb{R}}}w^{2}dx}+O((1-2s)^{1-\delta}).$

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

$\displaystyle|N_{1}(\phi)|\leq
C(w+\phi)^{2}T(\phi^{2})+C\left(2w\phi+\phi^{2}\right)T(w\phi)\leq
C\rho(x)\|\phi\|_{*}^{2}.$

A similar analysis yields

$\displaystyle|N_{2}(\phi)|\leq
C(1-2s)^{1-\delta}\left(|\phi|w+\rho^{2}\|\phi\|_{*}\right)+C|\phi|^{2},$

and therefore

$\displaystyle\|N(\phi)\|_{*}\leq
C(\|\phi\|_{*}^{2}+\sigma(1-2s)\|\phi\|_{*})$

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

$\|N(\phi_{1})-N(\phi_{2})\|_{*}\leq\kappa\|\phi_{1}-\phi_{2}\|_{*},$ (6.42)

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,

$(-\Delta)^{\frac{1}{2}}\phi+\phi-2w\phi+2w^{2}\frac{\int_{{\mathbb{R}}}w\phi
dx}{\int_{{\mathbb{R}}}w^{2}dx}=S(w)+N(\phi)+c\frac{dw}{dx}.$ (6.43)

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

$\displaystyle[left=\empheqlbrace]$
$\displaystyle(-\Delta)^{\frac{1}{2}}\phi+\phi-2V^{-1}U\phi+V^{-2}U^{2}\psi+\lambda_{s}\phi=0,$
$\displaystyle-\infty<x<\infty,$ (6.44a)
$\displaystyle(-\Delta)^{s}\psi-2U\phi+\tau\lambda_{s}\psi=0,$
$\displaystyle-\infty<x<\infty,$ (6.44b)

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

$\hat{U}=\tau_{s}^{-1}U,\quad\hat{V}=\tau_{s}^{-1}V.$

Then (6.44) can be rewritten as

$\displaystyle[left=\empheqlbrace]$
$\displaystyle(-\Delta)^{\frac{1}{2}}\phi+\phi-2\hat{V}^{-1}\hat{U}\phi+\hat{V}^{-2}\hat{U}^{2}\psi+\lambda_{s}\phi=0,$
$\displaystyle-\infty<x<\infty,$ (6.45a)
$\displaystyle(-\Delta)^{s}\psi-2\tau_{s}\hat{U}\phi+\tau\lambda_{s}\psi=0,$
$\displaystyle-\infty<x<\infty.$ (6.45b)

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

$\psi(x)=2\tau_{s}\int_{{\mathbb{R}}}{G(x-y)\hat{U}(y)\phi(y)dy}=2\frac{\int_{{\mathbb{R}}}w\phi
dy}{\int_{{\mathbb{R}}}w^{2}dy}+O((1-2s)^{1-\delta})\|\phi\|_{H^{1}({\mathbb{R}})}.$
(6.46)

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

$(-\Delta)^{\frac{1}{2}}\phi+\phi-2w\phi+2\frac{\int_{{\mathbb{R}}}w\phi
dx}{\int_{{\mathbb{R}}}w^{2}dx}w^{2}+\lambda_{0}\phi=0.$ (6.47)

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

$\displaystyle[left=\empheqlbrace]$
$\displaystyle(-\Delta)^{\frac{1}{2}}\frac{d\hat{U}}{dx}+\frac{d\hat{U}}{dx}-2\frac{\hat{U}}{\hat{V}}\frac{d\hat{U}}{dx}+\frac{\hat{U}^{2}}{\hat{V}^{2}}\frac{d\hat{V}}{dx}=0,$
$\displaystyle-\infty<x<\infty,$ (6.48a)
$\displaystyle(-\Delta)^{s}\frac{d\hat{V}}{dx}-2\tau_{s}\hat{U}\frac{d\hat{U}}{dx}=0,$
$\displaystyle-\infty<x<\infty.$ (6.48b)

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

$\phi=a\frac{d\hat{U}}{dx}+\phi^{\perp},\quad\mbox{and}\quad\psi=a\frac{d\hat{V}}{dx}+\psi^{\perp},$
(6.49)

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

$\displaystyle[left=\empheqlbrace]$
$\displaystyle(-\Delta)^{\frac{1}{2}}\phi^{\perp}+\phi^{\perp}-2\hat{V}^{-1}\hat{U}\phi^{\perp}+\hat{V}^{-2}\hat{U}^{2}\psi^{\perp}+\lambda_{s}\frac{d\hat{U}}{dx}+\lambda_{s}\phi^{\perp}=0,$
$\displaystyle-\infty<x<\infty,$ (6.50a)
$\displaystyle(-\Delta)^{s}\psi^{\perp}-2\tau_{s}\hat{U}\psi^{\perp}=0,$
$\displaystyle-\infty<x<\infty.$ (6.50b)

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

$\frac{\hat{U}}{\hat{V}}\to
w\quad\mbox{and}\quad\frac{\hat{U}^{2}}{\hat{V}^{2}}\psi^{\perp}\to
2\frac{\int_{{\mathbb{R}}}w\phi^{\perp}dy}{\int_{{\mathbb{R}}}w^{2}dy}w^{2}.$

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

$\lambda_{s}\int_{{\mathbb{R}}}|\phi^{\perp}|^{2}dx=-\int_{{\mathbb{R}}}\left((-\Delta)^{\frac{1}{2}}\phi^{\perp}+\phi^{\perp}-2\frac{\hat{U}}{\hat{V}}\phi^{\perp}+\frac{\hat{U}^{2}}{\hat{V}^{2}}\psi^{\perp}\right)\phi^{\perp}dx.$
(6.51)

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

$\displaystyle L_{1}(\phi^{\perp},\phi^{\perp})={}$
$\displaystyle\int_{{\mathbb{R}}}\left(|(-\Delta)^{\frac{1}{4}}\phi^{\perp}|^{2}+|\phi^{\perp}|^{2}-2w|\phi^{\perp}|^{2}+2\frac{\int_{{\mathbb{R}}}w\phi^{\perp}dx\int_{{\mathbb{R}}}w^{2}\phi^{\perp}dx}{\int_{{\mathbb{R}}}w^{2}dx}\right)$
$\displaystyle\geq{}$
$\displaystyle\frac{\int_{{\mathbb{R}}}w^{3}dx\left(\int_{{\mathbb{R}}}w\phi^{\perp}dx\right)^{2}}{\left(\int_{{\mathbb{R}}}w^{2}dx\right)^{2}}+a_{1}{\inf_{\psi\in
X_{1}}\|\phi^{\perp}-\psi\|_{L^{2}(\mathbb{R})}},$

where $a_{1}>0$ and

$X_{1}=\mbox{Span}\left\\{w,\frac{dw}{dx}\right\\}.$

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

$\lambda_{s}\int_{{\mathbb{R}}}|\phi^{\perp}|^{2}dx\leq 0.$ (6.52)

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.

## References

  * [1] P. Bressloff.  Stochastic Processes in Cell Biology.  Interdisciplinary Applied Mathematics. Springer International Publishing, 2014.
  * [2] NIST Digital Library of Mathematical Functions.  http://dlmf.nist.gov/, Release 1.0.28 of 2020-09-15.  F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.
  * [3] R. L. Frank and E. Lenzmann.  Uniqueness of non-linear ground states for fractional Laplacians in $\mathbb{R}$.  Acta Math., 210(2):261–318, 2013.
  * [4] A. Gierer and H. Meinhardt.  A theory of biological pattern formation.  Kybernetik, 12(1):30–39, Dec 1972.
  * [5] A. A. Golovin, B. J. Matkowsky, and V. A. Volpert.  Turing pattern formation in the Brusselator model with superdiffusion.  SIAM J. Appl. Math., 69(1):251–272, 2008.
  * [6] D. Gomez, M. J. Ward, and J. Wei.  An asymptotic analysis of localized three-dimensional spot patterns for the Gierer-Meinhardt model: Existence, linear stability, and slow dynamics.  SIAM Journal on Applied Mathematics, 81(2):378–406, 2021.
  * [7] D. Gomez, J. Wei, and W. Yang.  Multi-spike solutions to the fractional Gierer-Meinhardt system in a one-dimensional domain.  To appear in Numerical Mathematics: Theory, Methods and Applications.
  * [8] I. S. Gradshteyn and I. M. Ryzhik.  Table of integrals, series, and products.  Elsevier/Academic Press, Amsterdam, eighth edition, 2015.  Translated from the Russian, Translation edited and with a preface by Daniel Zwillinger and Victor Moll, Revised from the seventh edition [MR2360010].
  * [9] B. Henry, T. Langlands, and S. Wearne.  Turing pattern formation in fractional activator-inhibitor systems.  Physical review. E, Statistical, nonlinear, and soft matter physics, 72:026101, 09 2005.
  * [10] Y. Huang and A. Oberman.  Numerical methods for the fractional laplacian: A finite difference-quadrature approach.  SIAM Journal on Numerical Analysis, 52, 11 2013.
  * [11] D. Iron and M. J. Ward.  The dynamics of multispike solutions to the one-dimensional Gierer-Meinhardt model.  SIAM J. Appl. Math., 62(6):1924–1951, 2002.
  * [12] D. Iron, M. J. Ward, and J. Wei.  The stability of spike solutions to the one-dimensional Gierer-Meinhardt model.  Phys. D, 150(1-2):25–62, 2001.
  * [13] H. K. Khudhair, Y. Zhang, and N. Fukawa.  Pattern selection in the schnakenberg equations: From normal to anomalous diffusion.  Numerical Methods for Partial Differential Equations, 2021.
  * [14] A. Lischke, G. Pang, M. Gulian, F. Song, C. Glusa, X. Zheng, Z. Mao, W. Cai, M. M. Meerschaert, M. Ainsworth, and G. E. Karniadakis.  What is the fractional laplacian? a comparative review with new results.  Journal of Computational Physics, 404:109009, 2020.
  * [15] M. Medeiros, J. Wei, and W. Yang.  Existence and stability of symmetric and asymmetric patterns for the half-laplacian Gierer-Meinhardt system in one-dimensional domain.  To appear in Math. Models Methods Appl. Sci.
  * [16] J. D. Murray.  Mathematical biology. II, volume 18 of Interdisciplinary Applied Mathematics.  Springer-Verlag, New York, third edition, 2003.  Spatial models and biomedical applications.
  * [17] Y. Nec.  Spike-type solutions to one dimensional Gierer–Meinhardt model with Lévy flights.  Studies in Applied Mathematics, 129(3):272–299, 2012.
  * [18] Y. Nec and M. Ward.  Dynamics and stability of spike-type solutions to a one dimensional Gierer–Meinhardt model with sub-diffusion.  Physica D: Nonlinear Phenomena, 241:947–963, 05 2012.
  * [19] F. Oliveira, R. Ferreira, L. Lapas, and M. Vainstein.  Anomalous Diffusion: A Basic Mechanism for the Evolution of Inhomogeneous Systems.  Frontiers in Physics, 7, 02 2019.
  * [20] C. Pozrikidis.  The Fractional Laplacian.  Chapman and Hall/CRC, 2018.
  * [21] S. J. Ruuth.  Implicit-explicit methods for reaction-diffusion problems in pattern formation.  J. Math. Biol., 34(2):148–176, 1995.
  * [22] W. H. Tse, J. Wei, and M. Winter.  The Gierer-Meinhardt system on a compact two-dimensional Riemannian manifold: interaction of Gaussian curvature and Green’s function.  J. Math. Pures Appl. (9), 94(4):366–397, 2010.
  * [23] A. M. Turing.  The chemical basis of morphogenesis.  Philos. Trans. Roy. Soc. London Ser. B, 237(641):37–72, 1952.
  * [24] J. C. Tzou, S. Xie, T. Kolokolnikov, and M. J. Ward.  The stability and slow dynamics of localized spot patterns for the 3-D Schnakenberg reaction-diffusion model.  SIAM J. Appl. Dyn. Syst., 16(1):294–336, 2017.
  * [25] M. J. Ward and J. Wei.  Asymmetric spike patterns for the one-dimensional Gierer-Meinhardt model: equilibria and stability.  European J. Appl. Math., 13(3):283–320, 2002.
  * [26] M. J. Ward and J. Wei.  Hopf bifurcation of spike solutions for the shadow Gierer-Meinhardt model.  European J. Appl. Math., 14(6):677–711, 2003.
  * [27] J. Wei and M. Winter.  Spikes for the two-dimensional gierer-meinhardt system: the weak coupling case.  J. Nonlinear Sci., 11(6):415–458, 2001.
  * [28] J. Wei and M. Winter.  Spikes for the Gierer-Meinhardt system in two dimensions: the strong coupling case.  J. Differential Equations, 178(2):478–518, 2002.
  * [29] J. Wei and M. Winter.  Mathematial aspects of pattern formation in biological systems, volume 189.  Applied Mathematical Sciences Series, Springer, 2014.
  * [30] J. Wei and W. Yang.  Multi-bump ground states of the fractional Gierer-Meinhardt system on the real line.  J. Dynam. Differential Equations, 31(1):385–417, 2019.
  * [31] L. Zhang and C. Tian.  Turing pattern dynamics in an activator-inhibitor system with superdiffusion.  Physical Review E, 90, 12 2014.

##  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

$G_{D}(x)=\frac{1}{2}D+\sum_{n=1}^{\infty}\biggl{(}1+\frac{1}{D(\pi
n)^{2s_{2}}}\biggr{)}^{-1}\frac{\cos\pi nx}{(\pi n)^{2s_{2}}},$ (A.1)

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

$\bigl{(}1+\tfrac{1}{D(\pi n)^{2s_{2}}}\bigr{)}^{-1}=\bigl{(}1+\tfrac{1}{D(\pi
n)^{2s_{2}}}\bigr{)}^{-1}\bigl{(}\tfrac{-1}{D(\pi
n)^{2s_{2}}}\bigr{)}^{l}-D(\pi
n)^{2s_{2}}\sum_{k=1}^{l}\bigl{(}\tfrac{-1}{D(\pi n)^{2s_{2}}}\bigr{)}^{k}.$

Substituting into (A.1) then gives

$G_{D}(x)=\tfrac{D}{2}+\tfrac{(-1)^{k_{\max}}}{D^{k_{\max}}}\sum_{n=1}^{\infty}\bigl{(}1+\tfrac{1}{D(\pi
n)^{2s_{2}}}\bigr{)}^{-1}\tfrac{\cos\pi nx}{(\pi
n)^{2(1+k_{\max})s_{2}}}-D\sum_{k=1}^{k_{\max}}\tfrac{(-1)^{k}}{D^{k}}\sum_{n=1}^{\infty}\tfrac{\cos\pi
nx}{(\pi n)^{2ks_{2}}},$ (A.2)

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

$|x|^{\beta-1}=\tfrac{1}{\beta}+2\sum_{n=1}^{\infty}c_{\beta,n}\tfrac{\cos\pi
nx}{(\pi n)^{\beta}},\qquad c_{\beta,n}=\int_{0}^{\pi n}x^{\beta-1}\cos xdx$
(A.3)

Repeated integration by parts yields the identities

$\begin{cases}c_{\beta,n}=(-1)^{n}(\beta-1)(n\pi)^{\beta-2}-(\beta-1)(\beta-2)c_{\beta-2,n},&\beta>2,\\\
c_{\beta,n}=-(\beta-1)\int_{0}^{n\pi}x^{\beta-2}\sin xdx,&\beta>0.\end{cases}$

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

$c_{\beta,n}=(-1)^{n+1}\sum_{r=1}^{q-1}(-1)^{r}(\pi
n)^{\beta-2r}\prod_{l=1}^{2r-1}(\beta-l)+(-1)^{q}\prod_{l=1}^{2q-1}(\beta-l)\int_{0}^{n\pi}x^{\beta-2q}\sin
xdx.$ (A.4)

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

$|x|^{2}=\frac{1}{3}+4\sum_{n=1}^{\infty}(-1)^{n}\frac{\cos\pi nx}{(\pi
n)^{2}}.$ (A.5)

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

$\int_{0}^{\infty}x^{\beta-2q}\sin xdx=\int_{0}^{\infty}x^{\beta-2q}\sin
xdx-\int_{n\pi}^{\infty}x^{\beta-2q}\sin xdx.$

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

$\int_{0}^{\infty}x^{\beta-2q}\sin
xdx=\frac{\pi}{2\Gamma(2q-\beta)\cos(\pi(\beta-2q+1)/2)}=\frac{(-1)^{q}}{2\prod_{l=1}^{2q-1}(\beta-l)}\mathfrak{a}_{\beta/2}^{-1},$

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

$\int_{n\pi}^{\infty}x^{\beta-2q}\sin
xdx=(-1)^{n}(n\pi)^{\beta-2q}+(\beta-2q)\int_{n\pi}^{\infty}x^{\beta-2q-1}\cos
xdx,$

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

$c_{\beta,n}=\frac{1}{2\mathfrak{a}_{\beta/2}}+(-1)^{n}(\beta-1)(\pi
n)^{\beta-2}+a_{\beta,n},$ (A.6a) where
$a_{\beta,n}\equiv(-1)^{n+1}\sum_{r=2}^{q}(-1)^{r}(\pi
n)^{\beta-2r}\prod_{l=1}^{2r-1}(\beta-l)+(-1)^{q+1}\prod_{l=1}^{2q}(\beta-l)\int_{n\pi}^{\infty}x^{\beta-2q-1}\cos
xdx.$ (A.6b)

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

$\sum_{n=1}^{\infty}\frac{\cos\pi nx}{(\pi
n)^{2ks_{2}}}=\mathfrak{a}_{ks_{2}}\biggl{(}|x|^{2ks_{2}-1}-\frac{1}{2ks_{2}}-\frac{2ks_{2}-1}{2}\biggl{(}x^{2}-\frac{1}{3}\biggr{)}-2\sum_{n=1}^{\infty}a_{2ks_{2},n}\frac{\cos\pi
nx}{(\pi n)^{2ks_{2}}}\biggr{)}.$

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

$G_{D}(x)=\sum_{k=1}^{k_{\max}}\frac{(-1)^{k-1}\mathfrak{a}_{ks_{2}}}{D^{k-1}}|x|^{2ks_{2}-1}+R_{D}(x)$
(A.7a) where $k_{\max}=\lceil\tfrac{1}{2s_{2}}-1\rceil$ and the regular part
$R_{D}(x)$ is given by
$\begin{split}R_{D}(x)=&\frac{1}{2}D+\frac{(-1)^{k_{\max}}}{D^{k_{\max}}}\sum_{n=1}^{\infty}\biggl{(}1+\frac{1}{D(\pi
n)^{2s_{2}}}\biggr{)}^{-1}\frac{\cos\pi nx}{(\pi
n)^{2(1+{k_{\max}})s_{2}}}+\frac{1}{2s_{2}}\sum_{k=1}^{k_{\max}}\frac{(-1)^{k}\mathfrak{a}_{ks_{2}}}{kD^{k-1}}\\\
&+\sum_{k=1}^{k_{\max}}\frac{(-1)^{k}(2ks_{2}-1)\mathfrak{a}_{ks_{2}}}{2D^{k-1}}\biggl{(}x^{2}-\frac{1}{3}\biggr{)}+2\sum_{n=1}^{\infty}\biggl{(}\sum_{k=1}^{k_{\max}}\frac{(-1)^{k}\mathfrak{a}_{ks_{2}}a_{2ks_{2},n}}{D^{k-1}(\pi
n)^{2ks_{2}}}\biggr{)}\cos\pi nx,\end{split}$ (A.7b) where $a_{2ks_{2},n}$ is
given by (A.6b) with $\beta=2ks_{2}$ and $q=\lceil ks_{2}\rceil$, while
$\mathfrak{a}_{ks_{2}}$ is given by
$\mathfrak{a}_{ks_{2}}=-\frac{2ks_{2}}{\pi}\Gamma(-2ks_{2})\sin(\pi ks_{2}).$
(A.7c)

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

$\sum_{n=1}^{\infty}\frac{\cos n\pi
x}{n\pi}=-\frac{1}{\pi}\log|x|-\frac{1}{\pi}-\frac{2}{\pi}\sum_{n=1}^{\infty}\biggl{(}\text{Si}(n\pi)-\frac{\pi}{2}\biggr{)}\frac{\cos
n\pi x}{n\pi}$ (A.8)

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

$G_{D}(x)=\sum_{k=1}^{r-1}\frac{(-1)^{k-1}\mathfrak{a}_{ks_{2}}}{D^{k-1}}|x|^{2ks_{2}-1}+\frac{(-1)^{r}}{\pi
D^{r-1}}\log|x|+R_{D}(x)$ (A.9a) where $R_{D}(x)$ is given by
$\begin{split}R_{D}(x)=&\frac{1}{2}D+\frac{(-1)^{r}}{D^{r}}\sum_{n=1}^{\infty}\biggl{(}1+\frac{1}{D(\pi
n)^{1/r}}\biggr{)}^{-1}\frac{\cos\pi nx}{(\pi
n)^{1+1/r}}+r\sum_{k=1}^{r-1}\frac{(-1)^{k}\mathfrak{a}_{k/(2r)}}{kD^{k-1}}\\\
&+\sum_{k=1}^{r-1}\frac{(-1)^{k}(k/r-1)\mathfrak{a}_{k/(2r)}}{2D^{k-1}}\biggl{(}x^{2}-\frac{1}{3}\biggr{)}+2\sum_{n=1}^{\infty}\biggl{(}\sum_{k=1}^{r-1}\frac{(-1)^{k}\mathfrak{a}_{k/(2r)}a_{k/r,n}}{D^{k-1}(\pi
n)^{k/r}}\biggr{)}\cos\pi nx\\\ &+\frac{(-1)^{r}}{\pi
D^{r-1}}\biggl{(}1+2\sum_{n=1}^{\infty}\bigl{(}\text{Si}(n\pi)-\tfrac{\pi}{2}\bigr{)}\frac{\cos
n\pi x}{n\pi}\biggr{)},\end{split}$ (A.9b)

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

$\begin{cases}\text{if $a_{\pm}=0$ then $b_{\pm}=\varphi_{\pm N}L^{\beta}$ and
$\varphi_{\pm i}=\varphi_{\pm N}|y_{\pm i}/L|^{-\beta}$},&\text{(Case 1)},\\\
\text{if $b_{\pm}$ is given then $a_{\pm}=\varphi_{\pm N}-b_{\pm}L^{-\beta}$
and $\varphi_{\pm i}=\varphi_{\pm N}+b_{\pm}(|y_{\pm
i}|^{-\beta}-L^{-\beta})$},&\text{(Case 2)}.\end{cases}$ (B.1)

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