1 Introduction↩︎

Reaction–diffusion (RD) systems have been used as a fundamental framework to model pattern formation in biological, chemical, and physical systems [1]. A particularly influential example is the Gierer-Meinhardt (GM) model, which is originally proposed to describe the morphogen interactions underlying biological pattern formation [2]. The classical GM model involves two interacting chemical species: a short-range self-activator and a long-range inhibitor. Through the mechanism of local self-activation and lateral inhibition, this system captures the emergence of stationary spatial structures such as spots, stripes, and labyrinthine patterns [3], [4].

Over the past decades, the two-component GM model has been extensively analyzed both analytically and numerically. A focus of these work has been on identifying the conditions for Turing instability, the bifurcation structure of patterned solutions [5], and the construction of spatially localized patterns with their stability properties [6]–[8]. These studies have provided detailed asymptotic frameworks for understanding how spike equilibria form, interact, and evolve in one- and two-dimensional domains [9]–[11], and have provided a cornerstone for understanding the formation of self-organized biological patterns.

Despite its success, the two-component formulation can be restrictive for certain biological or chemical systems, where additional species often play important roles as secondary inhibitors or facilitators. To capture such complexity, extensions have been proposed to three or more-component GM-type systems [12]–[14]. The inclusion of an additional component can fundamentally alter dynamical behavior, leading to richer bifurcation structures, oscillatory instabilities, and novel patterns formation regimes that lie beyond the scope of the classical two-component model [15]–[17].

Motivated by understanding how extra components influence the existence and stability of localized structures, in this paper, we analyze a three-component GM model, which is given as follows:

\[\tag{1} \begin{alignat}{2} \frac{\partial u}{\partial t} & = a-u+\frac{u^3}{w v}+\delta_1 \frac{\partial^2 u}{\partial x^2}, \,\, x \in (-l,l), t>0, \tag{2}\\ \theta\frac{\partial v}{\partial t} & = u^2-b v+D_v\frac{\partial^2 v}{\partial x^2},\,\, x\in (-l,l), t>0, \tag{3}\\ \tau\frac{\partial w}{\partial t} & = u -cw +\delta_2 \frac{\partial^2 w}{\partial x^2}, \,\, x \in (-l,l), t>0, \tag{4}\\ \intertext{with Neumann boundary conditions} & \left. u_x\right|_{x=\pm l} = \left. v_x\right|_{x=\pm l} = \left. w_x\right|_{x=\pm l} =0. \nonumber \end{alignat}\] This model contains an additional reactant \(w\) that acts as an inhibitor of \(u\) and interacts with \(v\) indirectly through its dependence on the activator \(u\). The parameters \(a, b, c > 0\), and we focus on the semi-strong regime where \(\delta_1\ll 1, \delta_2=\delta_1^2\) and \(D_v=\mathcal{O}(1)\). In this regime, the activator \(u\) diffuses more slowly than \(v\), and \(w\) exhibits essentially the same behavior as \(u\). To account for the possibility that different species evolve on distinct intrinsic timescales, we introduce the time-scaling parameters \(\theta\) and \(\tau\). Our goal is to investigate the new dynamics introduced by the extra component \(w\). Moreover, we assume a nontrivial background \(a>0\) for the activator, which, as we will show below, is essential to generate new spikes.

1.1 Main Results↩︎

For the three-component GM model (1 ), we construct localized spike equilibria in the limit \(\delta_1 \to 0\) using asymptotic analysis and numerical-path following methods. In the presence of the nontrivial activator background \(a>0\), we show that spike nucleation (or insertion) can occur as the inhibitor diffusivity \(D_v\) decreases, whereby a new spike is generated either at the domain boundaries or from the quiescent background between adjacent spikes. This process is illustrated in Figure 1, where initially a interior spike is located at the domain center, then decreasing \(D_v\) first triggers the generation of two boundary spikes. With further decrease, new spikes nucleate in between the interior spike and each boundary spike.

By relaxing the two time-scaling parameters \(\theta\) and \(\tau\) in front of the equations of \(v\) and \(w\), respectively, we investigate the stability of the interior spike equilibria. In particular, two classes of eigenvalues are considered. The large eigenvalues that correspond to the structural stability of the pattern are studied by deriving a novel nonlocal eigenvalue problem (NLEP). The other type are the small eigenvalues that are associated with slow spike motion dynamics. For this instability, we derive asymptotic instability thresholds explicitly, which mark the onset of oscillatory instabilities in the spike motion.

Figure 2 and 3 reveal the two qualitatively distinct dynamical behaviors. As shown in Figure 2, as \(\theta\) increases beyond a critical value, the spike undergoes oscillations in amplitude, leading to instability of the localized structure. In contrast, in Figure 3, by increasing \(\tau\) and setting \(\theta=0\), the spike experiences an oscillatory instability in position, resulting in periodic motion of the spike across the domain. The oscillatory behavior decays as \(\tau\) increases, and and for sufficiently large \(\tau,\) the spike ultimately drifts to the boundary and remains there, as shown in Figure 3 (b).

1.2 Outline↩︎

The rest of this paper is organized as follows. In Section 2 we use the method of matched asymptotic expansions in the limit \(\delta_1\to 0\) to construct quasi-steady state spike solutions for (1 ). We will derive critical values for inhibitor diffusivity \(D_v\) at which spike nucleation will occur through slow passage beyond the saddle-node of a nonlinear boundary value problem defined in the outer region away from the core of a spike. In addition, for the regime \(a\ll 1\), we obtain an explicit analytical result for the spike profile.

Section 3 \(\&\) 4 are devoted to the stability of these spike equilibria in the regime \(a\ll 1\). In Section 3, we analyze the spectrum of large eigenvalues by deriving a novel Nonlocal eigenvalue problem (NLEP). Here, three scenarios are considered depending on the value of \(\theta\) and \(\tau\). We find that the large eigenvalue instability associated with \(\tau\) alone is not observed in full numerical simulations of an interior-spike pattern. Instead, before the large eigenvalue instability is reached, the small eigenvalue becomes unstable and triggers oscillation in slow spike dynamics.

Finally, Section 5 discusses the implications of our results in relation to existing theory and concludes with directions for future work.

2 Semi-strong interaction asymptotic analysis↩︎

In this section, we apply an asymptotic analysis of semi-strong interactions to construct a one-spike solution centered at \(x=0\) to the GM model (1 ) defined on the domain \(|x|\leq l\), which satisfies the following steady-state problem: \[\label{Gmhom} \begin{alignat}{2} a-u+\frac{u^3}{w v}+\delta_1 \frac{\partial^2 u}{\partial x^2}&=0,\\ u^2-b v+D_v\frac{\partial^2 v}{\partial x^2}&=0,\\ u -cw +\delta_1^2 \frac{\partial^2 w}{\partial x^2}&=0, \end{alignat}\tag{5}\] which are subject to Neumann boundary conditions.

We will show a novel type of instability that can arise when the system has a nonzero background, that is, \(a\neq 0\). This instability is known as spike insertion or nucleation, which involves the emergence of new spikes either near the domain boundaries or at the midpoint between adjacent spikes. For certain ranges of the parameters \(a, b,\) and \(c\), we show that this instability first occurs as \(D_v\) decreases below a saddle-node point associated with a nonlinear outer problem.

Our analysis follows a similar framework to the recent study [19], which couples a nonlinear inner problem for the spike profile to a nonlinear reduced scalar BVP defined in the outer region away from the spike.

In the mean time, we are interested in the case which arrive at a solution for equation 5 that exhibits homoclinic properties in space, before eventually reaching a homogeneous steady state \[u \to bc + a, \qquad v\to \frac{u^2}{b}, \qquad w\to \frac{u}{c} \quad asx \to \pm l.\] In this case, however, we will show that, in the presence of the homogeneous steady state, nucleation instabilities do not occur.

2.1 Asymptotic construction of quasi-equilibria: The inner solution↩︎

To construct a quasi steady-state spike solution for 5 , in the inner region near \(x=0\), we introduce the inner variables \(y, U(y), V(y)\) and \(W(y)\) by

\[\label{sec295inner95var} y= \frac{x}{\sqrt{\delta_1}}, \quad u=\frac{U}{\sqrt{\delta_1}}, \quad v=\frac{V}{\sqrt{\delta_1}}, \quad w=\frac{W}{\sqrt{\delta_1}},\tag{6}\] so that the steady-state problem (5 ) transforms to \[\tag{7} \begin{align} &U_{yy} -U+ \frac{U^3}{WV}+a\sqrt{\delta_1}=0,\tag{8}\\ &D_v V_{yy}-b\delta_1V+\sqrt{\delta_1} U^2=0,\tag{9}\\ &\delta_1W_{yy}-cW+U=0\tag{10} \end{align}\] on \(y \geq 0\), with \(U_y(0)=0\), \(V_y(0)=0\), and \(W_y(0)=0\). Upon expanding \(U(y), V(y), W(y)\) in the power of \(\sqrt{\delta_1}\) \[\label{sec295expansion} U(y)=U_0+\sqrt{\delta_1} U_1+O(\delta_1),\quad V(y)=V_0+\sqrt{\delta_1} V_1+O(\delta_1),\quad W(y)=W_0+\sqrt{\delta_1} W_1+O(\delta_1),\tag{11}\] and substituting them into (7 ) we find that \(W=c^{-1}U\) and \(V_0\) is a constant to be determined. Then \(V_1\) satisfies \[\label{sec295V1eqn} V_{1yy}=-\frac{U_0^2}{D_v}.\tag{12}\] To reflect the non-zero far-field of the activator, we take \(U_0\) to be the homoclinic solution of \[\label{sec295U095eqn} U_{0yy}-U_0+c\frac{U_0^2}{V_0}+a\sqrt{\delta_1}=0;\quad y\geq0, \quad U_{0y}(0)=0, \quad U_0(0)>0.\tag{13}\] Solving (13 ) yields the homoclinic solution \[\label{sec295U0} U_0=\frac{V_0}{c}\left(w_c(y)+\gamma\right),\tag{14}\] where \(w_c(y)\) is the unique solution to \[\label{sec295w95eqn} w_{cyy}-(1-2\gamma)w_c+w_c^2=0;\quad w_{cy}(0)=0, \quad w_c(0)>0 \quad \lim_{y\to\infty}w_c=0.\tag{15}\] Here, \(\gamma\) satisfies the quadratic equation \[\label{sec295gam95eqn} \gamma^2-\gamma+\frac{ac\sqrt{\delta_1}}{V_0}=0.\tag{16}\] Since we require \(\gamma < \frac{1}{2}\) in (15 ), we must take \(\gamma\) as the smaller root of (16 ), which is given by \[\label{sec295gam} \gamma=\frac{1-\sqrt{1-\frac{4ac\sqrt{\delta_1}}{V_0}}}{2}.\tag{17}\]

