Multi-Spike Patterns in the Gierer-Meinhardt System with a Non-Zero Activator Boundary Flux
MMULTI-SPIKE PATTERNS IN THE GIERER-MEINHARDT SYSTEM WITH ANON-ZERO ACTIVATOR BOUNDARY FLUX
DANIEL GOMEZ AND JUNCHENG WEI
Abstract.
The structure, linear stability, and dynamics of localized solutions to singularly per-turbed reaction-diffusion equations has been the focus of numerous rigorous, asymptotic, and nu-merical studies in the last few decades. However, with a few exceptions, these studies have oftenassumed homogeneous boundary conditions. Motivated by the recent focus on the analysis of bulk-surface coupled problems we consider the effect of inhomogeneous Neumann boundary conditionsfor the activator in the singularly perturbed one-dimensional Gierer-Meinhardt reaction-diffusionsystem. We show that these boundary conditions necessitate the formation of spikes that concen-trate in a boundary layer near the domain boundaries. Using the method of matched asymptoticexpansions we construct boundary layer spikes and derive a new class of shifted Nonlocal EigenvalueProblems for which we rigorously prove partial stability results. Moreover by using a combinationof asymptotic, rigorous, and numerical methods we investigate the structure and linear stability ofselected one- and two-spike patterns. In particular we find that inhomogeneous Neumann boundaryconditions increase both the range of parameter values over which asymmetric two-spike patternsexist and are stable. Introduction
Originally proposed by Gierer and Meinhardt in 1972, the Gierer-Meinhardt (GM) model illus-trates the pattern forming potential of two mechanisms: local self-activation and lateral inhibition[4]. In the original GM model two diffusing chemical species, an activator and inhibitor , inter-act through prescribed reaction-kinetics. Specifically, letting u and v denote the concentrationsof the activator and inhibitor respectively, the GM model takes the form of the two-componentreaction-diffusion-system u t = D u ∆ u − u + u p v − q , v t = D v ∆ v − v + u r v − s , x ∈ Ω , (1.1a)where it is assumed that the exponents ( p, q, m, s ) satisfy p > , q > , r > , s ≥ , and p − − ( s + 1) − qr < , (1.1b)and homogeneous Neumann boundary conditions are typically imposed. As is evident from thechoice of reaction kinetics, the GM model describes a process of self activation, or auto-catalysis, bythe activator and lateral inhibition by the inhibitor. In addition, the spatially homogeneous steadystate ( u, v ) = (1 ,
1) undergoes a spontaneous symmetry breaking bifurcation into a patterned statethrough a Turing instability when D u /D v is sufficiently small. While one of the early applicationsof the GM model was specifically to head generation and regeneration in Hydra [4], self-activationand lateral inhibition is considered to be an important pattern formation mechanism in generalbiological systems and this has been supported by the identification of several molecular candidatesfor the activator and inhibitor in activator-inhibitor systems [12].In the singularly perturbed limit of an asymptotically small diffusivity ratio D u /D v = O ( ε ) (cid:28) Date : August 12, 2020.
Key words and phrases.
Gierer-Meinhardt system, singular perturbation, matched asymptotic , nonlocal eigenvalueproblem (NLEP), spiky solution.J. Wei is partially supported by NSERC of Canada. D. Gomez is supported by NSERC of Canada. a r X i v : . [ n li n . PS ] A ug D. GOMEZ AND J. WEI spatial modes. In this singularly perturbed limit the GM system, as well as numerous other two-component reaction-diffusion systems, are known to exhibits multi-spike patterns in which theactivator concentrates at a discrete collection of asymptotically small regions. The rich structuraland stability properties of these multi-spike patterns have been the focus of extensive asymptotic,rigorous, and numerical studies over the past two decades [18, 6, 3, 20, 16]. While boundaryconditions have been identified as playing an important role in pattern formation [10, 2], relativelyfew studies have investigated the role of boundary conditions on the structure and stability of multi-spike solutions to singularly perturbed reaction diffusion systems. Instead most such studies haveassumed either homogeneous Neumann or homogeneous Dirichlet boundary conditions. Notableexceptions include the investigation of homogeneous Robin boundary conditions for the activator inthe GM model [1, 11] and inhomogeneous Robin boundary conditions for the inhibitor in the two-dimensional Brusselator model [15]. These two studies and their illustration of the effect of boundaryconditions on the structure and stability of multi-spike patterns serve as the primary motivationfor the present paper in which we consider inhomogeneous Neumann boundary conditions for theactivator in the one-dimensional singularly perturbed GM model. Additionally, this paper aims toaddress some of the technical issues that arise in bulk-surface coupled reaction diffusion systems forwhich inhomogeneous boundary conditions naturally arise [8, 14, 9, 5].To simplify the presentation in this paper we consider the singularly perturbed one-dimensionalGM model with exponents ( p, q, r, s ) = (2 , , ,
0) though we remark that different exponents resultonly in quantitative rather than qualitative differences. After possibly rescaling the spatial variablein (1.1a) we assume that the domain is the unit interval Ω = (0 ,
1) and that D u = ε while D v = O (1) where ε (cid:28) O ( ε ) length and by integrating the inhibitor equationin (1.1a) it is easy to see that if v = O (1) then we must have u = O ( ε − / ) in each interval onwhich it is concentrated. This motivates our choice of rescaling u = ε − ˜ u and v = ε − ˜ v which whensubstituted into (1.1a) and dropping the tildes gives u t = ε u xx − u + u v − , < x < , (1.2a) τ v t = Dv xx − v + ε − u , < x < , (1.2b)and for which we observe both u and v will be O (1) in each interval where u is concentrated.In addition we impose inhomogeneous and homogeneous Neumann boundary conditions for theactivator and inhibitor respectively which are given by − εu x (0) = A, εu x (1) = B, v x (0) = 0 , v x (1) = 0 , (1.2c)where we assume that A, B ≥ not have a spatially homogeneous steady state when A >
B >
A > B = 0. As we demonstrate in § shifted spike (see Figure 1a) that is analogous to the near-boundary spike solution of [1, 11] though it has different stability properties which we investigatein §
4. Furthermore, by considering examples of one- and two-spike patterns in § A = B = 0 [17], the non-zero boundary fluxes x One-Spike Solutions for A B = 0 A (a) x t u , v u ( x , t ) v ( x , t ) (b) Figure 1. (A) Examples of shifted one-spike solution concentrated at x = 0 forvarious values of A ≥ ε = 0 .
05 and D = 5. (B) Evolution of solutionto GM problem with D = 0 . ε = . τ = 0 .
1, and A = B = 0 .
08. The initialcondition is an unstable two-spike equilibrium where both spikes concentrate at theboundaries. A competition instability predicted by our asymptotic results in Figure5a is triggered and leads to the solution settling at an asymmetric pattern.force the solution to settle to an asymmetric pattern. This illustrates that one of the key features ofintroducing non-zero boundary fluxes is that it leads to a kind of robustness of asymmetric solutionssimilar to that observed in the presence of an inhomogeneous precursor gradient [7].The remainder of this paper is organized as follows. In § § O (1) time scale. Then, in § shifted nonlocal eigenvalue problems analogous to those studied in [11]. In § § Quasi-Equilibrium Multi-Spike Solutions and their Slow Dynamics
In this section we use the method of matched asymptotic expansions to derive an algebraic systemand an ordinary differential equation that determine the profile and slow dynamics of a multi-spikequasi-equilibrium solution to (1.2). The derivation uses techniques that are now common in thestudy of localized patterns in one dimension. Our presentation will therefore be brief, highlightingonly the novel aspects introduced by the inhomogeneous Neumann boundary conditions for theactivator. We begin by supposing that there are two spikes concentrated at the boundaries x L = 0and x R = 1 as well as N spikes concentrated in the interior at 0 < x < ... < x N <
1. In additionwe assume that the spikes are well separated in the sense that | x i − x j | = O (1) as ε → + for all i (cid:54) = j ∈ { L, , ..., N, R } . This last assumption is key for effectively applying the method of matchedasymptotic expansions. D. GOMEZ AND J. WEI
We first construct an asymptotic approximation for the solution near x = 0 by letting x = εy where y = O (1) and expanding u ∼ u L, ( y ) + O ( ε ) , v ∼ v L, ( y ) + εv L, ( y ) + O ( ε ) . (2.1)It is easy to see that v L = ξ L where ξ L is is an undetermined constant, and u L, ( y ) = ξ L w c ( y + y L )where w c ( y ) is the unique homoclinic solution to w (cid:48)(cid:48) c − w c + w c = 0 , < y < ∞ , w (cid:48) c (0) = 0 , w c ( y ) → y → ∞ , (2.2a)given explicitly by w c ( y ) = 32 sech y . (2.2b)Moreover, the undetermined shift parameter y L is chosen to satisfy the inhomogeneous Neumannboundary condition w (cid:48) c ( y L ) = − A/ξ L . (2.3)The unknown constants ξ L and y L are found by matching with the outer solution. To determine theappropriate Neumann boundary conditions for the outer problem we must first calculate v (cid:48) L, ( y ) as y → ∞ . This is done by integrating the O ( ε ) equation Dv (cid:48)(cid:48) L, = − ξ L w c ( y + y L ) , < y < ∞ ; v (cid:48) L, (0) = 0 . (2.4)over 0 < y < ∞ to obtain the limit lim y → + ∞ v (cid:48) L, ( y ) = − ξ L D η ( y L ) , (2.5)where η ( y ) ≡ (cid:90) ∞ w c ( y + y ) dy = 6 e − y (3 + e − y )(1 + e − y ) (2.6)Note that η (0) = 3 and η → + monotonically as z → ∞ . In a similar way we obtain the innersolution near x = 1 by letting x = 1 − εy and finding that u ∼ ξ R w c ( y + y R ) + O ( ε ) , v ∼ ξ R + εv R, ( y ) + O ( ε ) , where y R is determined by solving w (cid:48) c ( y R ) = − B/ξ R , (2.7)and for which we calculate the limit lim y → + ∞ v (cid:48) R, ( y ) = − ξ R D η ( y R ) . (2.8)We now consider the inner solution at each interior spike location. By balancing dominant termsin a higher order asymptotic expansion, it can be shown that the interior spike locations var on an O ( ε − ) timescale. Therefore we let x i = x i ( ε t ) for each i = 1 , ..., N and with x = x i ( ε t ) + y wecalculate the inner asymptotic expansions u ∼ ξ i w c ( y ) + O ( ε ) , v ∼ ξ i + εv i ( y ) + O ( ε ) , i = 1 , ..., N. (2.9)Furthermore, we must impose a solvability condition on the v i problem which gives1 ε dx i dt = − ξ i (cid:18) lim y → + ∞ v (cid:48) i ( y ) + lim y →−∞ v (cid:48) i ( y ) (cid:19) . (2.10)To determine the 2( N + 2) undetermined constants ξ i and y i where i ∈ { L, , ..., N, R } we mustnow calculate the outer solution, defined for | x − x i | = O (1) for each boundary and interior spikelocation, and match it with each of the inner solutions. Since w c decays to zero exponentially as y → ±∞ we determine that the activator is asymptotically small in the outer region. On the otherhand (2.5) and (2.8) imply the boundary conditions v x (0) ∼ − ξ L D η ( y L ) , v x (1) ∼ ξ R D η ( y R ) , while the exponential decay of w c implies that the following limits hold (in the sense of distributions) ε − u −→ N (cid:88) j =1 ξ j δ ( x − x j ) , ( ε → + ) , Thus, to leading order in ε (cid:28)
1, the outer problem for the inhibitor is given by Dv xx − v = − N (cid:88) j =1 ξ j δ ( x − x j ) , < x < , (2.11) Dv x (0) = − ξ η ( y L ) , Dv x (1) = ξ N +1 η ( y R ) . (2.12)This boundary value problem can be solved explicitly by letting G ω be the Green’s function satisfying G ω,xx − ω G ω = − δ ( x − ξ ) , < x < G ω,x (0 , ξ ) = 0 , G ω,x (1 , ξ ) = 0 , ω > G ω ( x, ξ ) = 1 ω sinh ω (cid:40) cosh ωx cosh ω (1 − ξ ) , < x < ξ, cosh ω (1 − x ) cosh ωξ, ξ < x < . (2.13b)Formally substituting ξ = 0 or ξ = 1 into the above expression gives G ω ( x,
0) = cosh ω (1 − x ) ω sinh ω , G ω ( x,
1) = cosh ωxω sinh ω , (2.13c)which is readily seen to satisfy G ω,xx − ω G ω = 0 , < x < , with boundary conditions G ω,x (0 ,
0) = − , G ω,x (1 ,
0) = 0 , G ω,x (0 ,
1) = 0 , G ω,x (1 ,
1) = 1 . Letting ω ≡ D − / we obtain the following leading order asymptotic expansion for the quasi-equilibrium solution to (1.2) u e ( x ) ∼ ξ L w c (cid:0) xε + y L (cid:1) + ξ R w c (cid:0) − xε + y R (cid:1) + N (cid:88) j =1 ξ j w c (cid:0) x − x j ε (cid:1) , (2.14a) v e ( x ) ∼ ω (cid:18) ξ L η ( y L ) G ω ( x,
0) + ξ R η ( y R ) G ω ( x,
1) + 6 N (cid:88) j =1 ξ j G ω ( x, x j ) (cid:19) . (2.14b)Furthermore, by imposing the consistency condition v e ( x i ) = ξ i for each i ∈ { L, , ..., N, R } weobtain the system of N + 2 nonlinear equations B ≡ ξ − ω G ω N ξ = 0 , (2.15a)where G ω and N are the ( N + 2) × ( N + 2) matrices given by( G ω ) ij = G ω ( x i , x j ) i, j = L, R, , ..., N, N ≡ diag( η ( y L ) , η ( y R ) , , ..., , (2.15b)and ξ ≡ ( ξ L , ξ R , ξ , ..., ξ N ) T , ξ = ( ξ L , ξ R , ξ , ..., ξ N ) T . (2.15c) D. GOMEZ AND J. WEI
Thus, for given spike configuration 0 < x < ... < x N <
1, the system (2.15a) together with (2.3)and (2.7) can be solved for the unknown spike heights ξ L , ξ , ..., ξ N , ξ R and boundary shifts y L and y R . Summarizing, we have the following proposition. Proposition 2.1.
In the limit ε → + and for t (cid:28) O ( ε − ) an N + 2 spike quasi-equilibriumsolution to (1.2) consisting of two boundary spikes and N well separated interior spikes concentratedat specified locations < x < ... < x N < is given asymptotically by (2.14) where G ω ( x, ξ ) isgiven explicitly by (2.13b) and ω = D − / . The boundary shifts, y L and y R , and spike heights, ξ L , ξ , ..., ξ N , ξ R , are found by solving the system of N + 4 equations (2.3) , (2.7) , and (2.15a) . The asymptotic solution constructed in the above proposition will not generally be an equilibriumof (1.2) due to the slow, O ( ε − ), drift motion of the interior spikes described by (2.10). However,this solution can be made into an equilibrium by choosing the interior spike locations x , ..., x N appropriately. Proposition 2.2.
The interior spike locations of a multi-spike pattern consisting of two boundaryspikes and N interior spikes vary on an O ( ε − ) time scale according to the differential equation ε dx i dt = − ξ i D (cid:10) ∂ x G ω ( x, x i ) (cid:11) x = x i − ξ i D (cid:88) j (cid:54) = i ξ j G x ( ξ i , ξ j ) − ξ i D (cid:20) ξ L η ( y R ) G x ( x i ,
0) + ξ R η ( y R ) G x ( x i , (cid:21) , ( i = 1 , ..., N ) , (2.16) where (cid:10) f ( x ) (cid:11) x = lim x → x +0 f ( x ) + lim x → x − f ( x ) , (2.17) which is to be solved together with (2.15a) , (2.3) , and (2.7) for the spike heights ξ L , ξ , ..., ξ N , ξ R andshifts y L and y R . In particular, if the configuration x , ..., x N is stationary with respect to the ODE (2.16) , then to leading order the quasi-equilibrium solution of Proposition 2.1 is an equilibrium forall t ≥ . Equilibrium Multi-Spike Solutions by the Gluing Method.
We now use an alternativemethod for constructing asymmetric multi-spike equilibrium solutions to (1.2). This method extendsthat of Ward and Wei [17] to account for inhomogeneous Neumann boundary conditions. The keyidea is to construct a single boundary spike solution in an interval of variable length and use thissolution to glue together a multi-spike solution. In particular, we begin considering the problem ε u xx − u + v − u = 0 , Dv xx − v + ε − u = 0 , < x < l, (2.18a) εu x (0) = − A, εu x ( l ) = 0 , Dv x (0) = 0 , Dv x ( l ) = 0 , (2.18b)where l > x = 0. Proceeding as in § x = O (1)) is given by u ( x ; l, A ) ∼ ξ w c (cid:0) ε − x + y (cid:1) , v ( x ; l, A ) ∼ ξ cosh ω ( l − x )cosh ω l , (2.19)where ω ≡ D − / while the shift parameter y and spike height ξ satisfy w (cid:48) c ( y ) = − Aξ , (2.20)and for which, using (2.2b) and (2.6), we explicitly calculate ξ = tanh ω lω η ( y ) , y = log (cid:18) q + (cid:112) q + 10 q + 12 (cid:19) , q ≡ ω A tanh ω l , (2.21) for which we remark that y ∼ q as q → ξ ∼ (3 ω ) − tanh ω l as A → + . Finally,we note that y is monotone increasing in A and monotone decreasing in D and l when A > < x < N + 2subintervals defined by x L = 0 , x i = l L + 2 i (cid:88) j =1 l j + l i ( i = 1 , ..., N ) , x R = 1 ,I L = [0 , l L ) , I i = [ x i − l i , x i + l i ) ( i = 1 , ..., N ) , I R = [1 − l R , , where l L , l , ..., l N , l R are chosen to satisfy the N + 2 constraints l L + 2 l + · · · l N + l R = 1 . (2.22a) v ( l L ; l L , A ) = v ( l ; l ,
0) = · · · = v ( l N ; l N ,
0) = v ( l R ; l R , B ) , . (2.22b)The first constraint guarantees that the intervals are mutually disjoint, while the second set of N + 1constraints guarantees the continuity of the multi-spike equilibrium solution u e ( x ) = u ( x ; l L , A ) , x ∈ I L u ( | x − x i | ; l i , , x ∈ I i u (1 − x, l R , B ) , x ∈ I R , v e ( x ) = v ( x ; l L , A ) , x ∈ I L v ( | x − x i | ; l i , , x ∈ I i v (1 − x, l R , B ) , x ∈ I R . (2.22c)We remark that the local symmetry of each interior spike implies that the interior spikes are sta-tionary with respect to the slow dynamics found in (2.10), and therefore the multi-spike solutionconstructed above is an equilibrium of (1.2).3. Linear Stability of Multi-Spike Pattern
In this section we derive a nonlocal eigenvalue problem (NLEP) that, to leading order in ε (cid:28) O (1)timescale. Letting u e and v e be the quasi-equilibrium solution from Proposition 2.1, we considerthe perturbations u = u e + e λt Φ and v = v e + e λt Ψ with which (1.2) becomes ε Φ xx − Φ + 2 u e v e Φ − u e v e Ψ = λ Φ , < x < , (3.1a) D Ψ xx − Ψ + 2 ε − u e Φ = τ λ Ψ , < x < . (3.1b)This problem admits both large and small eigenvalues characterized by λ = O (1) and O ( ε ) respec-tively. The small eigenvalues are closely related to the linearization of the slow-dynamics (2.16) andthe resulting instabilities therefore take place over a O ( ε − ) timescale [19]. In contrast, the largeeigenvalues lead to amplitude instabilities over a O (1) timescale. In this section we focus exclusivelyon the large eigenvalues and limit our discussion of the small eigenvalues to the specific examplegiven in § § ε (cid:28)
1, the inhibitor perturbation Ψ satisfies D Ψ xx − (1 + τ λ )Ψ = − N (cid:88) j =1 ξ j (cid:90) ∞−∞ w c ( y ) φ j ( y ) dyδ ( x − x j ) , < x < , (3.2a) D Ψ x (0) = − ξ L (cid:90) ∞ w c ( y + y L ) φ L ( y ) dy, D Ψ x (1) = 2 ξ R (cid:90) ∞ w c ( y + y R ) φ R ( y ) dy, (3.2b) D. GOMEZ AND J. WEI where φ L and φ R are the leading order inner expansions of the activator perturbation Φ at theboundaries satisfying L y i φ i − w c ( y + y i ) Ψ( x i ) = λφ i , < y < ∞ , φ (cid:48) i (0) = 0 , φ i → y → ∞ , (3.3a)for i = L, R respectively, while φ , ..., φ N are the leading order inner expansions of Φ at each of theinterior spike locations x , .., x N satisfying L φ i − w c ( y ) Ψ( x i ) = λφ i , −∞ < y < ∞ , φ → y → ±∞ , (3.3b)for each i = 1 , ..., N respectively. The linear differential operator L y parametrized by y ≥ L y φ ≡ φ (cid:48)(cid:48) − φ + 2 w c ( y + y ) φ. (3.4)Note that by decomposing each φ i = φ even i + φ odd i ( i = 1 , ..., N ) where φ even i and φ odd i are evenand odd about y = 0 respectively, we find that either φ odd i = 0 or else λ ≤
0. In particular, theodd components of each φ i ( i = 1 , ..., N ) do not contribute to any instabilities and without lossof generality we may therefore assume that each φ i is even about y = 0. Hence it suffices to pose(3.3b) on the half line with the same homogeneous Neumann boundary conditions used in (3.3a).Letting ω λ ≡ (cid:112) (1 + τ λ ) /D and recalling the definition of G ω in (2.13) we readily find that thesolution to (3.2) is explicitly given byΨ( x ) = 2 ω N (cid:88) j = L,R, ˆ ξ j G ω λ ( x, x j ) (cid:90) ∞ w c ( y + y j ) φ j ( y ) dy. where we let y = ... = y N = 0 , ˆ ξ i ≡ (cid:40) ξ i , i = L, R, ξ i , i = 1 , ..., N, . (3.5)Evaluating Ψ( x ) at each x = x i and substituting into (3.3) yields the system of NLEPs L y i φ i − ω w c ( y + y i ) N (cid:88) j = L,R, ˆ ξ j G ω λ ( x i , x j ) (cid:90) ∞ w c ( y + y j ) φ j ( y ) dy = λφ i , y > φ (cid:48) i (0) = 0 , φ i → y → + ∞ . (3.6b)for each i = L, R, , ..., N where L y i is defined by (3.4).The NLEP system (3.6a) has two key features that distinguish it from analogous NLEPs insingularly perturbed reaction diffusion systems [6, 17, 18, 19] and are explored in the rigorousstability results of § §
5. First, it considers both boundary-bound and interior-bound spikes. As explored in Examples 2 to 4 this has immediate consequencesfor both the existence and stability of asymmetric patterns even in the zero-flux case where A = B = 0. The second distinguishing feature of (3.6a) is the introduction of the shift parameters y L ≥ y R ≥
0. We remark that an analogous negative shift parameter has been examined inthe context of near-boundary spike solutions for homogeneous
Robin boundary conditions [1, 11].However, as highlighted in the stability results of §
4, the positive shift parameter plays a key role inthe stability properties of boundary-bound spikes. An important critical value of the shift parameteris the unique value y c > w (cid:48)(cid:48) c ( y c ) = 0 and which is explicitly given by y c = log(2 + √ . (3.7)In particular, it can be shown that if y ≶ y c then L y has an unstable and stable spectrumrespectively (see Lemma 4.1 below). Moreover the operator L y c has a one-dimensional kernelspanned by w (cid:48) c ( y + y c ). Reduction of NLEP to an Algebraic System.
It is particularly useful to rewrite (3.6a) asan algebraic system as follows. Assuming that λ is not an eigenvalue of L y i for all i = L, R, , ..., N we let φ i = c i ( L y i − λ ) − w c ( y + y i ) , i ∈ { L, R, , ..., N } , (3.8)where the coefficients c L , c R , c , ..., c N are undetermined. Note that in (3.8) the homogeneous Neu-mann boundary condition φ (cid:48) i (0) = 0 is assumed. In addition, note that if λ = 0 then (3.8) is onlyvalid if y L , y R (cid:54) = y c .Substituting into (3.6a) then yields the linear homogeneous system for c ≡ ( c L , c R , c , ..., c N ) T G ω λ D λ ccc = (2 ω ) − ccc, (3.9)where G ω λ is the ( N + 2) × ( N + 2) matrix with entries( G ω λ ) ij = G ω λ ( x i , x j ) , ( i, j = L, R, , ..., N ) , (3.10)while D λ is the diagonal ( N + 2) × ( N + 2) matrix given by( D λ ) ij = 1 ω η ( y i ) ξ i F y L ( λ ) , i = j = L, R, ξ i F ( λ ) , i = j = 1 , ..., N, , i (cid:54) = j, (3.11)where F y ( λ ) ≡ (cid:82) ∞ w c ( y + y )( L y − λ ) − w c ( y + y ) dy (cid:82) ∞ w c ( y + y ) dy , (3.12)and for which (4.3) and (4.4) below imply that for all y (cid:54) = y c F y (0) = 1 + w (cid:48) c ( y ) w c ( y ) w (cid:48)(cid:48) c ( y ) η ( y ) . (3.13)Comparing (3.9) and (3.6a), it follows that λ is an eigenvalue of (3.6a) if and only if (2 ω ) − is aneigenvalue of G ω λ D λ . In particular, when λ is not an eigenvalue of L y L , L y R , and L then it is aneigenvalue of the NLEP (3.6a) if and only if it satisfies the algebraic equationdet (cid:0) I N +2 − ω G ω λ D λ (cid:1) = 0 , (3.14)where I N +2 is the ( N + 2) × ( N + 2) identity matrix.We conclude by noting that if either y L = y c and/or y R = y c then the algebraic reduction failswhen searching for a zero eigenvalue λ = 0 since L y c is not invertible. However, in this case we candeduce an analogous system. In particular letting λ = 0 and assuming that y L = y c and y R (cid:54) = y c ,we multiply the i = L NLEP in (3.6a) by w (cid:48) c ( y + y c ) and integrate over 0 < y < ∞ to get ξ L (cid:90) ∞ w c ( y + y L ) φ L dy = − G ω (0 , N (cid:88) j = R, ˜ ξ j G ω (0 , x j ) (cid:90) ∞ w c ( y + y j ) φ j dy. (3.15)Proceeding as above we then deduce that the NLEP (3.6a) with λ = 0 is then equivalent to thealgebraic equation det( I N +1 − ω ˜ G ω ˜ D ) = 0 , (3.16)where ˜ G ω and ˜ D are the ( N + 1) × ( N + 1) matrices with entries( ˜ G ω ) ij = G ij − G ω (0 , G ω ( x i , G ω (0 , x j ) , ( ˜ D ) ij = D ij , i, j = R, , ..., N. (3.17)The same approach can likewise be used if y L = y R = y c . Zero-Eigenvalues of the NLEP and the Consistency Condition.
The conditions underwhich λ = 0 is an eigenvalue of the NLEP (3.6a) can be directly linked to the system (2.15a) ashighlighted in [19]. Specifically, assume that x , ..., x N are fixed (not necessarily at an equilibriumconfiguration of the slow dynamics ODE (2.16)) and let ξ L , ξ R , ξ , ..., ξ N together with y L and y R solve (2.15a), (2.3), and (2.7) with the additional assumption that y L , y R (cid:54) = y c . From the definitionof η in (2.6) and from (2.3) and (2.7) we calculate ∂η ( y i ) ∂ξ i = w c ( y i ) w (cid:48) c ( y i ) ξ i w (cid:48)(cid:48) c ( y i ) , (3.18)for i = L, R . Taking the Jacobian of the quasi-equilibrium system (2.15a) and recalling the definitionof D λ given in (3.11) we deduce that ∇ ξ B = I − ω G ω D . (3.19)Together with the discussion of § y L , y R (cid:54) = 0 and each x , ..., x N is independentof ξ L , ξ R , ξ , ..., ξ N , then λ = 0 is an eigenvalue of the NLEP (3.6a) if and only if the Jacobian ∇ ξ B is singular.4. Rigorous Stability and Instability Results for the Shifted NLEP
In this section we rigorously prove instability and stability results for the shifted
NLEP L y φ − µ (cid:82) ∞ wφ (cid:82) ∞ w w = λφ, < y < ∞ ; φ (cid:48) (0) = 0; φ → y → ∞ , (4.1)where µ is a real constant and for a fixed value of y ≥ L y φ ≡ φ (cid:48)(cid:48) − φ + 2 wφ, w ( y ) ≡ w c ( y + y ) , (4.2)and where w c is the unique solution to (2.2a). When y = 0 the NLEP (4.1) is stable if µ > µ < L y and its spectrum.First, we calculate L − y w = w − w (cid:48) (0) w (cid:48)(cid:48) (0) w (cid:48) , L − y w = w + 12 yw (cid:48) − w (cid:48) (0)2 w (cid:48)(cid:48) (0) w (cid:48) (4.3)where the additional terms are chosen so that homogeneous Neumann boundary conditions at y = 0are satisfied and which we use to compute (cid:90) ∞ w L − y w = (cid:90) ∞ w + w (cid:48) (0) w (0) w (cid:48)(cid:48) (0) , (cid:90) ∞ w L − y w = 34 (cid:90) ∞ w + 3 w (cid:48) (0) w (0) w (cid:48)(cid:48) (0) , (4.4a) (cid:90) ∞ w L − y w = (cid:90) ∞ w + w (cid:48) (0) w (0) w (cid:48)(cid:48) (0) , (cid:90) ∞ w = 65 (cid:90) ∞ w + 3 w (0) w (cid:48) (0)5 . (4.4b)In the next two lemmas, we describe some key properties of the eigenvalue problem L y Φ = ΛΦ , < y < ∞ ; Φ (cid:48) (0) = 0; Φ → , as y → + ∞ . (4.5) Lemma 4.1.
Let y ≥ an let Λ be the principal eigenvalue of (4.5) . Then Λ = 0 if y = y c and Λ ≶ if y ≷ y c . Furthermore, the eigenfunction corresponding to the principal eigenvalue isof one sign.Proof. Since L y is self-adjoint, the variational characterization − Λ = inf Φ ∈ H ([0 , ∞ )) (cid:82) ∞ | Φ (cid:48) | + | Φ | − w | Φ | (cid:82) ∞ | Φ | , (4.6) implies that the principal eigenfunction Φ is of one sign. Since Φ (0) (cid:54) = 0 we may, without loss ofgenerality, assume that Φ (0) = 1 and Φ >
0. Now we multiply (4.5) by w (cid:48) and integrate by partsto get Λ = w (cid:48)(cid:48) (0) (cid:82) ∞ w (cid:48) Φ dy , (4.7)where we remark that the denominator is negative since w (cid:48) ≤ y ≥
0. The claim follows bynoting that w (cid:48)(cid:48) (0) = 0 when y = 0 and w (cid:48)(cid:48) (0) ≷ y ≷ y c . (cid:3) Lemma 4.2.
Let Λ be the second eigenvalue of L y . Then Λ < for all y ≥ .Proof. First note that the second eigenfunction Φ must cross zero at least once since (cid:82) ∞ Φ Φ dy =0 and Φ is of one sign. Next we assume toward a contradiction that Λ ≥
0. We begin by showingthat Φ has exactly one zero in 0 < y < ∞ . Assume that Φ has more than one zero and choose0 < a < b < ∞ such that Φ ( a ) = Φ ( b ) = 0 and Φ > a < y < b . Then Φ (cid:48) ( a ) > (cid:48) ( b ) < ≥ Λ (cid:90) ba w (cid:48) Φ dy = (cid:90) ba w (cid:48) L y Φ dy = w (cid:48) ( b )Φ (cid:48) ( b ) − w (cid:48) ( a )Φ (cid:48) ( a ) > , (4.8)where we have used L y w (cid:48) = 0 and w (cid:48) < y >
0. Thus Φ has a unique zero 0 < a < ∞ andwe may assume that Φ ≶ y ≶ a . Setting b = ∞ in (4.8) we get a contradiction. (cid:3) In Figure 2a we plot the principal and second eigenvalues of the operator L y which we calculatednumerically (see Appendix B for details on the numerical method).Lemma 4.1 implies that the NLEP will have different stability properties depending on whether y is greater than or smaller than y c . We will henceforth refer to 0 ≤ y < y c and y > y c as the small-shift and large-shift cases respectively. When y = 0 it is known that for µ > L y . Since all the eigenvalues of L y are negative in thelarge-shift case we expect the spectrum of the NLEP (4.1) to remain stable for all µ ≥
0. Restrictingour attention to real eigenvalues, we have the following stability result for the large-shift case.
Theorem 4.1.
All real eigenvalues of the NLEP (4.1) are negative when y > y c . To prove this, we first prove the following lemma.
Lemma 4.3.
Let y > y c and suppose that φ satisfies L y φ − λφ ≥ , < y < ∞ ; φ (cid:48) (0) ≥ φ → , as y → + ∞ , (4.9) where λ ≥ . Then φ < for all y ≥ .Proof. Assume toward a contradiction that φ > ≤ a < y < b ≤ ∞ . Without loss of generalitywe may assume that φ ( a ) = 0 if a > φ (0) > a = 0. Then, for any such 0 ≤ a < b ≤ ∞ we have φ ( a ) ≥ , φ (cid:48) ( a ) ≥ , φ ( b ) = 0 , φ (cid:48) ( b ) ≤ . (4.10)Let g ( y ) ≡ w (cid:48)(cid:48) ( y ) − βw (cid:48) ( y ) where β ≡ max y ≥ | w (cid:48)(cid:48)(cid:48) ( y ) | w (cid:48)(cid:48) ( a ) is well-defined and positive. Then g > y ≥ g (cid:48) ( a ) ≤
0, and moreover ( L y − λ ) g = L y w (cid:48)(cid:48) − λg = − ( w (cid:48) ) − λg <
0. Integrating by partswe obtain the contradiction0 < (cid:90) ba g ( L y − λ ) φdy − (cid:90) ba φ ( L y − λ ) gdy = g ( b ) φ (cid:48) ( b ) − g ( a ) φ (cid:48) ( a ) − g (cid:48) ( b ) φ ( b )+ g (cid:48) ( a ) φ ( a ) ≤ . (4.11) (cid:3) Proof [Theorem 4.1. ] Suppose that λ ≥ L y − λ is invertible and from (4.1) we calculate φ = µ (cid:82) ∞ wφ (cid:82) ∞ w ( L y − λ ) − w . But w > µ (cid:82) ∞ w ( L y − λ ) − w (cid:82) ∞ w < . (4.12) (cid:3) From Theorem 4.1 we immediately deduce that the NLEP does not admit a zero eigenvalue forany µ ≥ y > y c . On the other hand, when 0 ≤ y ≤ y c we suspect that the NLEP admitsa zero eigenvalue for an appropriate choice of µ ≥
0. When y = y c this is the case for µ = 0.When 0 ≤ y < y c we set λ = 0 in (4.1) and obtain L y φ = µ (cid:82) ∞ wφ (cid:82) ∞ w w . (4.13)Using (4.3) we calculate φ = L − y w and substitute back into (4.13) to deduce that λ = 0 is aneigenvalue if and only if µ = µ c ( y ) ≡ (cid:82) ∞ w (cid:82) ∞ w L − y w = (cid:82) ∞ w (cid:82) ∞ w + w (cid:48) (0) w (0) w (cid:48)(cid:48) (0) . (4.14)Note that µ c ( y ) ≶ y ≷ y c . In terms of the critical value µ c we have the following instabilityresult for the small-shift case. Theorem 4.2.
Let ≤ y < y c and ≤ µ < µ c where critical value µ c is defined in (4.14) . Thenthe NLEP (4.1) admits a positive real eigenvalue.Proof. Let Λ be the principal eigenvalue of L y . First note that by Lemma 4.1 and 4.2 the principaland second eigenvalues of L y satisfy Λ < < Λ . Moreover the corresponding eigenfunction Φ is of one sign and we may assume that Φ > (cid:82) ∞ Φ = 1. Observe that if λ (cid:54) = Λ is a positiveeigenvalue of the NLEP (4.1) then φ = µ (cid:82) ∞ wφ (cid:82) ∞ w ( L y − λ ) − w , and since (cid:82) ∞ wφ (cid:54) = 0 the above equation is equivalent to h ( λ ) = 0 where h ( λ ) ≡ (cid:90) ∞ w ( L y − λ ) − w − (cid:82) ∞ w µ . (4.15)We now show that such a λ can always be found in 0 < λ < Λ for 0 ≤ µ < µ c .First we calculate h (0) = (cid:82) ∞ w ( µ − c − µ − ) <
0. Next we we let ψ be the unique solution to( L y − λ ) ψ = w , < y < ∞ ; ψ (cid:48) (0) = 0 . Decomposing ψ = c Φ + ψ ⊥ where (cid:82) ∞ Φ ψ ⊥ = 0 we find that ψ ⊥ satisfies( L y − λ ) ψ ⊥ = w − c (Λ − λ )Φ , < y < ∞ ; ( ψ ⊥ ) (cid:48) (0) = 0 . (4.16)Multiplying by Φ and integrating by parts we obtain c = (Λ − λ ) − (cid:82) ∞ w Φ and therefore h ( λ ) = (cid:82) ∞ w Φ (cid:82) ∞ w Φ Λ − λ + (cid:90) ∞ wψ ⊥ − (cid:82) ∞ w µ . (4.17) On the other hand, if we multiply (4.16) by ψ ⊥ and integrate then we obtain − (cid:90) ∞ | ψ ⊥ | (cid:18) λ + − (cid:82) ∞ ψ ⊥ L y ψ ⊥ (cid:82) ∞ | ψ ⊥ | (cid:19) = (cid:90) ∞ w ψ ⊥ . (4.18)By Lemma 4.2 and the variational characterization of the second eigenvalue of L y we obtain0 < − Λ = inf Φ ∈ H ([0 , ∞ )) (cid:82) ∞ ΦΦ =0 (cid:82) ∞ | Φ (cid:48) | + | Φ | − w | Φ | | Φ | ≤ − (cid:82) ∞ ψ ⊥ L y ψ ⊥ (cid:82) ∞ | ψ ⊥ | . Substituting into (4.18) we calculate || ψ ⊥ || L ([0 , ∞ )) ≤ λ − || w || L ([0 , ∞ )) || ψ ⊥ || L ([0 , ∞ )) so that || ψ ⊥ || L ([0 , ∞ )) and hence also (cid:82) ∞ wψ ⊥ are bounded as λ → Λ >
0. Therefore, from (4.17) we deduce h ( λ ) → + ∞ as λ → Λ − . By a continuity argument we deduce the existence of a λ ∈ (0 , Λ ) such that h ( λ ) = 0. (cid:3) We conclude this section by establishing sufficient conditions for the stability of the NLEP (4.1)in both the small- and large-shift cases. Suppose that φ = φ R + iφ I and λ = λ R + iλ I satisfies theNLEP. Separating real and imaginary parts in (4.1) then yields the system L y φ R − µ (cid:82) ∞ wφ R (cid:82) ∞ w w = λ R φ R − λ I φ I , L y φ I − µ (cid:82) ∞ wφ I (cid:82) ∞ w w = λ R φ I + λ I φ R . (4.19a)Multiplying the first and second equations by φ R and φ I respectively, integrating, and then addingthem together gives λ R (cid:90) ∞ | φ | = − L ( φ R , φ R ) − L ( φ I , φ I ) , (4.20)where we define L (Φ , Φ) ≡ (cid:90) ∞ | Φ (cid:48) | + Φ − w Φ + µ (cid:82) ∞ w Φ (cid:82) ∞ w Φ (cid:82) ∞ w . (4.21)It is clear that if L (Φ , Φ) > ∈ H ([0 , ∞ )) then the NLEP (4.1) will be linearly stable.In the next theorem we determine sufficient conditions on µ ≥ y ≥ Theorem 4.3. If ≤ y < y c and µ ( y ) < µ < µ ( y ) , or y > y c and ≤ µ < µ ( y ) where µ ( y ) ≡ (cid:82) ∞ w (cid:82) ∞ w L − y w + (cid:113)(cid:82) ∞ w L − y w (cid:82) ∞ w L − y w , (4.22) µ ( y ) ≡ (cid:82) ∞ w (cid:82) ∞ w L − y w − (cid:113)(cid:82) ∞ w L − y w (cid:82) ∞ w L − y w , (4.23) then Re λ < for all eigenvalues of the NLEP (4.1) .Proof. We first prove the result for 0 ≤ y < y c . When µ = 2 and y = 0, Lemma 5.1 (2) in [18]implies that L (Φ , Φ) > ∈ H (( ∞ , ∞ )) and hence, by restricting to even functions, alsofor all Φ ∈ H ([0 , ∞ )). In particular, by the variational characterization of the principal eigenvalue,this implies that the principal eigenvalue of the self-adjoint operator L Φ ≡ L y Φ − µ (cid:82) ∞ w Φ (cid:82) ∞ w w − µ (cid:82) ∞ w Φ (cid:82) ∞ w w, (4.24)must be negative. We then perturb y ≥ µ until L has a zero eigenvalue and for which wemay solve Φ = c L − y w + c L − y w. (4.25) y Eigenvalues of y (a) y Stability Thresholds (b) y ( y ) c ( y ) for 0 y y c (c) Figure 2. (A) Plot of the numerically computed principal and second eigenvaluesof the operator L y . The dashed vertical line corresponds to y = y c . (B) Plotof the stability thresholds µ and µ as functions of y . The dashed vertical andhorizontal lines correspond to y = y c and µ = 2 respectively. The NLEP (4.1) hasbeen rigorously demonstrated to be stable in the region bounded by the curves µ and µ . Note that µ and µ are interchanged as y passes through y c . (C) Plotof µ ( y ) − µ c ( y ) for 0 ≤ y < y c . The NLEP is unstable for µ < µ c and stablefor µ < µ < µ when 0 ≤ y < y c . It is conjectured that the NLEP is stable for µ > µ c .Substituting back into L Φ = 0 we obtain the system (cid:18) µ (cid:82) ∞ w L − y w (cid:82) ∞ w − (cid:19) c + µ (cid:82) ∞ w L − y w (cid:82) ∞ w c = 0 , µ (cid:82) ∞ w L − y w (cid:82) ∞ w c + (cid:18) µ (cid:82) ∞ w L − y w (cid:82) ∞ w − (cid:19) c = 0 . (4.26)Since (cid:82) ∞ w L − y w = (cid:82) ∞ w L − y w a nontrivial solutions exists if and only if (cid:18) µ (cid:82) ∞ w L − y w (cid:82) ∞ w − (cid:19) − µ (cid:82) ∞ w L − y w (cid:82) ∞ w L − y w (cid:0)(cid:82) ∞ w (cid:1) = 0 , (4.27)where explicit formulae for each integral can be found in (4.4). When µ = 2 and y = 0 the lefthand side of (4.27) equals − / < ≤ y < y c if (cid:18) µ (cid:82) ∞ w L − y w (cid:82) ∞ w − (cid:19) − µ (cid:82) ∞ w L − y w (cid:82) ∞ w L − y w (cid:0)(cid:82) ∞ w (cid:1) < , (4.28)which is easily seen to be equivalent to µ < µ < µ .The thresholds µ and µ are singular as y → y c and therefore the continuity argument fromabove does not extend to y > y c . However L (Φ , Φ) > y > y c and µ = 0.Therefore we proceed with the same continuity argument as above, but starting from µ = 0.This yields the same criteria, but since L − y w < µ ( y ) < ≤ µ < µ ( y ). (cid:3) Both of the stability thresholds µ and µ defined in (4.22) as well as the instability threshold µ c defined in (4.14) are easily computed using (4.4). In Figures 2b and 2c we plot the stabilitythresholds and the difference µ − µ c respectively. In particular, from the plot in 2c we see that µ > µ c . We conjecture, that as in the y = 0 case, the NLEP is stable for all µ > µ c . In Appendix A , I Hopf Thresholds when D hh (a) D h ( D , A ) A (b) D h h ( D , A ) A (c) Figure 3.
Hopf bifurcation threshold and accompanying eigenvalue for a singleboundary-spike solution with one-sided boundary flux A ≥ D → ∞ , and (B and C) for finite D > ≤ A < q c . In(A) the dashed vertical line corresponds to the threshold A = q c beyond which noHopf instabilities occur.B we provide numerical support for this conjecture by plotting Re λ versus µ and y in Figure 12a.In addition, we plot Λ − Re( λ ) in Figure 12b which suggest that Re( λ ) ≤ Λ .5. Examples
In this section we illustrate the effect of introducing a nonzero boundary flux for the activator byconsidering three distinct examples. Specifically, we first study the stability of a single boundaryspike concentrated at x = 0 when A ≥ B = 0. Using a winding number argument we illustratethat the stability of the single spike is improved by increasing the boundary flux A . Moreover, weillustrate that if A exceed a threshold, then the spike is stable independently of the parameters τ ≥ D >
0. We then consider the structure and stability of a two-boundary-spike patternwhen the boundary fluxes are equal, A = B ≥
0. One of the key findings is that if
A >
D > τ (cid:28) A ≥ B = 0. We demonstrate the existenceof several asymmetric patterns, with a certain branch of these patterns always being stable. Foreach example we include full numerical simulations of the GM system (1.2) using the finite elementsoftware FlexPDE [13].5.1. Example 1: One Boundary Spike at x = 0 with A > and B = 0 . In this example weassume that B = 0 and investigate the role of a non-negative flux, A ≥
0, on the stability of a singleboundary spike concentrated at x = 0. We denote the left boundary shift parameter by y = y L sothat using (2.13b) and (2.21) we reduce (3.6a) to (4.1) where µ ( λ ) = 2 ω tanh ω ω λ tanh ω λ . (5.1) t u ( , t ) D = 5, = 1 A (a) t u ( , t ) D = 5, = 2 A (b) t u ( , t ) D = 5, = 7 A (c) Figure 4.
Plots of u (0 , t ) for a one boundary-spike solution with one-sided bound-ary flux x = 0 (i.e. A ≥ B = 0) with ε = 0 . A stabilizes the single boundary-spike solution for fixed values of D and τ .Recalling that ω λ = (cid:112) (1 + τ λ ) /D we first observe that 0 < µ ( λ ) ≤ λ ≥ λ is not in the spectrum of L y we let φ = ( L y − λ ) − w c ( y + y ) so that as in § A y ( λ ) ≡ µ ( λ ) − F y ( λ ) = 0 , (5.2)where F y ( λ ) is defined in (3.12). In Figure 11 we plot F y (0) versus y , as well as the real andimaginary parts of F y ( iλ I ) versus λ I for select values of y . We integrate in λ over a closedcounter-clockwise contour consisting of the imaginary axis and a large semicircle in the right half-plane. Since F y ( λ ) = O ( | λ | − ) and µ ( λ ) = O ( λ − / ) as | λ | → ∞ , Re λ > , (5.3)the change in argument of A y ( λ ) over the large semicircle is π . Moreover, in Re λ > µ ( λ ) (cid:54) = 0whereas by Lemmas 4.1 and 4.2 we deduce that F y ( λ ) has one (resp. zero) pole(s) if y < y c (resp. y ≥ y c ). Letting Z denote the number of zeros of A y ( λ ) in Re λ > Z = 1 π ∆arg A y ( iλ I ) (cid:12)(cid:12) ∞ + (cid:40) / , y < y c , / , y ≥ y c , (5.4)where the first term on the right hand side denotes the change in argument of A y ( λ ) as λ followsthe imaginary axis from λ = + i ∞ to λ = 0. Note in addition that we have used A y (¯ λ ) = A y ( λ ) to obtain ∆ arg A y ( iλ I ) (cid:12)(cid:12) −∞ + ∞ = 2∆ arg A y ( iλ I ) (cid:12)(cid:12) ∞ . From (5.3) we immediately deducethat arg A y (+ i ∞ ) = π/
4. On the other hand, using (4.4), we evaluate A y (0) = − F y (0) ≶ y ≶ y c (see also Figure 11a). We will consider the cases y ≥ y c and y < y c separately below.If y ≥ y c then Re F y ( iλ I ) < λ I > A y ( iλ I ) never crosses theimaginary axis for all λ I >
0. As a result ∆ arg A y ( iλ I ) | ∞ = − π/ Z = 0 if y > y c .Since y is monotone decreasing in D when A >
0, we deduce that there is a threshold D c ( A ) such that the single spike pattern is stable for all τ ≥ D ≥ D c ( A ). Substituting y = y c and using(2.21) and (3.7) yields the threshold parameter q c = 4 + 3 √ ≈ . , (5.5)with which D c ( A ) is found by solving the transcendental equationtanh D − / c = q − c AD − / c . (5.6)It is easy to see that if A ≥ q c then D c = ∞ is the only positive solution. On the other hand, if0 < A < q c then this equation has a unique positive solution that is monotone increasing in A andsatisfies D c → + as A → + and D c → ∞ as A → q − c . In summary we deduce that the singlespike pattern is stable for all D > τ ≥ A ≥ q c , or for all τ ≥ < D ≤ D c ( A ) and0 < A < q c . To determine the stability when 0 ≤ A < q c and D > D c ( A ) we must consider thecase 0 ≤ y < y c .Next we assume that 0 ≤ y < y c . We begin by considering the shadow limit , defined by D → ∞ ,for which µ ( iλ I ) ∼ iτ λ I ) − and hence alsoRe A y ( iλ I ) ∼ − Re F y ( iλ I ) , Im A y ( iλ I ) ∼ τ λ I − Im F y ( iλ I ) , D → ∞ . Since F y (0) > / F y ( iλ I ) → λ I → ∞ (see Appendix A and accompanying Figure 11)we deduce that there exists a solution to Re A y ( iλ I ) = 0. Moreover, in Figure 11b we observe thatwhen Re F y ( iλ I ) is positive it is also monotone decreasing in λ I . Therefore there exists a uniqueeigenvalue λ ∞ h and time constant τ ∞ h = 2Im F y ( iλ I ) /λ ∞ h such that A y ( iλ ∞ h ) = 0. Furthermore,since Im A y ( iλ ∞ h ) ≶ τ ≶ τ ∞ h we get∆ arg A y ( iλ I ) | ∞ = (cid:40) − π/ , τ < τ ∞ h , π/ , τ > τ ∞ h . The single boundary spike solution therefore undergoes a Hopf bifurcation as τ exceeds the Hopfbifurcation threshold τ ∞ h . Using the shadow limit threshold as an initial guess, we numericallycontinue the Hopf bifurcation threshold for finite values of D > τ h ( D, A ) and accompanying critical eigenvalue λ = iλ h ( D, A ) shown in Figure 3b and 3crespectively.The above analysis, together with the plots of τ ∞ h ( A ) and τ h ( D, A ) in Figures 3a and 3b re-spectively, indicate that the single boundary spike solution is stabilized as
A > A exceeds the threshold q c given in (5.5), then the single boundary spike is stableindependently of the parameters τ ≥ D >
0. We illustrate the onset of Hopf instabilitieswhen D = 5 for τ = 1 , , A = 0 , . , . , . u (0 , t ) in Figure 4. In particular we observe that the single spikepattern is stabilized by increasing the boundary flux A . Additionally, our numerical simulationsshow good qualitative agreement with the Hopf bifurcation thresholds plotted in Figure 3b.5.2. Example 2: Two Boundary Spikes with A = B ≥ . In this example we investigate therole of equal boundary fluxes on the structure and stability of a two-boundary-spike pattern. Usingthe method of § l L = l and l R = 1 − l andsolving (2.22b), which is explicitly given bytanh ω lη ( y L ) cosh ω l − tanh ω (1 − l ) η ( y R ) cosh ω (1 − l ) = 0 , (5.7) for 0 < l < η is given by (2.6) and y L and y R are given by (2.21). Note that by (2.19) thealgebraic equation (5.7) is equivalent to ξ L ξ R = cosh ω l cosh ω (1 − l ) , (5.8)from which we deduce that l ≶ / ξ L ≶ ξ R . In particular l = 1 / A ≥ ξ L = ξ R we refer to it as the symmetric solution. For the remainder of thisexample we will construct asymmetric two-spike patterns for which 0 < l < / l > / A ≥ A = 0. Specifically, we let z = ω l L and ˜ z = ω l R so thatwhen A = 0 the system (2.22) (and hence also (5.7)) is equivalent to z + ˜ z = ω , b ( z ) = b (˜ z ) , b ( z ) ≡ tanh z cosh z . (5.9)It follows from Result 2.3 (with k = k = 1, µ = 1, and r = 1) of [17] that (5.7) has a uniquesolution 0 < l < < D < D c ≡ [2 log(1 + √ − ≈ . . (5.10)When A > D and A > < l < /
2. Rather than solving (5.7) numerically for l as a function of A and D , we found it more convenient to solve for A = A ( D, l ). The results ofour numerical calculations are shown in Figure 5a where we plot A = A ( D, l ) as well as the curve l = l max ( D ) along which A ( D, l ) is maximized for a fixed value of D , and the curve l = l c ( D ) alongwhich A ( D, l c ( D )) = A ( D, /
2) for D c < D < D c ≈ .
660 and l c ( D ) = 0 for 0 < D < D c .Consequently, (5.7) has zero solutions in 0 < l < / A > A max ( D ) ≡ A ( D, l max ( D )), whereas ithas two solutions, one with l < l c ( D ) and the other with l > l c ( D ), if A max ( D ) > A > A c ( D ) ≡ A ( D, l c ( D )) for 0 < D < D c . For all other values of D >
A > < l < /
2. We summarize these existence thresholds in Figure 5b wherewe note in particular that
A > D values over which asymmetrictwo-boundary-spike patterns exist.We now consider the linear stability of both the symmetric and asymmetric two-boundary-spikepatterns constructed above. Note that since x L = 0 and x R = 1 are fixed there are no slow driftdynamics and the linear stability of the two-spike patterns is completely determined by the O (1)eigenvalues calculated from the NLEP (3.6a). Moreover, we assume that τ = 0 so that no Hopfinstabilities arise (see Example 1) and for which we can exclusively focus on zero eigenvalue crossing,or competition, instabilities. We proceed by first using the rigorous results of § § L y φφφ − ω tanh (cid:0) ω (cid:1) (cid:82) ∞ w c ( y + y ) G ω φφφdy (cid:82) ∞ w c ( y + y ) dy w c ( y + y ) = λφφφ, < y < ∞ ; φφφ (cid:48) (0) = 0 , where φφφ = (cid:18) φ φ (cid:19) , G ω = 1 ω sinh ω (cid:18) cosh ω
11 cosh ω (cid:19) . (5.11) (a) (b) Figure 5. (A) Plots of A = A ( D, l ) for the construction of asymmetric two-boundary-spikes in Example 2. The curve l = l max ( D ) indicates the values of l at which A ( D, l ) is maximized, while l c ( D ) < / l at which A ( D, l ) has the same value as A ( D, / l = l max ( D ) also corresponds tothe competition instability threshold for asymmetric two-boundary-spike patterns. If l max ( D ) < l < / l < l max ( D ). (B) Thresholds for the existence of zero, one, or two asym-metric two-boundary-spike solutions in the presence of equal boundary fluxes. Thevalues of A = A c ( D ) for D c < D < D c and A = A max ( D ) for D > D c correspondto the competition instability threshold for the symmetric two-boundary-spike pat-tern. It is linearly unstable for values of A below this threshold, and stable otherwise.Note that no asymmetric solutions exist for A = 0 and D > D c .The Green’s matrix is symmetric and of constant row sum and therefore has eigenvectors (1 , ± T .Substituting φφφ = ( φ, ± φ ) T into the NLEP therefore yields two uncoupled scalar NLEPs of the form(4.1) with µ = µ ± where µ + ≡ , µ − ≡ (cid:0) ω (cid:1) . (5.12)From Theorem 4.3 and accompanying Figure 2b we immediately deduce that the φφφ + mode islinearly stable. On the other hand, by Theorem 4.2 the φφφ − mode is unstable if µ − < µ c , where µ c is defined by (4.14). We therefore calculate the competition instability threshold by numericallysolving 2 tanh ω = µ c ( y ) where y = y L = y R is the shift parameter given by (2.21) with l = 1 / A ( D, /