UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

The stability of spot patterns for the Brusselator reaction-diffusion system in two space dimensions… Chang, Yifan 2014

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


24-ubc_2014_september_chang_yifan.pdf [ 599.77kB ]
JSON: 24-1.0167555.json
JSON-LD: 24-1.0167555-ld.json
RDF/XML (Pretty): 24-1.0167555-rdf.xml
RDF/JSON: 24-1.0167555-rdf.json
Turtle: 24-1.0167555-turtle.txt
N-Triples: 24-1.0167555-rdf-ntriples.txt
Original Record: 24-1.0167555-source.json
Full Text

Full Text

The Stability of Spot Patterns for theBrusselator Reaction-Diffusion Systemin Two Space Dimensions:Periodic and Finite Domain SettingsbyChang, YifanB.Sc., Peking University, 2012A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Mathematics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)July 2014c© Chang, Yifan 2014AbstractIn this thesis, we asymptotically construct steady-state localized spot solu-tions to the Brusselator reaction-diffusion system in the semi-strong inter-action regime characterized by an asymptotically large diffusivity ratio. Weconsider two distinct settings: a periodic pattern of localized spots in R2 con-centrating at lattice points of a Bravais lattice, and multi-spot solutions thatconcentrate around some discrete points inside a finite domain. We use themethod of matched asymptotic expansions, Floquet-Bloch theory, and thestudy of certain nonlocal eigenvalue problems to perform a linear stabilityanalysis of these patterns. This analysis leads to a two-term approximationfor a certain stability threshold characterized by a zero-eigenvalue crossing.Numerical results for the stability threshold are obtained, and comparedwith various approximations. For the periodic problem, a key feature forthe determination of the stability threshold is to use an Ewald summationmethod to derive an explicit expression for the regular part of the BlochGreen function. Moreover, such an expression allows for the identificationof the lattice that offers the optimum stability threshold. For the finitedomain problem, we implement our asymptotic theory by calculating thestability threshold for an N -spot pattern where the spots are equidistantlyspaced on a circular ring that is concentric within the unit disk.iiPrefaceThis thesis is original, unpublished, independent work by the author, Y.Chang under the supervision Dr. M. Ward and Dr. J. Wei.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . viii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 Preliminaries: The Bravais Lattice and the Bloch Green’sFunction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.1 Bravais Lattice . . . . . . . . . . . . . . . . . . . . . . . . . . 62.2 Reciprocal Lattice . . . . . . . . . . . . . . . . . . . . . . . . 82.3 Bloch Theorem and Bloch Green Function . . . . . . . . . . 123 Periodic Spot Patterns for the Brusselator . . . . . . . . . . 193.1 Periodic Spot Solutions . . . . . . . . . . . . . . . . . . . . . 193.2 Linear Stability Analysis . . . . . . . . . . . . . . . . . . . . 263.3 A Quick Derivation of the Stability Threshold . . . . . . . . 354 Spot Patterns for the Brusselator on a Finite Domain . . 384.1 The N -Spot Solutions . . . . . . . . . . . . . . . . . . . . . . 394.2 Linear Stability Analysis . . . . . . . . . . . . . . . . . . . . 43ivTable of Contents4.2.1 λ 6= 0 and λ ∼ O(1) . . . . . . . . . . . . . . . . . . . 444.2.2 λ ∼ O(ν) and λ 6= 0 . . . . . . . . . . . . . . . . . . . 484.2.3 λ=0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 525 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . 555.1 Small S Asymptotics of χ(S, f) . . . . . . . . . . . . . . . . 555.2 Stability Threshold and the Optimal Lattice Arrangement . 575.3 Case Study: N Peaks on a Ring . . . . . . . . . . . . . . . . 656 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71vList of Tables5.1 Source strength threshold and its asymptotic approximationfor a regular hexagonal lattice with |Ω| = 1 and f = 0.4. . . . 655.2 The stability threshold in terms of the source strength S andits one- and two-term asymptotic approximation for a 5 spotpattern on a ring of radius 0.5 concentric within the unit diskwith f = 0.4. . . . . . . . . . . . . . . . . . . . . . . . . . . . 66viList of Figures2.1 Wigner-Seitz primitive cell . . . . . . . . . . . . . . . . . . . . 83.1 The Continuous Band of Spectra . . . . . . . . . . . . . . . . 345.1 Numerical solution (bottom curves) and asymptotic results(top curves) for χ(S, f). In the left panel we fix f = 0.4,while f = 0.5 for the right panel. In both pictures, the blue(bottom) curve is the numerical solution while the green (top)one is the two term asymptotic expansion. . . . . . . . . . . . 575.2 Numerical solution to (5.5) and the two-term asymptotic ap-proximations for Sc with different c. Left panel: f = 0.4 and = 0.01. Right panel: f = 0.5 and  = 0.05. The blue (top)curve is the numerical solution while the green (bottom) oneis the asymptotic approximations in both cases. . . . . . . . . 595.3 5 localized spots on a ring concentric within the unit disk. . . 65viiAcknowledgementsI would like to thank Dr. Michael Ward and Dr. Juncheng Wei for theirgreat guidance and generous support. I also want to thank all my friendshere at UBC who have made these two years a wonderful experience for me.In the end I want to thank my parents for their love through all these years.viiiChapter 1IntroductionInspired by the work of Allan Turing [14] in 1952, there has been much ef-fort over the past five decades in trying to characterize various patterns thatappear in the physical world through the modeling and analysis of reaction-diffusion (RD) systems. In [14], Turing proposed a mechanism for biologicalmorphogenesis based on his analytical study of a coupled two-componentsystem of reaction-diffusion equations with very different diffusion coeffi-cients. Using a linear stability analysis, he found that a small perturbationto a spatially homogeneous initial data can develop into certain spatial pat-terns through an instability. In the current literature, this type of instabilityis now referred to as a Turing instability. Since then, various RD systemshave been proposed and analyzed to model both biological and chemicalpatterns. Such systems include the Gray-Scott model (cf. [7]), and theGierer-Meinhardt system (cf [6]).Spatially localized spot patterns have been observed both in chemicaland numerical experiments. A survey of such patterns is given in [15].Over the past decade there has been a considerable focus on developinga theoretical mathematical framework for the study of localized patterns forsingularly perturbed reaction diffusion systems for which the ratio O(−2) ofthe two diffusivities is asymptotically large. The work of [17] gives a reviewon the existence, classification and stability of multiple-peaked solutionsfor the Gierer-Meinhardt system on an interval I ⊂ R1. The study [16]generalizes the results to spike solutions in a finite domain Ω ⊂ R2. In bothpapers, the existence of the multi-peaked solutions is proved rigorously byusing the Lyapunov-Schmidt reduction method, while the stability results onthe so-called large eigenvalues with λ = O(1) as → 0 are based on the studyof certain nonlocal eigenvalue problems. In [3] formal singular perturbation1Chapter 1. Introductiontechniques, based on the method of matched asymptotic expansions, areused to analyze the stability of spot patterns for the 2-D Gray-Scott system.In this work, the slow dynamics of the spot patterns is also characterized.The study [10] analyzes the self-replication process of spot patterns for theSchnakenburg model. A formal asymptotic analysis is used to derive an ODEsystem describing the slow dynamics of the spot patterns, which has thesame effect of summing infinite logarithmic series in powers of ν = −1/ ln .In this thesis, we construct localized spot solutions to the 2-D Brusselatorin the semi-strong interaction regime characterized by an asymptoticallylarge diffusivity ratio. We then study the stability of these localized patterns.The Brusselator, proposed by Prigogine and co-workers in Brussels in 1960s,is a theoretical model for a hypothetical autocatalytic reaction (cf. [12]). Thestandard form of this model can be written in terms of non-dimensional spacevariables asUT = 20∆U + E − (B + 1)U + U2V , VT = D∆V +BU − U2V , (1.1)where 20 = Du/L2, D = Dv/L2, L is a characteristic length-scale, while Duand Dv are the diffusivities of U and V . Many different patterns have beenobserved for this model through full scale numerical studies of the PDE, andvia Turing-type stability analysis augmented by weakly nonlinear theoriesfor the evolution of small amplitude patterns.Our goal in this thesis is to construct localized spot solutions underthe singularly perturbed limit 0 → 0 when E = O(0), while maintainingan asymptotically large diffusivity ratio Dv/Du = O(−20 ), so that D =Dv/L2 = (Dv/Du)(Du/L2) = O(1). Then, upon writing E = 0E0 withE0 = O(1), T = t/(B + 1), U = Bu/0, and V = 0v as in [13], thesingularly perturbed RD system becomesut = 2∆u+ 2E − u+ fu2v , τvt = D∆v +12(u− u2v) . (1.2)2Chapter 1. Introductionwhere we have defined f , τ , E, D, and , byf ≡BB + 1, τ ≡1f2, E ≡E0B, D ≡D(B + 1)B2,  ≡0√B + 1,This system has three independent bifurcation parameters E, D, and f ,depending on D, B, and E0. From the definition of f , we observe thatf ∈ (0, 1).In order to solve the system (1.2), we still need some information aboutthe domain and boundary conditions. We will consider two cases in thisthesis; periodic solutions on R2 and multi-spot solutions on a finite domain.For the periodic case, we look for periodic solutions to (1.2) with respectto some Bravais lattices, i.e.u(x+ li) = u(x), v(x+ li) = v(x), i = 1, 2, (1.3)where li, i = 1, 2 are two Bravais vectors. The study [8] undertakes a similaranalysis for the Gierer-Meinhardt and Schnakenburg reaction-diffusion sys-tems. Some basic facts concerning Bravais lattices, the Wigner-Seitz prim-itive cell, reciprocal lattices and Brillouin zones are reviewed in Chapter 2.Due to the periodicity, we can avoid solving the system on the whole plane.Instead, to construct an equilibrium, or steady-state solution, we need onlyconsider a Wigner-Seitz primitive cell together with periodic boundary con-ditions on the boundary of the cell. After constructing periodic solutions inthe primitive cell using the method of matched asymptotic expansions, weperturb this solution and perform a linear stability analysis. The solutionswe construct are stable for small D, unstable for D large, and the stabilitythreshold occurs when D ∼ O(ln 1 ). To calculate the stability threshold, weneed to find the spectrum of the singularly perturbed eigenvalue problemon the whole plane, which is equivalent to finding the eigenfunction of theoperator ∆ + V (x), where V (x) is a periodic function with respect to theBravais lattices. In this way, the Floquet-Bloch theory arises naturally inthe formulation of the stability problem. We prove some basic facts regard-ing the Bloch Green function in Section 2.3. According to the Floquet-Bloch3Chapter 1. Introductiontheory, instead of solving the perturbed system on the whole plane, we onlyneed to solve it within the the Wigner-Seitz primitive cell with the Blochquasi-periodic boundary conditions, which involve a Bloch vector k (notunique) in the first Brillouin zone. By analyzing this problem using themethod of matched asymptotic expansions, we derive a nonlocal eigenvalueproblem. From this problem we obtain the leading order result for the sta-bility threshold, which is independent of the geometry of the lattice andthe Bloch vector. To determine a higher-order approximation for the stabil-ity threshold, we perform a more refined perturbation analysis in order tocalculate a real-valued band of continuous spectrum lying within an O(ν)neighbourhood of the origin in the spectral plane. This band does dependon the geometry of the lattice and the Bloch vector k. For a given lattice,we determine the next term in the stability threshold from the requirementthat the rightmost edge of the real-valued continuous band of spectrum liein the left plane {λ|Re(λ) ≤ 0} for any Bloch vector k in the first Brillouinzone. This calculation involves minimizing the regular part of the BlochGreen’s function. Then, we can determine the optimal lattice arrangementwhich allows for the largest stability threshold. Overall, the identificationof the optimal stability threshold is through a min-max argument. Finally,in addition to this detailed way to calculate the stability threshold, whichinvolves calculating a continuous band of spectrum near the origin of thespectral plane, we also give a quick, but formal, alternative derivation of thefirst two terms of the stability threshold in Section 3.3.Our second problem concerns the analysis of the multi-spot patterns to(1.2) on a finite domain Ω ⊂ R2 with the no-flux boundary conditions:∂nu(x) = 0 , ∂nv(x) = 0 x ∈ ∂Ω . (1.4)Here n is the outer normal vector on the boundary. We focus on both theexistence and linear stability of multi-spot patterns for this problem. Thestudy [13] analyzes a similar problem on a sphere for the Brusselator. Weconstruct asymptotic spot solutions whose u component concentrates on Ngiven points x1,x2, ...,xn ∈ Ω. For simplicity, we will only consider the case4Chapter 1. Introductionfor which these N spots have a common height. As shown in [16] for therelated Gierer-Meinhardt system, the true steady-state positions of theseN points are not arbitrary but rather are close to the critical point of acertain objective function related to the Neumann Green function. It isanticipated that a similar result should hold for the true steady-state spotlocations for the Brusselator. In our approach, we focus not only on steady-state solutions, but also on quasi steady-state solutions that can persistover very long time intervals. Our only key assumption is that the spotpattern has sufficient symmetry so that the vector e = (1, 1, ..., 1)T is aneigenvector of a certain Neumann Green’s matrix. When this requirementis met, there is a common local behavior near each of the spots. Afterconstructing asymptotic solutions that have this symmetry, we introduce aperturbation and perform a linear stability analysis of these patterns as forthe periodic case. The key difference between the periodic and finite domainproblems is that, instead of solving the system only within one primitivecell in the periodic case, the interaction of the spots for the finite domainproblem arises through a Neumann Green matrix. A detailed calculation ofthe stability threshold for this problem, based on the method of matchedasymptotic expansions, involves the eigenvalues and eigenfunctions of theNeumann Green matrix in an essential way.The last chapter of this thesis is concerned with performing a few nu-merical experiments. Firstly, we numerically calculate a key term χ(S, f),which appears in the boundary conditions when solving the inner core prob-lems near a spot. The numerical results for this quantity are comparedwith asymptotic approximations for it that are derived in the small sourcestrength limit S  1. Then we calculate the stability threshold directlyfrom a nonlinear algebraic equation, and we compare the results with one-and two-term asymptotic approximations. For the periodic problem, we alsoshow that a regular hexagonal lattice of spots offers the optimum stabilitythreshold. Finally, we consider a case study for the finite domain problem,in which N spots are equidistantly-spaced on a circular ring that is concen-tric with the unit disk. For this configuration, a refined approximation forthe stability threshold is calculated.5Chapter 2Preliminaries: The BravaisLattice and the BlochGreen’s Function2.1 Bravais LatticeIn this section we review some basic definitions and results regarding theBravais lattice. The Bravais lattice was first introduced to describe theperiodicity of crystalline solids. The lattice points of a Bravais lattice in Rdcan be represented asΛ = {d∑i=1nili|ni ∈ Z, li are linear independent in Rd, i = 1, 2, ..., n} ,(2.1)so that the lattice looks exactly the same when translated by any integer lin-ear combination of li, which are called lattice vectors. We remark here thatthe choice of lattice vectors is not unique. In fact, any linear transformation{∑nj=1 aijlj}ni=1 with det(aij) = ±1 will give the same lattice.Definition 2.1.1. A primitive cell of a Bravais lattice is the smallest regionwhich when translated by all different integer linear combinations of latticevectors can cover the entire space without overlapping. The Wigner-Seitzprimitive cell of a lattice point is a special primitive cell consisting of allpoints in space that are closer to this lattice point than to any other latticepoint.As we can see from the definition, the Wigner-Seitz primitive cell of a62.1. Bravais Latticecertain lattice point is unique while a primitive cell is not. We can determinethe Wigner-Seitz primitive cell of a certain lattice point by identifying thesmallest region enclosed by all hyperplanes which perpendicularly bisect theline segments between this point and any other lattice point (it is sufficientto only consider the nearby points around the chosen point).Now we consider the lattice on R2. A Bravais lattice Λ on R2 is generatedby two lattice vectors l1, l2 which are not parallel to each other. Withoutlose of generality, we may assume l1 is (1, 0), and we can always choosel2 = (a, b), such that a ∈ (0, 1], b > 0 by adding some kl1, for k ∈ Z, to theoriginal l2. Then we consider the Wigner-Seitz primitive cell of the origin Oand show that it is either a hexagon or a rectangle. First of all, due to thecentral symmetry of the lattice, we know that the boundaries of the Wigner-Seitz primitive cell will come in pairs and so we only need to look at theright-hand side. Next, after some easy observation, it follows that b1 and b2,which is the perpendicular bisector of ±l2, form part of the boundaries of theWigner-Seitz primitive cell. The remainder of the boundary arises from someof the perpendicular bisectors of the line segments between the origin and thepoints of the second column, i.e. l1 + kl2, k ∈ Z. It is useful to observe themid-points of these line segments lie on a line parallel to l2 passing throughthe mid-point of OA and the distance between two mid-points next to eachother is |l2|2 . Since b1, b2 are perpendicular to l2, and the distance betweenthem is |l2|, then for most of the cases only two of these mid-points will liebetween them and the corresponding perpendicular bisectors intersectingwith b1, b2 form the boundaries of the primitive cell. This is the genericcase when the Wigner-Seitz primitive cell is a hexagon. A special case iswhen one mid-point lies in the middle of b1, b2 and the two mid-points nextto it on b1, b2, Then, the perpendicular bisectors corresponding to thesethree mid-points is the same line which is perpendicular to b1, b2. Thus, inthis case the Wigner-Seitz primitive cell will be a rectangle. In conclusion,if we use the polar coordinate to denote the angle between l1 and l2 by θ,and the length of l2 by r, then by the assumption we made on l2, we willhave θ ∈ (0, pi2 ), r ∈ (0,1cos θ ) and the Wigner-Seitz primitive cell will be ahexagon unless r = cos θk , k ∈ Z, k > 0. In this latter, degenerate case, the72.2. Reciprocal Latticeprimitive cell will be a rectangle.−1.5 −1 −0.5 0 0.5 1 1.5 2−1−0.500.511.5Figure 2.1: Wigner-Seitz primitive cell2.2 Reciprocal LatticeNext, we review the concept of the reciprocal lattice Λ∗ of a certain Bravaislattice. This notion arises from Fourier analysis. Firstly, we define a periodicfunction with respect to a Bravais lattice as follows:Definition 2.2.1. Given a Bravais lattice Λ in Rd with lattice vectors{li}di=1, a function u(x) is periodic with respect to the lattice Λ if u(x+li) =u(x), i = 1, 2, ..., d.We know that a 2pi periodic function u(x) on R can be decomposed intoFourier series as einx, n ∈ Z. According to Definition 2.2.1, we can viewu(x) as a periodic function with respect to the 1-D lattice Λ = {nl1|n ∈Z, l1 = 2pi} and the Fourier basis einx, n ∈ Z can be viewed as eik·x, k ∈Λ∗ = {nl∗1|n ∈ Z, l∗1 = 1}, where Λ∗ is another Bravais lattice relatedto Λ, which will be defined as the reciprocal lattice of Λ later. In higherdimensions, we first consider a 2pi periodic (in each direction) function u(x),i.e. u(x + 2piei) = u(x), i = 1, 2, ..., d, where ei are the standard basis ofRd. We know that u(x) can be decomposed into Fourier series of ein·x,n ∈ Zd. As explained above, u(x) now can be viewed as a periodic functionwith respect to the Bravais lattice Λ = {∑di=1 nili|ni ∈ Z, li = 2piei, i =1, 2, ..., d}. In addition, the Fourier basis ein·x, n ∈ Zd can be viewed as82.2. Reciprocal Latticeeik·x, k ∈ Λ∗ = {∑di=1 nil∗i |ni ∈ Z, l∗i = ei}. For the general case, if u(x) isa periodic function with respect to a Bravais lattice Λ with lattice vectors{li}di=1, we want to decompose it as we did above. In order to convert to theprevious case, we use a linear coordinate change x = 12pi [l1|...|ld]y and denoteA = 12pi [l1|...|ld]. Then u˜(y) = u(Ay) is a 2pi periodic function in each yidirection and can be decomposed into Fourier series as ein·y = ein·(A−1x) =ei((A−1)Tn)·x, n ∈ Zd. If we denote (A−1)T = [b1|...|bd], then (A−1)Tn,n ∈ Zd, represents a Bravais lattice with lattice vector {bi}d1. And sinceAT (A−1)T = I, we have the relation ( 12pi li) · bj = δij , so that li · bj = 2piδij ,which leads to the definition of reciprocal lattice.Definition 2.2.2. The reciprocal lattice Λ∗ of a Bravais lattice Λ with latticevectors {li}di=1 is the Bravais lattice given by the lattice vectors {bi}di=1 wherebi are the vectors satisfy li · bj = 2piδij.We remark that { 12pibi}di=1 are the dual basis for {li}di=1 in Rd. Thus,{bi}di=1 are well defined and uniquely determined by {li}di=1. We furtherremark that other authors choose the reciprocal lattice vectors {bi}di=1 tosatisfy li · bj = δij .Some important properties concerning reciprocal lattices Λ∗ of a Bravaislattice Λ are listed here as follows:• Λ∗ = {b ∈ Rd|eib·l = 1, ∀l ∈ Λ}• The reciprocal lattice of Λ∗ is Λ.• Any periodic function u(x) with respect to Λ can be decomposed intoFourier series of eib·x, where b ∈ Λ∗. Using a change of variable, wehave the formula:u(x) =1|Ω|∑b∈Λ∗(∫Ωu(y)e−ib·y dy)eib·x . (2.2)• The Fourier transform of∑l∈Λ δ(x− l) is(2pi)d|Ω|∑b∈Λ∗ δ(ξ− b), whereΩ is the Wigner-Seitz primitive cell of Λ.92.2. Reciprocal Lattice• If u(x) ∈ L1(Rd), then∑l∈Λ u(x + l) converges absolutely almosteverywhere, and we have:∑l∈Λu(x+ l) =1|Ω|∑b∈Λ∗uˆ(b)eix·b . (2.3)In particular, upon replacing u(x) by u(x)eix·k, we have that∑l∈Λu(x+ l)eik·l =1|Ω|∑b∈Λ∗uˆ(b− k)eix·(b−k) . (2.4)We remark that the last two properties are called Poisson summation for-mulae. We prove these properties as follows:Proof. First we establish the basic result that the Fourier transform of∑l∈Λ δ(x − l) is(2pi)dV∑b∈Λ∗ δ(ξ − b). The key step in the proof is thefollowing equality:∞∑n=−∞einx = 2pi∞∑k=−∞δ(x− 2pik) . (2.5)The proof of (2.5) can be found in many books; we simply sketch the outlineof the proof here. First of all, the equality holds in the sense of distribution,and since both sides are 2pi periodic, we may prove it only on the interval[−pi, pi], i.e. for any test function φ(x) ∈ C∞0 [−pi, pi],limN→∞∫ pi−pi(N∑n=−Neinx)φ(x) dx =∫ pi−pi(2pi+∞∑k=−∞δ(x− 2pik))φ(x) dx = 2piφ(0) .(2.6)This is true sinceN∑n=−Neinx =sin (N + 12x)sin x2→ δ(x) , as N →∞ ,in the sense of distributions.102.2. Reciprocal LatticeNotice that the 1-D case followed directly from this equality,̂∞∑n=−∞δ(x− nl) =∞∑n=−∞ein·ξl = 2pi∞∑k=−∞δ(lξ−2pik) =2pil∞∑k=−∞δ(ξ−k2pil) .For the higher dimensional case, we havê∑l∈Λδ(x− l) =∑mi∈Ze−i(∑di=1mili)·ξ =d∏i=1(∞∑mi=−∞e−imi(li·ξ)) (2.7)=d∏i=1(2pi∞∑ki=−∞δ(li · ξ − 2piki)) =d∏i=1(∞∑ki=−∞δ(li · ξ2pi− ki)) .If we denote ξ =∑di=1 ηibi = [b1|b1|...|bd]η = Bη, where {bi}di=1 are thereciprocal vectors to {li}di=1, then (2.7) becomesd∏i=1(∞∑ki=−∞δ(li · ξ2pi− ki)) =d∏i=1(∞∑ki=−∞δ(ηi − ki)) =∑k∈Zdδ(η − k) (2.8)=∑k∈Zdδ(B−1(ξ −Bk)) =1det(B−1)∑b∈Λ∗δ(ξ − b) .SinceB−1 = 12pi [l1|l1|...|ld], then det(B−1) = V(2pi)d , where V = det([l1|l1|...|ld])is the volume of the primitive cell. In this way, we conclude that̂∑l∈Λδ(x− l) =(2pi)dV∑b∈Λ∗δ(ξ − b) . (2.9)Then we prove the last property. First of all since∫Ω∑l∈Λ|u(x+ l)| dx =∫R2|u(x)| dx <∞ , (2.10)the series∑l∈Λ u(x+ l) converges absolutely almost everywhere. Moreover,since the series is periodic with respect to the lattice Λ, then by using the112.3. Bloch Theorem and Bloch Green Functionproperty above we can decompose it into a Fourier series as∑l∈Λu(x+ l) =1|Ω|∑b∈Λ∗(∫Ω(∑l∈Λu(y + l))e−ib·y dy)eib·x . (2.11)The last step in the derivation is to calculate the Fourier coefficients as∫Ω(∑l∈Λu(y + l))e−ib·y dy =∑l∈Λ∫Ωu(y + l)e−ib·y dy (2.12)=∫R2u(y)e−ib·y dy = uˆ(b) , (2.13)where the the second to last equality holds since R2 is tiled when Ω trans-lated by all lattice vector, and since ∀l ∈ Λ and ∀b ∈ Λ∗, we have eil·b = 1.In particular when we replace u(x) by u(x)eix·k, the Fourier transformis translated by k.There are two further useful concepts that relate to Bravais lattices.Definition 2.2.3. Bragg planes are the hyperplanes which perpendicularlybisect any line segment between two lattice points of Λ∗. The first BrillouinZone is the Wigner-Seitz primitive cell of Λ∗2.3 Bloch Theorem and Bloch Green FunctionIn this section, we review Bloch theorem and some properties of the BlochGreen function in the Wigner-Seitz primitive cell Ω of some lattice Λ. Theseresults will be used later when we consider the stability of a periodic patternof spots for the reaction-diffusion system. The proof below is similar to thatin [8].The Bloch theorem originates from quantum mechanics and states thatthe eigenfunction φ(x) of the operator −∆ + V (x), where the potentialfunction V (x) is periodic with respect to a Bravais lattice Λ, must havethe form φ(x) = eik·xφp(x), where φp(x) is also periodic with respect tothe lattice Λ, and k can be chosen to lie in the first Brillouin Zone, orequivalently ∀l ∈ Λ, φ(x + l) = eik·lφ(x). The Bloch theorem allows us122.3. Bloch Theorem and Bloch Green Functionto solve the eigenvalue problem within the primitive cell together with theBloch boundary conditions, instead of on the whole space. The proof ofthe Bloch theorem can be found in many solid physics books and the keyidea is that if two operators commute, which in this case are −∆ + V (x)and translation by any lattice vector, they share common non-degenerateeigenvectors while the eigenvalues may be different. By using this idea, wecan prove the Bloch theorem for a system, which is the basis for our linearstability analysis. Since the equilibrium solutions we construct are periodicwith respect to some Bravais lattice Λ, finding the eigenspace of the linearperturbed system is equivalent to finding the eigenfunction for an operatorsimilar to −∆ + V (x). In this way, the Bloch theorem arises naturally inthe stability analysis of a periodic arrangement of spots.Next we consider some key properties of the Bloch Green function. Thisis the Green function in the fundamental Wigner-Seitz cell that satisfies theBloch boundary conditions. As we have shown above, the primitive cell Ωwill be either a hexagon or a rectangle. We may assume that the boundariesof Ω consist of d±i, with i ≤ 2 for a rectangle and i ≤ 3 for a hexagon, whered±i represents the perpendicular bisector of ±Li ∈ Λ which come in pairs.The Bloch Green function in Ω is the solution to∆G0,k(x) = −δ(x), x ∈ Ω, (2.14)and satisfies the following Bloch boundary conditions, also referred to asquasi-periodic boundary conditions:∀x ∈ d−i, G0,k(x+Li) = e−ik·LiG0,k(x), (2.15)∂n−G0,k(x+Li) = e−ik·Li∂−n+G0,k(x),where the ± behind n in the directional derivatives denote one-side deriva-tives, n is the outer unit normal vector parallel to Li and k is some non-zerovector in the first Brillouin Zone Ω∗, i.e. in the Wigner-Seitz primitive cellof Λ∗. We first remark here that we require k 6= 0 since there is no solutionto (2.14) if k = 0. This is shown by integrating ∆G0,k(x) over Ω, and us-132.3. Bloch Theorem and Bloch Green Functioning the divergence theorem which results in a contradiction. The boundaryconditions are well-defined since x + Li ∈ di if x ∈ d−i. Moreover, we useone-side normal derivative since we only solve the equation inside Ω. Withthis boundary conditions, we can extend the solution to the whole planewith continuous normal derivative between contiguous cells.It is also useful to analyze the quasi-periodic reduced-wave Green’s func-tion, which is the solution to∆Gσ,k(x)− σ2Gσ,k(x) = −δ(x) , x ∈ Ω , σ ∈ R , (2.16)with the boundary conditions of (2.15). We observe that the Bloch Greenfunction is simply the special case of the reduced-wave Green’s functionwhen σ = 0. The first key property is that the regular part Rσ,k of Gσ,k(x),which is defined by subtracting the free space Green function − 12pi ln |x| fromGσ,k(x), i.e.Rσ,k = limx→0Gσ,k(x) +12piln |x| ,is real for k 6= 0.Lemma 1. The regular part Rσ,k of Gσ,k(x) is real-valued for |k| 6= 0.Proof. We let  > 0 and eliminate the singularity by cutting a small ballB(0, ) around the origin. We denote Ω ≡ Ω − B(0, ). Then we use thedivergence theorem and the fact that ∆Gσ,k(x) = σ2Gσ,k(x) for x ∈ Ω tocalculate∫∂ΩGσ,k∂nGσ,k dl =∫∂ΩGσ,k(∇Gσ,k · n) dl =∫∂Ω((Gσ,k∇Gσ,k) · (n dl)=∫Ω∇ · (Gσ,k∇Gσ,k) dx =∫Ω(∇Gσ,k · ∇Gσ,k +Gσ,k∆Gσ,k)dx=∫Ω(|∇Gσ,k|2 + σ2|Gσ,k|2) dx ,where n is the outer normal vector as usual. Upon calculating the boundary142.3. Bloch Theorem and Bloch Green Functionintegral directly, we have∫∂ΩGσ,k∂nGσ,k dl =∫∂ΩGσ,k∂nGσ,k dl −∫∂B(0,)Gσ,k∂nGσ,k dl . (2.17)We claim that the first term vanishes due to the Bloch boundary conditions(2.15). Since ∂Ω has either 4 or 6 perpendicular bisectors, which come inpairs, then according to (2.15), on each pair d±i, Gσ,k(x) is different by afactor of e−ik·Li and ∂nGσ,k is different by −e−ik·Li owing to the fact thatthe outer normal vectors are opposite. Then, if we integrate Gσ,k∂nGσ,k ond±i, these two terms will be different by e−ik·Li(−e−ik·Li) = −1. Thus, thefirst integral on the right hand-side of (2.17) banishes. Next, we calculatethe second term on the right hand-side of (2.17) as∫∂B(0,)Gσ,k∂xGσ,k dx ∼ ∫ 2pi0(−12piln +Rσ,k + o(1))(−12pi+O(1)) dθ∼12piln −Rσ,k +O( ln ) .In this way, we let → 0 to obtainRσ,k = lim→0[ ∫Ω(|∇Gσ,k|2 + σ2|Gσ,k|2) dx+12piln ], (2.18)which proves that Rσ,k is real-valued.As we mentioned before, there is no solution to (2.14) if k = 0. So thenext lemma will discuss the asymptotic behaviour of R0,k when k tends to0. To analyze this limiting behavior, we may assume k = σκ, where σ  1and |κ| = O(1). This form suggests that we can calculate a solution to(2.14) by a singular perturbation technique. As such, we expand G0,k asG0,k(x) =U0(x)σ2+U1(x)σ+ U2(x) + · · · . (2.19)In addition, the expansion of the boundary conditions (2.15) yield for ∀x ∈152.3. Bloch Theorem and Bloch Green Functiond−i thatU0σ2+U1σ+ ...∣∣x+Li= [1− iσ(κ ·Li) + ...](U0σ2+U1σ+ ...)∣∣x∂n−U0σ2+∂n−U1σ+ ...∣∣x+Li= − [1− iσ(κ ·Li) + ...](∂−n+U0σ2+∂−n+U1σ+ ...)∣∣x,where n is the outer normal vector on the boundary, which is parallel to Li.Then by equating the same order of σ, we get:O(−2) : ∆U0(x) = 0 , ∀x ∈ Ω ,∀x ∈ d−i , U0(x+Li) = U0(x) , ∂n−U0(x+Li) = ∂−n+U0(x) .O(−1) : ∆U1(x) = 0 , ∀x ∈ Ω ,∀x ∈ d−i , U1∣∣x+Li= U1 − i (κ ·Li)U0∣∣x ,∂n−U1∣∣x+Li= ∂−n+U1 − i (κ ·Li) ∂−n+U0∣∣x ,O(1) : ∆U2(x) = −δ(x) , ∀x ∈ Ω ,∀x ∈ d−i , U2∣∣x+Li= U2 − i (κ ·Li)U1 −(κ ·Li)22U0∣∣x ,∂n−U2∣∣x+Li= ∂−n+U2 − i (κ ·Li) ∂−n+U1 −(κ ·Li)22∂−n+U0∣∣x .From the leading order equation we conclude that U0 = a for some constanta. Form the O(−1) equation, we get that U1 is a linear function of the formU1(x) = (−iaκ) ·x+ b, where b is another constant. Next we integrate ∆U2over Ω defined by the O(1) problem. Upon using the boundary conditionsand the expression for U0 and U1, as derived above, we obtain that∫Ω∆U2 dx =∫∂Ω∂nU2(x) dl =L∑j=1∫±dj∂nU2(x) dl =∫Ω−δ(x) dx , ⇒L∑j=1−i(κ ·Lj)n · (−iaκ)|di| = −aL∑j=1(κ ·Lj)κ ·Li|Li||di| = −1 ,where L = 2 if the primitive cell is a rectangle, and L = 3 if the primitive162.3. Bloch Theorem and Bloch Green Functioncell is hexagon. In this way, we determine a asa =1∑Lj=1(κ ·Lj)2 |di||Li|=1κT (∑Lj=1|di||Li|LjLTj )κ=1κTQκ, (2.20)where we denote Q =∑Lj=1|di||Li|LjLTj . We remark that Q is a positivedefinite matrix related only to the lattice. This leads to our second lemma:Lemma 2. For sufficiently small k, the regular part R0,k of the Bloch Greenfunction G0,k has the asymptotic behaviourR0,k ∼1kTQk, (2.21)where Q is a positive definite matrix related to the lattice Λ as shown above.We remark here that when we try to find the singular perturbed solutionG0,k for k ∼ O(), the reason we expand G0,k from O(−2) is because it can’tsatisfy the boundary conditions if we start from O(−1).Next we state two other similar results which will be useful later.Lemma 3. The regular part Rσ,k of the reduced-wave Bloch Green functionGσ,k has the asymptotic behaviour Rσ,k ∼ R0,k + O(σ2), as σ → 0 and|k| ∼ O(1), where R0,k is the same as above.Proof. Just expand Gσ,k = G0 + σ2G1 + .... Then, we have G0 is the BlochGreen function G0,k. And since G1 is bounded for |k| ∼ O(1), we have theasymptotic behaviour Rσ,k ∼ R0,k +O(σ2).Lemma 4. The regular part Rσ,k of the reduced-wave Bloch Green functionGσ,k has the asymptotic behaviour Rσ,k ∼ 1σ2[|Ω|+κTQκ] , as σ → 0 and |k| ∼O(σ), where Q is the same positive definite matrix determined by the latticeΛ.Proof. The proof is basically the same as we do for R0,k, the only differenceis the equation for U2 became ∆U2 = U0−δ(x). So when we integrate ∆U2,U0|Ω| will appear.172.3. Bloch Theorem and Bloch Green FunctionSimilar to the Bloch boundary conditions (2.15), we will also use theperiodic boundary conditions:G(x+ li) = G(x), ∂n−G(x+ li) = ∂−n+G(x), ∀x ∈ d−i, (2.22)18Chapter 3Periodic Spot Patterns forthe BrusselatorIn this chapter we construct periodic spot solutions to (1.2) by first con-structing a single spot solution inside the Wigner-Seitz primitive cell Ω sub-ject to the periodic boundary conditions (2.22). We then periodically extendthis solution to the whole plane. Next, we analyze the linear stability of theequilibrium solutions. We first perturb the steady state solution, to derive asingular perturbed eigenvalue problem governing the linear stability of theperiodic pattern. From this eigenvalue problem we then provide an accu-rate calculation of the stability threshold corresponding to a zero eigenvaluecrossing. We also provide a more expedient approach to derive the stabilitythreshold.3.1 Periodic Spot SolutionsThe Brusselator reaction-diffusion model has the formu˜t = 2∆u˜+ 2E − u˜+ fu˜2v˜ ,τ v˜t = D∆v˜ +12(u˜− u˜2v˜) .From this first equation we observe that there is a spatially homogeneousequilibrium solution for  1 with u˜ ∼ 2E, and so we make a substitutionof the form u˜ = u+ 2E and v˜ = v. Then the system becomesut = 2∆u− u+ f(u2v + 22Euv + 4Ev) , (3.1)τvt = D∆v + E + −2(u− u2v)− 2Euv − 2E2v .193.1. Periodic Spot SolutionsSince the steady-state problem is singularly perturbed, we need to constructa localized asymptotic expansion solutions in the inner region near the originof the fundamental Wigner-Seitz cell. Then we use the method of matchedasymptotic expansions to match the inner solution to an appropriate outersolution.In the inner region, we use an inner coordinate and perform a scaling toeliminate the D. The inner variables y, U , and V , are defined byy = −1x , U(y) =u(x)√D, V (y) =√Dv(x) . (3.2)To motivate this scaling of U and V , we remark that if we assumed u ∼ DαU ,v ∼ DβV , then from the first equation we have that α = −β = k, while fromthe second equation we conclude that D · O(D−k) ∼ O(Dk). This leads tok = 1/2.In the inner region, to leading order the solution is radially symmetric.As such, if we denote ρ = |y| then to leading order we get the core problem∆ρU − U + fU2V = 0 , ∆ρV + U − U2V = 0 , 0 < ρ <∞ , (3.3)where U and V are functions of ρ. Therefore, (3.3) is an ODE system.It will have solutions when equipped with initial conditions or appropriateboundary conditions. In this case, we add boundary conditions at ρ = 0and an asymptotic condition as ρ→∞:U ′(0) = V ′(0) = 0 ; U → 0 , V ∼ S ln ρ+ χ(S, f) , as ρ→∞ , (3.4)where S is some unknown source strength to be determined, while χ(S, f) issome quantity depending on S and f that is to be computed from the coresolution.We remark here that we choose U ′(0) = V ′(0) = 0 since we are lookingfor differentiable radically symmetric solutions, U → 0 at infinity since wewant localized spot solutions, and we allow V to have logarithmic growth atinfinity since the solution to ∆V = −δ(y) on R2 is − 12pi ln |y|. To derive an203.1. Periodic Spot Solutionsidentity for S, we integrate ∆V over R2 to obtain2piS =∫R2(U2V − U) dx =∫ 2pi0dθ∫ ∞0(U2V − U)ρ dρ ,so thatS =∫ ∞0(U2V − U)ρ dρ . (3.5)Next, we construct the outer solution. Since u is localized, the secondequation in (3.1) reduces asymptotically to0 = D∆v + E +(∫Ω[−2(u− u2v)− 2Euv − 2E2v] dx)δ(x) . (3.6)Upon using (3.5), this equation becomesD∆v + E = 2piS√Dδ(x) , x ∈ Ω . (3.7)Then, if we integrate (3.7) over Ω and use the periodic boundary conditionson ∂Ω, we can calculate the source strength S asS =E|Ω|2pi√D, (3.8)so that (3.7) becomes∆v +ED=E|Ω|Dδ(x) . (3.9)To represent the solution v, we introduce the periodic Green functionGp for Ω as∆Gp =1|Ω|− δ(x) , x ∈ Ω , (3.10)with the periodic boundary conditions (2.22). Since Gp is determined onlyup to a constant, we impose that∫ΩGp dx = 0. As x → 0, Gp has thesingular behaviour Gp ∼ − 12pi ln |x|+Rp, where the regular part Rp can becalculated explicitly for any oblique Bravais lattice, as was shown in [4].In terms of this Green’s function, we have v = −E|Ω|D Gp + c. As x→ 0,213.1. Periodic Spot Solutionsthe limiting behavior of v isv ∼E|Ω|2piDln |x| −E|Ω|DRp + c . (3.11)We now match this behavior with the inner solution. From (3.4), we obtainas |y| → ∞ thatv ∼1√D(S ln |y|) + χ(s, f)) =E|Ω|2piD(lnx+1ν) +χ√D, (3.12)where ν = −1ln   1. Upon matching, we identify the constant c asc =χ√D+1νE|Ω|2piD+E|Ω|DRp . (3.13)Since the stability threshold occurs when D ∼ O( 1ν ), then from ourformula (3.8) for S, we conclude that S ∼ O( 1√D) ∼ O(ν12 ). Therefore, forS = O(ν1/2), we must look for an asymptotic solution to the core problem(3.3). To motivate the scalings for U and V , we observe from (3.12) and(3.3) that V =√Dv ∼ O( 1√Dν) ∼ O(ν−12 ), UV ∼ O(1), which leads toU ∼ O(ν12 ). Since S = O(ν12 ), we need χ(S, f) to be order O(ν−12 ) to matchwith V . These formal scaling arguments motivate the following asymptoticexpansion for the solution to the core problem:U ∼ ν12 (U0 + νU1 + ν2U2 + · · · ) , V ∼ ν− 12 (V0 + νV1 + ν2V2 + · · · ) ,χ ∼ ν−12 (χ0 + νχ1 + ν2χ2 + · · · ) , S ∼ ν12 (S0 + νS1 + ν2S2 + · · · ) .(3.14)Next, we substitute (3.14) into (3.3) and try to construct radially sym-metric solutions at each order. At leading order, we have∆U0 − U0 + fU20V0 = 0 , U′0(0) = 0 , U0 → 0 , as |y| → ∞ ,∆V0 = 0 , V′0(0) = 0 , V0 → χ0 , as |y| → ∞ .(3.15)223.1. Periodic Spot SolutionsFrom this system we conclude thatU0 =ωfχ0, V0 = χ0 , (3.16)where ω is the unique positive radially symmetric solution (see [5]) to∆ω − ω + ω2 = 0 , (3.17)with ω(y) having exponential decay as |y| → ∞. In addition, we have∫∞0 ω(ρ)ρdρ =∫∞0 ω2(ρ)ρ dρ ≡ b upon integrating the equation for ω. Fur-ther properties of this ground-state solution are given in [5].At next order, we have that∆U1 − U1 + f(ω2f2χ20V1 + 2ωfU1) = L0U1 +ω2fχ20V1 = 0 , (3.18)∆V1 +ωfχ0−ω2χ0f2= 0 , (3.19)U ′1(0) = V′1(0) = 0 ; U1 → 0 , V1 → S0 ln |y|+ χ1 , as |y| → ∞ ,where the operator L0 is defined as L0u = ∆u− u+ 2ωu. By applying thedivergence theorem to the V1 equation we obtain thatS0 =bfχ0−bfχ20⇒ χ0 =1− ff2bS0. (3.20)Then the solutions to (3.18) and (3.19) can be decomposed asU1 = −χ1fχ20ω −1f3χ30U1P , V1 = χ1 +V1Pf2χ0, (3.21)where U1P , V1P are the unique radially symmetric solution toL0U1P − ω2V1P = 0 , ∆V1P = ω2 − fω ,U ′1P (0) = V′1P (0) = 0 ; V1P → (1− f)b ln |y|+ o(1) , as |y| → ∞ .(3.22)233.1. Periodic Spot SolutionsTo eliminate the f -dependence in V1P , we introduce V1Q satisfying∆V1Q = ω , V′1Q(0) = 0 ; V1Q → b ln |y|+ o(1) , as |y| → ∞ .(3.23)Then we have∆(V1Q − ω) = ω −∆ω = ω2 ⇒ V1P = (1− f)V1Q − ω . (3.24)After substituting V1Q into (3.22) we getL0U1P = (1− f)ω2V1Q − ω3 . (3.25)This suggests that we should decompose U1P by introducing U1QI and U1QII ,which are taken to satisfyL0U1QI = ω3 , U ′1QI(0) = 0 , U1QI → 0 , as |y| → ∞ , (3.26)L0U1QII = ω2V1Q , U′1QII(0) = 0 , U1QII → 0 , as |y| → ∞ , (3.27)so that U1P = (1− f)U1QII − U1QI .At second order we obtain thatL0U2 + f(ω2f2χ20V2 + χ0U21 + 2ωfχ0U1V1) = 0 ,U ′2(0) = 0 , U2 → 0 , as |y| → ∞ ,∆V2 + U1 − U20V1 −2ωfU1 = 0 ,V ′2(0) = 0 , V2 → S1 ln |y|+ χ2 , as |y| → ∞ .By using the divergence theorem on the V2 equation we calculate S1 asS1 =∫ ∞0U20V1ρ dρ+2f∫ ∞0ωU1ρ dρ−∫ ∞0U1ρ dρ . (3.28)243.1. Periodic Spot SolutionsNext, upon integrating (3.18) we get∫R2L0U1 dx+∫R2ω2fχ20V1 dx = 0 ,⇒ −∫ ∞0U1ρ dρ+∫ ∞02ωU1ρ dρ+∫ ∞0ω2fχ20V1ρ dρ = 0 ,where∫R2 ∆U1 dx vanishes due to the boundary conditions. Upon substi-tuting this back into (3.28), together with (3.26) and (3.27), we obtain thatf1− fS1 =∫ ∞0U1ρ dρ =∫ ∞0(−χ1fχ20ω −1f3χ30((1− f)U1QII − U1QI))ρ dρ= −χ1bfχ20−1− ff3χ30∫ ∞0U1QIIρ dρ+1f3χ30∫ ∞0U1QIρ dρ .We then combine this expression with (3.21) and (3.20) to calculate χ1 asχ1 =fχ20b(−1− ff3χ30∫ ∞0U1QIIρ dρ+1f3χ30∫ ∞0U1QIρ dρ− S1f1− f) ,= −(1− f)bf2S1S20−S0b2∫ ∞0U1QIIρ dρ+S0(1− f)b2∫ ∞0U1QIρ dρ .In conclusion, we have constructed a two-term asymptotic spot solution tothe core problem (3.3) in the limit S → 0. We summarize our result asfollows:Principal Result 1. For S ∼ ν12 (S0 + νS1 + ...), where ν ≡ −1ln  , the coreproblem (3.3) has an asymptotic solution in the formU ∼ ν12 (U0+νU1+· · · ), V ∼ ν− 12 (V0+νV1+· · · ), χ ∼ ν− 12 (χ0+νχ1+· · · ) ,(3.29)where U0(ρ), U1(ρ), V0(ρ) and V1(ρ) are defined by:U0 =ωfχ0, U1 = −χ1fχ20ω −1f3χ30((1− f)U1QII − U1QI) ,V0 = χ0 , V1 = χ1 +1f2χ0((1− f)V1Q − ω) .253.2. Linear Stability AnalysisHere U1QI , U1QII , and V1Q are the unique radially symmetric solutions to(3.26), (3.27), and (3.24), respectively. In addition, ω is the unique positiveradially symmetric solution to (3.17), while χ0 and χ1 are defined byχ0 =(1− f)f2bS0,χ1 = −(1− f)bf2S1S20−S0b2∫ ∞0U1QIIρ dρ+S0(1− f)b2∫ ∞0U1QIρ dρ .3.2 Linear Stability AnalysisIn this section, we perform a linear stability analysis for the equilibrium spotsolution that we have just constructed. We perturb the steady state solu-tions and derive the singularly perturbed eigenvalue problem. By analyzingthis eigenvalue problem we get the stability threshold. As we mentioned inSection 2.3, since the steady state solutions we have just derived is periodicwith respect to the lattice Λ, solving the perturbed system is equivalentto finding eigenfunctions of the operator ∆ + f(x), where f(x) is a 2 × 2periodic matrix with respect to the lattice Λ. As a result of the general-ized Bloch theory, instead of solving the perturbed system on R2, we needonly consider the eigenvalue problem within the primitive cell Ω with theBloch boundary conditions (2.15), which involves a Bloch vector in the FirstBrillouin Zone.We begin perturbing the equilibrium solutions asu = ue + eλtφ , v = ve + eλtη . (3.30)We linearize the equations around this steady state solution and obtain, toleading order, the singularly perturbed eigenvalue problemλ(φτη)=(2∆φD∆η)+(−1 + 2fueve fu2e−2 − 2−2ueve −2u2e)(φη), (3.31)263.2. Linear Stability Analysisfor x ∈ Ω, where φ and η satisfy the Bloch boundary conditions (2.15) forsome k in the first Brillouin Zone.Since the equilibrium solution is localized then, as similar to the con-struction of the equilibrium spot solution, we introduce an inner region nearthe core of the spot. As such, we introduce the following inner coordinatesand scale the system asy = −1x , Φ(y) =φ(x)D=φ(y)D, N(y) = η(x) = η(y) . (3.32)We still look for radially symmetric solutions Φ(y) = Φ(ρ), N(y) = N(ρ),where ρ = |y|. Then, to leading order, the inner system becomes∆ρΦ− (λ+ 1)Φ + 2fUV Φ + fU2N = 0 , Φ′(0) = 0 ; Φ→ 0 as ρ→∞ ,(3.33)∆ρN + (1− 2UV )Φ− U2N = 0 , N ′(0) = 0 ; N → C ln |y|+B(C, f) ,(3.34)as ρ → ∞. We remark here that since this is a linear system for a fixed f ,it follows that the ratio BC is a constant. In addition, we remark that theasymptotic boundary condition for Φ is appropriate provided that Re(λ +1) > 0, i.e. that λ is not in the continuous spectrum. Finally, we obtain thefollowing identity by integrating the equation for N to getC =∫ ∞0(U2N − (1− 2UV )Φ)ρ dρ . (3.35)In terms of C, in the outer region the second equation of (3.31) is ap-proximated as∆η −τλDη = 2piCδ(x) , (3.36)with the Bloch boundary conditions (2.15) for some Bloch vector k. In termsof the reduced-wave Bloch Green function (2.16), we can write η asη = −2piCGµ,k(x) , µ =√τλD∼ O(ν12 ) . (3.37)273.2. Linear Stability AnalysisTherefore, as x→ 0, η(x) has the asymptotic behaviourη(x) ∼ −2piC(−12piln |x|+Rµ,k) ∼ C ln |x| − 2piCRµ,k . (3.38)We know from Lemma 3 that if |k| is bounded away from the origin, thenwe have Rµ,k = R0,k +O(µ2) = R0,k +O(ν).Next, we asymptotically match (3.38) with the far-field behavior of theinner solution (3.34), which is given byη ∼ C ln |x|+Cν+B . (3.39)Upon comparing (3.38) and (3.39), we conclude thatC(1 + 2piνR0,k +O(ν2)) = −νB . (3.40)Next, we will calculate the asymptotic solution to (3.33) and (3.34) forν small. Since we have UV ∼ O(1), U2 ∼ O(ν) and C ∼ νB, this suggeststhat we expand the solutions asΦ ∼ ν(Φ0 + νΦ1 + · · · ) , N ∼ N0 + νN1 + · · · ,B ∼ B0 + νB1 + · · · , C ∼ ν(C0 + νC1 + · · · ) .Upon substituting these expansions into (3.33) and (3.34), we obtain toleading order that∆Φ0 − (λ+ 1)Φ0 + 2fU0V0Φ0 + fU20N0 = 0 , ∆N0 = 0 ,Φ′0(0) = N′0(0) = 0 ; Φ0 → 0 , N0 → B0 as |y| → ∞ .Then, upon substituting the expressions for U0 and V0, as given in PrincipalResult 1, we getL0Φ0 = λΦ0 −B0fχ20ω2 , N0 = B0 . (3.41)283.2. Linear Stability AnalysisFrom (3.40), we conclude thatB0 = −C0 = N0 , B1 = −C1 − 2piR0,kC0 . (3.42)At next order, we have that Φ1 satisfies∆Φ1 − (λ+ 1)Φ1 + 2f(U0V1 + U1V0)Φ0 + 2fU0V0Φ1+ fU20N1 + 2fU0U1N0 = 0 ; Φ′1(0) = 0 , Φ1 → 0 as |y| → ∞ ,(3.43)while N1 satisfies∆N1 + (1− 2U0V0)Φ0 − U20N0 = 0 ,N ′1(0) = 0 , N1 → C0 ln |y|+B1 as |y| → ∞ .Upon substituting the results for U0 and V0 from Principal Result 1, weobtain from the divergence theorem on the N1 equation thatC0 =(1 +bf2χ20)−1 ∫ ∞0(2fω − 1)Φ0ρ dρ . (3.44)With C0 now determined, we substitute it back into the equation for Φ0 toobtainL0Φ0 = λΦ0 +fb+ f2χ20(2f∫ ∞0ωΦ0ρ dρ−∫ ∞0Φ0ρ dρ)ω2 . (3.45)In contrast to the nonlocal eigenvalue problems (NLEP’s) analyzed in[8], this NLEP is more intricate as it involves two nonlocal terms. In order toobtain an NLEP with only one nonlocal term, we integrate (3.45) to derive∫R2L0Φ0dx = 0− 2pi∫ ∞0Φ0ρ dρ+ 4pi∫ ∞0ωΦ0ρ dρ= 2piλ∫ ∞0Φ0ρ dρ+2pibfb+ f2χ20(2f∫ ∞0ωΦ0ρ dρ−∫ ∞0Φ0ρ dρ),293.2. Linear Stability Analysiswhich leads to(λ+ 1−bfb+ f2χ20)∫ ∞0Φ0ρ dρ = (2−2bb+ f2χ20)∫ ∞0Φ0ωρ dρ . (3.46)We substitute this back into (3.45) to getL0Φ0 −2b(λ+ 1− f)(λ+ 1)(b+ f2χ20)− bf∫∞0 Φ0ωρ dρ∫∞0 ω2ρ dρω2= L0Φ0 −2(λ+ 1− f)(λ+ 1)(1 + (1−ff )2 bS20)− f∫∞0 Φ0ωρ dρ∫∞0 ω2ρ dρω2 = λΦ0 .To analyze this NLEP we use some previous results on NLEP’s similarto the one above, which can be found in [16]. From this previous theory,we know the stability threshold occurs when β(0) = 1, where β(λ) is themultiplier of the nonlocal term defined byβ(λ) =2(λ+ 1− f)(λ+ 1)(1 + (1−ff )2 bS20)− f.This condition determines a critical value Sc0 of S0 as S2c0 =(1−f)bf2 . Then,by (3.8), this determines the following leading order result for the stabilitythreshold:Dc =Dc0ν+Dc1 + · · · , Dc0 =f2(1− f)bE2|Ω|24pi2. (3.47)We remark that at this leading-order stability threshold we have Φ0 = ω,which then yields f2χ20 = (1− f)b and C0 = fχ20 =1−ff b = −B0.It is critical here to emphasize that the leading order stability thresholdis independent of the details of the lattice Λ, and does not depend on theBloch vector k. In order to determine the effect of the lattice on the stabilitythreshold we must proceed to one higher order.We now continue the calculation to one higher order. At the leadingorder stability threshold we can simplify the equation for Φ and N , and weexpand the eigenvalue as λ = νλ1 + · · · in order to calculate the spectrum303.2. Linear Stability Analysisnear the origin. At the leading order threshold, the equation for N1 is now∆N1 =1fω2 − ω , N ′1(0) = 0 , N1 → C0 ln |y|+B1 as |y| → ∞ .(3.48)Upon recalling the definition of V1P , it follows thatN1 =V1Pf+B1 . (3.49)Since we are seeking the second order term then, as similar to the construc-tion of the equilibrium spot solution, we will need to analyze the third orderequation and impose a solvability condition. At the leading order stabilitythreshold, the equation for the third-order term N2 becomes∆N2 + (1−2fω)Φ1 − 3ω2V1Pf3χ20−B1f2χ20ω2 − 2χ1fχ0ω2 = 0 ,N ′2(0) = 0 , N2 → C1 ln |y|+B2 , as |y| → ∞ .The solvability condition for this equation yields thatC1 = −∫ ∞0(1−2fω)Φ1ρ dρ+3f3χ20∫ ∞0ω2V1Pρ dρ+B11− f+2√b1− fχ1 .(3.50)Upon combining this equation with −2piR0,kC0 − C1 = B1, we derive thatB1 =(1− f2− f)(∫ ∞0(1−2fω)Φ1ρ dρ−3f3χ20∫ ∞0ω2V1Pρ dρ+2f3χ20∫ ∞0U1Pρ dρ− 2piR0,k1− ffb+ 2√b1− fS1) .(3.51)Once again, we want to eliminate the term∫∞0 Φ1ρ dρ. This can be doneby using the same method as done previously above. To do so, we substitute(3.51) into the equation of Φ1 to getL0Φ1 − λ1ω +3f2χ20ω2V1P +2χ1χ0ω2 +B1fχ20ω2 = 0 . (3.52)313.2. Linear Stability AnalysisWe then integrate both sides of this expression we obtain∫ ∞0Φ1ρ dρ =∫ ∞0Φ1ωρ dρ−2− f2− 2fλ1b+32f2χ20∫ ∞0ω2V1Pρ dρ+bχ1χ0− pibR0,k .(3.53)Then we substitute (3.53) and (3.51) back into (3.52) to concludeLΦ1 ≡ L0Φ1 −∫∞0 Φ1ωρ dρ∫∞0 ω2ρ dρω2 = λ1(ω +f2− 2fω2)−3f2χ20ω2V1P+(32f2χ20b∫ ∞0ω2V1Pρ dρ−χ1χ0+ piR0,k)ω2 .(3.54)Finally, λ1 is determined by imposing a solvability condition on Φ1 in (3.54).The adjoint operator of L is simplyL?Ψ ≡ L0Ψ− ω∫∞0 ω2Ψρ dρ∫∞0 ω2ρ dρ. (3.55)It is readily verified that if we define Ψ? = w+ρw′/2, then we have L?Ψ∗ = 0.The null space of L? was first identified in [8]. By imposing the Fredholmalternative on (3.54) we getλ1∫ ∞0(ω +f2− 2fω2)Ψ∗ρ dρ−32f2χ20∫ ∞0ω2V1PΨ∗ρ dρ (3.56)+ (32f2χ20b∫ ∞0ω2V1Pρ dρ−χ1χ0+ piR0,k)∫ ∞0ω2Ψ∗ρ dρ = 0 .This expression can be simplified considerably by using the following iden-tities:∫ ∞0ω2Ψ∗ρ dρ =∫ ∞0(L0ω)(L−10 ω)=∫ ∞0ω2ρ dρ = b ,∫ ∞0ωΨ∗ρ dρ =∫ ∞0ρω(ω +ρ2ω′) dρ =∫ ∞0ω2ρ dρ+14∫ ∞0[ω2]′ρ2 dρ =b2,∫ ∞0ω2V1PΨ∗ρ dρ =∫ ∞0(L0U1P )(L−10 ω)=∫ ∞0U1Pωρ dρ ,323.2. Linear Stability Analysis2∫ ∞0ωU1PΨ∗ρ dρ−∫ ∞0U1Pρ dρ =∫ ∞0ω2V1Pρ dρ .In this way, we determine the spectrum near the origin in the spectralplane asλ1 = 2(1− f)(χ1χ0− piR0,k +32f2χ20b∫ ∞0U1Pρ dρ). (3.57)Since we are seeking a two-term expansion for the stability threshold Dc ofD, we need to write χ0 and χ1 in term of D. In Principal Result 1, χ waswritten in term of S, and S is in term of D asS = ν12 (S0 + νS1 + · · · ) =E|Ω|2pi√D=E|Ω|2pi[D0(1 + νD1D0+ · · ·)]− 12= ν−12E|Ω|D− 1202pi(1−ν2D1D0+ · · ·).Now with Dc0 =f2(1−f)bE2|Ω|24pi2 , we have Sc1 = −12S0Dc1Dc0, so thatχ1χ0=fχ0b(−1f3χ30∫ ∞0U1Pρ dρ− S1f1− f).Therefore, we haveλ1 = 2(1− f)(12(1− f)b2∫ ∞0U1Pρ dρ+2pi2(1− f)bf2E2|Ω|2D1 − piR0,k).(3.58)Notice that this is a continuous band of spectra parametrized by theBloch vector k. This is illustrated schematically in Figure 3.2. As proved inLemma 1 and Lemma 2, R0,k is real and tends to infinity as |k| → 0. Thisshows from (3.58) that λ1 is real, and leaves the ball of radius O(ν)  1near the origin along the negative real axis as |k| → 0.Therefore, in order to have stability, we need λ1 < 0, ∀k ∈ Ω∗. Wesummarize our result as follows:Principal Result 2. In the limit  → 0, D ∼ O( 1ν ), we have constructed333.2. Linear Stability AnalysisFigure 3.1: The Continuous Band of Spectraperiodic spot solutions(with respect to a fixed Bravais lattice Λ) in PrincipalResult 1. The linearized operator around this solution has a continuousspectrum within an O(ν) neighbourhood of the origin and is parametrized bya Bloch vector k ∈ Ω∗\0:λ(k) = 2ν(1− f)(12(1− f)b2∫ ∞0U1Pρ dρ+2pi2(1− f)bf2E2|Ω|2D1 − piR0,k).(3.59)To have linear stability, we need λ(k) < 0, ∀k ∈ Ω∗\0, which gives a twoterm asymptotic expansion for the stability threshold Dc:Dc =Dc0ν+Dc1 + · · · , where Dc0 =f2E2|Ω|24pi2(1− f)b, andDc1 = mink∈Ω∗\0{f2E2|Ω|22pi2(1− f)b(piR0,k −12(1− f)b2∫ ∞0U1Pρ dρ)} ,343.3. A Quick Derivation of the Stability Threshold= mink∈Ω∗\0{Dc0(2piR0,k −1(1− f)b2∫ ∞0U1Pρ dρ)}= mink∈Ω∗\0{Dc0(2piR0,k +1(1− f)b2∫ ∞0U1QIρ dρ−1b2∫ ∞0U1QIIρ dρ)} .We remark here that since R0,k is real-valued, then so is the thresholdDc,which is what we should expect. As shown above, the second order term inDc depends on the lattice Λ. In order to compare the stability threshold ondifferent lattices, we will fix |Ω| = 1. This leads to the following optimalitycriterion:Principal Result 3. With fixed primitive cell of area unity, the optimallattice arrangement Λop is the one with the largest stability threshold:Λop = arg maxΛ, |Ω|=1{Dc(Λ)} = arg maxΛ, |Ω|=1{ mink∈Ω∗\0{R0,k}} . (3.60)Some numerical results to identify the optimal lattice arrangement isgiven in Section A Quick Derivation of the StabilityThresholdIn this section, we give another much more expedient way to derive thestability threshold. Recall that in the inner core problem (3.3), S is a pa-rameter in the asymptotic boundary condition. More specifically, if we viewS as a parameter of the solution, then∆ρU(|y|, S)− U(|y|, S) + fU2(|y|, S)V (|y|, S) = 0 ,∂U∂ρ(0, S) = 0 , U → 0 as |y| → ∞ ,∆ρV (|y|, S) + U(|y|, S)− U2(|y|, S)V (|y|, S) = 0 ,∂V∂ρ(0, S) = 0 , V → S ln |y|+ χ(S, f) as |y| → ∞ .353.3. A Quick Derivation of the Stability ThresholdUpon taking the partial derivative of U and V with respect to S, we get∆(∂U∂S∂V∂S)+(−1 + 2fUV fU21− 2UV −U2)(∂U∂S∂V∂S)=(00), (3.61)with the boundary conditions∂2U∂S∂ρ(0, S) = 0 ,∂U∂S→ 0 , as |y| → ∞ ,∂2V∂S∂ρ(0, S) = 0 ,∂V∂S→ ln |y|+∂χ∂S(S, f) , as |y| → ∞ .Then, we observe that ∂U∂S and∂V∂S are, up to a scalar multiple, thesolution to the perturbed core problem (3.33) and (3.34) when λ = 0. Tofix the boundary conditions in (3.33) and (3.34) for N(y), we must chooseS appropriately. This constraint on S to hold when λ = 0 is the stabilitythreshold we are seeking. Since the solution to (3.33) and (3.34) is unique upto a constant scaling, then upon comparing the boundary conditions with(3.33) and (3.34), it follows that the stability threshold Sc occurs whenBC=∂χ∂S(Sc, f) . (3.62)Then, together with (3.40), we have∂∂Sχ(Sc, f) = −1ν− 2piR0,k +O(ν) . (3.63)We then use the expansion for χ from Principal Result 1, which we write asχ(S, f) = ν−12(b(1− f)f21S0+ ν(−(1− f)bf2S1S20−S0(1− f)b2∫ ∞0U1Pρ dρ) + · · ·),S = ν12 (S0 + νS1 + · · · ) ,to deriveχ(S, f) =b(1− f)f21S−∫∞0 U1Pρ dρ(1− f)b2S +O(ν) , (3.64)363.3. A Quick Derivation of the Stability Threshold=b(1− f)f21S+(1b2(1− f)∫ ∞0U1QIρ dρ−1b2∫ ∞0U1QIIρ dρ)S +O(ν) .Upon taking the partial derivative we obtain∂∂Sχ(S, f) = −b(1− f)f21S2−∫∞0 U1Pρ dρ(1− f)b2+O(ν) , (3.65)= −1νb(1− f)f21S20(1− 2S1S0ν + · · · )−∫∞0 U1Pρ dρ(1− f)b2+O(ν) ,= −1νb(1− f)f21S20+ (2b(1− f)f2S1S30−∫∞0 U1Pρ dρ(1− f)b2) +O(ν) .Upon comparing this expression with (3.63), we obtain from equating powersof ν that−1 = −b(1− f)f21S2c0,−2piR0,k = 2b(1− f)f2Sc1S3c0−∫∞0 U1Pρ dρ(1− f)b2,= 2b(1− f)f2Sc1S3c0+1b2(1− f)∫ ∞0U1QIρ dρ−1b2∫ ∞0U1QIIρ dρ .In this way, we can solve for Sc0 and Sc1, then obtain the corrections Dc0and Dc1 to the stability threshold from expanding the relation (3.8) betweenS and D. This yields the same result for the stability threshold as derivedin Principal Result 2. We remark that although this simple derivation isable to quickly isolate the stability threshold, it is unable to give any preciseaccount of the nature of the spectrum of the linearized operator near theorigin when D is near the leading-order stability threshold.37Chapter 4Spot Patterns for theBrusselator on a FiniteDomainIn this chapter, we first construct N -spot solutions to the Brusselator model,formulated asut = 2∆u− u+ f(u2v + 22Euv + 4Ev) ,τvt = D∆v + E + −2(u− u2v)− 2Euv − 2E2v ,(4.1)on a finite domain x ∈ Ω ⊂ R2, with no-flux boundary conditions∂nu(x) = ∂nv(x) = 0 , x ∈ ∂Ω , (4.2)where n is the outer normal vector on ∂Ω. After constructing such multi-spot patterns, we then perform a linear stability analysis to calculate astability threshold corresponding to a zero-eigenvalue crossing.The asymptotic construction of the N -spot pattern and the linear sta-bility analysis is very similar to that for the periodic case. One of the keydifferences between the periodic and finite-domain problems, is that for theperiodic case we need only construct a one-spot solution in one primitivecell. In contrast, for the finite domain case, spots interact with each otherin the domain through the Green’s function. As a result, a Neumann Greenmatrix together with its eigenvalues and eigenvectors play a key role in theanalysis.384.1. The N -Spot Solutions4.1 The N-Spot SolutionsIn this section, our goal is to construct N -spot solutions where u(x) concen-trates around N distinct O() balls centred at N given points x1,x2, ...,xn ∈Ω. The position of these N points are not arbitrary, but satisfy some con-ditions to be derived and explained below.We first introduce local coordinates around each of these N points inthe formy = −1(x− xj) , Uj(y) =u(x)√D, Vj(y) =√Dv(x) , j = 1, 2, ..., N .(4.3)Then, in an inner region around xj , we look for a radially symmetric solutionin the form Uj(y) = Uj(ρ), Vj(y) = Vj(ρ), where ρ = |y|. To leading order,we get the core problem∆ρUj − Uj + fU2j Vj = 0 , ∆ρVj + Uj − U2j Vj = 0 , 0 < ρ <∞ , (4.4)with the boundary conditionsU ′j(0) = V′j (0) = 0 , Uj → 0 , Vj → Sj ln ρ+ χ(Sj , f) , as ρ→∞ .(4.5)Here the introduction of the source strength Sj and the function χ(Sj , f)is the same as in (3.4). Upon integrating ∆Vj and using the divergencetheorem we get2piSj =∫R2(U2j Vj − Uj) dx =∫ 2pi0dθ∫ ∞0(U2j Vj − Uj)ρ dρ ,which yieldsSj =∫ ∞0(U2j Vj − Uj)ρ dρ . (4.6)In the outer region, given that u is localized around {xi}Ni=1, then byusing the relation (4.6) we obtain that the leading order outer solution v394.1. The N -Spot SolutionssatisfiesD∆v + E = 2pi√DN∑i=1Siδ(x− xi) , x ∈ Ω , (4.7)with ∂nv = 0 for x ∈ ∂Ω. Upon integrating (4.7) over Ω, and by using theno-flux boundary conditions (4.2), we get the solvability conditionE =2pi√D|Ω|N∑j=1Sj . (4.8)To represent the solution v, it is convenient to introduce the Neumann Greenfunction G(x, ξ), which satisfies∆xG0(x, ξ) =1|Ω|− δ(x− ξ) , (4.9)∇xG0(x, ξ) · n = 0 , ∀x ∈ ∂Ω ,∫ΩG0(x, ξ) dx = 0 , (4.10)G0(x, ξ)→ −12piln |x− ξ|+R0(ξ) , as x→ ξ , (4.11)where n is the outer normal vector to ∂Ω. We remark here that R(ξ) is theregular part of the Green function while − 12pi ln |x− ξ| is the singular part.The right hand side of the equation 1|Ω| − δ(x − ξ) is consistent with theno-flux boundary conditions ∇xG0(x, ξ) ·n = 0, ∀x ∈ ∂Ω. We also imposethe uniqueness condition,∫ΩG0(x, ξ) dx = 0, since the solution to the PDEwith the no-flux boundary conditions is only unique up to a constant.In terms of G0, the solution v to (4.7) can be represented asv(x) = −N∑i=12piSi√DG0(x,xi) + c , (4.12)where c is some constant. Then, by letting x→ xj , and by matching to theinner solution near each spot, we obtain the following nonlinear system of404.1. The N -Spot SolutionsN equations for Sj , j = 1, . . . , N , and for c:Sjν+ χ(Sj) = −2piR0(xj)Sj − 2pi∑i 6=jG0(xi,xj)Si +√Dc , 1 ≤ j ≤ N .(4.13)The remaining equation to complete the system is (4.8). It is convenient torewrite this system in matrix form by introducing the following notation:S ≡S1...SN , e ≡1...1 , χ(S, f) ≡χ(S1, f)...χ(SN , f) ,G ≡R0(x1) G0(x1,x2) · · · G0(x1,xN )G0(x2,x1) R0(x2) · · ·.... . ....G0(xN ,x1) · · · R0(xN ).We shall refer to G as the Neumann Green matrix of x1,x2, ...,xN . Then,the N matching conditions (4.13) and the solvability condition (4.8) can bewritten in matrix form asS + νχ(S, f) = −ν2piGS + ν√Dce ,eTS =|Ω|E2pi√D.(4.14)For simplicity we will assume that the N spots have a common sourcestrength S, i.e. for Sj = S, j = 1, 2, ..., N . For such a pattern, we have thatS = Se, χ(S, f) = χ(S, f)e, and that (4.14) reduces toGS = SGe = −12piν(S + νχ(S, f)− ν√Dc)e , S =E|Ω|2piN√D. (4.15)Therefore, it follows that for such a common source strength pattern toexist we must have that e is an eigenvector of G(x1,x2, ...,xN ). We remarkthat the existence of such a special eigenvalue does not generally occur fora pattern of N arbitrarily-located spots. Since the Green function satisfies414.1. The N -Spot Solutionsthe usual reciprocity condition, it follows that G is a symmetric matrix andcan be diagonalized with an orthonormal basis as 1√Ne, {qj}Nj=2, i.e.G(1√Ne) = k1(1√Ne) , Gqj = kjqj , j = 2, 3, ..., N, (4.16)qTj1√Ne = 0 , j = 2, 3, ..., N , qTj qi = 0 , ∀ 2 ≤ i, j ≤ N .Upon using this relation, we obtain for a common source strength patternthat (4.14) reduces to12piν(S + νχ(S, f)− ν√Dc) = −Sk1 . (4.17)Next, since the stability threshold occurs in the regime where D ∼O(ν−1), it follows that S ∼ O(ν12 ) at a zero eigenvalue crossing. By fol-lowing the same procedure as for the periodic problem, we can calculate thesmall S asymptotics of the solution to the core problem (4.4) as follows:Principal Result 4. For S ∼ ν12 (S0 + νS1 + · · · ), where ν ≡ − 1ln  , theradially symmetric asymptotic solutions to the core problem (4.4) in an O()ball centred at xj is given by:Uj ∼ ν12 (Uj0 + νUj1 + · · · ) , Vj ∼ ν− 12 (Vj0 + νVj1 + · · · ) ,χ ∼ ν−12 (χ0 + νχ1 + · · · ) , (4.18)where Uj0(ρ), Uj1(ρ), Vj0(ρ) and Vj1(ρ) are defined byUj0 =ωfχ0, Uj1 = −χ1fχ20ω −1f3χ30((1− f)U1QII − U1QI) , (4.19)Vj0 = χ0 , Vj1 = χ1 +1f2χ0((1− f)V1Q − ω) . (4.20)Here U1QI , U1QII , V1Q are the unique radially symmetric solutions to (3.26),(3.26), (3.24) as before, ω is the unique radially symmetric solution to424.2. Linear Stability Analysis(3.17), while χ0 and χ1 are defined asχ0 =(1− f)f2bS0, (4.21)χ1 = −(1− f)bf2S1S20−S0b2∫ ∞0U1QIIρ dρ+S0(1− f)b2∫ ∞0U1QIρ dρ . (4.22)We remark here that since the source strength is the same for each j,then so are the inner solutions. Moreover, we remark that such patterns caneither be true steady-state solutions if the spot locations satisfy some furthercondition, or simply quasi-steady patterns that can persist for very long timeintervals provided that there is no unstable O(1) eigenvalue in the spectrumof the linearization. The analysis of the spectrum of the linearization isanalyzed in the next section. This completes our construction of an N -spotsolution where the spots have a common source strength.4.2 Linear Stability AnalysisWe denote the N -spot solution constructed above as ue(x), ve(x) and weintroduce the perturbationu(x) = ue(x) + eλtφ , v(x) = ve(x) + eλtη . (4.23)Upon substituting (4.23) into (4.1), we linearize around the N -spot solutionto obtain the singularly perturbed problem (3.31), written again asλ(φτη)=(2∆φD∆η)+(−1 + 2fueve fu2e−2 − 2−2ueve −2u2e)(φη).In the inner region around each xj , we introduce the local variables asy = −1(x− xj) , Φj(y) =φ(y + xj)D, Nj(y) = η(y + xj) . (4.24)We look for radially symmetric solution of the form Φj(y) = Φj(ρ), Nj(y) =434.2. Linear Stability AnalysisNj(ρ), where ρ = |y|. Then, the inner problem near the j-th spot at xjreduces asymptotically to∆ρΦj − (λ+ 1)Φj + 2fUjVjΦ + fU2jNj = 0 , (4.25)∆ρNj + (1− 2UjVj)Φj − U2jNj = 0 ,Φ′j(0) = N′j(0) = 0 ; Φj → 0 , Nj → Cj ln |y|+Bj(Cj , f) , as ρ→∞ .We remark here that for a fixed f , the ratio BjCj is a constant since the systemis linear. The solvability condition for the Nj equation yieldsCj =∫ ∞0(U2jNj − (1− 2UjVj)Φj)ρ dρ . (4.26)Since both ue(x) and φ(x) are localized near {xj}Nj=1, the outer approx-imation for the eigenfunction component η(x) satisfies the leading orderproblem∆η −τλDη = 2piN∑i=1Ciδ(x− xi) . (4.27)Notice that when λ = 0, corresponding to the stability threshold, then ifwe integrate (4.27) over the domain and use the no-flux boundary conditionwe conclude that∑Nj=1Cj = 0. However, we do not have this constraintwhen λ 6= 0. This observation suggest that we need to split our analysis intoseveral cases.4.2.1 λ 6= 0 and λ ∼ O(1)First we introduce the reduced-wave Green function, which satisfies∆xGσ(x, ξ)− σ2Gσ(x, ξ) = −δ(x− ξ) , (4.28)∇xGσ(x, ξ) · n = 0 , ∀x ∈ ∂Ω , (4.29)Gσ(x, ξ)→ −12piln |x− ξ|+Rσ(ξ) , as x→ ξ . (4.30)We need an approximation to this Green’s function when σ  1. Sincethere is no solution when σ = 0, then as we did in (2.3) we must seek Gσ444.2. Linear Stability Analysisfor σ  1 in the form of a singular asymptotic expansion asGσ(x, ξ) =1σ2|Ω|+G0(x, ξ) +O(σ2) , (4.31)where G0(x, ξ) is the Neumann Green function defined in (4.9). Then if wedenote σ = τλD , we can represent the solution to (4.27) in the formη(x) = −2piN∑i=1CiGσ(x,xi) = 2piN∑i=1Ci(−Dτλ|Ω|−G0(x,xi) +O(σ2)).(4.32)By expanding this solution as x → xj and comparing the resulting ex-pression with the far-field behavior of the corresponding inner solution, weobtain the following matching condition near each spot:Cjν+Bj = −2piDτλ|Ω|N∑i=1Ci−2piR(xj)Cj−∑i 6=j2piG0(xi,xj)Ci , 1 ≤ j ≤ N .(4.33)We then rewrite this system in matrix form by first introducing thenotationC ≡C1...CN , B ≡B1...BN , Φ ≡Φ1...ΦN , N ≡N1...NN , (4.34)E ≡1NeeT =1N1 1 · · · 1.... . ....1 · · · 1 . (4.35)In terms of these new variables, (4.33) can be written in matrix form as(I +2piνDNτλ|Ω|E + 2piνG)C = −νB . (4.36)From this system it follows that B is one order higher in ν than Cwhen λ ∼ O(1). Moreover, since the system (4.4) is linear, we may assume454.2. Linear Stability AnalysisN(y) ∼ O(1). This suggests that we introduce the asymptotic expansion asΦj ∼ ν(Φj0 + νΦj1 + · · · ) , Nj ∼ Nj0 + νNj1 + · · · , (4.37)Bj ∼ Bj0 + νBj1 + · · · , Cj ∼ ν(Cj0 + νCj1 + · · · ) .Upon substituting this expansion into (4.25), and by using the results fromPrincipal Result 4, we obtain at leading order that∆ρΦj0 − (λ+ 1)Φj0 + 2fUj0Vj0Φj0 + fU2j0Nj0 = 0 , ∆ρNj0 = 0 , (4.38)Φ′j0(0) = N′j0(0) = 0 ; Φj0 → 0, Nj0 → Bj0 as ρ→∞ . (4.39)The solution of this system in terms of the linear operator L0 of Chapter 3,defined by L0φ = ∆φ− φ+ 2ωφ, is simplyL0Φj0 = λΦj0 −Bj0fχ20ω2 , Nj0 = Bj0 , (4.40)⇒ L0Φ0 = λΦ0 −ω2fχ20B0 , N0 = B0 .Moreover, to leading order from the matching condition (4.36), we concludethat(I + µE)C0 = −B0 , µ ≡2piD0Nτλ|Ω|. (4.41)At the next order, we obtain from (4.25) that the equation for Nj1 is∆ρNj1 + (1− 2Uj0Vj0)Φj0 − U2j0Nj0 = 0 ,N ′j1(0) = 0 ; Nj1 → Cj0 ln ρ+Bj1 , as ρ→∞ .Upon integrating this equation, and using the divergence theorem, we obtainthe consistency condition thatCj0 =bf2χ20Bj0 −∫ ∞0(1− 2ωf)Φj0ρ dρ , (4.42)⇒ C0 =bf2χ20B0 −∫ ∞0(1− 2ωf)Φ0ρ dρ .464.2. Linear Stability AnalysisWe want to eliminate the∫∞0 Φ0ρ dρ term as before by integrating the equa-tion for Φ0 in (4.40). This yields that(λ+ 1)∫ ∞0Φ0ρdρ = 2∫ ∞0ωΦ0ρ dρ+bfχ20B0 . (4.43)Then, upon combining (4.41), (4.42) and (4.43), we conclude that[(1 + a)I + aµE ]B0 = −2(λ+ 1− f)f(λ+ 1)(I + µE)∫ ∞0ωΦ0ρ dρ , (4.44)where a ≡ b(λ+1−f)f2χ20(λ+1). Using the fact that (I + kE)−1 = I − kk+1E , we getB0 = −2m(1 +bmfχ20)−1I +µτλ11 + bmfχ20(1 + µτλ)E∫ ∞0ωΦ0ρ dρ ,(4.45)where m is defined as m ≡ 1f −1λ+1 . Upon substituting this expression backinto L0Φ0 = λΦ0 − ω2fχ20B0, we get a vector nonlocal eigenvalue problem(NLEP).Next, we will decompose Φ0 into certain directions, most of which areunaffected by the matrix E . We observe that the geometric meaning of thematrix E is that of projecting a vector into the direction e = (1, 1, ..., 1)T .This suggests that we decompose Φ0 into ke+ r, where k is some constantand r⊥e. Recall from (4.16) that the eigenvectors of the Neumann Greenmatrix G are such that 1√Ne, {qj}Nj=2 forms a orthonormal basis of RN andqj⊥e for j = 2, 3, ..., N . We decompose Φ0 into this basis by writingΦ0(|x|) = k(|x|)e+N∑j=2uj(|x|)qj , (4.46)where k(|x|), uj(|x|) are radially symmetric coefficient functions. Uponsubstituting this into the equation for Φ0, and by using the formula for B0,474.2. Linear Stability Analysiswe obtain the two distinct NLEP’sL0k(|x|)−2b(λ+ 1− f)(τλ+ µ)(λ+ 1)f2χ20τλ+ b(1− f + λ)(τλ+ µ)∫∞0 kωρ dρ∫∞0 ω2ρ dρω2 = λk(|x|) ,L0uj −2b(λ+ 1− f)(λ+ 1)f2χ20 + b(1− f + λ)∫∞0 ujωρ dρ∫∞0 ω2ρ dρω2 = λuj . (4.47)These two distinct NLEP’s are similar to those considered previously. Assuch we conclude that the component k(|x|) is linearly stable, but that thereis a stability threshold for uj(|x|) whenβ(λ)|λ=0 =2b(λ+ 1− f)(λ+ 1)f2χ20 + b(1− f + λ)|λ=0 = 1 ,which yields thatχ20 =(1− f)bf2, Sc0 =√b(1− f)f. (4.48)We remark here that at the stability threshold we have λ = 0, whichseems to contradict our starting assumption λ 6= 0. However, the previ-ous theorem of [16], states that when β(0) < 1 the NLEP has a positivereal eigenvalue. Another difficulty is that near the stability threshold theeigenvalue can be very small. In particular, if λ = O(ν), then the 2piνDNτλ|Ω| Eterm in the matching condition (4.36) becomes the dominate term and theasymptotic expansions need a little modification. We will handle these twosituations in Section 4.2.2 and Section 4.2.3, respectively.4.2.2 λ ∼ O(ν) and λ 6= 0Next we treat the case where λ ∼ O(ν) and λ 6= 0, and we expand λ asλ = νλ1 + ν2λ2 + .... For this case the results (4.27), (4.28), and (4.31)still hold. However, the only difference is that the Dτλ|Ω| term is no longerO(ν−1), but instead is O(ν−2). This implies that the matching condition ismodified as (2piD0Nτλ1|Ω|E + νI + 2piν2G)C = −ν2B . (4.49)484.2. Linear Stability AnalysisSince the system is linear, then without loss of generality we may assume thatB ∼ O(ν) and from (4.49) it seems that C ∼ O(ν2). However, this scalingassumption would be inconsistent with the logarithmic growth condition inthe equation for Nj1. In fact, the leading-order solution to (4.49 is thatEC = 0, which is equivalent to C⊥e. We then expand the solutions as in(4.37) and obtain that the matching condition becomes(I + 2piνG)C = −νB . (4.50)Due to the properties of G in (4.16), it follows that B⊥e. Therefore, atleading order, we haveB0 = −C0⊥e, L0Φ0 = −ω2fχ20B0, N0 = B0, ⇒ Φ0 = −ωfχ20B0 .(4.51)We may assume that C0 =∑Nj=2 djqj = −B0, where qj are other eigenvec-tors of G that are orthogonal to e and dj are some constant coefficients. Atnext order, the equation for N1 is:∆N1 =ωfχ20B0 −ω2f2χ20B0 , (4.52)N ′1(0) = 0 , N1 → C0 ln |y|+B1 , as |y| → ∞ .The solvability condition for this equation gives:C0 =bfχ20(1−1f)B0 = −B0, ⇒ χ0 =√b(1− f)f2, S0 =√b(1− f)f,(4.53)which is precisely the same threshold we obtained from (4.47). The solutionto (4.52) can then be written asN1 = B1 −V1Pb(1− f)B0 = B1 +V1Pb(1− f)N∑j=2djqj . (4.54)494.2. Linear Stability AnalysisUpon substituting this expression into the equation for N2 we get∆N2 + (1−2ωf)Φ1 − 2(ωfχ0(χ1 +V1Pf2χ0) + χ0(−χ1fχ20ω −U1Pf3χ30))Φ0−ω2f2χ20N1 − 2ωfχ20(−χ1fχ20ω −U1Pf3χ30)N0 = 0 ,N ′2(0) = 0 , N2 → C1 ln |y|+B2 , as |y| → ∞ .Upon using the divergence theorem on this equation for N2, and by using(4.51) and (4.54), we obtain that the following consistency condition musthold:C1 +∫ ∞0Φ1ρdρ−2ωf∫ ∞0Φ1ρ dρ−11− fB1 −3b2(1− f)2∫ ∞0ω2V1Pρ dρN∑j=2djqj − 2fχ1b12 (1− f)32N∑j=2djqj = 0 . (4.55)Similarly, the matching condition (4.50) givesC1 + 2piGC0 = −B1, ⇒ C1 + 2piN∑j=2kjdjqj = −B1. (4.56)The equation for Φ1 then becomesL0Φ1 +3fω2V1pb2(1− f)2N∑j=2djqj +2f2χ1ω2(b(1− f))32N∑j=2djqj (4.57)+fω2b(1− f)B1 = λ1fb(1− f)ωN∑j=2djqj .Upon integrating this equation and then substituting into (4.55) and (4.56),we get thatB1 = piN∑j=2djkjqj −1− ff∫ ∞0ωΦ1ρ dρ504.2. Linear Stability Analysis−fb(1− f)N∑j=2djqj(bχ1(1− f)fχ0+32bf∫ ∞0ω2V1Pρ dρ+b2λ1).Finally, we substitute this expression back into (4.57) to obtain the vectorNLEPLΦ1 ≡ L0Φ1 −∫∞0 Φ1ωρ dρ∫∞0 ω2ρ dρω2 = λ1(ω +f2− 2fω2)fb(1− f)N∑j=2djqj−3fω2V1Pb2(1− f)2N∑j=2djqj −fb(1− f)N∑j=2djqj(−pikj +χ1χ0−3∫∞0 ω2V1Pρ dρ2b2(1− f))ω2 .Upon decomposing Φ1(x) = r(x)e+∑Nj=2 tj(x)qj as before, we obtain thatthe coefficient functions tj(x) satisfyLtj(x) =fdjb(1− f){λ1(ω +f2− 2fω2)−3b(1− f)ω2V1P− ω2(−pikj +χ1χ0−32b2(1− f)∫ ∞0ω2V1Pρ dρ)}.Finally, we use a solvability condition on this problem to calculate λ1. FromSection 3.2 the adjoint operator L∗ has a one-dimensional nullspace Ψ∗(x)in the class of radially symmetric functions, where Ψ? = w + ρw′/2 wasgiven in Section 3.2. We then impose the solvability condition as similar tothat done in Section 3.2 to conclude that∫ ∞0RHS ·Ψ∗ρ dρ = 0 ⇒ λ1 = 2(1− f)(−pikj +χ1χ0+3∫∞0 U1Pρ dρ2(1− f)b2)= 2(1− f)(−pikj −S1S0+∫∞0 U1Pρ dρ2(1− f)b2).With this expression we can calculate the next order term in the stabilitythreshold that makes λ1 = 0 asSc1 = max2≤j≤N{S02(−2pikj +1(1− f)b2∫ ∞0U1Pρ dρ)} = max2≤j≤N{√b(1− f)2f514.2. Linear Stability Analysis(−2pikj +1(1− f)b2((1− f)∫ ∞0U1QIIρ dρ−∫ ∞0U1QIρ dρ))} . (4.58)We remark here that since G is real symmetric, all of its eigenvalues kj arereal. Therefore, as expected, the stability threshold is real-valued.4.2.3 λ=0Finally we consider the case where λ = 0, which corresponding to the sta-bility threshold. The equation for η now becomes∆η = 2piN∑j=1Cjδ(x− xj) . (4.59)As we mentioned before, if we integrate this PDE with the no-flux boundaryconditions on ∂Ω, we obtain thatN∑j=1Cj = 0 , η(x) = −2piN∑j=1CjG0(x,xj) , (4.60)where G0(x,xj) is the Neumann Green function defined as before. Usingthe same method in Section 3.3, we derive that ∂∂SUj(y, Sj),∂∂SVj(y, Sj) aresolutions to (4.25) when λ = 0. This leads to the relationBjCj=∂∂Sχ(Sj , f) . (4.61)Since we are assuming a common source strength Sj = S = Sc, (4.61) givesB = ∂∂Sχ(Sc, f)C. We substitute this expression back into the matchingcondition (4.36) and, upon noticing that EC = 0, we have(I + ν∂χ∂S(Sc, f))C = −2piνGC . (4.62)This implies that C is an eigenvector of G, and we obtain the relation1 + ν∂χ∂S(Sc, f) = −2piνkj . (4.63)524.2. Linear Stability AnalysisAs we derived in (3.64) and (3.65), we calculate∂∂Sχ(Sc, f) = −1νb(1− f)f21S2c0+ (2b(1− f)f2Sc1S3c0−∫∞0 U1Pρ dρ(1− f)b2) +O(ν) .(4.64)We substitute this expression back into (3.63) and equate terms of a commonorder in ν. In this way, we derive the stability threshold resultsSc0 =√b(1− f)f, (4.65)Sc1 = max2≤j≤N{√b(1− f)2f(−2pikj +1(1− f)b2((1− f)∫ ∞0U1QIIρ dρ−∫ ∞0U1QIρ dρ))} .Finally, by using (4.15) which relates S to D, we calculate the stabilitythreshold in terms of D. The results are summarized as follows:Principal Result 5. In the limit  → 0, and on the range D ∼ O( 1ν ), themulti-spot patterns constructed in Principal Result 4 are linearly stable ifD < Dc =Dc0ν+Dc1 + ... , (4.66)where Dc0 and Dc1 are defined byDc0 =f2E2|Ω|24pi2N2(1− f)b, (4.67)Dc1 = min2≤j≤N{Dc0(2pikj −1(1− f)b2∫ ∞0U1Pρ dρ)} (4.68)= min2≤j≤N{Dc0(2pikj +1(1− f)b2∫ ∞0U1QIρ dρ−1b2∫ ∞0U1QIIρ dρ)} .(4.69)We remark here that since we are solving for eigenvalues of a self-adjointoperator, all the eigenvalues are real-valued, and so the stability threshold isreal. Although we have discussed the three cases of λ separately, the analysisindeed provides a uniform transition between the ranges of λ. More specif-534.2. Linear Stability Analysisically, we obtain the same leading order results for the stability thresholdfrom (4.48), (4.53) and (4.65). Moreover, (4.58) also agrees with (4.65) atsecond order. Finally, in Section 4.2.2, we must have∑Nj=1Cj = 0 to havea solution, which also agrees with the solvability condition in Section 5Numerical ResultsIn this chapter, we perform some numerical experiments and compare theresults with the two term asymptotic approximations for the stability thresh-old derived in the previous chapter. For the periodic case, we will identifythe optimal lattice arrangement of spots. For the finite domain problem weillustrate our theory for the case of equally-spaced spots on a ring that isconcentric with the unit disk. For this finite domain problem, there is anexplicit formula for the Neumann Green function that will be used.5.1 Small S Asymptotics of χ(S, f)For both the periodic and the finite domain problems, the same core prob-lems (3.3) and (4.4) arise in the asymptotic construction of the spot pattern.In these common inner problems, a key quantity is χ(S, f), which appearsin the asymptotic boundary condition. Two-term asymptotic expansions forχ(S, f) have been derived previously in Principal Result 1 and Principal Re-sult 4. We now solve the core problem numerically to compute the χ(S, f),and we compare the numerical results for χ(S, f) with the correspondingtwo-term asymptotic results derived in the small S limit.Since we are seeking radially symmetric solutions, solving the core prob-lem is actually solving an ODE system. We use the ODE boundary valueproblem solver BVP4C in Matlab. We now remark on a few details of thenumerical implementation.• The Laplace operator in R2 polar coordinates is expressed as:∆ρ =d2dρ2+1ρddρ.555.1. Small S Asymptotics of χ(S, f)• Instead of solving the ODE systems on the whole interval [0,∞), weuse the interval [µ,R], where R is a sufficiently large number so as toapproximate the infinite domain, while µ is a sufficiently small numberto avoid the singularity at r = 0. We chose R = 15 and µ = 0.005 inour computations.• Instead of using the boundary conditions for U(ρ) and V (ρ) as ρ→∞in (3.4) directly, we set U(R) = 0 and V′(R) ∼ SR . Then, after solvingthe core problem numerically, we define χ(S, f) by χ(S, f) = V (R)−S lnR.• The boundary value solver requires a good initial guess consistent withthe boundary conditions. First we may want to use our small S asymp-totic approximation in Principal Result 1 as an initial guess. To dothis, we in principle need to know the radially symmetric solution of∆ω−ω+ω2 = 0 in R2. However, this explicit solution is not availablein R2, and is only available in R1. As such, we adapt a homotopyalgorithm to find the initial guess. For a fixed f , we first start with asmall S0, and use the asymptotic approximation for the core problemin R. We then slowly increase the dimension from 1 to 2 and solve thecore problem using the previous step as an initial guess. After hav-ing obtained the core solution for S = S0 in R2 using this homotopystrategy, we then increase S and solve the core problem based on thepreviously computed solution.After computing the χ(S, f) in this way, we compare the results withthe two term asymptotic expansions in (3.64). In the asymptotic approx-imation, we require numerical values for a few integrals. We obtain thatb ≈∫∞0 ωρdρ = 4.9343,∫∞0 U1QIρdρ ≈ 11.9131 and∫∞0 U1QIρ dρ ≈ 11.4384.Figure 5.1 shows results for f = 0.4 and f = 0.5, where the green (top)curves are the asymptotic approximations and the blue (bottom) curves arethe full numerical results. We observe that the two curves are rather closefor small S, which is what we should expect.565.2. Stability Threshold and the Optimal Lattice Arrangement0 1 2 3 4 5 6020406080100120140160180200χsf=0.40 1 2 3 4 50102030405060708090100χsf=0.5Figure 5.1: Numerical solution (bottom curves) and asymptotic results (topcurves) for χ(S, f). In the left panel we fix f = 0.4, while f = 0.5 forthe right panel. In both pictures, the blue (bottom) curve is the numericalsolution while the green (top) one is the two term asymptotic expansion.5.2 Stability Threshold and the Optimal LatticeArrangementIn this section, we compute the stability threshold numerically and comparethe results with the two term asymptotic approximation.As derived in (3.63) and (3.65), the stability threshold for S, labeledby Sc, for the periodic spot problem is the largest S, corresponding to thesmallest D, that solves the transcendental equation∂∂Sχ(Sc, f) = −1ν− 2piR0,k +O(ν) , (5.1)for some Bloch vector k in the first Brillouin zone. For the correspondingfinite domain problem, the stability threshold is the largest S that solves∂∂Sχ(Sc, f) = −1ν− 2piki , (5.2)for certain eigenvalues ki of the Neumann Green matrix.Since we have already computed χ(S, f) above, we can use a cubic spline575.2. Stability Threshold and the Optimal Lattice Arrangementand a numerical derivative to get ∂∂Sχ(S, f). Then we use a nonlinear equa-tion solver to compute the threshold directly.On the other hand, we have derived the two term asymptotic approxi-mation for the stability threshold in both cases. For the periodic case, wederived previously thatSc = ν12√b(1− f)f(1 + ν(12b2∫ ∞0U1QIIρ dρ−12b2(1− f)∫ ∞0U1QIρ dρ− pi mink∈Ω∗\0R0,k)) , (5.3)while for the finite domain problem we derived thatSc = ν12√b(1− f)f(1 + ν(12b2∫ ∞0U1QIIρ dρ−12b2(1− f)∫ ∞0U1QIρ dρ− pi min2≤i≤Nki)) . (5.4)Notice that the only key difference between these two expressions isthat mink∈Ω∗\0R0,k is replaced by min2≤i≤N ki. Therefore, we introduce aparameter c and our goal is to compare the solution to∂∂Sχ(Sc, f) = −1ν− c , (5.5)with the expressionSc = ν12√b(1− f)f(1 + ν(12b2∫ ∞0U1QIIρ dρ−12b2(1− f)∫ ∞0U1QIρ dρ−c2)),for some small . In Figure 5.2 we show numerical results that confirm thatthe full numerical results and asymptotic results agree rather well.Next, we identify the optimal lattice arrangement for the periodic case.As stated in Principal Result 3, the optimal lattice arrangement Λop withfixed primitive cell of area unity is the one which solve the following max-minproblem:arg maxΛ, |Ω|=1{ mink∈Ω∗\0{R0,k}} (5.6)585.2. Stability Threshold and the Optimal Lattice Arrangement−1 −0.5 0 0.5 10.4270.4280.4290.430.4310.4320.433cS c−1 −0.5 0 0.5 10.6750.680.6850.690.6950.70.7050.710.7150.72cS cFigure 5.2: Numerical solution to (5.5) and the two-term asymptotic ap-proximations for Sc with different c. Left panel: f = 0.4 and  = 0.01.Right panel: f = 0.5 and  = 0.05. The blue (top) curve is the numericalsolution while the green (bottom) one is the asymptotic approximations inboth cases.Therefore, if we want to find the optimal lattice numerically, we need toknow how to calculate R0,k. First we follow the process in [2] and [8], toderive an explicit expression for the Bloch Green function and its regularpart R0,k. Recall that the Bloch Green function satisfies:∆G0,k(x) = −δ(x), ∀x ∈ Ω ,G0,k(x+Li) = e−ik·LiG0,k(x) , ∀x ∈ d−i ,∂n−G0,k(x+Li) = e−ik·Li∂n+G0,k(x) , ∀x ∈ d−i .The free space Green’s function in the absence of any boundary conditionsis Gfree(x) = − 12pi ln |x|. We then observe that the infinite sumG(x) =∑l∈ΛGfree(x+ l)eik·l ,satisfies the PDE together with the quasi-periodic boundary conditions. To595.2. Stability Threshold and the Optimal Lattice Arrangementverify that it satisfies these boundary conditions we calculateG(x+Li) =∑l∈ΛGfree(x+Li + l)eik·l ,=∑Li+l∈ΛGfree(x+ (Li + l))eik·(Li+l)eik·(−Li) ,= e−ik·Li∑l′∈ΛGfree(x+ l′)eik·l′= e−ik·LiG(x) .The second line above follows since Li ∈ Λ and l′= Li + l.By the Poisson summation formula proved in (2.4) in Chapter 2, andthe fact that Ĝfree(ξ) = 1|ξ|2 and |Ω| = 1, we haveG(x) =∑l∈ΛGfree(x+ l)eik·l =∑d∈Λ∗Ĝfree(d− k)eix·(d−k) =∑d∈Λ∗eix·(d−k)|d− k|2.It is easy to prove from an integral test that the last series is not absolutelyconvergent. However, we can show it is actually conditionally convergent forx 6= 0 by decomposing it into two parts as was done in [2]. We pick some ηin η ∈ (0, 1), and rewrite the infinite series as∑d∈Λ∗eix·(d−k)|d− k|2=∑d∈Λ∗eix·(d−k)|d− k|2(1− e− |d−k|24η2 + e− |d−k|24η2 ) ,=∑d∈Λ∗eix·(d−k)|d− k|2e− |d−k|24η2 +∑d∈Λ∗eix·(d−k)|d− k|2(1− e− |d−k|24η2 ) .The first term is an absolutely convergent series, which we denote asGFourier(x) =∑d∈Λ∗eix·(d−k)|d− k|2e− |d−k|24η2 . (5.7)We claim that the second term equals another absolutely convergent seriesover the original lattice Λ. We can write this series in a convenient formusing the following lemma:605.2. Stability Threshold and the Optimal Lattice ArrangementLemma 5.∑d∈Λ∗eix·(d−k)|d− k|2(1− e− |d−k|24η2 ) =∑l∈ΛFsing(x+ l)eik·l , ∀x ∈ Ω\0 , (5.8)where Fsing(x) ≡ 12pi∫∞ln(2η) e− |x|24 e2sds = 14piE1(η2|x|2), and E1(z) is theexponential integral defined by E1(z) =∫∞z t−1e−t dt.Proof. Firstly, we observe that 12pi∫∞ln(2η) e− |x|24 e2sds = 14piE1(η2|x|2) by us-ing a simple change of variables. Then, according to [1], the exponentialintegral E1(z) has the decay property, E1(z) < e−z ln(1 + 1z ), so that the se-ries over Λ given by the right hand-side of (5.8) converges absolutely. Then,by using the Poisson summation formula as proved in (2.4), we get∑l∈ΛFsing(x+ l)eik·l =∑d∈Λ∗F̂sing(d− k)eix·(d−k) . (5.9)Upon comparing this result with the statement that we want to prove, weneed only show that F̂sing(ξ) = 1|ξ|2 (1 − e− |ξ|24η2 ). To prove this we showthat the inverse Fourier transform of the right hand-side is Fsing(x). No-tice that both Fsing(x) and the right hand-side are radially symmetric,and that the inverse Fourier transform of a radially symmetric functionis the inverse Hankel transform of order zero (cf. [11]), so that f(r) =(2pi)−1∫∞0 fˆ(ρ)J0(ρr)ρ dρ. Upon using the well-known inverse Hankel trans-form (cf. [11]) ∫ ∞0e−ρ2e−2sρJ0(ρr) dρ =12e2s−r2 e2s4 ,and the fact that1ρ2(1− e− ρ24η2 ) = 2∫ ∞ln(2η)e−ρ2e−2s−2sds ,we calculate the inverse Fourier transform of 1ρ2 (1− e− ρ24η2 ) as12pi∫ ∞0(1ρ2(1− e− ρ24η2 ))J0(ρr)ρ dρ =12pi∫ ∞0(∫ ∞ln(2η)2e−ρ2e−2s−2s ds)J0(ρr)ρ dρ615.2. Stability Threshold and the Optimal Lattice Arrangement=1pi∫ ∞ln(2η)e−2s(∫ ∞0e−ρ2e−2sρJ0(ρr) dρ)ds =12pi∫ ∞ln(2η)e−2se2s−r24 e2sds=12pi∫ ∞ln(2η))e−r24 e2sds .Thus, we conclude that F−1( 1|ξ|2 (1 − e− |ξ|24η2 )) = Fsing(x), which completesthe proof of the lemma.We remark here that the series over the reciprocal lattice Λ∗, given bythe left hand-side of (5.8), is only conditionally convergent, while the seriesover the original lattice Λ, given by the right hand-side of (5.8), convergesabsolutely and we denote it byGspatial(x) =∑l∈ΛFsing(x+ l)eik·l . (5.10)In this way, we have an explicit expression for the Bloch Green func-tion G0,k = G(x), and have separated it into the sum of two absolutelyconvergent series asG(x) = GFourier(x) +Gspatial(x) . (5.11)We remark here that strictly speaking the demonstration above is not com-pletely rigorous. The Poisson summation formula proved previously requiresthat the function to be in L1, but Gfree(x) is not. The way to circumventthis technical difficulty is to first define the two absolutely convergent se-ries GFourier(x) and Gspatial(x) as in (5.7) and (5.10). Then, we defineG(x) = GFourier(x) + Gspatial(x) and simply prove it satisfies the differen-tial equation and the quasi periodic boundary conditions. Notice that G(x)is independent of the choice of η since ∀x ∈ Ω\0,GFourier(x) +Gspatial(x) =∑d∈Λ∗eix·(d−k)|d− k|2e− |d−k|24η2 +eix·(d−k)|d− k|2(1− e− |d−k|24η2 ) ,=∑d∈Λ∗eix·(d−k)|d− k|2,625.2. Stability Threshold and the Optimal Lattice Arrangementwhich is independent of η.Next we calculate the singular behaviour of G(x) as x → 0. The termGFourier(0) is finite so this term is readily calculated. However, in the seriesGspatial(x), there is a singularity as x → 0 for the term corresponding tol = 0, owing to the fact that the exponential integral E1(z) has a singularityat 0. Upon using the well-known series expansion of E1(z)E1(z) = −γ − ln(z)−∞∑n=1(−1)nznnn!, | arg z| < pi , (5.12)as given in §5.1.11 of [1], where γ = 0.57721 · · · is Euler’s constant, we derivethatFsing(x) ∼ −γ4pi−ln η2pi−ln |x|2pi+ o(1), as x→ 0. (5.13)This shows that G(x) has the expected logarithmic singularity as x → 0,and that the regular part of the Bloch Green’s function isR0,k = limx→0(G(x) +12piln |x|), (5.14)=∑d∈Λ∗1|d− k|2e− |d−k|24η2 +∑l∈Λ\0eik·lFsing(l)−γ4pi−ln η2pi.We remark here that if we take conjugate of this expression, we get thesame quantity due to the symmetry of the lattice. This gives an alternativeproof that R0,k is real-valued. In addition, since R0,k only depends on Ωand k, the expression above should be independent of the choice of η. Toestablish this result, we take the derivative of (5.14) with respect to η andprove it vanishes. Upon differentiating (5.14), we obtain∂∂ηR0,k =∑d∈Λ∗12η3e− |d−k|24η2 −12piη∑l∈Λ\0e−|l|2η2eik·l −12piη. (5.15)635.2. Stability Threshold and the Optimal Lattice ArrangementTo show that this expression vanishes, it is equivalent to show that∑d∈Λ∗piη2e− |x−d|24η2 = 1 +∑l∈Λ\0e−|l|2η2eix·l =∑l∈Λe−|l|2η2eix·l. (5.16)Notice that the left hand-side is an absolutely convergent series and anintegrable function due to the exponential decay. Thus, as we have shown in(2.2), it can be decomposed into a Fourier series of eix·l, where l ∈ (Λ∗)∗ = Λand the coefficient of eix·l is calculated as1|Ω∗|∫Ω∗f(y)e−iy·l dy =1|Ω∗|∫Ω∗(∑d∈Λ∗piη2e− |y−d|24η2 )e−iy·l dy(∀d ∈ Ω∗, ∀l ∈ Ω, l · d = 2kpi, k ∈ Z) =1|Ω∗|∫R2piη2e− |y|24η2 e−iy·l dy(the Fourier transform of a Gaussian) =1|Ω∗|piη24piη2e−|l|2η2(|Ω| = 1, then |Ω∗| = 4pi2) = e−|l|2η2 .This establishes that ∂∂ηR0,k = 0, which yields that R0,k is independent ofη.The explicit expression (5.14) provides a way to calculate R0,k numer-ically. Since the two series converge absolutely, for a fixed lattice Λ with|Ω| = 1, we can truncate Λ, Λ∗ by a finite subset to get a good approxi-mation of R0,k. We then minimize it numerically over k ∈ Ω∗\0. Noticethat Lemma 2 is useful here since it tells us that R0,k blows up as k → 0,thus we can minimize it away from 0. Then, we change the lattice andmaximize R(Λ) ≡ mink∈Ω∗\0{R0,k} over different lattices with |Ω| = 1. Thenumerical results shown in [8] indicates that R(Λ) is maximized for a regularhexagonal lattice Λop and that R(Λ∗) = −0.079124. For a regular hexagonallattice, Table 5.1 compares the numerical results for the stability threshold,measured in terms of the source strength, and the corresponding one- andtwo-term asymptotic approximations.645.3. Case Study: N Peaks on a RingLattice Sc Leading order Two term approximation = 0.1 1.3854 1.3603 1.3706 = 0.01 0.4306 0.4302 0.430526Table 5.1: Source strength threshold and its asymptotic approximation fora regular hexagonal lattice with |Ω| = 1 and f = Case Study: N Peaks on a RingIn this section, we implement our stability theory for the finite domain fora particular arrangement of spots inside the unit disk Ω = D1. We take 5points {xi}5i=1 equally distributed on a circle of radius 0.5 concentric withinthe unit disk, as shown in Figure 5.3. The centers of the localized spotscorresponds to the locations of these points.Figure 5.3: 5 localized spots on a ring concentric within the unit disk.For this special symmetric configuration of 5 equally-spaced spots ona ring, the corresponding Neumann Green matrix G has a constant rowsum, which implies that e is an eigenvector of G. This spot configurationconsisting of equally-spaced spots is one of the simplest ways to ensure that655.3. Case Study: N Peaks on a Ringe is an eigenvector of G.For the unit disk, there is an explicit formula for the Neumann Green’sfunction G0(x, ξ) and its regular part. As derived in [9], we haveG0(x, ξ) =12pi(− ln |x− ξ| − ln ||ξ|x−1|ξ|ξ|+12(|x|2 + |ξ|2)−34) , (5.17)and thus the regular part R0(x) is given byR0(x) =12pi(− ln ||x|x−1|x|x|+ |x|2 −34) . (5.18)Without loss of generality we label the spot locations xi on the ringas xi = (12 cos2pi(i−1)5 ,12 sin2pi(i−1)5 )T , for i = 1, 2, ..., 5. We then substitutethis into (5.17) and (5.18) to obtain G. By using Matlab, we numericallycalculate all of the eigenvalues of G as k1 = −0.2126, k2 = k3 = 0.1392,and k4 = k5 = −0.1174. Next, we choose the smallest eigenvalue other thanthe one that corresponding to e. This is the eigenvalue k4 = k5 = −0.1174,which we then use to calculate the stability threshold in terms of the sourcestrength. This allows us to numerically evaluate the second order term inthe stability threshold. The full numerical results for the stability threshold,measured in terms of the source strength, are compared versus the one- andtwo-term asymptotic results in Table 5.2.Lattice Sc Leading order Two term approximation = 0.05 0.9742 0.9713 0.9619 = 0.02 0.6110 0.6083 0.6107Table 5.2: The stability threshold in terms of the source strength S and itsone- and two-term asymptotic approximation for a 5 spot pattern on a ringof radius 0.5 concentric within the unit disk with f = 0.4.We remark that in this case, there is an analytical way to determine theeigenvalues of the Neumann Green matrix G. Since this matrix is cyclic,its eigenvectors are qi = (1, ωi−1, (ωi−1)2, ..., (ωi−1)n−1)T , for i = 1, 2, ..., n,while the corresponding eigenvalue is f(ωi−1), where ω is the n-th root of665.3. Case Study: N Peaks on a Ringunity. Here f(x) =∑nk=1 ckxk−1 and (c1, c2, ..., cn) is the first row of G.We can calculate the eigenvalues in this way and obtain the same results asgiven above by a direct numerical calculation of the eigenvalues by Matlab.67Chapter 6SummaryIn this thesis we have studied the linear stability of steady-state localizedspot patterns for a singularly perturbed Brusselator reaction-diffusion sys-tem in both a periodic and finite domain setting. For both problems, thereis a stability threshold Dc ∼ O(− 1ln ) that characterizes a zero eigenvaluecrossing. We have calculated a two term asymptotic approximation forDc through an asymptotic solution of a singularly perturbed linear eigen-value problem. In the periodic setting, we first use Floquet-Bloch theory toconvert the whole plane problem into a problem posed on a primitive celltogether with the Bloch boundary conditions. Then we obtain the leadingorder approximation for Dc by analyzing a leading order nonlocal eigenvalueproblem (NLEP) derived using the method of matched asymptotic expan-sions. This leading order NLEP is independent of the geometry of the latticeΛ and the Bloch vector k. In order to characterize the effect of the latticeand the Bloch vector on the stability threshold, we calculated a higher or-der approximation for Dc by imposing a solvability condition to the nextorder equations. The calculation leads to a formula for a real-valued contin-uous band of spectra of the linearization that lies within a small ball nearthe origin in the spectral plane when D is near the leading order stabilitythreshold. The refined approximation to the stability threshold is obtainedfrom the requirement that this band of spectrum lies in the left half of thespectral plane. The correction to the leading order stability threshold ob-tained in this way depends on the regular part R0,k of the Bloch Greenfunction, which in turn is determined by the lattice and the Bloch vector k.An explicit formula for R0,k is also derived for numerical computation usingEwald summation methods. This formula is used to determine the optimallattice arrangement which allows for the largest stability threshold.68Chapter 6. SummaryThe analysis for the finite domain problem is similar, with the key dif-ference being that the N spots interact with each other through a NeumannGreen matrix G. For a pattern with arbitrarily-located spots, this leads usto analyze N distinct problems, one near each of the spots. For simplicity,we restrict the locations of the spots so that the spots have a common sourcestrength. In this way, the local problem near each of the spots is the same.By decomposing the solution to the linearized problem into the directionsof the eigenvectors of G, the analysis becomes very similar to that for theperiodic problem. More specifically, we obtain the leading order approxima-tion for the stability threshold through an NLEP that is the same in eachof N − 1 directions. We then calculate the second order approximation byimposing a solvability condition in each direction on the second order terms.This higher order approximation to the stability threshold depends on thematrix eigenvalues of G.For both the periodic and finite domain problems, we also provide a quickway to derive the stability threshold, which avoids any detailed calculationof spectra near the origin in the spectral plane. This simplified analysisshows that the stability threshold can be determined by solving a nonlinearequation. Numerical comparison between the two-term approximation forthe stability threshold in terms of the source strength Sc and the resultsobtained from solving the nonlinear equation is provided. For the finitedomain problem we illustrate our theory for a case study of N = 5 equally-spaced localized spots on a circular ring that is concentric with the unitdisk.There are some open problems suggested by this study. Firstly, forthe periodic case, numerical evidence obtained from computing the regularpart of the Bloch Green’s function indicates that it is the regular hexago-nal lattice that offers the optimum stability threshold. However, it wouldbe preferable to obtain a rigorous analytic proof of this result. Secondly,although we have employed a systematic asymptotic method to calculate arefined approximation to the stability threshold for the finite domain prob-lem, it would be interesting to extend the rigorous leading-order analysis in[16] to rigorously derive the second order term for the stability threshold.69Chapter 6. SummaryThirdly, it would be interesting to try to extend the rigorous frameworkof [16] to rigorously analyze the periodic problem. The technical difficultyhere is that, in contrast to the finite domain problems considered in [16] thathave discrete spectra, the periodic problem requires analyzing the edges ofa band of continuous spectra. Fourthly, it would be interesting to give aprecise relationship between the stability threshold for a multi-spot patternwith regularly spaced spots on a very large but finite domain and that forthe periodic problem. It is expected that the stability thresholds for thesetwo problems would be similar, with the only difference being essentiallythe perturbing effect of a distant domain boundary. More specifically, uponcomparing the stability thresholds for the periodic and finite domain prob-lems, we identify a formal correspondence that the regular part of the BlochGreen function R0,k is replaced by the eigenvalues ki of the Neumann Greenmatrix. Since we may view the periodic case as the limit of a truncatedlattice, i.e. ΛN = {n1l1 + n2l2||ni| ≤ N, i = 1, 2}, the question then ishow do the matrix eigenvalues ki approximate, or discretize, the continuousband R0,k? Finally, we remark that the analysis in this thesis has focusedon determining refined formulae for the stability thresholds associated withO(1) eigenvalues that result from zero eigenvalue crossings. However, it iswell-known that there are additional small eigenvalues of order λ ∼ O(2)that are associated with the translation modes. Unstable eigenvalues in thisclass lead to weak instabilities that are only manifested over very long timeintervals. It would be interesting to calculate the stability thresholds forthese eigenvalues for both the periodic and finite domain problems. For thefinite domain problem, a leading order analysis of these eigenvalues is givenin [16] for a related Gierer-Meinhardt reaction-diffusion system.70Bibliography[1] M. Abramowitz, I. A. Stegun, et al. Handbook of mathematical func-tions, volume 1. Dover New York, 1972.[2] G. Beylkin, C. Kurcz, and L. Monzo´n. Fast algorithms for helmholtzgreen’s functions. Proceedings of the Royal Society A: Mathematical,Physical and Engineering Science, 464(2100):3301–3326, 2008.[3] W. Chen and M. J. Ward. The stability and dynamics of localizedspot patterns in the two-dimensional gray-scott model. arXiv preprintarXiv:1009.2805, 2010.[4] X. Chen and Y. Oshita. An application of the modular function in non-local variational problems. Archive for Rational Mechanics and Analy-sis, 186(1):109–132, 2007.[5] B. Gidas, W. M. Ni, and L. Nirenberg. Symmetry of positive solutionsof nonlinear elliptic equations in Rn. Adv. Math. Suppl. Stud. A, 7:369–402, 1981.[6] A. Gierer and H. Meinhardt. A theory of biological pattern formation.Kybernetik, 12(1):30–39, 1972.[7] P. Gray and S. K. Scott. Sustained oscillations and other exotic pat-terns of behaviour in isothermal reactions. The Journal of PhysicalChemistry, 89(1):22–32, 1985.[8] D. Iron, John R., M. J. Ward, and J. Wei. Logarithmic expansions andthe stability of periodic patterns of localized spots for reaction-diffusionsystems in R2. Journal of nonlinear science, submitted.71Bibliography[9] T. Kolokolnikov, M. S. Titcombe, and M. J. Ward. Optimizing thefundamental neumann eigenvalue for the laplacian in a domain withsmall traps. European Journal of Applied Mathematics, 16(02):161–200, 2005.[10] T. Kolokolnikov, M. J. Ward, and J. Wei. Spot self-replication anddynamics for the schnakenburg model in a two-dimensional domain.Journal of nonlinear science, 19(1):1–56, 2009.[11] R. Piessens. The hankel transform. The Transforms and ApplicationsHandbook, 2:9–1, 2000.[12] I. Prigogine and R. Lefever. Symmetry breaking instabilities in dissi-pative systems. ii. The Journal of Chemical Physics, 48(4):1695–1700,1968.[13] I. Rozada, S. J. Ruuth, and M. J. Ward. The stability of localized spotpatterns for the brusselator on the sphere. SIAM Journal on AppliedDynamical Systems, 13(1):564–627, 2014.[14] A. M. Turing. The chemical basis of morphogenesis. PhilosophicalTransactions of the Royal Society of London. Series B, Biological Sci-ences, 237(641):37–72, 1952.[15] V. K. Vanag and I. R. Epstein. Localized patterns in reaction-diffusionsystems. Chaos: An Interdisciplinary Journal of Nonlinear Science,17(3):037110, 2007.[16] J. Wei and M. Winter. Spikes for the two-dimensional gierer-meinhardtsystem: the weak coupling case. Journal of Nonlinear Science,11(6):415–458, 2001.[17] M. Winter and J. Wei. Existence, classification and stability analysisof multiple-peaked solutions for the gierer-meinhardt system in R1.Methods and Applications of Analysis, 14(2):119–164, 2007.72


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items