Note that for \(\delta_1\ll 1\), we have from Taylor expansion that \(\gamma \sim \frac{ac\sqrt{\delta_1}}{V_0} +O(\delta_1)\). Then the explicit solution to (15 ) is calculated as \[\label{sec295w} w_c(y)=\frac{3}{2}\left(1-2\gamma\right) \text{sech}^2\left(\frac{\sqrt{1-2\gamma}}{2}y\right),\tag{18}\] which, from (14 ), determines the homoclinic solution to (13 ) up to a constant \(V_0\) to be found.

We now summarize the result in the inner region as follows:

Next, we determine the far-field behavior of \(U\) and \(V\). By letting \(y \to \infty\), we use (6 ), (11 ) and (14 ) to conclude that \[\label{sec295u95far} u\sim\frac{V_0\gamma}{c\sqrt{\delta_1}},\quad \text{as} \quad y\to \infty.\tag{19}\] In terms of the far-field behavior of \(V\), by substituting (14 ) into (12 ), we obtain that \[\begin{align} \label{sec295V1yy} V_{1yy}&=\frac{1}{D_v}\left(\frac{V_0}{c}\right)^2\left(w_c+\gamma\right)^2 \nonumber \\ &=-\frac{V_0^2}{D_vc^2}\gamma^2-\frac{V_0^2}{D_vc^2}\left(w_c^2+2\gamma w_c\right). \end{align}\tag{20}\] Upon integrating (20 ) using the boundary condition \(V_{1y}(0) = 0\), we obtain for any \(y > 0\) that \[\label{sec295V1y} V_{1y}=-\frac{V_0^2}{D_vc^2}\gamma^2 y-\frac{V_0^2}{D_vc^2}\left(\int_0^yw_c^2\,ds+2\gamma\int_0^yw_c \,ds\right).\tag{21}\] To determine the limiting behavior as \(y \to \infty\), we use (18 ) to calculate \[\label{sec295int95w} \int_0^{\infty} w_c \,dy=3\sqrt{1-2\gamma}, \quad \int_0^{\infty} w_c^2 \,dy=3\left(1-2\gamma\right)^{3/2}.\tag{22}\] Now using (22 ) in (21 ), we conclude that \[\label{sec295V1y95far} \lim_{y\to\infty} \left(V_{1y}+\frac{V_0^2}{D_vc^2}\gamma^2 y\right) = -\frac{V_0^2}{D_vc^2}\left( 3\left(1-2\gamma\right)^{3/2}+6\gamma\sqrt{1-2\gamma} \right) = -3\frac{V_0^2}{D_vc^2}\sqrt{1-2\gamma}.\tag{23}\] In this way, by using (23 ) together with (6 ) and (11 ) we obtain that \(v\) has far-field behavior \[\label{sec295v95far} v\sim \frac{V_0}{\sqrt{\delta_1}}-3\frac{V_0^2}{D_vc^2}\sqrt{1-2\gamma}y-\frac{V_02^2}{2D_vc^2}\gamma^2y^2, \quad \text{as} \quad y\to \infty.\tag{24}\]

2.2 Asymptotic construction of quasi-equilibria: The outer solution↩︎

In this section, we match the far-field behavior of the inner solution with an outer solution valid on \(0^+ < |x| < l\), to determine \(V_0\) and complete the construction of the spike solution. As this outer solution is even, we focus on half of the domain, where it satisfies \[\tag{25} \begin{alignat}{2} a-u+\frac{u^3}{w v}&=0, \quad u_x(l)=0, \tag{26}\\ D_vv_{xx}+u^2-b v&=0, \quad v_x(l)=0,\tag{27}\\ u -cw &=0, \quad w_x(l)=0. \tag{28} \end{alignat}\] Solving \(w\) in (28 ) and substitute it into (26 ), we obtain that \[\label{sec295v} v=\frac{cu^2}{u-a} \quad \text{for} \quad u>a,\tag{29}\] which implies \[\label{sec295v95x} v_x=-\frac{cu(2a-u)}{(u-a)^2}u_x \quad \text{for} \quad u>a.\tag{30}\] Upon substituting (29 ) and (30 ) into (27 ) we obtain that the outer problem for \(u\) is

\[\label{sec295u95outer95problem} D_v\left(f(u)u_x\right)_x=R(u), \quad 0^+<x<l, \quad u_x(l)=0,\tag{31}\] where \[\label{sec295f44R} f(u)=\frac{cu(2a-u)}{(u-a)^2}, \quad\text{and} \quad R(u)=u^2-b\frac{cu^2}{u-a}.\tag{32}\]

The problem (31 ) is well-posed when \(u>0\) and \(u_x>0\) on \((0^+,l)\). This implies that we must have \(a<u<2a\) on \((0^+,l)\).

Next, we derive the matching conditions between the inner and outer solution. From (19 ), together with (17 ) for \(\gamma\), the first matching condition for the outer solution is \[\label{sec295u400p41} u(0^+)=\frac{V_0}{2c\sqrt{\delta_1}}\left(1-\sqrt{1-\frac{4ac\sqrt{\delta_1}}{V_0}}\right).\tag{33}\]

For \(\frac{4ac\sqrt{\delta_1}}{V_0}<1\), we claim that \(u(0^+)>a\). To establish this inequality, we introduce \(z=\frac{2ac\sqrt{\delta_1}}{V_0}\) and observe that \[\frac{u(0^+)}{a}=\frac{1-\sqrt{1-2z}}{z}.\] Since \(\sqrt{1-2z}<1-z\) on \(0<z<\frac{1}{2}\), the expression above yields \(u(0^+)>a\) whenever \(\frac{4ac\sqrt{\delta_1}}{V_0}<1\).

The second matching condition involves matching the flux \(u_x\) as \(x\to 0^+.\) We first observe that the \(\mathcal{O}(y^2)\) term in the far-field behavior (24 ) matched exactly with the quadratic term in the near-field behavior of \(u\) as \(x\to 0^+\) that arises from the \(u^2\) term in (27 ). From (24 ) we conclude that

\[\label{sec295v95x40041} v_x=-\frac{3V_0^2}{\sqrt{\delta_1}D_vc^2}\sqrt{1-2\gamma} \quad \text{as} \quad x\to 0^+.\tag{34}\]

Using (30 ) and (34 ), we obtain the second matching condition for the outer solution in terms of \(v\), which is given by \[\label{sec295matching} \lim_{x\to 0^+}-v_x=\lim_{x\to 0^+}f(u)u_x=\frac{3V_0^2}{\sqrt{\delta_1}D_vc^2}\sqrt{1-2\gamma}.\tag{35}\]

Next, we establish the following lemma.

The remaining steps in the analysis to construct the quasi-steady state is to solve (31 ). We define \(\mathcal{G}\) by

\[\label{sec295G95prime} \mathcal{G'(\xi)} \equiv -R(\xi) f(\xi)=c(bc+a-\xi)(2a-\xi)\frac{\xi^3}{(\xi-a)^3}>0 \quad \text{on} \quad a<\xi<2a.\tag{38}\] A first integral of (38 ) yields \[\label{sec295G} G(\xi)=c\left(\frac{1}{3}(\xi-a)^3+\frac{2a-bc}{2}(\xi-a)^2-2abc(\xi-a)-2a^3\ln(\xi-a)+\frac{a^4-2a^3bc}{\xi-a}-\frac{a^4bc}{2(\xi-a)^2}\right).\tag{39}\]

Upon first multiplying (31 ) by \(f(u)u_x\) and then integrating, we use the monotonicity of \(u(x)\) and equation (38 ) to obtain \[\label{sec295first32integral32of32outer} -\frac{1}{2}D_v\left(f(u)u_x\right)^2=\int_{x}^lR(u)f(u)u_x dx dx=-\int_{u(x)}^{\mu}\mathcal{G'(\xi)}d\xi =\mathcal{G}(\mu)-\mathcal{G}(u(x)).\tag{40}\]

Here \(\mathcal{G}(\xi)\) is given in (39 ) and \(\mu\equiv u(l)\) satisfies \(a < \mu \leq 2a\). By taking the positive square root in (40 ), we have \[\label{sec295sqrt95outer} f(u)u_x=\sqrt{\frac{2}{D_v}}\sqrt{\mathcal{G}(\mu)-\mathcal{G}(u)}.\tag{41}\]

Letting \(x\to 0^+\) in (41 ) and imposing the matching condition (35 ) we conclude that \(V_0\) is related to \(\mu\) by the following nonlinear algebraic equation \[\label{sec295Newton951} \frac{3V_0^2}{\sqrt{2\delta_1}\sqrt{D_v}c^2}\sqrt{1-2\gamma}=\sqrt{\mathcal{G}(\mu)-\mathcal{G}(u(0^+))},\tag{42}\] where \(u(0^+)\) and \(\gamma\) are given in terms of \(V_0\) by \[\label{sec295u400944341} u(0^+)=\frac{V_0}{c\sqrt{\delta_1}}\gamma, \quad \text{where}\quad \gamma=\frac{1-\sqrt{1-\frac{4ac\sqrt{\delta_1}}{V_0}}}{2}.\tag{43}\]

Next, upon integrating the separable ODE (41 ), we get an implicit relation for \(u(x)\) on \(0^+ < x < l\) given by \[\label{sec295out95solution} \chi(u(x))=\sqrt{\frac{2}{D_v}}x, \quad \text{where} \quad \chi(u(x))=\int_{u(0^+)}^{u(x)} \frac{f(\xi)}{\sqrt{\mathcal{G}(\mu)-\mathcal{G}(\xi)}}d \xi.\tag{44}\]

Then, by setting \(x = l\) and \(\mu\equiv u(l)\) in (44 ), we obtain an implicit equation for \(\mu\) given by

\[\label{sec295implicit95mu} \chi(\mu)=\int_{u(0^+)}^{\mu} \frac{f(\xi)}{\sqrt{\mathcal{G}(\mu)-\mathcal{G}(\xi)}}d \xi=\sqrt{\frac{2}{D_v}}l,\tag{45}\] with \(u(0^+) > a\) given by (43 ).

Since the integral in (44 ) is improper at \(\xi = \mu\), to obtain a more tractable formula for \(\chi(\mu)\) we integrate the expression in (45 ) by parts by using \(f(\xi) = -\mathcal{G}'(\xi)/R(\xi)\). This yields the proper integral

\[\label{sec295proper95chi} \chi(\mu)=-2\frac{\sqrt{\mathcal{G}(\mu)-\mathcal{G}(u(0^+))}}{R(u(0^+))}+2\int_{u(0^+)}^{\mu}\frac{\sqrt{\mathcal{G}(\mu)-\mathcal{G}(\xi)}}{R^2(\xi)}R'(\xi) d\xi,\tag{46}\] where \(R(\xi) = \xi^2 - bc\frac{\xi^2}{\xi-a}\) and \(\mathcal{G}(\xi)\) is given in (39 ). On the range \(\mu > u(0^+)\), we observe that \(\chi(\mu)\) is positive since \(R(\xi) < 0\) and \(R'(\xi) > 0\) on \(a < \xi < 2a\). Moreover, differentiating (46 ) yields \[\label{sec295chi95prime} \chi'(\mu)= \frac{f(\mu)}{\sqrt{\mathcal{G}(\mu)-\mathcal{G}(u(0^+))}}\frac{R(\mu)}{R(u(0^+))}+\mathcal{G}'(\mu)\int_{u(0^+)}^{\mu}\frac{R'(\xi)}{\sqrt{\mathcal{G}(\mu)-\mathcal{G}(u(0^+))}R^2(\xi)}d\xi.\tag{47}\] Since \(f(\xi) > 0, \mathcal{G}'(\xi) > 0\) and \(R'(\xi) > 0\) on \(a < \xi < 2a\), it follows that \(\chi(\mu)\) is a monotonically increasing function of \(\mu\) on \(u(0^+) < \mu < 2a\) that reaches its maximum value at \(\mu = 2a\). As such, recalling that \(\mu = u(l)\), we define

\[\label{sec295chi95max} \chi_{max} \equiv \chi(2a), \mu_{max} \equiv 2a.\tag{48}\] We now summarize our asymptotic construction of a one-spike solution to (5 ) and the mechanism of nucleation instability as follows:

To validate the analytical spike construction, Figure 4 compares the asymptotic spike profile with full numerical simulations. In particular, the predicted spike amplitude, obtained by solving the coupled nonlinear system (42 ) and (45 ) using Newton’s method, is shown with the numerically computed steady-state solution. The good agreement between theory and computation confirms the accuracy of the asymptotic construction.

In Figure 5 we compare the asymptotic prediction \(D_{nuc}\) in (?? ) for the critical diffusion threshold of spike nucleation with the numerically computed location of the saddle-node bifurcation point of the full system (1 ) using pde2path [20]. Given the condition that \(a<bc=1\), we fix \(b=1\) and \(c = 1\) in Figure 5, and plot \(D_{nuc}\) as a function of \(a\) for values satisfying \(a<1.\) The comparison shows a nice agreement, which verifies the accuracy of the asymptotic prediction.

2.3 Global bifurcation diagram and full PDE simulations↩︎

To illustrate the global bifurcation structure, in Figure 6, we consider the parameter set for \(a = 0.5, \delta_1 = 0.01^2, b=1, c=1\) and \(l=4\), and obtain the numerically computed global bifurcation diagram by path-following the single-spike branch for the GM model (1 ) using pde2path [18]. The saddle-node bifurcation point at point \((b)\), indicated by the red star, occurs at \(D_{nuc} \approx 1.06\) and marks the onset of spike nucleation. This value is well-approximated by the asymptotic prediction for \(D_{nuc}\) in (?? ) from our asymptotic theory, at which the outer solution ceases to exist (see Figure 5 for \(a = 0.5\)).

This global bifurcation diagram of single-spike steady-states shown in Figure 6 provides a detailed view on how new spikes are created near the domain boundaries as \(D_v\) approaches the saddle-node bifurcation point. In the right panel of Figure 6 the solution profiles \(u(x)\) corresponding to the points marked in the left panel are shown. Starting from the bottom branch at point \((a)\), a single interior spike is present; this state terminates as \(D\) decreases below \(D_{nuc}\). Traversing the middle branch from the saddle-node bifurcation, we observe at point \((c)\) that the nucleation of new boundary spikes is fully developed at the domain endpoints.

To further validate these bifurcation scenarios, Figure 1 presents full time-dependent PDE simulations of (1 ) computed with FlexPDE [18] as the diffusion rate \(D_v\) slowly decreases in time by \(D = e^{-\rho t}\) with \(\rho = 10^{-4}.\) for parameters \(a=0.5, \delta_1 = 0.01^2, b=1, c=1\), and \(l=4\). As \(D_v\) decreases, boundary spikes first appear when \(D_v\approx 1.07\). Upon further decrease, new spikes nucleate at the midpoint between the interior spike and each boundary spike. These transition values of \(D_v\), obtained from the time-dependent simulations, are well-approximated by the critical thresholds predicted asymptotically for the non-existence of the outer solution.

2.4 Asymptotics of the outer solution for \(a\) small↩︎

In this subsection, we approximate the outer problem (31 ) by a linear problem when \(a\ll 1\). In the limit \(a\ll 1\), no spike nucleation behavior occurs.

For \(a \ll 1\), we obtain from (5 c) that \(w=\frac{u}{c}\), and plug it in (5 a), we get \[\label{sec295u95out95a95small} u\sim a+O(a^2)\tag{49}\] in the outer region. As a result, from (5 b), we obtain that \(v(x)\) satisfies \[\label{sec295v95eqn95a95small} D_vv_{xx}-bv=-a^2+O(a^3), \quad 0^+<x<l, \quad v_x(l)=0,\tag{50}\] with the matching condition \(v(0^+)=\frac{V_0}{\sqrt{\delta_1}}\). Upon neglecting the \(O(a^3)\) term in (50 ), we calculate that \[\label{sec95295v95a95small} v(x)=\frac{a^2}{b}+\left(\frac{V_0}{\sqrt{\delta_1}}-\frac{a^2}{b}\right)\frac{\cosh\left(\sqrt{b}(l-|x|)/\sqrt{D_v}\right)}{\cosh(\sqrt{b}l/\sqrt{D_v})}.\tag{51}\]

By imposing the second matching condition given by (34 ), we obtain that \(V_0\) must satisfy \[\label{sec295V095match} \frac{3V_0^2}{\sqrt{\delta_1}D_vc^2}\sqrt{1-2\gamma}=\sqrt{\frac{b}{D_v}}\left(\frac{V_0}{\sqrt{\delta_1}}-\frac{a^2}{b}\right)\tanh\left(\sqrt{\frac{b}{D_v}}l\right).\tag{52}\]

For \(a\ll 1\), (43 ) yields \(\gamma\sim \frac{ac\sqrt{\delta_1}}{V_0} \ll 1\). By setting \(\gamma\ll 1\) in (52 ), we get \(\sqrt{1-2\gamma}\sim1-\gamma\) and that

\[\label{sec295V095quadratic} 3V_0^2-\left(3ac\sqrt{\delta_1}+\sqrt{bD_v}c^2\tanh\left(\sqrt{\frac{b}{D_v}}l\right)\right)V_0+\sqrt{\frac{D_v}{b}}a^2c^2\sqrt{\delta_1}\tanh\left(\sqrt{\frac{b}{D_v}}l\right)=0.\tag{53}\]

In the limit of \(\delta_1\to 0\) and \(a\ll 1\), we get two asymptotic roots of (53 )

\[\label{sec295V095roots95a95small} V_{0+}\sim \frac{\sqrt{bD_v}c^2}{3}\tanh\left(\sqrt{\frac{b}{D_v}}l\right), \quad V_{0-}\sim \frac{\sqrt{\delta_1}a^2}{b}.\tag{54}\]

In this way, for \(a\ll 1\) our asymptotic result for \(v(0) = V_0/\sqrt{\delta_1}\), after using the scaling relation (6 ), is that

\[\label{sec295v095a95small} v(0)_{+}\sim \frac{\sqrt{bD_v}c^2}{3\sqrt{\delta_1}}\tanh\left(\sqrt{\frac{b}{D_v}}l\right), \quad v(0)_{-}\sim \frac{a^2}{b}.\tag{55}\]

In contrast, for \(a = \mathcal{O}(1)\), we predict that \(v(0)\sim \frac{V_0}{\sqrt{\delta_1}}\), where \(V_0\) must be determined from the coupled nonlinear algebraic system (42 ) and (44 ). This highlights the qualitative difference between the small \(a\) and \(\mathcal{O}(1)\) regimes: in the former case a closed-form approximation is available, whereas in the latter case \(V_0\) must be computed numerically, for example using Newton’s method.

Figure 7 compares the asymptotic predictions with full numerical simulations of (1 ). In particular, we plot the the spike amplitude \(V_0\) by the simple closed-form result in terms of \(D_v\) for various values of \(a\). Although (54 ) and (55 ) were derived under the assumption \(a\ll 1\). The good agreements observed in Figure 7 show that the small-\(a\) asymptotics remain accurate even for moderately small \(a\).

2.5 Conditions for the absence of nucleation instability↩︎

Our analyses in Section 2.2 and 2.3 have shown that, on certain parameter ranges, spike nucleation behavior can occur as the diffusivity \(D_v\) decreases for the GM model in the limit \(\delta_1 \ll 1\). A key mechanism underlying this spike nucleation behavior was that there is a saddle-node bifurcation point at some finite critical value of \(D_{nuc}\), beyond which the outer solution ceases to exist.

One key condition for the existence of this saddle node point is that the parameters in the GM model (1 ) are such that there is no spatially homogeneous steady-state solution in the range where the outer problem is well defined. More specifically, as shown in the derivation in Lemmas 1, we require that in (32 and (36 ) \[R(u) < 0 \quad \text{and} \quad R'(u) < 0,\] in the range \(a<u<2a\) where the outer problem is well-posed. Under this condition, the integral \(\chi(\mu)\) defined from (46 ) increases monotonically in \(\mu\) but has finite values as \(\mu\to 2a\) below. This limiting value, based on the nonexistence of the outer asymptotic solution, was used in (?? ) to obtain a leading order prediction for the critical threshold of \(D_{nuc}\) where a saddle-node point must occur along the single-spike solution branch. The existence of such a saddle-node point is the signature of the onset of spike nucleation behavior.

On the other hand, if the GM system (1 ) admits a spatially homogeneous steady state, denoted by \(u_{\infty}\), within the well-posedness of the outer problem, it follows that \(R(u_{\infty})=0\), with \(R'(u_{\infty})>0\). This yields \[u_\infty= a+bc < 2a,\] and the saddle-node bifurcation for spike nucleation does not occur as the threshold \(2a\) is not reached.

In Figure 8, we show a global bifurcation diagram obtained from path-following a single-spike steady-state solution of the GM model (1 ) using pde2path [20] for the parameters \(a = 1.5, b = 1, c=1, \delta_1=0.01^2, D_v = 1.\) Since \(bc <a\), and \(\mu\) approaches to the spatially homogeneous steady state \(u_\infty=a+b=2.5\) but never reaches the limit \(2a=3\). Therefore, we observe that no spike nucleation events will occur.

In the meantime, we compute the single-spike steady state by modifying the calculation in Section 2.4 and replacing the outer approximation \(u\sim a+O(a^2)\) by \(u\sim a+bc\). this leads to the solution for \(v(x)\) as \[\label{sec95295v95no95nuc} v(x)=\frac{(a+bc)^2}{b}+\left(\frac{V_0}{\sqrt{\delta_1}}-\frac{(a+bc)^2}{b}\right)\frac{\cosh\left(\sqrt{b}(l-|x|)/\sqrt{D_v}\right)}{\cosh(\sqrt{b}l/\sqrt{D_v})}.\tag{56}\] Then by matching the condition (34 ) as in Section 2.4, we obtain two asymptotical roots for \(V_0\) that satisfy \[\label{sec295V095quadratic95nonuc} 3V_0^2-\left(3ac\sqrt{\delta_1}+\sqrt{bD_v}c^2\tanh\left(\sqrt{\frac{b}{D_v}}l\right)\right)V_0+\sqrt{\frac{D_v}{b}}(a+bc)^2c^2\sqrt{\delta_1}\tanh\left(\sqrt{\frac{b}{D_v}}l\right)=0.\tag{57}\] In Figure 9, we compare the asymptote solution with the full numerical simulation of model 1 obtained using Auto. The two asymptotic values, \(v_{0+}=1.47\) and \(v_{0-}=0.69\) are obtained from \(v_{0\pm}=\frac{V_{0\pm}}{\sqrt{\delta_1}}\) by solving the two roots of (56 ). In the following section, e will show that the upper branch with \(v=v_{0+}\) is stable, while the lower branch with \(v=v_{0-}\) is unstable, and the asymptotic predictions in Figure 9 align well with both the stable and unstable branches.

3 Stability of spike equilibrium: Large eigenvalues↩︎

In this section we study the linear stability of the one-spike equilibrium obtained in the previous section. We will consider the “extended" system obtained by adding two time-scaling parameters, \(\theta\) and \(\tau\), into system (1 ), which takes the form \[\tag{58} \begin{alignat}{2} \frac{\partial u}{\partial t} & = a-u+\frac{u^3}{w v}+\delta_1 \frac{\partial^2 u}{\partial x^2}, \,\, x \in (-l,l), t>0, \tag{59}\\ \theta \frac{\partial v}{\partial t} & = u^2-b v+D_v\frac{\partial^2 v}{\partial x^2},\,\, x\in (-l,l), t>0, \tag{60}\\ \tau \frac{\partial w}{\partial t} & = u -cw +\delta_2 \frac{\partial^2 w}{\partial x^2}, \,\, x \in (-l,l), t>0, \tag{61}\\ \intertext{with boundary conditions} & \left. u_x\right|_{\pm l} = \left. v_x\right|_{\pm l} = \left. w_x\right|_{x=\pm l} =0. \nonumber \end{alignat}\]

We are particularly interested in the new dynamics introduced by the additional variable \(w\), in contrast to previous studies of spike dynamics in two-component RD systems [16], where only one type of spike oscillation was observed. In this three-component setting, we demonstrate that two distinct types of spike oscillation can arise. The first corresponds to oscillations in the spike amplitude, associated with the destabilization of a “large” eigenvalue, and occurs when \(\theta\) increases beyond a critical threshold. The second corresponds to oscillations in the spike position, triggered when a “small” eigenvalue becomes unstable. This latter instability is induced by the slower timescale introduced through the variable \(w\), and emerges as \(\tau\) becomes sufficiently large. These results reveal a richer variety of instabilities than in the classical case of two components, leading to more complex dynamics and even the possibility of chaotic behavior.

In this section, we carry out the large eigenvalue analysis, in which the localized spike equilibrium, denoted by \(u_e(x), v_e(x)\) and \(w_e(x)\), takes the form given in Sec. 2. Since \(v_e\) is available in closed form only in the regime \(a\ll 1\) as shown in Sec. 2.4, we have for arbitrary parameter values of \(a\) that \[\label{sec395ss95u44w} u_e\sim \frac{V_{0}}{c\sqrt{\delta_1}}\left(w_c\left(\frac{x}{\sqrt{\delta_1}}\right)+\gamma\right), w_e=\frac{u_e}{c}.\tag{62}\]

By introducing a perturbation around the steady state as \[\label{sec395perturbations} u(x,t)=u_e(x)+\phi(x)e^{\lambda t}, v(x,t)=v_e(x)+\psi(x)e^{\lambda t}, w(x,t)=w_e(x)+\xi(x)e^{\lambda t},\tag{63}\] we obtain the following linearized eigenvalue problem: \[\tag{64} \begin{align} \lambda \phi &= \delta_1 \phi_{xx} -\phi+ \frac{3u_e^2}{w_e v_e}\phi-\frac{u_e^3}{w_ev_e^2}\psi-\frac{u_e^3}{w_e^2v_e}\xi,\tag{65}\\ \theta \lambda \psi &=D_v \psi_{xx}-b\psi+2u_e\phi,\tag{66}\\ \tau \lambda\xi &=\delta_2\xi_{xx}-c\xi+\phi\tag{67}. \end{align}\] As \(\delta_2=\delta_1^2\ll 1\), we get from (67 ) that \[\label{sec395xi} \xi=\frac{\phi}{\tau \lambda+c},\tag{68}\] and then (64 ) reduces to \[\tag{69} \begin{align} \lambda \phi &= \delta_1 \phi_{xx} -\phi+ \left(2+\frac{\tau\lambda}{\tau \lambda+c}\right)\frac{cu_e}{v_e}\phi-c\frac{u_e^2}{v_e^2}\psi,\tag{70}\\ \theta \lambda \psi &=D_v \psi_{xx}-b\psi+2u_e\phi.\tag{71} \end{align}\] Introducing the inner variable \(x=\sqrt{\delta_1}y\), and using the fact that \(\frac{u_e}{v_e}=\frac{w_c+\gamma}{c}\) in the inner region, we obtain \(\psi\sim \psi_0\), where \(\psi_0\) is a constant to be determined. Then (70 ) reduces to \[\label{sec395inner95phi} \lambda \left(1-\frac{\tau}{c+\tau\lambda}\gamma\right) \phi = \phi_{yy} -(1-2\gamma)\phi+ \left(2+\frac{\tau\lambda}{c+\tau \lambda}\right)w_c\phi-\frac{(w_c+\gamma)^2}{c}\psi_0.\tag{72}\] To determine the value of \(\psi_0\), we consider the outer region, where \(\phi\) is localized and can be approximated as a delta function. Therefore, \(\psi\) satisfies \[\begin{align} \label{sec395outer95psi} &D_v\psi_{xx}-(b+\theta\lambda)\psi= A \delta(x;0), \quad \psi_x(\pm l)=0, \quad \text{where} \\ &A=-2\int_{0^-}^{0^+}u_e\phi dx\sim -\frac{2V_{0}}{c}\int_{-\infty}^{\infty}w_c\phi dy. \end{align}\tag{73}\] This implies that \[\label{sec395psi} \psi =-\frac{A}{b+\theta\lambda}G(x;0),\tag{74}\] where \(G(x;0)\) satisfies \[\label{sec395Green39s32eqn} \frac{D_v}{b+\theta \lambda}G_{xx}-G= -\delta(x;0), \quad G_x(\pm l)=0.\tag{75}\] Solving (75 ) yields \[\label{sec395Green} G(x;0)=\frac{\sqrt{b+\theta \lambda}}{2\sqrt{D_v}}\frac{\cosh\left(\sqrt{\frac{b+\theta\lambda}{D_v}}(l-|x|)\right)}{\sinh\left(\sqrt{\frac{b+\theta\lambda}{D_v}}l\right)}.\tag{76}\] By the matching condition \(\psi(0)=\psi_0\), we obtain \[\label{sec395psi0} \psi_0=\frac{V_{0}}{c\sqrt{D_v(b+\theta\lambda)}}\frac{1}{\tanh\left(\sqrt{\frac{b+\theta\lambda}{D_v}}l\right)}\int_{-\infty}^{\infty}w_c\phi dy.\tag{77}\]

Now substituting (77 ) into (72 ) leads to the following non-local eigenvalue problem (NLEP) for arbitrary parameter values of \(a\):

\[\label{sec395NLEP95full} \lambda \left(1-\frac{\tau}{c+\tau\lambda}\gamma\right) \phi = \phi_{yy} -(1-2\gamma)\phi+\left(2+\frac{\tau\lambda}{c+\tau\lambda}\right) w_c\phi -\frac{V_{0}}{c^2\sqrt{D_v}\sqrt{b+\theta\lambda}}\frac{(w_c+\gamma)^2}{\tanh\left(\sqrt{\frac{b+\theta\lambda}{D_v}}l\right)}\int_{-\infty}^{\infty}w_c\phi dy.\tag{78}\]

Since the NLEP (78 ) is complicated to analyze, in this paper we will focus on the regime at the limit \(a \ll 1\), so that \(\gamma\sim \frac{ac\sqrt{\delta_1}}{V_0}\ll 1\), and (78 ) reduced to \[\label{sec395NLEP} \lambda \phi = \phi_{yy} -\phi+\left(2+\frac{\tau\lambda}{c+\tau\lambda}\right) w_c\phi -\frac{V_{0}}{c^2\sqrt{D_v}\sqrt{b+\theta\lambda}}\frac{w_c^2}{\tanh\left(\sqrt{\frac{b+\theta\lambda}{D_v}}l\right)}\int_{-\infty}^{\infty}w_c\phi dy.\tag{79}\]

To analyze the large eigenvalues in the NLEP (79 ), we treat \(\tau\) and \(\theta\) as bifurcation parameters and investigate their influence on stability. Specifically, we consider the following three distinct cases: (1) \(\tau = \theta = 0;\) (2) \(\theta > 0, \tau = 0;\) (3) \(\theta = 0, \tau > 0\). The stability results in case (1) and (2) in terms of large eigenvalues exhibit a structure similar to that of the NLEP analyzed in [16]. While in case (3), they extend to a novel eigenvalue problem, which leads to new stability conditions not present in the previous study.

Moreover, in contrast to the results of [16], which show that both \(\theta\) and \(\tau\) can induce Hopf bifurcations in the spike amplitude through large eigenvalue instabilities, we find that the instability associated with \(\tau\) in the third component could not be observed in full numerical simulations. Instead, before the large eigenvalue loses stability, a different mechanism is triggered arising from small eigenvalue instabilities. This small-eigenvalue instability and its role in spike dynamics will be studied in detail in the next section.

3.1 Case 1: \(\theta=0, \tau=0\).↩︎

We begin with the simplest case by setting \(\theta=0, \tau=0\), so that equation (79 ) reduces to the following NLEP

\[\label{sec395NLEP95c1} \lambda \phi(y)=L_0{\phi}-\frac{w_c^2}{A}\int_{-\infty}^{\infty}w_c\phi dy, \quad \text{where} \quad \frac{1}{A}=\frac{V_{0}}{c^2\sqrt{bD_v}}\frac{1}{\tanh\left(\sqrt{\frac{b}{D_v}}l\right)},\tag{80}\] in which \(L_0\phi\equiv \phi_{yy}-\phi+2w_c\phi\). This well-known NLEP was first studied in [21]. It has the following basic result:

Now we determine the stability of the two branches of one-spike solution with amplitudes \(V_{0-}\) and \(V_{0+}\) shown in (54 ). We substitute \(V_{0\pm}\) into (79 ) to obtain \[A(V_{0-})=\frac{c^4b\sqrt{bD_v}}{a^2\sqrt{\delta_1}}.\\\] Since \(\delta_1\ll1\), it follows that \(A(V_{0-})\gg 6\). By Theorem 1 we conclude that the lower branch of the one-spike steady state with \(V_0=V_{0-}\) is unstable. In contrast, for the upper branch with \(V_0=V_{0+}\), we have \(A(V_{0+})=3<6\), so Theorem 1 implies that this branch is stable. This stability result also explains Figure 4 and Figure 7, where the full simulations agree with the stable branch \(V_{0+}\).

3.2 Case 2: \(\theta>0, \tau=0\).↩︎

When choosing \(\theta>0\) and \(\tau=0\), the nonlocal term of (79 ), with \(V_0=V_{0-}\) takes the form \[A(V_{0-})=\frac{bc^4\sqrt{D_v(b+\theta\lambda)}}{a^2\sqrt{\delta_1}}\frac{\tanh\left(l\sqrt{(\theta\lambda+b)/D_v}\right)}{\tanh\left(l\sqrt{b/D_v} \right)}\gg 1 >6, \quad \text{as} \quad \delta_1 \to 0.\] Therefore, according to Theorem 1, the lower branch of the one-spike steady state with \(V_0 = V_{0-}\) is unstable.

Now we focus on the stability of the branch of solutions with \(V_0=V_{0+}\), and we are interested to investigate how relaxing \(\theta\) can destabilize the one-spike equilibrium. With \(V_0=V_{0+}\), the NLEP (79 ) simplifies to \[\label{sec395NLEP95c2} \lambda \phi=L_0 \phi-w_c^2\frac{\int w_c \phi dy}{A(\lambda;\theta)}, \quad\text{where} \quad A(\lambda;\theta)=3\sqrt{1+\frac{\theta \lambda}{b}} \frac{\tanh\left(l\sqrt{(\theta\lambda+b)/D_v}\right)}{\tanh\left(l\sqrt{b/D_v} \right)}.\\\tag{81}\]

Equation (81 ) has a structure similar to the NLEP studied in [22] and [2023 royal A]. Here we use the same idea to derive the stability result of (81 ). First, we rewrite (81 ) in the following form \[(L_0-\lambda)\phi =w_c^2, \qquad \mathrm{where} \quad \int w_c \phi \: dy=A(\lambda;\theta)\,,\] or \[\label{sec395f} f(\lambda):=\int w_c(L_0-\lambda)^{-1}w_c^2 \: dy =A(\lambda;\theta).\tag{82}\] The global behavior of the same \(f(\lambda)\) was studied in [11], from which we obtain the following basic results:

The graph of \(f(\lambda)\) is shown in Figure 10 (a). We then study the stability of the solution branch with \(V_{0}=V_{0+}\). By introducing the rescaled parameter \(\hat{\theta}:=\frac{\theta}{b}\), the function \(A(\lambda;\theta)\) can be rewritten as

\[\label{sec395A2} A(\lambda)=3\sqrt{1+\hat{\theta}\lambda}\frac{\tanh\left(l\sqrt{\frac{b}{D_v}}\sqrt{1+\hat{\theta}\lambda}\right)}{\tanh\left(l\sqrt{\frac{b}{D_v}}\right)}. \\\tag{83}\]

When \(\hat{\theta}\) is sufficiently large, the system can be destabilized via a Hopf bifurcation. This result was first proved in [11]. Although there is no closed-form expression for the Hopf threshold \(\hat{\theta}_h\) at which \(Re(\lambda) = 0\), this critical value can be computed numerically by discretizing the NLEP (79 ) using finite differences. This result is illustrated in Figure 10 (b), where \(l, b\) and \(c\) are fixed, and the above method is applied to calculate \(\hat{\theta}_h\). Close agreement is observed between theoretical predictions and full numerical simulations, with errors of no more than \(4.5\%\). Moreover, we observe that as \(D_v\) decreases or \(l\) increases, (83 ) reduces to \(A(\lambda;\theta)\sim 3\sqrt{1+\hat{\theta} \lambda}\), in which case the Hopf threshold approaches the asymptotic value \(\hat{\theta}_h\to 2.7492\), as also shown in the figure.

3.3 Case 3: \(\theta=0, \tau>0\).↩︎

In the case where \(\theta=0, \tau>0\), upon substituting \(V_0=V_{0+}\), NLEP (79 ) reduces to \[\label{sec395NLEP95c3} \lambda \phi= L_\lambda \phi- \frac{1}{3}w_c^2\int w_c \phi dy,\\\tag{84}\] where \[\label{sec395L95lambda} L_\lambda \phi= \phi_{yy}-\phi+\left(2+\frac{\tau\lambda}{c+\tau\lambda}\right)w_c\phi.\tag{85}\] Note that in this case, we obtain a new operator \(L_\lambda\), with nonlinear dependence on \(\lambda\) and the parameters \(\tau\) and \(c\). In contrast to \(L_0\) defined in Sec. 3.1, which has a single positive eigenvalue at \(\lambda=\frac{5}{4}\), and corresponds to the singularity of \(f(\lambda)\) given in (82 ). Here we investigate the operator \(L_\lambda\), by solving the nonlinear local eigenvalue problem \[\label{sec395L95lambda95eig} L_\lambda\phi=\lambda \phi.\tag{86}\] The problem can be rewritten as a Pöschl-Teller equation, whose solutions are well studied, and we summarize the result as follows; further details are given in the Appendix.

Now we investigate the nonlocal eigenvalue problem (84 ) by rewriting it in the following form, \[(L_\lambda-\lambda)\phi =w_c^2, \qquad \mathrm{where} \quad \int w_c \phi \: dy=3\,,\] or \[\label{sec395g} g(\lambda;c, \tau):=\int w_c(L_\lambda-\lambda)^{-1}w_c^2 \: dy =3.\tag{87}\] As \(L_\lambda\) has a unique positive eigenvalue \(\lambda=\lambda_0\), which can be solved via (?? ), the global behavior of \(g(\lambda; c, \tau)\) has a structure similar to \(f(\lambda)\), but with singularity now depends on \(\tau\) and \(c\). In Figure 11 we plot \(g(\lambda;c, \tau)\) for fixed \(c=0.5\). The right panel shows the singularity curve as \(\tau\) increases, in contrast with theorem 1, where the singularity point is fixed. In particular, as \(\tau \to \infty\), the operator approaches \(L_\lambda \to \phi_{yy}-\phi+3w_c\phi\), which (see Appendix) has a unique positive eigenvalue \(\lambda_\infty=2.5643\). This value is shown as a horizontal asymptote in the figure.

When \({\tau}\) is sufficiently large, the single spike steady state loses stability through a Hopf bifurcation at \(\tau=\tau_{lh}\). Although no closed-form expression is available for the Hopf threshold \(\tau_{lh}\) at which \(Re(\lambda) = 0\), we apply similar methods as in Sec. 3.2 and compute this critical value numerically by discretizing the NLEP (84 ) using finite differences. The predicted threshold shows excellent agreement with the full time-dependent simulations, as shown in Figure 12 (a). In the panel (b), we plot the time series of the spike amplitude at the boundary for \(\tau\) just above the threshold \(\tau_{lh}\). We also notice that the result \(\tau_{lh}\) appears linearly on the parameter \(c.\)

Note that in the full time-dependent simulations to verify the large-eigenvalue Hopf threshold \(\tau_{lh}\) for a single interior spike, simply increasing \(\tau\)on the symmetric domain \([-l,l]\) does not show the expected oscillations in spike amplitude (as would be associated with a large-eigenvalue Hopf mode). Instead, we observe oscillations in spike position, indicating that a small-eigenvalue (translational) instability intervenes first. This effect will be analyzed in the next section. Here, to suppress this effect and isolate the large-eigenvalue mechanism, we simulate on the half-domain \([0,l]\) with Neumann boundary conditions, so that a boundary spike at \(x=0\) represents the even extension of an interior spike and the odd translational mode is removed. In this setting, time-dependent PDE simulations performed in FlexPDE yield the threshold shown in 12. Moreover, a color plot of the spatiotemporal evolution of the half-spike amplitude as \(\tau\) passes \(\tau_{lh}\) is shown in Figure 13.

4 Stability of spike equilibrium with \(a\ll 1\): Small eigenvalues↩︎

In this section, we investigate the stability of the one-spike solution with respect to small eigenvalues. We will show that the parameter \(\tau\) can trigger a small-eigenvalue instability, leading to a Hopf bifurcation associated with oscillatory spike motion. Numerical evidence suggests that this bifurcation is supercritical, giving rise to stable, time-periodic spike dynamics. For simplicity, we will focus on the regime where \(\theta=0, \tau>0\) and begin our analysis from the eigenvalue problem (64 ).

To distinguish the notation from that used in the large-eigenvalue analysis, in terms of \(y=\frac{x}{\sqrt{\delta_1}}\), we introduce the following inner variables for the one-spike equilibrium and perturbations: \[\label{sec495small95eig95var} \begin{align} u_e(x)&=U_e(y), v_e(x)=V_e(y), w_e(x)=W_e(y), \\ \phi(x)&=\Phi(y), \psi(x)=\Psi(y), \xi(x)=\Xi(y). \end{align}\tag{88}\] Then (64 ) becomes \[\tag{89} \begin{align} \lambda \Phi &= \Phi_{yy} -\Phi+ \frac{3U_e^2}{W_e V_e}\Phi-\frac{U_e^3}{W_eV_e^2}\Psi-\frac{U_e^3}{W_e^2V_e}\Xi,\tag{90}\\ 0&=D_v \Psi_{yy}-\delta_1 b\Psi+2U_e\Phi,\tag{91}\\ \tau \lambda\Xi &=\frac{\delta_2}{\delta_1}\Xi_{yy}-c\Xi+\Phi\tag{92}. \end{align}\] Since \(\frac{\delta_2}{\delta_1}\ll 1\), we get from (92 ) that \(\Xi=\frac{\Phi}{c+\tau\lambda}\), then (89 ) reduces to \[\tag{93} \begin{align} \lambda \Phi &= \Phi_{yy} -\Phi+ \frac{3U_e^2}{W_e V_e}\Phi-\frac{U_e^3}{W_eV_e^2}\Psi-\frac{U_e^3}{W_e^2V_e}\frac{\Phi}{c+\tau \lambda},\tag{94}\\ 0&=D_v \Psi_{yy}-\delta_1 b\Psi+2\delta_1 U_e\Phi,\tag{95} \end{align}\] We expand \[\label{sec495expand1} \lambda = \delta_1 \lambda+0 + \cdots,\quad \Phi(y)=\Phi_0(y)+\sqrt{\delta_1}\Phi_1(y)+ \cdots, \Psi(y)=\sqrt{\delta_1}\Psi_0(y)+ \cdots,\tag{96}\] and \[\label{sec495expand2} U_e(y)=U_{e0}+\sqrt{\delta_1} U_{e1}+\cdots, V_e(y)=V_{e0}+\sqrt{\delta_1} V_{e1}+\cdots, W_e(y)=W_{e0}+\sqrt{\delta_1} W_{e1}+\cdots.\tag{97}\] Substituting (96 ) and (97 ) into (93 ) and collecting \(\mathcal{O}(1)\) terms, we get \(\Phi_0\) satisfies \[\label{sec495Phi95eqn} \Phi_{0yy}-\Phi_0+2w_c\Phi_0=0.\tag{98}\] The solution to (98 ) is \[\label{sec495Phi0} \Phi_0=w_{cy}, \quad \text{where} \quad w_c=\frac{3}{2}\text{sech}^2\left(\frac{y}{2}\right).\tag{99}\]

Now collecting the next-order \(\mathcal{O}(\varepsilon)\) terms of (93 ), and using the fact \(\frac{U_{e0}}{V_{e0}}=\frac{w_c}{c}\), we obtain \[\tag{100} \begin{align} \Phi_{1yy} -\Phi_1+ 2w_c\Phi_1&= \frac{w_c^2}{c}\Phi_0+(k_2-k_1)\Phi_0,\tag{101}\\ \Psi_{0yy}&=-\frac{2V_0}{cD_v}w_c(y)w_{cy}(y)\tag{102}, \end{align}\] where \(k_1, k_2\) are the \(\mathcal{O}(\sqrt{\delta_1})\) coefficients of the two terms \(\frac{3U_e^2}{W_e^2 V_e}\) and \(\frac{U_e^3}{W_e V_e}\), separately. Recall the fact that \(\mathcal{L}_0w_c=w_c''-w_c+2w_c^2=w_c^2\), so \(\Phi_1\) could be written as \[\label{sec495Phi1} \Phi_1=\Phi_0(0)w_c+\Phi_{1,\text{odd}},\tag{103}\] where \(\Phi_{1,\text{odd}}\) is some odd function.

Next, multiply (93 ) by \(U_y\) and integrate over \((-\infty, \infty)\), we obtain by integration by parts that \[\label{sec495int} \lambda \int\Phi U_{ey} \left(1-w_c\frac{\tau}{c+\tau \lambda}\right) dy=\int \Phi \left( U_{eyyy}- U_{ey}+2w_c U_{ey}\right) dy-\int\frac{w_c^2}{c}\Psi U_{ey}dy.\tag{104}\] Here the shorthand notation \(\int f dy\) denotes \(\int_{-\infty}^{\infty} f(y) dy\). Moreover, using the fact that \(U_e\) satisfies (8 ), we reduce (104 ) to the following small eigenvalue problem \[\label{sec495int2} \lambda \int\Phi U_{ey} \left(1-\frac{\tau}{c+\tau\lambda}w_c\right) dy=\int \frac{w_c^2}{c}\Phi V_{ey} dy-\int\frac{w_c^2}{c}\Psi U_{ey}dy.\tag{105}\]

Now, to continue to reduce (105 ), we estimate the left-hand side of (105 ) as \[\begin{align} \label{sec495LHS1} \text{LHS}&=\lambda\frac{V_0}{c\sqrt{\delta_1}}\left(\int w_{cy}^2 dy-\frac{\tau}{c+\tau\lambda}\int w_cw_{cy}^2 dy\right), \end{align}\tag{106}\] and using the leading-order term \(\Phi\sim \Phi_0=w_{cy}\), the right-hand side of (105 ) becomes \[\begin{align} \label{sec495RHS1} \text{RHS}&=\frac{1}{c}\int w_c^2w_{cy}\left(V_{ey}-\frac{1}{c}\frac{V_0}{\sqrt{\delta_1}}\Psi\right) dy \nonumber\\ &=-\frac{1}{c}\int \frac{ w_c^3}{3}\left(V_{eyy}-\frac{1}{c}\frac{V_0}{\sqrt{\delta_1}}\Psi_y\right) dy, \end{align}\tag{107}\] where from (9 ) \(V_{eyy}\) can be expressed as \[\label{sec495v95eyy} V_{eyy}=\frac{\delta_1bV_e-\delta_1U_e^2}{D_v}.\tag{108}\] To estimate the term \(\int\frac{w_c^3}{3}\Psi_y\) in (107 ), we define \(F := \int_0^y\frac{w_c^3(s)}{3}ds\) and write \[\label{sec495rhsterm1} \int\frac{w_c^3}{3}\Psi_y=\int\Psi_ydF=\Psi_yF\vert_{-\infty}^\infty-\int F(y)\Psi_{yy}.\tag{109}\] Note that from (91 ), \(\Psi_{yy}\) can be approximated as \[\label{sec495Psi95yy} \Psi_{yy}\sim -2\frac{\sqrt{\delta_1V_0}}{cD_v}w_c\Phi_0\nonumber=-2\frac{\sqrt{\delta_1V_0}}{cD_v}w_cw_{cy},\tag{110}\] then applying integration by parts yields \[\begin{align} \label{sec495rhsterm2} \int F(y)\Psi_{yy} &=-2\int F(y)\frac{\sqrt{\delta_1}V_0}{cD_v}w_cw_{cy}\nonumber\\ &=\frac{\sqrt{\delta_1}V_0}{cD_v}\int \frac{w_c^5}{3}dy. \end{align}\tag{111}\] Substituting (111 ) into (109 ), we obtain \[\label{sec495rhsterm3} \;int\frac{w_c^3}{3}\Psi_y= \int\frac{w_c^3}{3}\left(\langle \Psi_y \rangle - \frac{\sqrt{\delta_1}V_0}{cD_v}w_c^2\right) dy,\tag{112}\] where \(\langle \Psi_y \rangle\) is defined as \(\langle \Psi_y \rangle := \frac{\Psi_y(\infty)+\Psi_y(-\infty)}{2}\). Now substituting (108 ) and (112 ) into (107 ), the right-hand side of (105 ) simplifies to \[\text{RHS}=-\frac{1}{c}\int \frac{w_c^3}{3}dy\left(\frac{\sqrt{\delta_1}b}{D_v}V_0-\frac{1}{c}\frac{V_0}{\sqrt{\delta_1}}\langle \Psi_y \rangle\right).\] Therefore, the small eigenvalue problem (105 ) reduces to \[\label{sec495int3} \lambda\left(\int w_{cy}^2 dy-\frac{\tau}{c+\tau\lambda}\int w_cw_{cy}^2 dy\right)=\int \frac{w_c^3}{3}dy\left(\frac{1}{c}\langle \Psi_y \rangle-\frac{\delta_1b}{D_v}\right),\tag{113}\] in which the only term to determine is \(\langle \Psi_y \rangle := \frac{\Psi_y(\infty)+\Psi_y(-\infty)}{2}\). To evaluate this term, we consider the outer problem (95 ) for \(\Psi\).

We have shown from (99 ) and (103 ) that \[\phi(x)=\Phi(y)\sim \Phi_0+\sqrt{\delta_1}\Phi_1=w_{cy}+\sqrt{\delta_1}\Psi_0(0)w_c+\sqrt{\delta_1}\Phi_{1,odd}.\] Note that \(w_{cy}\) is like a dipole (odd) and \(w_c\) behaves like a Dirac delta function, so taking \(\sqrt{\delta_1}\ll \delta\ll 1\) and integrating (95 ) on \((-\delta,\delta)\) yields the following condition \[\label{sec495outer95conditon1} \psi_x(\delta)- \psi_x(-\delta)\sim -\frac{2}{D_v} \int_{-\delta}^\delta U_e\phi dx\sim-\sqrt{\delta_1}\frac{2V_0}{cD_v}\Psi(0)\int w_c^2 dy.\tag{114}\]

To derive the second jump condition for \(\psi\), we multiply (95 ) by \(x\) and integrate over \((-\delta,\delta)\). Applying integration by parts yields \[\label{sec495outer95conditon2} \psi(\delta)-\psi(-\delta)=-\frac{\sqrt{\delta_1}V_0}{cD_v}\int w_c^2dy\tag{115}\] Therefore, the outer variable \(\psi(x)\) satisfies: \[\label{sec495psi95out} \begin{align} &D_v\psi_{xx}-b\psi=0, \quad x\neq 0, \\ &\psi(0^+)-\psi(0^-)=-\frac{\sqrt{\delta_1}V_0}{cD_v} \int w_c^2 dy,\\ &\psi_x(0^+)-\psi_x(0^-)=-\frac{2\sqrt{\delta_1}V_0}{cD_v}\Psi_0(0) \int w_c^2 dy. \end{align}\tag{116}\] Since (110 ) implies that \(\Psi_{yy}\) is an odd function, so \(\Psi=\Psi(0)+\text{odd function}\) and \[\Psi(0)=\frac{\Psi(\infty)+\Psi(-\infty)}{2}.\] By matching \(\psi(x)=\Psi(y)=\sqrt{\delta_1}\Psi_0(y)\) and introducing \(\psi(x)=\sqrt{\delta_1}\eta(x)\), we get \[\label{sec495eta95out} \begin{align} &D_v\eta_{xx}-b\eta=0, \quad x\neq 0, \\ &\eta(0^+)-\eta(0^-)=-\frac{V_{0}}{cD_v} \int w_c^2 dy,\\ &\eta_x(0^+)-\eta_x(0^-)=-\frac{2V_0}{cD_v}\langle \eta\rangle \int w_c^2 dy. \end{align}\tag{117}\] and the inner problem (113 ) written in terms of \(\eta\) as \[\label{sec495int4} \lambda\left(\int w_{cy}^2 dy-\frac{\tau}{c+\tau\lambda}\int w_cw_{cy}^2 dy\right)=\delta_1\int \frac{w_c^3}{3}dy\left(\frac{1}{c}\langle \eta_x \rangle-\frac{b}{D_v}\right),\tag{118}\] where \(\langle \eta_x \rangle:=\frac{\eta_x(0^-)+\eta_x(0^+)}{2}\). Now to solve \(\eta\) from (118 ) we specify the boundary conditions of \(\eta\) for single spike equilibrium as \(\eta'(\pm L)=0\) and \(\eta\) is odd. Therefore, \(\eta_x(0^+)-\eta_x(0^-)=0\) and solving (117 ) yields \[\label{sec495eta} \eta(x)=-\frac{V_0}{cD_v}\frac{\int w_c^2 dy}{2\cosh\left(\sqrt{\frac{b}{D_v}}l\right)}\left\{ \begin{array} [c]{ll}\cosh\left(\sqrt{\frac{b}{D_v}}(x+L)\right), & -l<x<0 \\ -\cosh\left(\sqrt{\frac{b}{D_v}}(x-L)\right), & 0<x<l. \\ \end{array} \right.\tag{119}\] We then calculate \(\langle\eta_x \rangle=\frac{\eta_x(0^-)+\eta_x(0^+)}{2}=\frac{V_0}{2cD_v}\sqrt{\frac{b}{D_v}}\int w_c^2 dy\tanh\left(\sqrt{\frac{b}{D_v}}l\right)\). With \(V_0=V_{0+}=\frac{\sqrt{bD_v}c^2}{3}\tanh\left(\sqrt{\frac{b}{D_v}}l\right)\), we obtain \[\label{sec495eta95x} \langle\eta_x \rangle=\frac{bc}{6D_v}\int w_c^2 dy\tanh^2\left(\sqrt{\frac{b}{D_v}}l\right).\tag{120}\] Finally, plug (120 ) into the eigenvalue problem (118 ) and using the fact that \(\int w_c^2 dy=6\), we get \[\label{sec495int5} \lambda\left(\int w_{cy}^2 dy-\frac{\tau}{c+\tau\lambda}\int w_cw_{cy}^2 dy\right)=-\delta_1\int \frac{w_c^3}{3}dy\frac{b}{D_v}\text{sech}^2\left(\sqrt{\frac{b}{D_v}}l\right)\tag{121}\]

As \(\lambda=\mathcal{O}(\delta_1)\ll 1\), we approximate \(\frac{\tau}{c+\tau \lambda}\sim \frac{\tau}{c} +\mathcal{O}(\delta_1)\) to get the leading order of \(\lambda\) \[\label{sec495small95eig1} \lambda\sim\frac{-\delta_1\int \frac{w_c^3}{3}dy\frac{b}{D_v}\text{sech}^2\left(\sqrt{\frac{b}{D_v}}l\right)}{\int w_{cy}^2 dy-\frac{\tau}{c}\int w_cw_{cy}^2 dy}.\tag{122}\] It is obvious that the top of (122 ) is always negative, therefore, using the fact that \(\int w_{cy}^2 dy=\frac{6}{5}\) and \(\int w_cw_{cy}^2 dy=\frac{36}{35}\), \(\lambda\) crosses \(0\) whenever \[\label{sec495tau95h} \tau>\tau_h=\frac{c\int w_{cy}^2 dy}{\int w_cw_{cy}^2 dy}=\frac{7}{6}c.\tag{123}\]

To obtain the full expression of \(\lambda\), we rewrite (121 ) as the following quadratic equation \[\label{sec495quadratic95lambda} \frac{7}{6}\tau\lambda^2-\left(\tau-\frac{7}{6}c-\frac{35}{36}\delta_1\tau k\right)\lambda+\frac{35}{36}\delta_1kc=0,\tag{124}\] where \(k=\int \frac{w_c^3}{3}dy\frac{b}{D_v}\text{sech}^2\left(\sqrt{\frac{b}{D_v}}l\right)\).

This implies that as \(\tau\) increases over \(\tau_h\), a pair of complex conjugate eigenvalues enter the unstable right half-plane, triggering an oscillatory instability in the motion of the spike, with \[\label{sec495real95eig} Re(\lambda_{\pm})=\frac{\tau-\frac{7}{6}c-\frac{35}{36}\delta_1\tau k}{\frac{7}{3}\tau}\sim\frac{3}{7\tau}\left(\tau-\frac{7}{6}c\right),\tag{125}\] and \[\label{sec495imag95eig} Im(\lambda_{\pm})=\pm \sqrt{\delta_1\frac{5}{6}\frac{kc}{\tau}}i.\tag{126}\] We now summarize the discussion in the following.

In Figure 14, we compare asymptotic predictions \(\tau_h\) with full numerical simulations of (58 ).In the simulations, we record the critical value of \(\tau\), beyond which the spike begins to oscillate periodically around the center of the domain. In terms of \(c\), \(\tau_h\) varies linearly and shows excellent agreement between analysis and simulations. Moreover, as shown in the Introduction 1, Figure 3 presents full numerical simulations with FlexPDE illustrating the spike dynamics for values of \(\tau\) exceeding the Hopf threshold \(\tau_h\). For \(\tau\) slightly larger than \(\tau_h\), the interior spike undergoes small-amplitude oscillations around its equilibrium location, which implies the onset of Hopf instability. As \(\tau\) is increased further, these oscillations grow in amplitude, and the spike begins to drift significantly. In this regime, the motion eventually drives the spike toward the boundary, indicating a transition from oscillatory dynamics to drift-dominated motion, which is consistent with the decrease of the imaginary part of the eigenvalue as \(\tau\) increases.

5 Discussion↩︎

In this paper, we present an extension of the GM model (1 ) in the semi-strong interaction regime, also characterized by an asymptotically large diffusivity ratio. Within this framework, we have constructed a single spike equilibrium for an arbitrary value of \(a>0\). From a mathematical perspective, the novelty of the analysis in contrast to previous studies [22], [23] is that we now must couple a nonlinear inner problem for the spike profile to a nonlinear reduced scalar boundary value problem (BVP) defined in the outer region away from the spike. This coupling leads to a more complex inner-outer interaction, but it applies to all values of \(a>0\).

For the non-trivial background state \(a>0\), we analyzed a global bifurcation mechanism that is responsible for the generation of spatial patterns as the inhibitor diffusivity \(D_v\) decreases. In parameter regimes where a one-spike solution exists on the infinite line, in Section 2.5 we showed that spike nucleation will not occur as \(D_v\) decreases in the semi-strong interaction regime.

By introducing the time-scaling parameter \(\theta\) and \(\tau\) in front of the equations for \(v\) and \(w\), respectively, we studied the novel behavior introduced by the third component \(w\). The dynamics now exhibits not only large-scale oscillatory motion in the amplitude, which is triggered by large-eigenvalue instabilities, but also oscillatory spike motion associated with small-eigenvalues crossing into the right half-plane. These two mechanisms highlight a key difference from classical two-component RD systems [11], [22], [24], as well as some three-component framework [16], where only large-scale oscillatory motion in the amplitude is observed. As a result, the extended model supports a richer variety of oscillatory dynamics.

There are numerous open questions for future study. In Section 3 we have derived a novel nonlocal eigenvalue problem (NLEP) due to the presence of the \(\tau\)-dependent term, it would be interesting to provide a rigorous study for the spectrum of NLEP. Previous work on NLEPs with eigenvalue dependence has focused primarily on cases where the eigenvalue enters the non-local term rather than the operator itself [11]. A systematic investigation of the spectrum in the present three-component setting could reveal new bifurcation structures.

In this paper, we have studied the effects of \(\theta\) and \(\tau\) separately under the regime \(a\ll 1\). It would be natural to extend the stability analysis to the case \(a=\mathcal{O}(1)\). Preliminary simulations suggest that this setting gives rise to a rich variety of spike dynamics that deserve further study, including spike motion, spike nucleation, and spike competition leading to spike death. It would therefore be interesting to investigate the detailed spike dynamics in this regime, as well as the possible interplay with multi-spike patterns.

An example is shown in Figure 15 for \(a>0\), and \(\tau\) is sufficiently large. In the left panel, the spike motion leads to the nucleation as the spike radius near the far boundary becomes large enough. However, as the spike oscillates back, this motion induces competition: the boundary spike is annihilated, while a new spike is nucleated on the opposite side. In the right panel, when \(\tau\) is increased further, the interior spike no longer oscillates but instead drifts toward the boundary, and triggers nucleation at the opposite boundary. This newly generated boundary spike then transitions into an interior spike, which in turn annihilates the previously existing interior spike once it becomes a boundary.

Acknowledgments↩︎

Chunyi Gai gratefully acknowledges the support of the NSERC Discovery Grant Program. The work of Fahad Al Saadi is funded by the Ministry of Higher Education, Research, and Innovation (MoHERI) under the Block Funding Program and conducted with support from the Military Technological College (MTC), Oman.

Appendix: The Local Eigenvalue Problem \(L_\lambda\phi=\lambda\phi\)↩︎

In this appendix, we solve the local eigenvalue problem (86 ) as follows. \[\label{app95local95problem} L_\lambda \phi:= \phi_{yy}-\phi+\left(2+\frac{\tau\lambda}{c+\tau\lambda}\right)w_c\phi=\lambda \phi,\tag{127}\] where \(w_c=\frac{3}{2}\text{sech}^2(\frac{y}{2}).\) We are interested in finding any positive solution \(\lambda>0\). Let \(z=\frac{x}{2}\), then (127 ) becomes \[\label{app95local95scaled} \phi_{zz}+\left[\left(12+6\frac{\tau\lambda}{c+\tau\lambda}\right)\text{sech}^2(z)-4(1+\lambda)\right]\phi=0,\tag{128}\] which is of the well-known Pöschl Teller type, and has known exact eigenvalues[25]. In particular, in the canonical form \[\label{app95Poschl--Teller32} \phi_{zz}+\left(p(p+1)\text{sech}^2(\frac{y}{2})-k^2\right)\phi=0, \quad v>0.\tag{129}\] The parameter \(p\) controls the strength of the potential bump, and there are \([p]\) localized eigenfunctions (here \([p]\) denotes the floor of \(p\)), which can be expressed in terms of associated Legendre or hypergeometric functions. Their discrete levels are \[\label{app95kn} k_n^2=(v-n)^2, \quad n=0,1,..[v]-1.\tag{130}\] Now we match the parameters in our local eigenvalue problem (127 ), which yields the following system. \[\Psi(0)=\frac{\Psi(\infty)+\Psi(-\infty)}{2}.\] By matching \(\psi(x)=\Psi(y)=\sqrt{\delta_1}\Psi_0(y)\) and introducing \(\psi(x)=\sqrt{\delta_1}\eta(x)\), we get \[\label{app95match} \begin{align} v(v+1)&=12+6\frac{\tau\lambda}{c+\tau\lambda}, \\ k^2&=4(1+\lambda),\\ k_n^2&=(v-n)^2. \end{align}\tag{131}\] Solving the system (131 ), we find that \(\lambda\) satisfies \[\label{app95lambda95eqn} F(\lambda):= \frac{\left(p(\lambda)-n\right)^2}{4}-1-\lambda=0,\quad n=0,1,..[p]-1,\tag{132}\] where \[\label{app95p40lambda41} p(\lambda)=\frac{-1+\sqrt{1+4\left(12+6\frac{\tau\lambda}{\tau\lambda+c}\right)}}{2}.\tag{133}\] 132 can be further reduced to the following equation, \[\label{app95lambda95eqn2} 12+6\frac{\tau\lambda}{\tau\lambda+c}=n^2+n+4+4\lambda+\sqrt{n^4-n^2+(8+8\lambda)\left(2n^2+2n+\frac{1}{2}\right)}.\tag{134}\] As we are only interested in positive solutions and \(\lambda\) decreases with the mode index \(n\). In fact, for \(n=1\), (134 ) reduces to \[1+\frac{\tau\lambda}{c+\tau\lambda}-\frac{2}{3}\lambda-\sqrt{1+\lambda}=0,\] which always admits the solution \(\lambda=0\) for an arbitrary value of \(c\) and \(\tau\), so we seek a positive solution for the lowest mode \(n=0\), which satisfies \[\label{app95thm3} 4+3\frac{\tau\lambda}{c+\tau\lambda}-2\lambda-\sqrt{1+\lambda}=0.\tag{135}\] Moreover, as \(\tau\to \infty\), equation (135 ) reduces to \[7-2\lambda-\sqrt{1+\lambda}=0,\] which has a solution \(\lambda\sim 2.56\).

References↩︎