Pattern Formation on CurvedSurfacesbyLaurent CharetteA THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Mathematics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)March 2019c© Laurent Charette 2019The following individuals certify that they have read, and recommend to the Facultyof Graduate and Postdoctoral Studies for acceptance, the dissertation entitled:Pattern Formation on Curved Surfacessubmitted by Laurent Charette in partial fulfillment of the requirements forthe degree of Doctor of Philosophyin MathematicsExamining Committee:Wayne Nagata, MathematicsCo-SupervisorColin B. Macdonald, MathematicsCo-SupervisorCarl Ollivier-Gooch, Mechanical EngineeringUniversity ExaminerBhushan Gopaluni, Chemical and Biological EngineeringUniversity ExaminerShigui Ruan, Mathematics at The University of MiamiExternal ExaminerAdditional Supervisory Committee Members:Michael Jeffrey Ward, MathematicsSupervisory Committee MemberBrian Wetton, MathematicsSupervisory Committee MemberiiAbstractPatterns emerge in various growing biological organisms, like conifer embryos, of-ten from an homogeneous preceding state. The embryo tip, initially hemisphericalshaped gradually flattens as cotyledons arranged in a roughly regular pattern emerge.A common way to model these patterns is through reaction-diffusion systems of partialdifferential equations. This thesis relays results obtained studying such systems andprovides various new results about pattern formation on curved surfaces via diffusiondriven bifurcations.We first describe situations where, for certain critical values of domain or differentialequation parameters, the unpatterned solution is unstable to two different linear normalmodes. We use centre manifold and normal form reductions to analyze the existenceand stability of pure and mixed modes of nonlinear patterned solutions of the reaction-diffusion system, for parameters near two cases of critical values. In one case, thesystem reduces to a well known example of mode interaction. In the other case, themode interaction is new, due to very small quadratic terms in the normal form.We then perform a reduction of a nonautonomous Brusselator reaction-diffusionsystem of partial differential equations on a spherical cap with time dependent curvatureusing an asymptotic series expansion on the centre manifold reduction. Parameter valuesare chosen such that the change in curvature would cross critical values which wouldchange the stability of the patternless solution in the constant domain case. The non-isotropic nature of the domain evolution insert a small patterned component to thepreviously patternless state, which we call the ‘quasi-patternless solution’. The evolvingdomain functions and quasi-patternless solutions are derived as well as a method toobtain this nonautonomous normal form. The obtained reduction solutions are thencompared to numerical solutions.Finally we provide an adaptation of the closet point method to evolving domains. Weperform several convergence analysis experiments of the heat equation on test surfacesand obtain quadratic convergence. Then, a few reaction diffusion simulations are shownon evolving surfaces, some featuring non-isotropic evolution. The previously shownapplication of the Brusselator system on an evolving spherical cap are compared to acentre manifold reduction and also show quadratic convergence.iiiLay SummaryMotivated by cotyledon patterns in conifer embryos, we study pattern emergence insolutions of systems reaction-diffusion partial differential equations on evolving curvedsurfaces. We start by reviewing previously obtained results on a static domain using anormal form differential equation. Then, we cover the case where two possible patternscan emerge in a constant spherical cap. Next, we go over the reduction to a non-autonomous normal form in the case of a flattening spherical cap. Finally, we shownumerical results, using the closest point method, for a variety of different surfaces.ivPrefaceThis thesis is original work by the author, Laurent Charette.Figure 2.1 in Chapter two is taken from [43] and has been used with authorizationfrom the author.A version of chapter 3 has been published: Charette, L., and Nagata, W. Bifurcationof mixed mode reaction-diffusion patterns in spherical caps, International Journal ofBifurcation and Chaos, 28(6), 15 June 2018, 1830017. I found all the analytical results,performed all the numerical simulations and wrote some of the manuscript. W. Nagataalso participated in the writing as supervisory author.A version of chapter 4 has been submitted: Charette, L., Macdonald, C. B., andNagata, W. Pattern formation in a slowly flattening spherical cap: delayed bifurcation.Ifound all the analytical results, performed all the numerical simulations and wrote mostof the manuscript. W. Nagata suggested edits on the article and acted as supervisoryauthor. C. B. Macdonald provided the numerical simulation files that were adapted insection 4.3.5 “Numerical simulations” and acted as supervisory author on the numericalaspects of the paper.A version of chapter 5 will soon be prepared as an article and be submitted toa journal to be determined. I developed the numerical method, based on a previousmethod provided by C. B Macdonald, who will also act as supervisory author. W.Nagata will also participate as supervisory author on the theoretical dynamical systemsand motivation sections.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 Model Description and Analytical Background . . . . . . . . . . . . . 32.1 Biological Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32.2 Spherical Cap Domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.3 Diffusion-Driven Instability . . . . . . . . . . . . . . . . . . . . . . . . . 72.4 Reaction-Diffusion and Linear Stability in a Fixed Spherical Cap . . . . 92.4.1 Marginal Stability Curves . . . . . . . . . . . . . . . . . . . . . . 112.5 Bifurcation Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.6 Delayed Bifurcation Analysis . . . . . . . . . . . . . . . . . . . . . . . . 152.7 Closest Point Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.7.1 Spatial Discretization . . . . . . . . . . . . . . . . . . . . . . . . 192.7.2 Dirichlet Boundary Conditions . . . . . . . . . . . . . . . . . . . 202.7.3 Explicit Time-Stepping Scheme Example . . . . . . . . . . . . . 213 Bifurcation with Two Competing Modes in a Constant Spherical Cap 233.1 Case 1: Non-Zero Orders . . . . . . . . . . . . . . . . . . . . . . . . . . 253.2 Case 2: Zero Order in One of the Modes . . . . . . . . . . . . . . . . . . 293.3 Numerical Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 323.3.1 Case 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 333.3.2 Case 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33vi4 Delayed Bifurcation in a Slowly Changing Spherical Cap . . . . . . . 464.1 The Quasi-Patternless Solution . . . . . . . . . . . . . . . . . . . . . . . 494.2 Deviation Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 504.3 Pattern Formation in a Slowly Changing Domain . . . . . . . . . . . . . 514.3.1 Linearization About the Quasi-Patternless Solution . . . . . . . . 514.3.2 The Centre Subspace . . . . . . . . . . . . . . . . . . . . . . . . 534.3.3 Centre Manifold Equations . . . . . . . . . . . . . . . . . . . . . 554.3.4 Normal Form . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 574.3.5 Numerical simulations . . . . . . . . . . . . . . . . . . . . . . . . 585 Closest Point Method on Evolving Domains . . . . . . . . . . . . . . . 645.1 Reaction-Diffusion Systems in a Time Dependent Domain . . . . . . . . 645.1.1 Curves in Two-Dimensional Space . . . . . . . . . . . . . . . . . 655.1.2 Surfaces in Three Dimensions . . . . . . . . . . . . . . . . . . . . 685.2 Pattern Formation for Reaction-Diffusion Systems . . . . . . . . . . . . 685.3 Adaptation of the Closest Point Method on an Evolving Domain . . . . 695.3.1 Time Derivatives of Closest Point Extensions . . . . . . . . . . . 705.3.2 Eulerian Approach . . . . . . . . . . . . . . . . . . . . . . . . . . 705.3.3 Lagrangian Approach . . . . . . . . . . . . . . . . . . . . . . . . 715.3.4 Time-Stepping in an Evolving Domain . . . . . . . . . . . . . . . 725.3.5 Time-Stepping on Moving Surfaces . . . . . . . . . . . . . . . . . 735.4 Convergence of the Heat Equation . . . . . . . . . . . . . . . . . . . . . 785.5 Schnakenberg RDE on Different Types of 2D Domain Evolution . . . . . 805.5.1 Exponential Growth . . . . . . . . . . . . . . . . . . . . . . . . . 805.5.2 Logistic Growth . . . . . . . . . . . . . . . . . . . . . . . . . . . 815.5.3 Exponential Growth Along Bipolar Trajectories . . . . . . . . . . 845.5.4 Logistic Growth Along Bipolar Trajectories . . . . . . . . . . . . 855.5.5 Stretching Ellipse . . . . . . . . . . . . . . . . . . . . . . . . . . 865.6 Numerical Examples in 3D . . . . . . . . . . . . . . . . . . . . . . . . . 895.6.1 Schnakenberg System . . . . . . . . . . . . . . . . . . . . . . . . 905.6.2 Brusselator System on a Linearly Growing Sphere . . . . . . . . 915.6.3 Unevenly Stretching Ellipsoid . . . . . . . . . . . . . . . . . . . . 935.7 Brusselator System on a Flattening Spherical Cap . . . . . . . . . . . . 946 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101AppendicesA Computation of Normal Form Coefficients for Codimension 2 Bifurca-tions in a Constant Spherical Cap . . . . . . . . . . . . . . . . . . . . . . 106viiA.1 Case 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107A.2 Case 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108B Computation of Bifurcation Curves, for Codimension 2 Bifurcationsin a Constant Spherical Cap, Case 2 . . . . . . . . . . . . . . . . . . . . 110B.1 The Quadratic Normal Form . . . . . . . . . . . . . . . . . . . . . . . . 110B.1.1 Bifurcations of the Origin . . . . . . . . . . . . . . . . . . . . . . 111B.1.2 Bifurcations of x -Modes . . . . . . . . . . . . . . . . . . . . . . . 111B.1.3 Bifurcations of Mixed Modes . . . . . . . . . . . . . . . . . . . . 112C Time Derivative of the Laplace–Beltrami Eigenfunction . . . . . . . . 113D Exact Solutions of the Heat Equation on an Evolving Domain . . . . 115D.1 Exponentially Growing Circle . . . . . . . . . . . . . . . . . . . . . . . . 115D.2 Exponentially Growing Sphere . . . . . . . . . . . . . . . . . . . . . . . 116viiiList of Tables5.1 Convergence of the heat equation on a uniformly, exponentially growingcircle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 785.2 Convergence of the heat equation on a uniformly, exponentially growingcircle with smaller time step. . . . . . . . . . . . . . . . . . . . . . . . . . 795.3 Convergence of the heat equation on a uniformly, exponentially growingcircle with even smaller time step. . . . . . . . . . . . . . . . . . . . . . . 795.4 Convergence of the heat equation on a uniformly rotating circle. . . . . . 795.5 Convergence of the heat equation on a simultaneously growing and rota-ting circle using the Lagrangian approach. . . . . . . . . . . . . . . . . . 795.6 Convergence of the heat equation on a simultaneously growing and rota-ting circle using the Eulerian approach. . . . . . . . . . . . . . . . . . . . 795.7 Convergence of the heat equation on a uniformly, exponentially growingsemicircle with Dirichlet boundary conditions. . . . . . . . . . . . . . . . 795.8 Convergence of the heat equation on a uniformly, exponentially growingsphere. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 805.9 Convergence of the heat equation on a uniformly, exponentially growingsphere over a longer time period. . . . . . . . . . . . . . . . . . . . . . . 805.10 Comparison of heat equation solutions on a stretching ellipse with aninitial circle of radius one shape and final semi-major and semi-minoraxes of lengths 2 and 0.75 respectively. . . . . . . . . . . . . . . . . . . . 89ixList of Figures2.1 Scanning electron micrographs showing conifer embryo cotyledon deve-lopment. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.2 Cross section of a spherical cap. . . . . . . . . . . . . . . . . . . . . . . . 62.3 Curves of constant bipolar coordinates η and ξ. . . . . . . . . . . . . . . 72.4 Marginal stability curves for four different (m,n) modes. . . . . . . . . . 132.5 Sketch of a delayed pitchfork bifurcation solution along with the constantparameter bifurcation diagram. . . . . . . . . . . . . . . . . . . . . . . . 172.6 Sketch of a forward Euler closest point method time step. . . . . . . . . . 223.1 Marginal stability curves with highlighted intersections. . . . . . . . . . . 243.2 Hopf-Hopf parametric portrait for Case 1 amplitude equations. . . . . . . 283.3 Hopf-Hopf phase portraits for Case 1 amplitude equations. . . . . . . . . 293.4 Bifurcation diagram for the Case 2 normal form for a large σ1, σ2 range. 313.5 Bifurcation diagram for the Case 2 normal form for a decreasing σ1, σ2range. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 353.6 Bifurcation diagram for the Case 2 normal form for a very small σ1, σ2range. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 363.7 Phase portraits of the Case 2 normal form equations. . . . . . . . . . . . 373.8 Phase portraits of the Case 2 normal form equations (continued). . . . . 383.9 Phase portraits of the Case 2 normal form equations (end). . . . . . . . . 393.10 A bifurcation diagram of a Z2-symmetric pitchfork bifurcation and twopossible results of a small symmetry-breaking perturbation. . . . . . . . . 403.11 Effects of changing the curvature γ on the critical eigenvalues and ampli-tude equation parameters σ1, σ2. . . . . . . . . . . . . . . . . . . . . . . 413.12 Effects of changing the cap radius R on the critical eigenvalues and am-plitude equation parameters σ1, σ2. . . . . . . . . . . . . . . . . . . . . . 413.13 Effects of changing the parameter A on the critical eigenvalues and am-plitude equation parameters σ1, σ2. . . . . . . . . . . . . . . . . . . . . . 423.14 A finite element simulation showing a stable pure (5, 1) mode steady stateof the Brusselator in a spherical cap, viewed from above. . . . . . . . . . 433.15 A finite element simulation showing a stable pure (3, 2) mode steady stateof the Brusselator in a spherical cap, viewed from above. . . . . . . . . . 443.16 A finite element simulation showing a stable mixed mode steady state ofthe Brusselator in a spherical cap, viewed from above. . . . . . . . . . . . 45x4.1 Flattening spherical cap normal form solutions for different ε values andthe constant domain bifurcation diagram. . . . . . . . . . . . . . . . . . . 584.2 Quasi-patternless solution of the Brusselator system on a flattening spher-ical cap with modified affine operator, simulated using the closest pointmethod. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 594.3 Quasi-patternless solution simulation computed using the closest pointmethod compared to a five-term truncated series approximation. . . . . . 604.4 Numerical simulations compared to the normal form solution as a functionof curvature index γ by extracting (5, 1) coefficients using a projection. . 624.5 Snapshots of numerical solutions for the U component of the Brusselatorsystem at four different curvature values. . . . . . . . . . . . . . . . . . . 635.1 Diagram showing the evolution of the closest point function on a surfacerotating with respect to the closest point. . . . . . . . . . . . . . . . . . . 715.2 Comparison of the Eulerian and Lagrangian approaches to backtracking. 745.3 Schnakenberg simulation on an exponentially growing circle. . . . . . . . 815.4 Band solution of the Schnakenberg system on an exponentially growingcircle for four different time. . . . . . . . . . . . . . . . . . . . . . . . . . 825.5 Circle solutions of each snapshot in Figure 5.4, obtained by interpolatingthe band solution onto regularly angled points on the circle. . . . . . . . 835.6 Schnakenberg Simulation on a logistically growing circle. . . . . . . . . . 845.7 Time plot of the solution on the circle as a function of the angle aroundits centre with divergence term precomputed. . . . . . . . . . . . . . . . 855.8 Time plot of the solution on the circle as a function of the angle aroundits centre with divergence term computed numerically. . . . . . . . . . . . 865.9 Schnakenberg Simulation on a logistically growing circle along bipolartrajectories. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 875.10 Schnakenberg Simulation on a logistically growing circle along bipolartrajectories with smaller final size. . . . . . . . . . . . . . . . . . . . . . . 875.11 Schnakenberg Simulation on a logistically growing circle along bipolartrajectories with even smaller final size. . . . . . . . . . . . . . . . . . . . 885.12 Three snapshots of Schnakenberg solutions on a stretching ellipse. . . . . 895.13 Three snapshots of Schnakenberg solutions on an ellipse expanding onone direction and contracting in the other. . . . . . . . . . . . . . . . . . 905.14 Three snapshots of Schnakenberg solutions on an exponentially growingsphere. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 915.15 Three snapshots of Schnakenberg solutions on an exponentially shrinkingsphere. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 925.16 Six snapshots of Schnakenberg solutions on an exponentially growing sp-here along toroidal coordinates. . . . . . . . . . . . . . . . . . . . . . . . 925.17 Six snapshots of Brusselator solutions on a linearly growing sphere. . . . 935.18 Snapshots of Schnakenberg solutions on an unevenly growing shape. . . . 94xi5.19 Convergence analysis plot for the closest point method Brusselator simu-lations on a flattening spherical cap. . . . . . . . . . . . . . . . . . . . . . 95xiiAcknowledgementsI want to thank the graduate students, faculty and staff of the UBC mathematics andthe Institute of Applied Mathematics department for contributing everyday to a friendlyand motivating environment.In particular I thank profs Wayne Nagata and Colin B. Macdonald for the help theyprovided by supervising this research.I would also like to thank profs Michael J. Ward and Brian Wetton for providingsome clarifications and feedback on more specific parts of the research subject.This research project would not have been possible without the financial support ofthe UBC department of Mathematics, the Institute of Applied Mathematics and theNatural Sciences and Engineering Research Council of Canada.Finally, I want to thank my family for being very supportive throughout my studies.xiiiChapter 1IntroductionPattern formation occurs in a wide array of phenomena. Some notable examples includecoloration of animals, swarming phenomena, and vortex arrangements in fluid dynamics[40, 9, 12]. There are also many potential applications in very active research fields suchas tumour research [18] and nanoparticle assembly [55]. Plants also exhibit differentpatterns through their branches, roots, leaves or petals [20, 7]. These plant developmentevents are often modelled by a set of partial differential equations and the emergenceof patterns would be bifurcations of the solutions of the equations, where a constantand uniform patternless state loses stability and the solution transitions to anotherpatterned, stable state. The different types of patterns can sometimes be controlled bya few parameters (size, curvature of domain, chemcal concentrations, temperature forexample).One theory of pattern formation, formulated by Turing [56], postulates that there arechemical agents called morphogens that react among each other in a reaction-diffusionprocess on plant or animal tissue. The tissue then becomes differentiated and growth isthen encouraged or suppressed according to the chemical pattern template. The growthhormone auxin is known to behave in many ways similarly to a morphogen and is anactive agent in plant growth [57].The goal of this project is not to observe or detect the morphogens, but ratherto provide mathematical models of the transition between a patternless state and dif-ferent patterns. The models obtain patterns that can be similar to those observed inexperiments. These patterns often are similar to harmonic functions. We use a reaction-diffusion system to model pattern formation, more specifically we use Brusselator [50]and Schnakenberg [54] reaction kinetics. Both these reaction-diffusion equations havebeen studied extensively and yield bifurcation solutions that can often be predictedusing linear analysis. In particular, the Brusselator system has been used previously inother pattern formation studies as well as in plant growth models [28, 29].Previous research on bifurcations of the Brusselator system in a constant sphericalcap by Nagata et al. [43] described a family of marginal stability curves in parameterspace. The curves are the interface between the patternless and patterned regimesand can also serve to predict the type of spotted or ringed pattern emerging from thepatternless solution. They also established a center manifold reduction of the Brusselatorpartial differential equation to a single ordinary differential equation representing thetime evolution of a coefficient to a harmonic spherical cap function.In order to approximate the solutions of the nonlinear reaction-diffusion systemswe can employ numerical methods such as finite elements and finite differences. The1former is well established and can be adapted easily to match any surface for whichwe can produce a mesh. The latter, however, is much simpler to use and has betterpointwise convergence than the finite elements method, but is usually restricted tosimpler domains. In this thesis we use the closest point method, which uses finitedifference concepts on more complicated surfaces. Introduced by Ruuth and Merriman[53], the closest point method is an embedding method whose main concept is to solvethe solution of the partial differential equations on a band surrounding the surface ofcurve with the restriction that solutions must be constant on any line normal to thedomain.The closest point method has been used on evolving surfaces before in computergraphics applications [30], especially for fluid effects on moving interfaces [2], [31]. Themethod has also been used paired with a level set method on an evolving surface, inorder to generate new closest point functions [44]. The closest point method has alsobeen used on nonautonomous systems of reaction-diffusion that have been rescaled fromevolving domain simulations [60], [52].The aim of this thesis is to further the research on pattern formation in curvedsurfaces and the numerical methods used to simulate these pattern phenomena. Inchapter 2, following this Introduction, we present the spherical cap model along withan outline of some of the single mode bifurcation results already mentioned here. Then,in Chapter 3 we consider the special cases when two possible pattern conditions aresatisfied simultaneously, giving rise to a larger number of equations and the possibilityof mixed patterns. Following that is Chapter 4 where we document the centre manifoldreduction method used in the context of a slowly flattening spherical cap, the first suchreduction performed on a curved surface evolving non-uniformly. Both this chapter andthe preceding one will feature numerical simulations supporting the analytical results.Chapter 5 will introduce and study an adaptation of the closest point method to evolvingsurfaces, where convergence tests are performed as well as several patterned reaction-diffusion experiments in curves and surfaces. The thesis concludes with Chapter 6, wherewe summarize and discuss our results, as well as suggest future research directions.2Chapter 2Model Description and AnalyticalBackgroundWe start by giving a brief description of developing conifer embryo tips, of Turingbifurcations and of the reaction-diffusion plant tip model used in much of this document.We also give some background information on the analytical and numerical methods usedin this thesis. The plant tip model is the same one used in [43], with the only differencebeing a time dependent curvature parameter in Chapters 4 and 5. We give an overviewof marginal stability curves and centre manifold reduction in the spherical cap for singlemode solutions of the constant curvature domain and give a short introduction on theconcept of delayed bifurcations in a simple ordinary differential equation system with atime dependent parameter. We finish this section by briefly explaining the principles ofthe closest point method, which we use in numerical simulations in Chapters 4 and 5.2.1 Biological MotivationPlants exhibit a variety of different shapes and forms during growth. One may onlythink of the various branch, root, leaf and petal arrangements displayed in specimensfor examples. These deformation of the plant tissue are controlled by chemical agentsthat soften the cell walls, which then allows these outgrowth patterns to emerge [46].In order to further understand these mechanisms there have been numerous studies onnonlinear chemical kinetic dynamics underlying pattern formation [57], [7], [20].One specific pattern formation event of interest here is the branching pattern inplant embryo. Previous studies performed on conifer embryos found that the formationof cotyledons, or embryonic leaves, occurs in regular patterns that can be mapped toharmonic functions of the tip surfaces, such as Bessel functions for the disk or spher-ical harmonic functions with Legendre polynomials for hemispheres. Figure 2.1 showsscanning electron micrographs of the emergence of a common pattern of five cotyledons.The number of cotyledons was demonstrated to depend on the embryo size [27], andthus have fixed spacing. Other parameters that can influence the pattern number aredomain curvature and boundary conditions. Another observation is that shortly prior tocotyledon emergence there is an observable flattening of the embryo tip. This is shownin panel B of Figure 2.1. The flattening embryo tip motivates us to study the effect ofthe tip curvature on pattern formation.3Figure 2.1: Scanning electron micrographs showing conifer embryo cotyledon develop-ment. The shape of the initial embryo tip is roughly hemispherical, as shown in panelA, then flattens prior to the emergence cotyledons as ‘bumps’, as shown on panels B, Cand D. The final cotyledon pattern, containing 5 of roughly similar size arranged in apattern similar to a (5, 1) mode referred in the later sections. Taken from [43].42.2 Spherical Cap DomainWe model the tip of a plant embryo by a spherical cap at the end of a cylindrical stalk.A spherical cap is a part of a sphere with a circular base. One may visualize it as thesurface obtained if a sphere was sliced using a horizontal plane, intersecting the plane atthe base. With a fixed radius of the base R and radius of the sphere % = R/γ, dependingon a positive curvature parameter γ, the caps can range over values 0 < γ ≤ 1, whereγ = 1 corresponds to a hemisphere and the limit γ → 0 corresponds to a flat disk. Inrectangular coordinates (x1, x2, x3) we align the centre of the cap on the x3 axis with itscircular boundary in the x1x2-plane.With this orientation we can use spherical coordinates to parametrize the surfacex1 =Rγsin θ cosϕ, x2 =Rγsin θ sinϕ, x3 =Rγ(cos θ −√1− γ2),θ ∈ [0, θmax), ϕ ∈ [0, 2pi], γ ∈ (0, 1],(2.1)where the maximal co-latitude angle θmax is dependent on the curvature of the cap, andthe longitudinal angle ϕ is the rotation angle around the x3-axis (also called the polaraxis) computed from the x1-axis. We restrict ourselves to a maximal co-latitude angleup to pi/2 to account for plant tip observations. In spherical coordinates the sphericalcap can be described asΩ = {(θ, ϕ)| 0 ≤ θ < θmax, 0 ≤ ϕ ≤ 2pi} . (2.2)The variable R is also the radius of the cylindrical stalk and γ is the sine of the anglethe cap makes with the x1x2-plane at the boundary. We can relate that parameter withthe maximal angle usingθmax = arcsin γ. (2.3)Figure 2.2 illustrates how the different parameters describe the cap.Another useful set of coordinates to describe spherical caps are bipolar and toroidalcoordinates. Bipolar coordinates are defined by the following formulasx1 =R sinh ηcosh η − cos ξ , x2 =R sin ξcosh η − cos ξη ∈ [0,∞), ξ ∈ [pi/2, pi).(2.4)The curves of constant ξ will create a circle arc with endpoints at (−R, 0) and (R, 0)and centre at (0, R cot(ξ)). We may extend the arc into a full circle by adding the arcof ξ + pi or equivalently extending the η values with η + pii. Figure 2.3 shows curves ofconstant η (in blue) and ξ (in red).Toroidal coordinates are the extension of bipolar coordinates into three-space, usinga longitudinal angle ϕx1 =R sinh η cosϕcosh η − cos ξ , x2 =R sinh η sinϕcosh η − cos ξ , x3 =R sin ξcosh η − cos ξη ∈ [0,∞), ϕ ∈ [0, 2pi], ξ ∈ [pi/2, pi).(2.5)5x1, x2x3γ = 1γ = γo2RC̺ = Rγθmaxγ = 0Figure 2.2: Cross section of a spherical cap along its diameter and the x3 axis withcurvature parameter γ, cap radius R and spherical radius %. The spherical centre of thecap is noted by C. Dashed curves describe hemispherical cap (γ = 0, C at the origin)and disk (γ = 0, C at infinity).In toroidal coordinates (2.5), a spherical cap is a surface with a constant ξ value,Ω = {(η, ϕ)|0 ≤ η <∞, 0 ≤ ϕ ≤ 2pi} . (2.6)The variable ξ measures the role of curvature here and is related to γ byξ = pi − arcsin γ. (2.7)Also worth noting is that the bounds of the toroidal coordinates (η, ϕ) on the sphericalcap do not depend on the value of ξ, unlike θmax in spherical coordinates, which dependson γ.6Figure 2.3: Curves of constant bipolar coordinates η and ξ. Circle arcs correspondingto a constant value of ξ are shown in red and circles of constant η are shown in blue. Inthis case we selected the focal distance R = 12.3 Diffusion-Driven InstabilityIn this section we undergo the general conditions for a general reaction-diffusion systemof two equations on a surface Ω∂X∂t= D∇2ΩX + f(X), (2.8)where ∇2Ω is the Laplace–Beltrami operator for surface Ω, D is a diagonal matrix of thediffusion coefficients and X = (X, Y )ᵀ is the solution in vector form, more details inSection 2.4. This was initially done by Turing [56] and restated in many reference works(e.g. [42], Chapter 2).Let X00 = (X00, Y00)ᵀ be a homogeneous steady state of system (2.8). Its sta-bility can be determined by the eigenvalues of the Jacobian of the function f(X) =(f(X, Y ), g(X, Y ))ᵀ, which represents chemical reaction kinetics, evaluated at X00.A0 =(fX fYgX gY .)(2.9)We need the patternless homogeneous solution to be stable, therefore both eigenvaluesneed to be negative, which is equivalent to the trace being negative and the determinantpositiveTr(A0) = fX + gY < 0 (2.10)det(A0) = fXgY − fY gX > 0. (2.11)7Now assume that the Laplace–Beltrami operator has a family of eigenfunctions Φj or-dered in eigenvalues∇2ΩΦj = −µjΦj, (2.12)for µj ≥ 0. For a small patterned eigenfunction u1 = u10 + u1jΦj, and similarly for u2,added to the patternless solution the Jacobian for the jth mode is, to leading orderA0j =(−DXµj + fX fYgX −DY µj + gY), (2.13)and its eigenvalues σ1j, σ2j must solve the characteristic equationσ2j − σj(Tr(A0)− (DX +DY )µj) +R(µj) = 0, (2.14)R(µj) = DXDY µ2j +DXµjgY +DY µjfX + det(A0), (2.15)where the partial derivatives of f and g are evaluated for the patternless solution values.Assuming that only the patternless solution is present, then the eigenvalues of everyeigenmode should all be negative, thereforeTr(A0j) = Tr(A0)− (DX +DY )µj < 0 (2.16)det(A0j) = R(µj) > 0. (2.17)For a Turing instability to occur we need one of the eigenmodes jc to become unstable,that is the sum of the eigenvalues of the Jacobian stay negative, but their productbecomes negative as well,Tr(A0j) = Tr(A0)− (DX +DY )µjc < 0 (2.18)det(A0j) = R(µjc) < 0. (2.19)Since DX , DY and µjc are all non negative the trace inequality will be preserved. Thedeterminant can be written as a quadratic function of µj that is concave up, so it canonly be negative if its discriminant is positive, or(DXgY +DY fX)2 > 4DXDY det(A0), (2.20)and we know that if µj = 0 then R(µj) = det(J0) > 0, therefore for R(µj) to be negativeboth its roots must be positive, thereforeDXgY +DY fX > 0. (2.21)Note that Turing instabilities cannot occur if the diffusion coefficients are equal, asconditions (2.10) and (2.21) cannot both be satisfied.82.4 Reaction-Diffusion and Linear Stability in aFixed Spherical CapIn this section we give a brief overview of simple (codimension one) bifurcations onspherical caps in a constant domain, a fixed spherical cap. We will also introduce thestructure and some notation that will be followed in the slowly changing domain case.In the constant domain case, the curvature γ is a constant parameter and will define aunique spherical cap domain Ω. We solve a reaction-diffusion equation on this sphericalcap domain (2.8), where X = (X, Y )ᵀ are the surface concentrations of the activechemical species, f(X) = (f(X), g(X))ᵀ give the chemical reaction kinetics, D is adiagonal matrix with positive constants DX and DY on the diagonal, and ∇2Ω is theLaplace–Beltrami operator on the surface Ω. The Laplace–Beltrami operator can beexpressed in spherical coordinates as∇2Ω =γR2 sin θ[∂∂θ(sin θ∂∂θ)+1sin2 θ∂2∂ϕ2], (2.22)or in toroidal coordinates as∇2Ω =(cosh η − cos ξ)2R2 sinh η[∂∂η(sinh η∂∂η)+1sinh η∂2∂ϕ2], (2.23)with periodic boundary conditions for ϕ at ϕ = 0, 2pi. The functions f and g give thechemical kinetics, for which we use the Brusselator:f(X) =(aA− dX − bBX + cX2YbBX − cX2Y), (2.24)where a, b, c, d are positive reaction rate constants and A, B are abundant chemical con-centrations, also considered to be positive constants. The algebraic system of equationsf(X, Y ) = g(X, Y ) = 0 in the Brusselator has a unique solutionX00 = (X00, Y00)ᵀ =(aAd,bBdaAc)ᵀ. (2.25)We take Dirichlet boundary conditions, setting X on the boundary to be equal to thispatternless solutionX = X00 on ∂Ω, (2.26)to match the concentrations on the cylindrical stalk, where they are assumed to beconstant.The reaction function f(X) has the Taylor expansion about the constant patternlesssolution X00,f(X) = f(X00) + fX(X00)(X−X00) + 12 fXX(X00)(X−X00,X−X00)+ 16fXXX(X00)(X−X00,X−X00,X−X00),(2.27)9wheref(X00) = 0, (2.28)and the different derivatives of f represent linear, then symmetric bilinear and trilinearfunctionsfX(X00) = K0,12fXX(X00) = B0,16fXXX(X00) = C0, (2.29)withK0 =(k1 k2k3 k4), (2.30)k1 = bB − d, k2 = a2A2cd2, k3 = −bB, k4 = −a2A2cd2, (2.31)andB0(U1,U2) =(1−1)[bBdaAU1U2 +aAcd(U1V2 + V1U2)], (2.32)C0(U1,U2,U3) =(1−1)c3[U1U2V3 + U1V2U3 + V1U2U3] , (2.33)forUj =(UjVj), j = 1, 2, 3. (2.34)We use the deviations U from the patternless solution X00X = X00 + U, U = (U, V )ᵀ (2.35)in order to study the stability of X00 and the bifurcating patterned solutions that emergefrom the patternless solution, now corresponding to U = 0. The reaction-diffusionsystem (2.8) becomes∂U∂t= D∇2ΩU + K0U + B0(U,U) + C0(U,U,U), (2.36)and the Dirichlet boundary conditions for U become homogeneous,U = 0 on ∂Ω. (2.37)The spherical cap domain and Brusselator system both have symmetries under rota-tion of any angle ϕ0 around the x3-axis and reflections through vertical planes containingthe x3-axis, generated by the two transformationsϕ→ ϕ+ ϕ0, ϕ→ −ϕ. (2.38)This implies that for any existing solution of the differential equation (2.8), rotationsand reflections of that solution will also be solutions.102.4.1 Marginal Stability CurvesTo determine the stability of the patternless solution, we linearize (2.36) about U = 0to obtain∂Uˆ∂t= A0Uˆ, Uˆ = 0 on ∂Ω, (2.39)whereA0 = D∇2Ω + K0. (2.40)We then use the linear stability ansatz,Uˆ = eσtU0 (2.41)to get the eigenvalue problem for A0,A0U0 = σU0. (2.42)To solve for the eigenvalues σ we putU0 = u0,mnΦmn, u0,mn ∈ R2 (2.43)where, for m = 0,±1,±2, · · · , n = 1, 2, 3, · · · the scalar function Φmn is an eigenfunctionof the Laplace–Beltrami operator with homogeneous boundary conditions,∇2ΩΦmn = −µmnΦmn, Φmn = 0 on ∂Ω. (2.44)The eigenfunctions Φmn are given in spherical coordinates byΦmn(θ, ϕ) = eimϕP|m|λmn(γ)(cos θ), (2.45)or in toroidal coordinates byΦmn(η, ϕ) = eimϕP|m|λmn(γ)(1− cosh η cos ξcosh η − cos ξ), (2.46)with ξ = pi− arcsin γ. The function P |m|λ is the associated Legendre function of the firstkind of integer order |m| and real degree λ. The value of λ = λmn(γ) is the nth positiveroot of the equationP|m|λ(√1− γ2)= 0. (2.47)The eigenvalues −µmn of the Laplace–Beltrami operator corresponding to the eigen-functions Φmn are given by−µmn = −λmn(γ)[λmn(γ) + 1]γ2R−2. (2.48)11Then for each m, n we obtain the algebraic eigenvalue problem in R2,A0,mnu0,mn = σu0,mn, (2.49)whereA0,mn =(−DXµmn + k1 k2k3 −DY µmn + k4)(2.50)We can find the eigenvalues σ as the roots of the characteristic equation,σ2 + σ[(DX +DY )µmn − (k1 + k4)] +DXDY µ2mn− (DXk4 +DY k1)µmn + k1k4 − k2k3 = 0.(2.51)We consider parameter values where the solutions of (2.51) are real and distinct σ = σ±mnwith σ+mn > σ−mn. If σ+mn < 0 for all values of m, n then the patternless solutionis asymptotically stable because all eigenvalues of (2.42) are negative. In this case,solutions will approach the patternless solution exponentially in time. If σ+m,n > 0 forsome m, n then the patternless solution is unstable.By fixing all parameters in the model except for two we can use (2.51) to obtain amarginal stability curve, a curve where the eigenvalue σ+mn = 0. We call the eigenvalueσ+mn the critical eigenvalue when it is close to this marginal curve and the mode (m,n)is called the critical mode. For example if we fix all parameters except for A and γ weobtain a marginal stability curve A = Amn(γ). We can compute these curves for any(m,n) pair to obtain a collection of curves such as in Figure 2.4 for R = 0.96747 andthe parameter valuesDX = 0.005, DY = 0.1,a = 0.01, bB = 1.5, c = 1.8, d = 0.375.(2.52)In this case for parameter values (A, γ) with A above all the curves the patternlesssolution is stable, while if A is below any one curve then the patternless solution isunstable.2.5 Bifurcation AnalysisIn this section we introduce the concepts of centre manifold reduction and equivariantnormal forms and as an example we walk through the single mode fixed spherical cappitchfork bifurcation analysis.The principle of centre manifold reduction is to reduce an ordinary or partial diffe-rential equation system onto a manifold of finite dimensionality without losing any ofthe important information about the local dynamics of the system.Letx˙ = f(x) (2.53)120.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.876.176.276.376.476.576.676.7γA(0, 3) (5, 1)(3, 2)(2,2)Figure 2.4: Marginal stability curves for four different (m,n) modes using A and γ asparameters. Other parameter values are given in equation (2.52) and R = 0.96747.be a system of ordinary differential equations for some nonlinear, continuously differen-tiable vector field f . We may represent partial differential equations, such as reaction-diffusion systems like equation (2.36) as an infinite system of ordinary differential equa-tions over all eigenfunctions, such as those described in equations (2.45) and (2.46).Assume the origin is an equilibrium, f(0) = 0, then local differential equation theorystates that the dynamics close to the origin will often be determined by the linearizationof the systemx˙ = A0x, (2.54)where A0 is the Jacobian matrix for the function f evaluated at the equilibrium X = 0.If A0 has k eigenvalues (counting multiplicities) with negative real part, l eigenva-lues with positive real part, and m eigenvalues with zero real part, then there willbe a k-dimensional stable subspace Es, a l-dimensional unstable subspace Eu, and anm-dimensional centre subspace Ec where solutions of the linearized system (2.54) areexponentially attracted, repulsed or marginal with respect to the origin.The Stable Manifold and Centre Manifold Theorems (e.g. [23]) state that there is ak-dimensional invariant local stable manifold W s, where any solution xs(t) ∈ W s to thenonlinear system (2.53) near the origin has the limitlimt→∞xs(t) = 0, (2.55)a l-dimensional invariant local unstable manifold W u, where any solution xu ∈ W u near13the origin has the limitlimt→−∞xu(t) = 0, (2.56)and an m-dimensional invariant local centre manifold W c. Each manifold is tangent toits associated subspace and is locally invariant under the flow of (2.53), which meansthat solutions contained in each manifold stay within the manifold with time evolution.If the system (2.53) is finite dimensional with dimension r, then all three manifolds arealso finite dimensional and m = r − k − l.For centre manifold reduction to be most useful in local bifurcation analysis, thecentre manifold should have a finite dimension m and the unstable manifold should beabsent (l = 0). In this case, we know that any solution starting near the origin expo-nentially decays to the local centre manifold and any interesting bifurcation dynamicsoccur within W c.We can then use the differential equation (2.36) to find an equation for the centremanifold, which we can approximate using a Taylor expansion (e.g. Section 3.2 of [23])as a function of a finite number of manifold coordinates. Expressing the differentialequation restricted to the centre manifold will result in an m-dimensional system ofordinary equations. In a bifurcating system the dynamics near the origin will changequalitatively as we change one or more system parameters. This system can then betransformed using locally defined phase space coordinate and parameter homeomor-phisms into a locally topologically equivalent normal form system (e.g. Chapter 2 of[32]). Two systems are locally topologically equivalent if there is a homeomorphism (acontinuously differentiable coordinate change) between the two systems that locally pre-serves their dynamics. Normal forms are the fundamental systems on which bifurcationanalysis is performed. Sometimes prior knowledge of the system such as symmetries inthe domain and the system of differential equations (2.36) help us finding the normalform.As an example, we use the marginal stability curves on the spherical cap describedin section 2.4.1. Restricting parameter changes to within the neighbourhood of a singlemarginal stability curve we can apply a simple bifurcation analysis. We use projectionmethods on the nonlinear system to build an invariant centre manifold for the reaction-diffusion equations and obtain a normal form. This section briefly outlines this processfor m 6= 0.Setting our parameters so that the largest eigenvalue of A0 is critical, achieved withσ+mn = 0 for some fixed (m,n), m 6= 0, we find an associated eigenfunctionU(0) =(u(0)mnv(0)mn)Φmn (2.57)for the critical (m,n) values. At these parameter values, all the other eigenvalues of A0are negative and bounded away below zero. This yields a two-dimensional (if m 6= 0)14centre subspace for the linearization, which we can describe using a complex variable zasEc ={zU(0) + z¯U¯(0)∣∣ z ∈ C} , (2.58)where the overhead bar designates the complex conjugate and U(0) is any fixed choiceof an eigenfunction (2.57). This space denotes the dominant mode that will emerge atthe bifurcation.Because the largest eigenvalue is close to zero in the neighbourhood of the parameterspace, there will be, for |z| sufficiently small, a locally invariant centre manifoldW c ={zU(0) + z¯U¯(0) +O (|z|2)∣∣ z ∈ C, |z| sufficiently small} (2.59)for the nonlinear equation. All nearby solutions of (2.36) will decay exponentially tothis manifold and the general long term local behaviour of the system can essentially bedescribed by the system restricted to W c. When solving for the normal form, it mustbe equivariant under the transformationsz → eimϕ0z, z → z¯, (2.60)associated to the rotation by an angle ϕ0 and reflection symmetries (2.38) respectively.This leads us to a normal form for a pitchfork bifurcationz˙ = σz + C|z|2z, z ∈ C (2.61)where σ is the critical eigenvalue, which is real and close to zero, and acts as thebifurcation parameter and a cubic coefficient C, which is real because of the modelsymmetries (2.60).Because of the equivariance under rotations, we may only consider the real part ofthe normal formx˙ = σx+ Cx3, x ∈ R (2.62)where σ and C are the same as in (2.61). The nontrivial equilibrium x =√−σ/C, withσ and C of opposite sign, for this normal form equation gives us the magnitude of thebifurcating mode so under appropriate scaling we can obtain a good approximation ofa patterned solution of the partial differential equation (2.36) asU =√−σC(U(0) + U¯(0)). (2.63)2.6 Delayed Bifurcation AnalysisIn Chapter 4, we make the curvature index γ for the spherical cap domain slowly decre-ase in time and choose parameters such that parameter values cross a marginal stability15curve from the patternless region. This added feature will give rise to delayed bifurcati-ons, which show a time delay in the emergence of the patterned solution. This sectiongoes through the process of estimating the delay time to leading order in a simple delayedpitchfork bifurcation in a similar fashion as explained by Lobry [34].The delayed pitchfork bifurcation occurs when we replace the constant bifurcationparameter in the pitchfork bifurcation normal form. As a simple example we can modifythe pitchfork normal form (2.62) with C = −1 tox˙ = σ(t)x− x3, x(0) = x0,σ(t) = σ0 + εt, σ(0) = σ0,(2.64)where 0 < ε 1 a small parameter. This is a nonautonomous system, though it canbe made autonomous if we include an separate equation for σ:x˙ = σx− x3, x(0) = x0,σ˙ = ε, σ(0) = σ0.(2.65)In the situation of the delayed pitchfork bifurcation, we put σ0 < 0 and make theparameter σ increase through the value σ = 0 (where the bifurcation occurs in theconstant parameter case ε = 0). We start our solution with x0 6= 0, as otherwise this isan equilibrium for all parameter values and no bifurcation would be observed.In systems with slowly evolving quantities, such as σ in this situation we define aslow timescale τ = εt, which changes the system (2.65) toεx′ = σˆx− x3, x(0) = x0,σˆ′ = 1, σˆ(0) = σ0,(2.66)where the prime notation denotes a derivative with respect to τ , thus σˆ(τ) = σ0 + τ .Initially, because σ(0) < 0 the solution is quickly attracted to the origin. Then as theparameter changes sign at τ = −σ0 the solution starts to grow again, but very slowly,as we are still in proximity of an unstable equilibrium, with respect to x at the origin,now unstable. Once the solution grows enough, it will quickly grow and then approachthe stable equilibrium at x =√σˆ(τ). The delay in the bifurcation is measured by thetime elapsed between the change of sign of σˆ and when the solution is close enough tothe constant σ equilibrium. In the limit where ε is very small, the quick growth intervalis short enough that a suitable measure of the delay is when the solution growth is at amaximum.In order to quantify this delay, we assume that the solution x is very close to zero,as it is the majority of the time until that jump. Because of this the cubic term in thedifferential equation is negligible and we may solve the linear equation exactly:x(t) = x(0)eσ0τ+τ22ε ,16which is a symmetric function centered at τ = −σ0: decreasing for σˆ < 0 and increasingfor σˆ > 0, as expected. Due to this symmetry the time the solution takes to reach theconstant σ curve is the same time the solution took to get from a value O(1) to σ = 0and that is approximately t0 =−σ0ε. If we choose instead to start exactly at σ0 = 0 orwith a very small initial condition on x we can extrapolate the delay by using the valueof the solution at σˆ = 0 and finding where it is equal to some value of order 1.σ−σ0xx1x0σ0Figure 2.5: Sketch of a delayed pitchfork bifurcation solution (solid line) of an equationof the type (2.64) along with the constant parameter bifurcation diagram (dashed). Westart at time t = 0 at the constant parameter bifurcation point σ = 0 and with a verysmall initial condition x0, relative to the value x1, when the solution approaches closelythe constant parameter curve. In worked examples we found that for ε = 0.001 we getthat σ0 ≈ −0.1.Figure 2.5 shows a sketch of a typical delayed pitchfork solution. One can see thatit is required for the solution to be close to the origin for the delay to be observable. Ifwe choose to start at the constant domain bifurcation point, that is σ0 = 0 and measurethe delay, we need to refer to the solution close to the origin, which is now:x(t) = x(0)eτ22ε .One way to know we are in the transition from one equilibrium to the other is whenthe solution is no longer close to the origin, so at some value x1. Normally we wouldwant this value to be either the expected value of the solution at the stable equilibriumor the value when it is increasing at a maximal rate, but the solution is increasing soquickly under these circumstances that any value in the same order of magnitude as theexpected value will work, to leading order in ε. In this example the choice of x1 = 1 isadequate.We need to solve the following for the delay time τd:x1 = x(0)eτ2d2ε17and the result isτd =√2ε[ln(x1)− ln(x(0))],which is order√ε in scaled time. Using similar methods we find that if σ0 < 0 we getthe following delay time resultτd = σ0 +√σ20 + 2ε(ln(x1)− ln(x0)) (2.67)and that if x0 is large enough, say equal to x1, then our time delay is t = 2σ0/ε,that is twice the time from the initial σ0 until σ = 0. In the end the function willapproximately rejoin the non-zero stable branch at a σ value of −σ0. This result is agood approximation as long as x0 remains in the same order of magnitude as x1.Finally, the estimation of the delay time only depends on the linear term of thedifferential equation, so the results will still be applicable for a slowly varying cubiccoefficient, such as with the following equation:x˙ = σ(t)x− (1 + aσ(t))x3, σ(t) = σ0 + εt.However if one includes other slowly varying terms, the zero solution x = 0 might nolonger be an equilibrium for all values of σ and in this case a computation of the slowlyvarying equilibrium will be necessary [24].2.7 Closest Point MethodTo validate our results performed on ordinary differential equations in Chapters 3 and 4,we performed numerical simulations of the different partial differential equations usingeither the finite element method or the closest point method. In this section we give ashort description of the latter method.Let Ω be an (n−1)-dimensional continuously differentiable surface in Rn. Let∇Ω and∇2Ω be the surface gradient and Laplace–Beltrami operator associated to the surface Ω.We may then use these operators to formulate reaction-diffusion differential equationson the surface, for example in a system of two equations∂u1∂t= DX∇2Ωu1 + f(u1, u2),∂u2∂t= DY∇2Ωu2 + f(u1, u2),(2.68)where x ∈ Rn. The aim of the closest point method is to solve such equations numericallyby associating points on a band surrounding the surface to their closest points on thesurface.The closest point function cp(x) map each point in Rn to a point on Ω such that theEuclidian distance between the two points is minimizedcp(x) = argminy∈Ω‖x− y‖. (2.69)18In a neighbourhood near enough to the smooth surface, this point is unique. Thisfunction allows for a surface function u(x) on Ω to be extended to the surroundingspace using the extension operator Ev(x) = Eu(x) = u(cp(x)). (2.70)The extended function v is constant in a normal direction to the surface, due to the factthat points on a line perpendicular to a surface all share the same closest point, at leastin a sufficiently small neighbourhood of a smooth surface.Because extended functions are constant in a normal direction to Ω there is an equiva-lence between the surface differential operators and the Cartesian differential operatorsin the surrounding space of the extended function evaluated on the surface. This is notedas three principles: the Gradient Principle, the Divergence Principle and the LaplacianPrinciple∇(Eu)(y) = ∇Ωu(y), y ∈ Ω,∇ · (Eg)(y) = ∇Ω · g(y), y ∈ Ω, (2.71)∇2(Eu)(y) = ∇2Ωu(y), y ∈ Ω,respectively, for any smooth scalar function u : Ω→ R and smooth vector field g : Ω→Rn. These principles have been developed rigourously in previous work [53], [39], [59].The main corollary of the principles (2.71) is that we can evaluate surface derivativesby evaluating the corresponding operators in Rn, provided we keep the solution constantin the normal direction to the surface Ω. In other words we reuse algorithms and methodsfor solving the non-surface PDE∂u1∂t= DX∇2u1 + f(u1, u2),∂u2∂t= DY∇2u2 + f(u1, u2),(2.72)and then restrict our solutions to points lying in Ω.2.7.1 Spatial DiscretizationThe closest point extension onto the surrounding space to the surface Ω allows us todiscretize the extended system (2.72), then solve that discretized system using a finitedifference method on a band B(Ω) around Ω with regular grid spacing h.This discretized version of the PDE has the following form∂ts = As, (2.73)where s is the discretization of the solution on a regular cartesian grid and A is thediscretization of the differential operator. We describe this in more detail in 2D. We use19a 2D five-point (seven-point in 3D) stencil for the discrete LaplacianLvi,j =1h2(vi+1,j + vi−1,j + vi,j+1 + vi,j−1 − 4vi,j), (2.74)and two-point centered finite differences for any first derivatives (such as in advectionor divergence operators),D1vi,j =12h(vi+1,j − vi−1,j),D2vi,j = 12h(vi,j+1 − vi,j−1), (2.75)where the indices i, j denote the ith column and jth row in the grid space.In order to force the system to be constant in a normal direction, when solvingusing an explicit time stepping method we need to apply an extension from the closestpoints on the surface to the points on the band. We do this numerically by applying aninterpolation of the solution on the band to the surface then reassigning the values ontothe band. This operator can be written as a matrix E [36]. In this document we useBarycentric Lagrange interpolation [6], which uses polynomials at points on the grid ina close neighbourhood of the interpolated point.The spatial discretization and extension steps may be recombined together into asingle operator that uses a penalty term [59]Apen = EA− γpen(I− E), (2.76)where γpen is a penalty factor. In the case of the finite difference discrete Laplacianwe follow [59] and take γpen = 2 · dim/h2, the value of the diagonal elements of theoperator with dim being the dimension of the embedding space (for example dim = 3for a surface embedded in R3).2.7.2 Dirichlet Boundary ConditionsFor surfaces with a boundary the closest point method automatically applies homoge-neous Neumann boundary conditions. As previously done by Macdonald et al. [35], inorder to use Dirichlet boundary conditions we can modify the closest point functionc¯p(x) = cp (x + 2(cp(x)− x)) = cp (2cp(x)− x) . (2.77)For points whose closest points lie in the interior of the surface this new function takesthe same value as cp(x). For closest points xg on the boundary the function creates aghost point on the surface to which we can apply a boundary conditionu(xg) = 2u(cp(xg))− u(c¯p(xg)), (2.78)where u(cp(xg)) is the value of the solution on the boundary. In the case of homogeneousDirichlet boundary conditions we have u(cp(xg)) = 0 and thusu(xg) = −u(c¯p(xg)). (2.79)20To apply this to the numerical simulations we include it in the interpolation operatorEDir(boundary) = −E(boundary). (2.80)In practice, this involves changing the sign of certain rows of E2.7.3 Explicit Time-Stepping Scheme ExampleIn order to perform numerical simulations of partial differential equations using theclosest point method, we can use a finite difference scheme on a regularly spaced griddiscretization of size h of the embedding space in each time step ti. We can then findthe closest point of each grid point to Ω to get our approximate solution on Ω.However, the Laplacian discretization will no longer satisfy the three principles(2.71), which means that after a few time steps solutions would probably no longerbe constant on normals of Ω. To remedy this we perform an interpolation onto Ω andextension back into the embedding space grid using the closest point function at everytime step in order to force the solution to remain constant in the normal directions.Because we perform an interpolation and an extension at every time step, the solutionat all grid points will be determined from the interpolated solution, which is determinedwith grid points sufficiently close to Ω. As a result we may perform the time stepson a narrow band around Ω in the embedding space to significantly reduce computingtime and avoid grid points with multiple closest points. The band should include allgrid points used in interpolation plus the numbers of neighbours needed to compute thediscretized Laplacian. Figure 2.6 shows an explicit forward Euler time step using theclosest point method.This concludes the introduction of basic concepts and tools which we will expandupon in this document in order to study different pattern formation problems on curvedsurfaces. More specifically, in Chapter 3 we will extend the Brusselator reaction-diffusionsystem bifurcation diagram shown in Section 2.4 by describing the bifurcations occuringnear the intersection to two marginal stability curves. In Chapter 4 we will extend thissection by studying delayed bifurcations resulting from a slowly decreasing curvatureparameter σ. Then, in Chapter 5, we will adapt the closest point method in order tocompute numerical simulations on evolving curved domains, including our flatteningspherical cap model.21(a) (b) (c)(d) (e) (f)Figure 2.6: Sketch of a forward Euler closest point method time step. In panel a, wehave the solution on the curve Ω (in red). In panels b and c, the solution is extendedonto the grid points using the closest point method. In panel d, the explicit time step isapplied to the point grid. Points on the edge of the band (squares) will be defective (ingreen) because they are missing neighbours for a valid Laplacian computation and innerpoints (circles) will have a solution value for the next time step (in blue). In panels eand f, barycentric Lagrangian interpolation is applied on points along Ω and we continueto panel a for the next time step. Note that only inner band points (circles) should beused for interpolation and that after extension (panel b) the outer points (squares) willbecome valid again.22Chapter 3Bifurcation with Two CompetingModes in a Constant Spherical CapIn this chapter we study the bifurcations for the Brusselator reaction-diffusion system ina constant spherical cap domain near the intersection of two marginal stability curves.In this situation, there will be two distinct possible marginally stable linear modes, andnonlinear “mode competition” that involves multiple bifurcations. We use centre ma-nifold analysis to reduce the Brusselator system of partial differential equations (2.36),(2.37) to a small system of ordinary differential equations, which we call a normal form.Two specific cases are considered that depend on the symmetries (2.38) of the margi-nally stable modes. Analysis of the normal forms is then performed to determine theemerging patterns in the Brusselator.We use a similar setup as in Sections 2.2 and 2.4, using the same spherical capdomain with base radius R and (constant) curvature index γ. We express the Brusselatorequations (2.8) as an evolution equation as in (2.36) with the same operators K0, B0,C0. The same marginal stability curves as described in Section 2.4.1 apply here. Figure3.1 highlights two such intersections in a marginal stability curve diagram, which definecritical parameter values, for parameter values (2.52) and R = 0.96747.In observations of conifer embryos [58] a five-cotyledon pattern is commonly obser-ved, although the number of cotyledons can range from two to twelve. In this chapter,we turn our focus to stable patterns emerging near the (m,n) = (5, 1) marginal stabilitycurves. It has been previously determined that there is a simple bifurcation of a pure(5, 1) mode if the curve is traversed in parameter space from the region where the pat-ternless solution is stable (above each curve in Figure 3.1), without traversing any othercurve. We extend those results to the intersections of the (5, 1) marginal stability curvewith the (3, 2) and (0, 3) curves, labeled as Case 1 and Case 2 respectively in Figure 3.1.These two curves are the closest to the (5, 1) curve for our choice of parameter values(2.52) and are the only possible transitions from the patternless solution region into aregion where there are two possible competing patterned modes, including (5, 1).We begin, however, with a general approach. We wish to study the bifurcation ofsolutions to (2.36) from the patternless steady state U = (0, 0)ᵀ to patterned steadystates when the parameters are close to the critical values. For two intersecting marginalstability curves there will be a pair of modes (m1, n2) and (m2, n2) and two linearlyindependent normal modes U(0)1 and U(0)2 associated to the zero eigenvalues of A0.Both have the same structure as in the single mode case (2.57), but for different (m,n)23Figure 3.1: Marginal stability curves as in Figure 2.4 with highlighted intersections. Case1 corresponds to the case where both modes have m 6= 0 and Case 2 corresponds to thecase where, for one mode, m = 0. We used parameter values (2.52) and R = 0.96747.valuesU(0)j =(u(0)mj ,njv(0)mj ,nj)eimjϕP|mj |λmj,nj (γ)(cos θ), j = 1, 2, (3.1)where u(0)mj ,nj and v(0)mj ,nj solve the eigenvector problem (2.49) and are normalized asu(0)mj ,nj = k2/N(0)j , v(0)mj ,nj= k5,j/N(0)j , k5,j = DXµmj ,nj − k1, (3.2)N(0)j =√Mj(k22 + k25,j), Mj =∫ θmax0[P|mj |λmj,nj (γ)(cos θ)]2sin θdθ, (3.3)where µmj ,nj is the appropriate Laplace–Beltrami operator eigenvalue. We point outthat u(0)mj ,nj , v(0)mj ,nj and P|mj |λmj,nj (γ)(cos θ) in (3.1) are all real and independent of the signof mj. For each curve j = 1, 2 the real eigenspace associated to U(0)j is one-dimensionalif mj = 0 or two-dimensional if mj 6= 0. The centre subspace Ec at the intersectionbetween the curves is the direct sum of the eigenspaces corresponding to the doublezero eigenvalue, and is spanned by U(0)1 , U(0)2 and their complex conjugates. The stablesubspace Es, complementary to Ec, is the sum of the eigenspaces corresponding to allthe other possible eigenvalues and has infinite dimension.By centre manifold theory on reaction-diffusion systems [11] there exists a locallyinvariant, exponentially attracting centre manifold W c for the system (2.36) at the24critical parameter values. This manifold has the same dimension as the centre subspaceand is tangent to Ec at the patternless equilibrium U = 0.In the neighbourhood of the marginal stability curve intersection in parameter spacethere are two real eigenvalues σ1, σ2 of (2.49) that will remain close to 0 and dependsmoothly on the sytem parameters and σ1 = σ2 = 0 occurs at the critical parametervalues. All other eigenvalues are real, negative and bounded away from zero. The centremanifold W c also depends smoothly on the parameters and will retain its local invarianceand exponential attraction. This signifies that solutions to (2.36) will quickly decay tothis centre manifold, where any local bifurcation may occur. To find patterns emergingfrom the patternless state we can thus reduce the dyamics of the infinite dimensionalspace of functions in Ω to a manifold of a finite dimension. The reduction process isdescribed in Appendix A.We study two cases here that feature different dynamics, as shown in Figure 3.1. Thefirst case has m1 6= 0 and m2 6= 0, which we call Case 1, has a four-dimensional centresubspace and manifold. Case 2 has m1 6= 0 and m2 = 0 and has a three-dimensionalcentre subspace and manifold.3.1 Case 1: m1 6= 0 and m2 6= 0We express the four-dimensional centre subspace using complex notationEc = {z1U(0)1 + z2U(0)2 + z¯1U¯(0)1 + z¯2U¯(0)2 |z1, z2 ∈ C}. (3.4)We can then construct the four-dimensional centre manifold, tangent to the centresubspace, in terms of z1, z2 asW c = {z1U(0)1 + z2U(0)2 + z¯1U¯(0)1 + z¯2U¯(0)2 +O(‖(z1, z2)‖2) |z1, z2 ∈ C}, (3.5)with ‖(z1, z2)‖ sufficiently small. Some of the higher order terms O (‖(z1, z2)‖2) in theTaylor series of the centre manifold are needed explicitly for computing normal formcoefficients and their derivations are described in Appendix A.For parameter values in a neighbourhood of the critical values, the dynamics of(2.36) reduce to a four-dimensional system of ordinary differential equations on thecentre manifold. The reduced system has the following general formz˙1 = f1(z1, z2, z¯1, z¯2), (3.6)z˙2 = f2(z1, z2, z¯1, z¯2), (3.7)where the complex f1 and f2 also depend on the many parameters of the Brusselatorsystem. The symmetries (2.60) of the Brusselator system act, according to (3.2) and(3.4), on the centre subspace Ec aszj → eimjφzj, zj → z¯j, (3.8)25for j = 1, 2, for any angle rotation φ. This symmetry carries through the centre manifoldreduction and implies that the manifold is invariant under (3.8), which in turn impliesthat the reduced system (3.6) must be equivariant under the transformations (3.8), thatisfj(eim1φz1, eim2φz2, e−im1φz¯1, e−im2φz¯2) = eimjφfj(z1, z2, z¯1, z¯2), (3.9)for j = 1, 2 and any real φ.The equivariance equation (3.9) implies that the cubic Taylor polynomial of thereduced system (3.6) has no quadratic terms and the following linear and cubic termsz˙1 = σ1z1 + C11z21 z¯1 + C12z1z2z¯2,z˙2 = σ2z2 + C21z1z2z¯1 + C22z22 z¯2,(3.10)with higher order terms sized at least O (‖z1, z2‖5) and where the coefficients σj andCjk are real numbers depending on the parameter values near the critical values. Infact, the σj are the two eigenvalues of the linearized system (2.39) closest to zero forparameter values near the critical values, as mentioned above. We refer to equation(3.10) as the normal form, mj 6= 0 case. In the case where we use parameter values(2.52) to numerically compute the different coefficients the normal form turns out tobe generic, thus predictions of stability and local bifurcation of steady states obtainedfrom it are valid and accurate to leading order for the full system (2.36) (or equivalently(2.8)) as long as parameter values stay near the critical values, despite neglecting higherorder terms.The equivariance (3.9) of the reduced system (3.6) also implies that the argumentsarg(zj) remain constant in time t, and thus it suffices to consider only the amplitudesrj = |zj| ≥ 0, j = 1, 2. At leading order, these satisfy the amplitude equationsr˙1 = σ1r1 + C11r31 + C12r1r22,r˙2 = σ2r2 + C21r21r2 + C22r32.(3.11)The amplitude equations (3.11) can have equilibria of different types. The origin(r1, r2) = (0, 0) is always an equilibrium, corresponding of course to the equilibrium(z1, z2) = (0, 0) in the normal form (3.10), and to the patternless solution U = 0 ofthe Brusselator (2.36). Thus we often refer to this equilibrium at the origin as thepatternless state. We call equilibria of the type (r1, 0) or (0, r2), with rj > 0, (non-linear) r1-modes or r2-modes, respectively. For the complex normal form (3.10) thesecorrespond to circles of equilibria (z1, 0) or (0, z2), respectively, with fixed amplitudes(r1 and r2) and arbitrary arguments. For the Brusselator (2.36) these correspond to‘pure patterned modes’, steady state solutions U = zjU(0)j + z¯jU¯(0)j + O (|zj|2) of thenonlinear Brusselator system closely resembling a constant multiple zjU(0)j + z¯jU¯(0)j ofone of the critical normal modes or eigenvectors (3.1), j = 1, 2. There are continuousfamilies of these pure patterned modes, corresponding to arbitrary rotations of a specific26patterned mode. Finally we call equilibria of the type (r1, r2), with both r1 > 0 andr2 > 0, (nonlinear) mixed modes. For the complex normal form (3.10) these mixed mo-des correspond to 2-tori of equilibria, and for the Brusselator (2.36) these correspond to‘mixed patterned modes’, steady state solutions of (2.36) that closely ressemble linearcombinations of the two critical eigenfunctions (3.1) of the linearization (2.39). Theamplitude equations (3.11) have no periodic solutions in the region of interest.The system of amplitude equations (3.11) is the normal form for a double pitchforkbifurcation with Z2 ⊕ Z2 symmetry, which has been well studied (e.g. [22], ChapterX). For (σ1, σ2) in a neighbourhood of (0, 0) the bifurcation and stability of equilibrianear the origin for (3.11) and (3.10) (asymptotic stability of equilibria for the amplitudeequations corresponds to the asymptotic orbital stability of the families (circles or 2-tori) of equilibria in the complex normal form) are qualitatively valid and accurate toleading order for the reduced system (3.6), not depending on the neglected remainderterms O (‖(z1, z2)‖5), and thus are also valid for the full Brusselator system (2.36). Theamplitude equations (3.11) also correspond to the amplitude equations of the doubleHopf, or Hopf-Hopf bifurcation normal form (e.g. [32], Section 8.6).The possible bifurcations and stability configurations of equilibria for the amplitudeequations (3.11) depend on relationships between the cubic coefficients Cjk. We givesome details in Appendix A on the calculation of these cubic coefficients. To accuratelydescribe the bifurcation and stability of equilibria near the critical parameter values, itis sufficient to evaluate the cubic coefficients at the critical parameter values themselves.For the choice of parameter values (2.52) andγ = 0.5, A = 76.344, R = 1.03348, (3.12)corresponding to the (m1, n1) = (5, 1) and (m2, n2) = (3, 2) marginal stability curvesintersection, we obtain the following numerical values for the cubic coefficientsC11 = −2.1509, C12 = −3.2905, C22 = −2.5687, C21 = −3.7082. (3.13)Different values of the critical eigenvalues σ1, σ2 can give qualitatively differentphase portraits. To use the terminology of [32], a parametric portrait is a partition ofparameter space into subsets of parameter values that give qualitatively the same phaseportraits. For the values of the coefficients obtained in (3.13), the (σ1, σ2)-parametricportrait corresponds to Case I in Figure 8.25 of [32], Section 8.6, which we reproducehere in Figure 3.2, with the phase portrait corresponding to each region sketched inFigure 3.3. All curves in Figure 3.2 correspond to pitchfork bifurcations, and thesecurves divide the parameter plane into six open regions shown that are all contained inany open neighbourhood of the origin.We describe the relationship between the parametric portrait and the phase portraitsof the amplitude equations (3.11) by taking a counter-clockwise path around the origin inthe parametric portrait, Figure 3.2. In region 1, the only equilibrium is the patternlessstate (r1, r2) = (0, 0) and it is asymptotically stable. When we cross the negative271 23456σ1σ2L2L1Figure 3.2: Parametric portrait for Case 1 m1 6= 0, m2 6= 0 amplitude equations (3.11)at parameter values (2.52). Curves of pitchfork bifurcations subdivide the (σ1, σ2)-plane into six regions. Phase portraits for the amplitude equations for parameter values(σ1, σ2) in each of the six numbered regions are shown in Figure 3.3 below.σ2-axis from region 1 to region 2, we get a pitchfork bifurcation at the origin andan asymptotically stable equilibrium on the r1-axis appears, an r1-mode. If we thencross the positive σ1-axis we find ourselves in region 3 and there is another pitchforkbifurcation at the origin that marks the appearance of an r2-mode. The r2-mode is asaddle, thus is unstable for parameter values in region 3, then we cross into region 4, anda pitchfork bifurcation at the r2-mode creates a mixed mode. The mixed mode is a saddlepoint and both pure modes (the r1-mode and the r2-mode) are asymptotically stable.Thus for parameter values in region 4 we have bistability, two distinct asymptoticallystable equilibria, the two pure modes: a typical trajectory is attracted to only one ofthe pure modes, depending on the initial condition. The mixed mode, a saddle, onlyattracts trajectories that start on the codimension 1 stable manifold of the mixed mode,an unlikely event if trajectories are chosen randomly. As parameter values move fromregion 4 to region 5, the mixed mode vanishes in a pitchfork bifurcation at the r1-modeleaving an unstable r1-mode and an asymptotically stable r2-mode. As parameter valuesare varied between regions 3, 4 and 5, we have hysteresis or ‘mode jumping’. As wecontinue our parameter path from region 5 to region 6, the unstable r1-mode vanishesin a pitchfork bifurcation at the origin while the asymptotically stable r2-mode and the28r1r21r1r22r1r23r1r24r1r25r1r26Figure 3.3: Phase portraits for Case 1 m1 6= 0, m2 6= 0 amplitude equations (3.11).Numbers on the phase portraits correspond to parameters (σ1, σ2) in the numberedopen regions in the parametric portrait in Figure 3.2 above.unstable origin remain. Finally from region 6 we cross the negative σ1-axis back toregion 1 where the r2-mode also vanishes in a pitchfork bifurcation at the origin andonly the asymptotically stable patternless state remains at the origin.3.2 Case 2: m1 6= 0 and m2 = 0In the case when one of the marginally stable modes (3.1) has mj = 0 the action of therotation symmetry (2.38) is trivial. If m1 6= 0 and m2 = 0 the center subspaceEc = {zU(0)1 + z¯U¯(0)1 + xU(0)2 |, z ∈ C, x ∈ R},and the corresponding center manifold W c is three-dimensional. The reduced three-dimensional system or ordinary differential equations has the general formz˙ = f(z, z¯, x),x˙ = g(z, z¯, x),(3.14)where f is complex and g is real. Equivariance under the remaining symmetriesz 7→ eim1ψz, z 7→ z¯. (3.15)29for arbitrary angle ψ restricts (3.14) so that the Case 2 normal form isz˙ = σ1z +B11xz + C11z2z¯ + C12zx2x˙ = σ2x+B21zz¯ +B22x2 + C11zz¯x+ C22x3,(3.16)with remainders O (‖(z, x)‖4).By equivariance all the coefficients in the normal form (3.16) are real and the argu-ment arg(z) remains constant in t, and thus may only consider the amplitude r = |z| ≥ 0in (3.16), which satisfiesr˙ = σ1r +B11xr + C11r3 + C12rx2x˙ = σ2x+B21r2 +B22x2 + C21r2x+ C22x3.(3.17)This normal form (3.17) coincides with the amplitude equations for a transcritical-Hopfbifurcation [33], but without periodic solutions. As in Case 1 the system may haveequilibria of different types. The equilibrium at the origin (r, x) = (0, 0) correspondsto the equilibrium (z, x) = (0, 0) in (3.16), and to the patternless steady state U = 0of the Brusselator (2.8). There are also x-mode equilibria for the amplitude equations(0, x) with x 6= 0, which corresponds to equilibrium of (3.16) with z = 0, x 6= 0. Theseequilibria correspond to pure single patterned modes for the Brusselator (2.8) of the formU = xU(0)2 +O (x2) that are rotationally invariant (independent of ϕ). There may alsobe mixed modes corresponding to circles of equilibria (reiφ, x) of the normal form (3.16),and to mixed patterned modes of the form U = zU(0)1 +z¯U¯(0)1 +xU(0)2 +O (‖z, x‖). Thesemixed mode patterned solutions have m1-fold rotational symmetry about the sphericalcap centre.Derivations of the quadratic and cubic coefficients Bjk, Cjk in the normal forms(3.16) or (3.17) are explained in Appendix A. Using parameter values (2.52) andγ = 0.5, A = 76.531, R = 0.96747, (3.18)with the critical values at the intersection between the (m1, n1) = (5, 1) and (m2, n2) =(0, 3) mode marginal stability curves we obtain the following values for the quadraticcoefficientsB11 = 0.0108× 10−2, B21 = 1.4821× 10−2, B22 = −0.9413× 10−4. (3.19)Since the values of these coefficients are numerically small, truncating the normal form(3.17) after the quadratic term will yield the behaviour of the solutions in a region toosmall around the patternless solution to be useful. We thus also compute the cubiccoefficients and after some more work we obtainedC11 = −3.5712, C12 = −1.5367, C21 = −3.0553, C22 = −2.0732. (3.20)These values are all at least two orders of magnitude larger than the quadratic coeffi-cients, supporting truncation of the normal form at cubic order.30Using the numerical values (3.19), (3.20) for the normal form coefficient we find the(σ1, σ2)-parametric portrait for the amplitude equations (3.17), shown in Figures 3.4 to3.6. The portrait contains fourteen different open regions of the parameter plane forwhich solutions have a distinct stability configuration, the phase diagram of those areshown in Figures 3.7 to 3.9. Details on the calculations behind the regions’ boundariesare in Appendix B.Figure 3.4: Bifurcation diagram for the Case 2 normal form (3.17) for a large σ1, σ2range. The diagram is very similar to the double pitchfork diagram shown in Figure 3.2The parametric portrait has different features depending on the scales used for theparameters σ1, σ2. At the large scale |σ1| < 1, |σ2| < 1 in Figure 3.4 the portrait lookssimilar to the double pitchfork normal form portrait shown in Figure 3.2, where allquadratic coefficients Bjk = 0.As we zoom in on the origin, more structure is revealed. Figures 3.5a to 3.5c showthe parametric portrait for successively smaller scales. For example what appeared tobe a semi line separating regions 4 and 5 in the large scale portrait in Figure 3.4 isactually two nearby curves (dashed red) of pitchfork bifurcations from equilibria on thex-axis, with a region between them, where one pitchfork bifurcation happened beforethe other. There is also a curve of fold bifurcations (dotted blue) between regions 3 and4. Similarly, the apparent pitchfork bifurcation lines in Figure 3.4 on the σ1 and σ2-axesare actually perturbed pitchfork bifurcations consisting of a line of fold bifurcations(cyan dash-dotted line) and a line of transcritical bifurcations (the σ1 axis, solid black)31and a line of pitchfork bifurcations (the σ2 axis, solid black) with regions 8, 9, 10, and11 between the fold and transcritical lines, as shown in Figure 3.5b. Getting even closerto the origin reveals a short curve of fold bifurcations (dotted blue) separating regions6 and 12 and regions 8 and 13. We also observe that the pitchfork curves (dashed red)encroach to the left side of the σ2-axis which creates region 14. Regions 12, 13 and 14are pictured in Figure 3.5c.Extremely near the origin, at the very small scale |σ1| < 10−11, |σ2| < 10−9 in Figure3.6, the parametric portrait consists of three curves intersecting at the origin (the σ1-axis, the σ2-axis and a curve of tangent of slope B22/B11 at the origin), with six openregions labeled 7, 9, 10, 12, 13 and 14 closest to the origin. This is the portrait onewould obtain by truncating the normal form (3.17) to quadratic order. The region ofvalidity of this truncation is very small, as expected.One may think of the small quadratic terms in the normal form equations (3.17) as asmall symmetry-breaking perturbation of the double pitchfork bifurcation with Z2⊕Z2symmetry (3.11). This is a direct effect of the absent symmetry of x. In the perturbedequations (3.17) only one of the reflection (Z2) symmetries (r, x)→ (−r, x) remains andthe other (r, x) → (r,−x) is broken, so the x-axis r = 0 remains invariant while ther-axis x = 0 is no longer invariant. The perturbed pitchfork bifurcation is split intofold and transcritical bifurcations (Figures 3.10a and 3.5b) or simply fold bifurcations(Figures 3.10b and 3.5a).To translate predictions from the normal form equation (3.17) to predictions forthe Brusselator (2.36), we study the effects of changing one of the model parametersγ, R or A on the critical eigenvalues and normal form linear coefficients σ1 and σ2.We use equation (2.51) to find the two roots σj closest to zero for model parametervalues, with other parameter values fixed by (2.52). Figures 3.11, 3.12 and 3.13 illustraterelationships between the model parameters γ, R and A and the normal form coefficientsσ1 and σ2. We see that qualitatively, increasing γ or R has roughly the same effect ofclockwise rotation around the origin in the (σ1, σ2)-plane, from parameter region 5 inFigure 3.5c, across region 7 and into region 3. Increasing A moves the point (σ1, σ2) fromregion 7 near the boundary with region 4, towards the lower left passing close to theorigin, into region 1. We note stable mixed mode solutions are present in regions 4 and 7(with phase portraits described in Figures 3.7d and 3.8a). These solutions are mixturesof (5, 1) and (0, 3) pure modes. Region 5 features stable pure (0, 3) mode solutions andregion 3 features stable pure (5, 1) mode solutions. By varying model parameter valuesnear the mj = 0 case intersection point, stable mixed mode appear for a range of values.3.3 Numerical SimulationsIn this section, to support our analytical results for the normal forms, we show somenumerical simulation results. While the analytical results are valid in some neighbour-hood in parameter space and in some neighbourhood of the patternless steady state,32it is not clear how large these neighbourhoods are. We explore parameter space nearcritical values using finite elements simulations of the Brusselator system on sphericalcap domains using an implicit-explicit time-stepping scheme. The finite element methodis a well established numerical simulation tool that can discretize a surface into smallflat elements and approximate solutions using simpler polynomials on each element, forfurther reference, see [3]. The simulations in this section are based on similar simulationsperformed in [43] and [28].3.3.1 Case 1: m1 6= 0 and m2 6= 0For the choice of model parameter values (2.52), there is an intersection of the (m1, n1) =(5, 1) and (m2, n2) = (3, 2) marginal stability curves. In the corresponding (σ1, σ2)-parametric portrait Figure 3.2, in region 4 of the parameter plane we have bistability,parameter values where two distinct pure patterned modes are both asymptoticallystable. For these parameter values, there is also a mixed mode steady state (phaseportrait 4 of Figure 3.3), but it is unstable and therefore would not normally be directlyobserved in simulations.To illustrate the bistability we choose model parameter values (2.52) andγ = 0.5, A = 76.344, R = 1.031, (3.21)so that the corresponding (σ1, σ2) values lie in the region 4 of Figure 3.2, where there isbistability. We can observe evidence of bistability in the Brusselator system by selectingan initial condition near one of the predicted pure patterned modes and then near theother. Since we are so close to the multiple bifurcation point, for initial conditions we uselinear pure mode eigenfunctions with coefficients obtained from the pure mode equilibrain the nonlinear amplitude equations (3.11), plus some small randomly generated noise.Figures 3.14 and 3.15 show the results of two simulations with the same parametervalues, but with two different initial conditions, after evolving to time t = 5000 withstep size 0.1, where we found numerical evidence for exponential convergence to a stablepure mode steady state. The former shows a stable pure (5, 1) mode and the latter astable pure (3, 2) mode at the same model parameter values.We found that it is important for the spatial mesh on the spherical cap to be quitefine, as the region of interest in parameter space is quite small and the stability ofpatterns quite weak near the multiple bifurcation point. As we refined the mesh, weobserved a roughly order h convergence between the numerical steady state and theanalytical predictions.3.3.2 Case 2: m1 6= 0 and m2 = 0For model parameter values (2.52), we have an intersection of the (m1, n1) = (5, 1) and(m2, n2) = (0, 3) marginal stability curves. In contrast with the previous case, there isan orbitally asymptotically stable mixed mode in parameter regions near the (5, 1) and33(0, 3) intersection. In Figures 3.4 to 3.6, for (σ1, σ2) in parameter regions 4 and 7 thereis a stable mixed mode equilibrium in the phase portrait.Figure 3.16 shows the results of a simulation of the Brusselator system with modelparameter values (2.52) andγ = 0.5, A = 76.513885, R = 0.96608, (3.22)so that the corresponding (σ1, σ2) values lie in region 7 of Figures 3.4 to 3.6, whichgives a stable mixed mode equilibrium. The initial conditions for the simulations werea linear combination of the (5, 1) and (0, 3) eigenfunctions, with coefficients obtainedfrom the stable mixed mode equilibrium in the nonlinear amplitude equation (3.17)together with small random noise terms. The simulations were run on successively finermeshes for various times, with 5000 time units showing numerical evidence of exponentialconvergence in time to a stable mixed mode steady state.As parameter values are changed, the stable mixed mode countinuously changes fromone pure mode to the other. However, in simulations it is very difficult to distinguishbetween a stable pure mode and a stable mixed mode that is nearly pure, and so wecould not use simulations to numerically test with much accuracy the boundaries in theparameter space that divide stable pure mode phase portraits from stable mixed modephase portraits. Moreover, because the (r, x)→ (r,−x) equivariance is (slightly) brokenin Case 2, there are no pure r-mode equilibria corresponding to pure (5, 1) mode steadystates: generically the x-components are nonzero, even if they are small.What we have shown here is that we can extend the possible stable patterned statesby including the intersections of two marginal stability curves at the boundary of thepatternless solution parameter region. In the two specific cases studied here we haveobserved the mode competition dynamics of bistability in Case 1 and stable coexistencein the Case 2. In the latter case we observed the new (5, 1)− (0, 3) stable mixed mode,which along with the pure (5, 1) mode, could correspond to observed five-cotyledonpatterns [10]. We also developed the codimension 2 bifurcation diagrams close to theintersection for both cases. In Case 2 we observed 14 distinct stability regions correspon-ding to a perturbed pitchfork bifurcation. We finally performed finite element methodsimulations that show the different possible patterns.34Figure 3.5: Bifurcation diagram for the Case 2 normal form (3.17) for a decreasing σ1,σ2 range. as we restrict the parameter range closer to the origin, additional structureof the bifurcations are revealed, totaling 14 possible equilibrium distributions portrayedin the 14 parametric zones labeled in the figures.35Figure 3.6: Bifurcation diagram for the Case 2 normal form (3.17) for a very small σ1,σ2 range. In this range the quadratic terms of the normal form (3.17) dominate and wesee a pitchfork-transcritical bifurcation diagram.36rx 1O rx 2OMrx 3OPPMrx 4MMMOPPrx 5POPM rx 6POPFigure 3.7: Phase portraits of the normal form equations (3.17) using coefficient values(3.19), (3.20) corresponding to parameter values in the regions described in Figures 3.4to 3.6.37rx 7OPPMMrx 8OPPrx 9OPPMrx 10OPPMMrx 11OPPMrx 12MMPOPFigure 3.8: Phase portraits of the normal form equations (3.17) using coefficient values(3.19), (3.20) corresponding to parameter values in the regions described in Figures 3.4to 3.6 (continued).38rx 13MMOPPrx 14POPMFigure 3.9: Phase portraits of the normal form equations (3.17) using coefficient values(3.19), (3.20) corresponding to parameter values in the regions described in Figures 3.4to 3.6 (end).39σσ σa bFigure 3.10: A bifurcation diagram of a Z2-symmetric pitchfork bifurcation (top) andtwo possible results of a small symmetry-breaking perturbation. In Case a, the pitchforkbifurcation is perturbed into a fold bifurcation and a transcritical bifurcation. In Caseb, there is only a fold bifurcation.40Figure 3.11: Effects of changing the curvature γ on the critical eigenvalues and amplitudeequation parameters σ1, σ2. In this figure γ increases from 0.4 to 0.6, with fixed R =0.96608 and A = 76.508.Figure 3.12: Effects of changing the cap radius R on the critical eigenvalues and ampli-tude equation parameters σ1, σ2. In this figure R increases from 0.962 to 0.970, withfixed γ = 0.5 and A = 76.508.41Figure 3.13: Effects of changing the parameter A on the critical eigenvalues and ampli-tude equation parameters σ1, σ2. In this figure A increases from 76.46 to 76.54, withfixed γ = 0.5 and R = 0.96608.42Figure 3.14: A finite element simulation showing a stable pure (5, 1) mode steady stateof the Brusselator in a spherical cap, viewed from above. The parameters used areγ = 0.5, A = 76.344, R = 1.031, placing (σ1, σ2) in region 4 of Figure 3.2. The initialcondition is obtained from the stable r1-mode equilibrium for the Case 1 amplitudeequations (3.11), with small added random noise. The mesh used had 25095 nodes and50527 elements distributed across the spherical cap surface, and the simulation was runto t = 5000 time units. Concentrations of the X morphogen is color coded, with redhigh and blue low, relative to the X00 value on the boundary.43Figure 3.15: A finite element simulation showing a stable pure (3, 2) mode steady stateof the Brusselator in a spherical cap, viewed from above. The parameters used are thesame as in Figure 3.14, only the initial condition has changed, obtained from the stabler2-mode equilibrium for the Case 1 amplitude equations(3.11), with small added randomnoise. The mesh and color-coding used are the same as in Figure 3.14 and t = 5000.44Figure 3.16: A finite element simulation showing a stable mixed mode steady state ofthe Brusselator in a spherical cap, viewed from above. The parameters used are γ = 0.5,A = 76.513885, R = 0.96608, placing (σ1, σ2) in region 7 of Figures 3.4-3.6. The initialcondition is obtained from the stable mixed mode equilibrium of the amplitude equations(3.17), with small added noise. The mesh used and the color-coding is the same as inFigure 3.14 and t = 5000.45Chapter 4Delayed Bifurcation in a SlowlyChanging Spherical CapA way to produce a more realistic model is to have the curvature parameter γ forthe spherical cap domain as a slowly changing function of time, making the governingreaction-diffusion system non-autonomous by virtue of a slowly changing domain. Inthis section, we have this slowly changing γ cross a marginal stability curve for a con-stant domain. The result features a delayed bifurcation, where the bifurcating solutionconverges with the expected patterned equilibrium after a time delay. Other numericalstudies (e.g. [5], [21], [38], [41], [48]) have shown that slowly changing domains haveeffects on pattern formation, and it would be interesting to determine these effects inthe geometric context of our spherical cap model.In order to model the flattening tip in plant embryos we make the curvature para-meter γ a function of time. This allows the spherical cap to change curvature whileretaining a fixed radius R at its boundary. For a slowly changing domain we will wantthis curvature parameter to change at a rate of order ε, where ε is a small parameter.We can then define a slow time τ = εt and setγ = γ(τ) = γ(εt). (4.1)As a result, the time derivative becomes ∂∂t= ε ∂∂τ.As we change curvature over time the surface area of our domain also changes. Thischange of surface area adds extra terms in the system of equations in order to satisfythe Reynolds transport theorem [51] and maintain conservation of total chemical agents.Nonautonomous reaction-diffusion equations for chemical reactions on a growing surfaceΩ = Ω(τ) = Ω(εt) are∂X∂t= DX∇2Ω(εt)X − (∇Ω(εt) · v)X + f(X, Y ) (4.2)∂Y∂t= DY∇2Ω(εt)Y − (∇Ω(εt) · v)Y + g(X, Y ), (4.3)where f and g are the Brusselator kinetics (2.24), v is the velocity vector of a point onour surface, ∇Ω(εt) is the tangential gradient to the surface. In the constant domain,where γ is constant, we have v = 0 in (4.2), (4.3) and in this case the solution (2.25)represents a homogeneous equilibrium solution of the reaction-diffusion system. The46Laplace–Beltrami operator ∇2Ω(τ) is now dependent on the slow time τ = εt,∇2Ω(τ) =γ(τ)R2 sin θ[∂∂θ(sin θ∂∂θ)+1sin2 θ∂2∂ϕ2](4.4)in spherical coordinates, or∇2Ω(τ) =(cosh η − cos ξ(τ))2R2 sinh η[∂∂η(sinh η∂∂η)+1sinh η∂2∂ϕ2](4.5)in toroidal coordinates. We use the same equilibrium boundary values as in the constantdomain case(X, Y ) = (X00, Y00) on ∂Ω(τ). (4.6)The reasoning for this boundary condition is that the base of the spherical cap remainsconstant over time and thus concentrations should have the same value as in the constantcase.The form of the divergence term will depend on how exactly we want the surface tochange curvature. Here we choose to have points on the surface travel through pointson a circle where a coordinate in toroidal coordinates is constant (η in this case). Thishas the advantages of being able to be expressed by a fixed domain, having the surfacevelocity always normal to the surface and preserving the ratio of distances between apoint and the two foci, located at each end of the bisection of the spherical cap thepoint in question. In general, a time-dependent surface domain can be described as acontinuous mapping ψt of the initial surfaceΩ(τ) = ψt(Ω(0)) = Ω(εt).Under these circumstances this divergence factor can be expressed with the use of scalefactors, as done by Plaza et al. [48], for example. If our surface is described by twovariables, associated with the scale factors h1 and h2, then∇Ω(εt) · v = ∂∂t[ln(h1h2)].Using toroidal coordinates to compute the scale factors can greatly simplify thecomputations. If we keep using the spherical coordinates as before, the bounds onthe polar angle θ will change as we change the curvature and finding an appropriatetrajectory for a point on the moving surface becomes very cumbersome. Consequentlywe use toroidal coordinates and fix our parametrization of the spherical cap to ϕ and η,both having the same bounds for any given curvature.If one expresses the position x = (x1, x2, x3)ᵀ on the surface in toroidal coordinates(2.5), the scale factors are given byh1 =∣∣∣∣∂x∂η∣∣∣∣ = Rcosh η − cos ξ(εt) (4.7)h2 =∣∣∣∣∂x∂ϕ∣∣∣∣ = R sinh ηcosh η − cos ξ(εt) . (4.8)47The divergence term may thus be computed:∇Ω(τ) · v = ∂∂t[ln(h1h2)] =−2εξ′ sin ξcosh η − cos ξ =−2εzξ′R, (4.9)where the prime represents the τ derivative. We can then use the formulation of the x3coordinate to transform back to spherical coordinates:∇Ω(τ) · v = 2εγ′γ(cos θ√1− γ2 − 1). (4.10)We note here that when using spherical coordinates the θ-coordinate will have a timedependent upper bound.We find that the divergence term is proportional to the x3-coordinate, which willdepend both explicitly on time and on the distance from the point on the cap to thex3-axis. We thus obtain the following equations in terms of the slow time variable τ :ε∂X∂τ= DX∇2Ω(τ)X − εγ′(τ) Q(τ)X + f(X, Y )ε∂Y∂τ= DY∇2Ω(τ)Y − εγ′(τ) Q(τ)Y + g(X, Y ),(4.11)where Q(τ) depends on the position on the cap and may be expressed either in sphericalcoordinates byQ(τ)(θ) =2γ(τ)(cos θ√1− γ(τ)2 − 1), 0 6 θ < θmax(τ), (4.12)with θmax(τ) = arcsin γ(τ), or in toroidal coordinates byQ(τ)(η) =2ε tan ξ(τ)cosh η − cos ξ(τ) , 0 6 η <∞, (4.13)where ξ(τ) = pi − arcsin γ(τ). We observe that the function Q(τ) does not depend onthe longitudinal angle ϕ.We write the system of equations (4.11) in vector formε∂X∂τ= D∇2Ω(τ)X− εγ′(τ)Q(τ)X + f(X), (4.14)with boundary conditionsX = X00 on ∂Ω(τ), (4.15)Because the nonautonomous modifications to the equations are independent of thelongitudinal angle ϕ, the reaction-diffusion system with our choice of domain evolutionpreserves the rotation and reflection symmetries (2.38) and these will feature in thederivation of the reduced normal forms.484.1 The Quasi-Patternless SolutionUnder the slow domain evolution the concentrations X, Y change slightly so that theconstant patternless solution of the autonomous constant domain system is no longera solution of the nonautonomous system (4.14). This introduces the idea of a “quasi-patternless” solution that will depend on time. Due to the slow curvature change, weuse the slow timescale τ .We assume there is a slowly evolving, radially symmetric quasi-patternless solutionX0(τ, ε) that is O(ε)-close to the constant domain patternless solution, and coincideswith it when ε = 0. Such a solution would have an asymptotic seriesX0(τ, ε) = X00 + εX01(τ) +O(ε2), (4.16)where each X0j vector has an X and Y component, and the leading order term X00 isthe constant domain patternless solution (2.25).We reduce the problem to solving for the higher order correction terms in ε for thesolution. We start by substituting expansion (4.16) into (4.14, 4.15) to obtainε2∂X01∂τ+O (ε3) = εD∇2Ω(τ)X01 − εγ′(τ)Q(τ)X00 + f (X00 + εX01 +O (ε2))+O (ε2) .(4.17)The order ε terms must satisfy0 = D∇2Ω(τ)X01(τ)− γ′(τ)Q(τ)X00 + K0X01(τ), X01 = 0 on ∂Ω(τ). (4.18)We can then express the solution to (4.18) as a series of the now τ dependent eigen-functions of the Laplace–Beltrami operator Φmn(τ), where∇2Ω(τ)Φmn(τ) = −µmn(τ)Φmn(τ), Φmn(τ) = 0 on ∂Ω(τ) for all τ. (4.19)We use a series expression for the first correction termX01(τ) =∞∑n=1(α0n(τ)β0n(τ))Φ0n(τ). (4.20)Since (4.18) is independent of ϕ, we only write the series for order m = 0. This makesthe quasi-patternless solution invariant under rotations of the angle ϕ and reflections.To find an approximation of the quasi-patternless solution we solve for the firstseveral coefficients α0n and β0n, and again for an increasing number of n to check forconvergence. We found that computing modes from n = 1 to a few after the mostunstable m = 0 mode to be sufficient. We computed coefficients in the case when the(0, 3) mode is the most unstable and the difference between computing 5 modes andcomputing 12 modes is very small. We compare our truncated series approximation ofthe quasi-patternless solution with results of a numerical computation in Section 4.3.5.49We also expand the scalar function Q(τ) in an eigenfunction seriesQ(τ) =∞∑n=1qn(τ)Φ0n(τ). (4.21)We find the coefficients qn by applying a projection with a Legendre function in sphericalor toroidal coordinates, for example in spherical coordinatesqn(τ) =∫ arcsin γ(τ)0Q(τ)(θ)P 0λ0,n(τ)(cos θ) sin θ dθ∫ arcsin γ(τ)0(P 0λ0,n(τ)(cos θ))2sin θ dθ. (4.22)In the end, to solve for the coefficients α0n and β0n, we need to find the solution of thelinear equation(00)= −γ′(τ)qn(τ)(X00Y00)+(−DXµ0n(τ) + k1 k2k3 −DY µ0n(τ) + k4)(α0n(τ)β0n(τ)),(4.23)which is (α0n(τ)β0n(τ))= γ′(τ)qn(τ)A−1(τ, 0)(X00Y00), (4.24)where A(τ, 0) is the matrix in the right side of (4.23). We assume here that a criticaleigenvalue 0 of A(τ, 0) does not occur for any (0, n) mode for n ≥ 1, for any τ , so theequations are all solvable. If required, we could solve for correction terms of higher orderin ε.4.2 Deviation EquationsIn this section we re-express the slowly changing system (4.14)–(4.15) in terms of thedeviations from the quasi-patternless solution using the same methodology as in theconstant domain case.We define the deviation U from the quasi-patternless solution, byX = X0(τ, ε) + U, (4.25)where X0 = (X0, Y0)ᵀ is the quasi-patternless solution and U = (U, V )ᵀ. Substitutingthe expansion (4.16) in the system (4.14)–(4.15) with the use of expressions (2.30) to(2.34) and we get∂U∂t= A(τ, ε)U + B(τ, ε)(U,U) + C0(U,U,U), U = 0 on ∂Ω(τ), (4.26)where A, B are respectively nonautonomous linear and bilinear operators defined byA(τ, ε)U = D∇2Ω(τ)U + K(τ, ε)U− εγ′(τ)Q(τ)U, (4.27)50B(τ, ε)(U1,U2) = B0(U1,U2) + 3C0(X0(τ, ε)−X00,U1,U2), (4.28)and keeping terms up to first order in ε we haveX0(τ, ε) = X00 + εX01(τ) +O(ε2), (4.29)K(τ, ε) = K0 + 2B0(X0(τ, ε)−X00, ·) + 3C0(X0(τ, ε)−X00,X0(τ, ε)−X00, ·)= K0 + εK1(τ) +O(ε2),(4.30)K1(τ) =(2bBdaAX01(τ) +2aAcdY01(τ)2aAcdX01(τ)−2bBdaAX01(τ)− 2aAcd Y01(τ) −2aAcd X01(τ)), (4.31)3C0(X0(τ, ε)−X00,U1,U2) =(1−1)[(bBdaA+ εcY01(τ) +O (ε2))U1U2+(aAcd+ εcX01(τ) +O (ε2))](U1V2 + U2V1).(4.32)4.3 Pattern Formation in a Slowly ChangingDomainIn this section we describe emerging patterned solutions using a bifurcation theory met-hod adapted to this slowly evolving domain case. We start by analyzing the linearizationof the system (4.26) to determine the stability of the quasi-patternless solution. We thenperform a nonautonomous centre manifold reduction of the nonlinear system [49] witha projection method to build our nonautonomous normal form. Our analysis includesthe effects of the non-isotropically evolving domain.4.3.1 Linearization About the Quasi-Patternless SolutionThe linearization of (4.26) is∂Uˆ∂t= A(εt, ε)Uˆ, Uˆ = 0 on ∂Ω(εt). (4.33)Following WKB theory, e.g. [19], we assume the following ansatzUˆ(τ, ε) = eΨ(τ,ε)/ε U(τ, ε), (4.34)to obtainεU′(τ, ε) = A(τ, ε)U(τ, ε)−Ψ′(τ, ε)U(τ, ε), U = 0 on ∂Ω(τ). (4.35)Expanding in asymptotic seriesA(τ, ε) = A0(τ) + εA1(τ) +O(ε2),Ψ(τ, ε) = Ψ0(τ) + εΨ1(τ) +O(ε2),U(τ, ε) = U0(τ) + εU1(τ) +O(ε2),(4.36)51then substituting back into (4.35), we obtain that the leading order terms must solveA0(τ)U0(τ) = Ψ′0(τ)U0(τ), U0(τ) = 0 on ∂Ω(τ). (4.37)For each τ , this is an eigenvalue problem for the leading order term U0. If we setU0(τ) = u0,mn(τ)Φmn(τ), (4.38)where u0,mn(τ) ∈ R2, and Φmn(τ) are as defined in (4.19), we obtain an algebraiceigenvalue problem for u0,mn(τ)A0,mn(τ)u0,mn(τ) = Ψ′0(τ)u0,mn(τ), (4.39)whereA0,mn(τ) =(−DXµmn(τ) + k1 k2k3 −DY µmn(τ) + k4). (4.40)We obtain that Ψ′0(τ) = σ±mn(τ) is one of the two (now τ -dependent) eigenvalues of thematrix A0,mn(τ), so it solves the characteristic polynomial (2.51), with µmn = µmn(τ)depending on the slow time variable. Then u0,mn(τ) is the associated eigenvectoru0,mn(τ) =(k2DXµmn(τ)− k1 + σ±mn(τ)). (4.41)The order ε terms in (4.35) give us[A0(τ)−Ψ′0(τ)I]U1(τ) = Ψ′1(τ)U0(τ) + U′0(τ)−A1(τ)U0(τ),U1(τ) = 0 on ∂Ω(τ).(4.42)with Ψ′0(τ) = σ±mn(τ). We assume U1 has a series expansionU1(τ) =∞∑i=1u1,mi(τ)Φmi(τ), (4.43)then (4.42) becomes[A0,mi(τ)−Ψ′0(τ)I] u1,mi(τ) = δni[Ψ′1(τ)u0(τ) + u′0,mn(τ)]+ di(τ)u0,mn(τ)−A1,mi(τ)u0,mn(τ),(4.44)where δij is the Kronecker delta function, andΦ′mn(τ) =∞∑i=1di(τ)Φmi(τ), (4.45)A1(τ)u0,mn(τ)Φmn(τ) =∞∑i=1A1,mi(τ)u0,mn(τ)Φmi(τ), (4.46)52with Ψ′0(τ) = σ±mn(τ). The coefficients di are obtained by computing a time derivativeof the Laplace–Beltrami eigenfunction Φmi, considering trajectories with constant ηtoroidal coordinate. More details about this computation are provided in Appendix C.There will be two different cases whether i = n or not. In the latter case A0,mi(τ) −Ψ′0(τ)I is non-singular and we may solve (4.44) for u1,mi. In the case where i = n,A0(τ)−Ψ′0(τ)I is singular. We can still solve (4.44) if the right side of (4.46) satisfies thesolvability condition of being orthogonal to any null eigenvector u(∗)0 (τ) of the transposeof the matrix A0,mn(τ)−Ψ′0(τ)I. We obtain the following expressionΨ′1(τ) =u(∗)0 (τ)ᵀ[u′0(τ)− dn(τ)u0,mn(τ) + A1,mn(τ)u0,mn(τ)]u(∗)0 (τ)ᵀu0,mn(τ). (4.47)4.3.2 The Centre SubspaceWe choose parameter values in our model (4.2) - (4.6) so that (as in Section 2.4.1) forthe curvature parameter γ = γ0, the eigenvalues of the constant domain autonomouslinearization (2.39) are all real and negative, except for one isolated zero eigenvalueσ±m0n0 = 0. Furthermore for γ ∈ [γ1, γ2], with γ1 < γ < γ2, the eigenvalue σ±m0n0 is realand near zero with σ±m0n0 < 0 for γ = γ2 and σ±m0n0> 0 for γ = γ1, and the remainingeigenvalues are real and negative, uniformly bounded away from zero.Then we take γ = γ(τ) in the slow timescale, decreasing from γ2, through γ0, to γ1.For example, we may takeγ(τ) = γ0 − τ, (4.48)for τ belonging to a suitable interval.For parameter values corresponding to Figure 2.4, the parameter point (A, γ(τ))moves, as τ increases, to the left along the dotted arrow that crosses the (m0, n0) = (5, 1)marginal stability curve from the regime of stability into the region of instability.For the slowly changing domain, the leading order τ -dependent linearization (4.39)has eigenvalues Ψ0′(τ) = σ±mn(τ) with1. σ(0)0 (τ) = σ+m0n0(τ) is the only such value of all the σ±mn(τ) that is near 0, itsabsolute value remains small, for all τ , is strictly increasing with respect to τ andgoes from negative to positive as τ is increased.2. σ−m0n0(τ) and all other σ±mn(τ), (m,n) 6= (m0, n0) are strictly negative and largerin absolute value than |σ(0)0 (τ)| for all τ .Then from the WKB approximation (4.36), (4.38), for small ε the linearization (4.34)about the quasi-patternless solution makes a transition, as the slow time τ increases,from the stable to unstable regimes. We assume for small ε that the solution space of(4.34) splits into a direct sum of two invariant subspaces: a finite dimensional centresubspace associated with σ(0)0 (τ), and a complementary infinite dimensional stable sub-space. Solutions of (4.34) in the centre subspace decay or grow at a slow rate, while53solutions in the stable subspace decay much faster, exponentially, throughout the timeinterval.There are two linearly independent critical solutions to the linearized system (4.35).One of them isUˆ(0)(τ, ε) = eΨ(0)(τ,ε)/ε U(0)(τ, ε), (4.49)where ∂∂τΨ(0)(τ, ε) = σ(0)0 (τ) + O(ε), and the other critical solution is the complexconjugate of (4.49).Let X c(τ, ε) be the span of Uˆ(0)(τ, ε) and its complex conjugate, then this space istwo-dimensional and may be expressed asX c(τ, ε) = {zU(0)(τ, ε) + z¯U¯(0)(τ, ε)|z ∈ C} . (4.50)The action of the symmetries (2.38) on the space X c(τ, ε) has the following effects on z:z → eiϕ0z, z → z¯. (4.51)Let Pc(τ, ε) be a projection onto X c(τ, ε), which we assume has the formPc(τ, ε) = Pc0(τ) + εPc1(τ) +O(ε2). (4.52)We will only need the leading order term Pc0(τ), given byPc0(τ) = 〈U(∗)0 (τ),U〉U(0)0 (τ) + 〈U¯(∗)0 (τ),U〉U¯(0)0 (τ), (4.53)where the angled brackets denote the inner product〈U1,U2〉 =∫ 2pi0∫ arcsin γ(τ)0(U¯1U2 + V¯1V2) sin θ dθ dϕ=∫ 2pi0∫ ∞0(U¯1U2 + V¯1V2)sinh η sin2 ξ(τ)[cosh η − cos ξ(τ)]2 dη dϕ,(4.54)for Uj = (Uj, Vj)ᵀ, j = 1, 2, andU(∗)0 (τ) = N(∗)(τ)(k3DXµm0n0(τ)− k1 + σ(0)0 (τ))Φm0n0(τ) (4.55)is the solution of the adjoint problem to (4.39) with m = m0, n = n0, Ψ′0(τ) = σ(0)0 (τ)and the normalization coefficient N (∗)(τ) selected so that〈U(∗)0 (τ),U(0)0 (τ)〉 = 1 (4.56)for all τ . Let X s(τ, ε) be the kernel of P c(ε, τ), then the projection onto X s(τ, ε) isPs(τ, ε) = I−Pc(τ, ε). (4.57)544.3.3 Centre Manifold EquationsFollowing the analysis in the constant domain case, for the slowly changing domainwe reduce the nonlinear partial differential equation system (4.26) into an ordinarydifferential equation on a centre manifold. The solutions of the reduced system give agood approximation of the transition to the patterned state. In this case, however, theequations are nonautonomous, so the normal form coefficients we obtain after reductionwill depend on time.We split the function space X of all sufficiently regular functions on the sphericalcap satisfying homogeneous Dirichlet boundary conditions into the two subspacesX = X c(τ, ε)⊕X s(τ, ε). (4.58)We assume the splitting (4.58) is invariant under the linearization (4.26) for all τ , ε.Each function U ∈ X may then be written uniquely as a sum of two componentsbelonging to each subspace and each component may be retrieved using projections,U = Uc(τ, ε) + Us(τ, ε), (4.59)whereUc(τ, ε) = Pc(τ, ε)U ∈ X c(τ, ε), Us(τ, ε) = Ps(τ, ε)U ∈ X s(τ, ε). (4.60)Thus, equation (4.26) may be rewritten as the system∂Uc∂t= A(εt, ε)Uc + Pc(εt, ε)F(Uc(εt, ε) + Us(εt, ε), εt, ε) (4.61)∂Us∂t= A(εt, ε)Us + Ps(εt, ε)F(Uc(εt, ε) + Us(εt, ε), εt, ε), (4.62)with F(U, εt, ε) = B(εt, ε)(U,U) + C(U,U,U). We assume that there is an exponen-tially attracting centre manifold tangent to X c(εt, ε)Us = W(Uc, εt, ε) (4.63)for all Uc in the subspace X c(εt, ε) with sufficiently small magnitude. Equation (4.62)is then written as (cf. [49])∂∂tW(Uc, εt, ε) +∂∂UcW(Uc, εt, ε) [A(εt, ε)Uc + Pc(εt, ε)F(Uc + W(Uc, εt, ε), εt, ε)]= A(εt, ε)W(Uc, εt, ε) + Ps(εt, ε)F(Uc + W(Uc, εt, ε), εt, ε),(4.64)after substituting for ∂Uc∂t. The Taylor expansion of W around Uc = 0 starts at thesecond degree due to the tangency conditionW(Uc, εt, ε) = W(1)(εt, ε)(Uc,Uc) +O (‖Uc‖3) , (4.65)55where W(1)(εt, ε) is a symmetric bilinear form. Applying the product rule to the Ucderivative yields∂∂UcW(Uc, εt, ε)V = 2W(1)(εt, ε)(Uc,V) +O(‖Uc‖2‖V‖) (4.66)for any vector V.As a result of equivariance of the system of differential equations (2.36) under ro-tations and reflections (4.51) we may only study the real part of z, labeled x and thenapply a phase φ to get the other solutions z = xeiφ. Due to this we will continue withusing x in our derivations as a solution representative of the family of solutions.If we choose Uc = xU˜(0), where U˜(0) = Re U(0) for x real then equation (4.64) canbe used to find the order x2 terms with∂∂tW(1)(εt, ε)(U(0)(εt, ε),U(0)(εt, ε)) + 2W(1)(U(0)(εt, ε),A(εt, ε)U(0)(εt, ε))= A(εt, ε)W(1)(εt, ε)(U(0)(εt, ε),U(0)(εt, ε)) + Ps(εt, ε)B(U0(εt, ε),U(0)(εt, ε)),(4.67)dropping the tilde for the real part of U(0). By definingU(1)(εt, ε) = W(1)(εt, ε)(U(0)(εt, ε),U(0)(εt, ε)), (4.68)we can write the solution U as a series in xU = xU(0)(εt, ε) + x2U(1)(εt, ε) +O(x3). (4.69)We use εt = τ as well asA(τ, ε)U(0)(τ, ε) =∂Ψ(0)∂τ(τ, ε)U(0)(τ, ε) + ε∂∂τU(0)(τ, ε) (4.70)and∂∂τU(1)(τ, ε) =∂∂τW(1)(τ, ε)(U(0)(τ, ε),U(0)(τ, ε))+ 2W(1)(τ, ε)(U(0)(τ, ε),∂∂τU(0)(τ, ε)),(4.71)to get the expression[ε∂∂τ−A(τ, ε) + 2∂Ψ(0)∂τ(τ, ε)]U(1)(τ, ε) = Ps(τ, ε)B(τ, ε)(U(0)(τ, ε),U(0)(τ, ε)),(4.72)which we can use to solve for U(1)(εt, ε) = U(1)0 (εt) +O(ε), assuming that U(1)0 (εt) hasa series expressionU(1)0 (εt) =∞∑j=1[u(1)0j (εt)Φ0j(εt) + u(1)2mj(εt)Φ2mj(εt)]. (4.73)564.3.4 Normal FormHaving found the leading terms in the centre manifold, we can use a projection on thecentre subspace X c(εt, ε). This is given by equation (4.61)∂∂tUc = A(εt, ε)Uc + Pc(εt, ε)F(Uc + Us, εt, ε). (4.74)We proceed here as if Wc(εt, ε) is such thatUc = xU(0)(εt, ε), Us = x2U(1)(εt, ε) +O (x3) . (4.75)Using this in (4.74) yieldsdxdtU(0)(εt, ε) + xε∂∂τU(0)(εt, ε) = xA(εt, ε)U(0)(εt, ε)+ Pc(εt, ε)F(xU(0)(εt, ε) + x2U(1)(εt, ε) +O (x3) , εt, ε). (4.76)We then use (4.70) and carry out a projection in the direction of U(0)(εt, ε) to obtainx˙ = σ(0)(εt, ε)x+ C(εt, ε)x3 +O (x5) , (4.77)withσ(0)(εt, ε) =∂∂τΨ(0)(εt, ε) = σ(0)0 (εt) + εσ(0)1 (εt) +O(ε2), (4.78)C(εt, ε) = C0(εt) +O(ε) (4.79)C0(εt) =〈U(∗)0 (εt), 2B0(U(0)0 (εt),U(1)0 (εt))+ C0(U(0)0 (εt),U(0)0 (εt),U(0)0 (εt))〉.(4.80)Equations (4.39), (4.47) and (4.80) allow us to find approximate normal form coeffi-cients for any time t. In order to use these equations we need a good approximation forthe infinite sums of eigenfunctions of modes 0, m0 or 2m0 present, for example, in thequadratic component of the centre manifold U(1) or the derivative of the eigenfunctionΦ′m0n(τ), both of which are found through equations (4.73) and (4.45). For m0 = 5we found that truncating the series after the first five terms for each order yield similarresults to computing the first eight or twelve terms, which suggests that truncated seriesapproximations are sufficiently accurate.Solving the Normal Form EquationsTo solve the normal form ordinary differential equation, we need the coefficient valuesthroughout the time interval. As the coefficients do not have an algebraic form thatcan easily solved using symbolic computing tools, we cannot solve at all time values.To get a good approximation, we collect the coefficient values at a sample of equally57Figure 4.1: Normal form solutions and the constant domain bifurcation diagram forparameter values (2.52), A = 76.5198 and different ε values. The values of ε are 1×10−6for the long-dashed orange curve, 3 × 10−7 for the dotted gold curve, 1 × 10−7 for thedot-dashed green curve and 3 × 10−8 for the dashed blue curve. Each solution has thesame initial conditions γ(0) = 0.51, x(0) = 0.002305. The solid black curve is theconstant coefficient pitchfork bifurcation diagram.spaced times and interpolate between them. We used a sample of 36 time points andspline interpolation using the CurveFitting tool in the Maple computational softwarepackage, which uses piecewise cubic polynomials to create an interpolation that matchesall sample point values and is globally twice continuously differentiable.Figure 4.1 show solutions to the normal form equations for different ε values supe-rimposed with the constant domain pitchfork bifurcation nontrivial stable equilibriumbranch. We clearly see that a time dependent cubic coefficient C(τ, ε) is important toget a good prediction out of our normal form reduction.The normal form solutions were computed using the dsolve function in the Maplesoftware package using a fourth-and-fifth order Runge–Kutta method with adaptivetime stepping. We chose the initial condition x(0) = 2.3057 · 10−3 to match the initialcondition used in the numerical simulations.4.3.5 Numerical simulationsTo validate the normal form predictions we computed numerical simulations of solutionsto the nonautonomous Brusselator system (4.14) – (4.15) using the closest point met-hod [53] with SBDF2 implicit-explicit time stepping [1] and a five-point stencil spatial58discretization (2.74) for the Laplacian operator. We used parameter values (2.52) andA = 76.51981. (4.81)We first checked the accuracy of our truncated series approximations of the quasi-patternless solution, by computing simulations of a modified version of (4.26), wherethe linear operator A is replaced with the affine linear operator defined byA˜(τ, ε)U = D∇2Ω(τ)U + K(τ, ε)U− εγ′(τ)Q(τ) [U + X00] , (4.82)where U is the deviation from the patternless solution of the constant domain system.We used the same parameter values except that we had our curvature index γ(τ) go from0.5015 to 0.4915, where the quasi-patternless solution remains stable. The simulationsran for 10000 time units, with ε = 10−6, and with a spatial grid size h = 0.02 and timestep size of 0.1. We extracted the numerical solution of U over the polar angle θ for twolongitudinal angles ϕ separated by pi/4 and compared it with the five-term truncatedseries from our normal form calculations. The truncated series agrees well with thesimulations and the simulations converge with decreasing spatial step h. Truncatedseries with twelve terms were also compared and yielded very similar results. Figures4.2 and 4.3 show a numerical simulation that show the quasi-patternless solution on aspherical cap and a cross-section along one longitudinal angle respectively. In Figure4.3 the simulation is compared to a truncated series approximation.Figure 4.2: Quasi-patternless solution of (2.36) with modified affine operator (4.82),simulated using the closest point method, after 10000 time units from a homogeneousinitial condition with h = 0.02, ∆t = 0.1. Parameter values are ε = 10−6, γ = 0.4915,A = 76.51981, and the same as (2.52).We computed other numerical simulations of the full system (4.26) that show thedelayed bifurcation. We ran a typical simulation for 40000 time units, with ε = 10−6,59Figure 4.3: Quasi-patternless solution simulation computed using the closest point met-hod (blue dots) for h = 0.02 with the same parameter values as in Figure 4.2, comparedto a five-term truncated series approximation (red curve) outlined in (4.20) with coeffi-cients as computed in (4.24).and with a time step size of 0.1 and different grid sizes between 0.04 and 0.015. Due tothe slow domain evolution, we updated the domain geometry only every 50 time stepsinstead of every time step to save computation time. We used an initial conditionU(0) =0.021u0,51 · M˜51· u0,51 · Φ˜51, (4.83)where Φ˜51 is the approximate (5, 1) eigenfunction of the Laplace–Beltrami operatorcomputed using the closest point method [35] and the MATLAB eigs function, M˜51is its approximate maximal value, u0,mn is computed by (4.41) and u0,51 = k2 is itsfirst coordinate. This makes 0.021 the maximal value of the initial condition. We useε = 10−6 and the curvature starts at γ = 0.4915 and ends at γ = 0.4515. We startbelow the constant domain critical curvature γ0 = 0.5 because, as seen in Section 4.3.4the order ε terms shift the transition value for γ, where solutions go from decreasingto increasing in time, to a lower value. We want to start as close as possible to thistransition value in order to see solutions “level off” to the patterned relative equilibriumat later times, while also keeping computation times reasonable.To make the comparison between the numerical simulations and normal form solu-tions, we extract the (5, 1) part of the numerical simulation using the projection (4.53)via numerical integration. We then solve for the x value by dividing the maximum ofthe residue by u0,51(τ) times the maximum of the Legendre function at each time. Forexample consider the initial condition stated in (4.83) with an exact eigenfunction Φ51.It lies in the centre subspace, therefore the projection will yield the same functionPcU(0) = U(0). (4.84)60Assuming the centre subspace and centre manifold to be close at this value we use thereal version of the centre subspace definition (4.50) and comparingU(0) = x(0)U(0)(0, ε) = U0 =0.021u0,51 ·M51(0) ·U(0)(0, ε), (4.85)we obtainx(0) =0.021k2 ·M51(0) ≈ 2.3057 · 10−3. (4.86)For a later time t a similar computation is used to find x(t),x(t) = ‖U(t)‖51k2 ·M51(εt), (4.87)where ‖ · ‖51 denotes the L∞ norm of the (5, 1)-mode projection and M51(εt) is theabsolute maximum of the (5, 1) Legendre function for curvature γ(εt).Figure 4.4 shows the resulting equivalent x(t) values for the closest point methodsimulations for spatial grid sizes between h = 0.04 and h = 0.015 (dashed or dottedlines) compared with the equivalent normal form computation (solid line). The twoinitial conditions are chosen so that each numerical simulation match the normal formsolution to leading order in x and ε at the initial time. We see that as the spatialgrid is refined the simulations approach the results of the normal form computation.Snapshots of a h = 0.02 solution are provided in Figure 4.5 for a random uniform noiseof amplitude of about the same scale as the final solution at γ = 0.4315. We can seefrom panels (a) to (b) that the noise quickly decays, then the remaining (5, 1) dominantsolution increases in amplitude as x increases.For each simulation in Figure 4.4 we also computed a prediction from the normalform equations using the leading term of the centre manifold (4.69) as well as a truncatedseries of the quadratic terms for the dominant modes for the final curvature value ofγ = 0.4515. We then computed an error approximation by comparing the numericalsolution to this prediction for the different spatial grid sizes h. Figure 5.19 showsthe differences compared to the grid size h superimposed with a reference quadraticfunction to show that we have quadratic convergence in both L2 and L∞ norms, thelatter implying pointwise convergence. This is expected convergence behaviour of theclosest point method on constant surfaces [53], so this result validates both of numericalsimulations and normal form results. More details on the numerical method used maybe found in Section 5.7.61Figure 4.4: Numerical simulations compared to the normal form solution as a functionof curvature index γ by extracting (5, 1) coefficients using a projection. The spatial gridspacings are h = 0.04 for the red dash-dotted curve, 0.03 for the turquoise long dashedcurve, 0.02 for the black dotted curve and 0.015 for the purple dashed line and ε = 10−6for each curve. The normal form solution is the blue solid curve. The horizontal axisrepresents the curvature index γ so solutions go from right to left as time increases. Thevertical axis is the variable x for the normal form or the appropriately scaled L∞ normof the (5, 1)-mode projection of the numerical solutions.62(a) γ = 0.4915 (b) γ = 0.4875(c) γ = 0.4675 (d) γ = 0.4515Figure 4.5: Snapshots of numerical solutions for U using grid size h = 0.02 and timestep ∆t = 0.1 running for 60000 time units at four different curvature values. Parametervalues used are listed in (2.52) and the curvature of the domain linearly decreases fromγ = 0.4915 to γ = 0.4315, corresponding to a value of ε = 10−6. The initial conditionsused is uniformly random noise of amplitude 0.05 for U and size 0.005 for V . Thesolution will eventually approach the black dotted curve in Figure 4.4.63Chapter 5Closest Point Method on EvolvingDomainsIn this chapter we present an extension to the closest point method that allows for simu-lations on evolving domains. Such numerical simulations of partial differential equationsare usually performed using finite element methods [48], [38], where the domain meshmay be updated in between time steps according to the domain velocity field. Previ-ous attempts of simulating growing domains using the closest point method have beenperformed on similar reaction diffusion systems by substituting domain growth withequivalent diffusion coefficient modifications [60], [52]. Actual growing domains havebeen simulated using a closest point method by using a grid based particle method orradial basis functions approximations with the closest point method [47]. The methodpresented here that can simulate evolving domains while retaining the classical closestpoint method steps.5.1 Reaction-Diffusion Systems in a TimeDependent DomainLet Ω(t) be a domain in Rn space at time t with any point ~x(t) in the domain followinga fixed trajectory in time, which will generate the time-dependent velocity field~s(~x, t) =d~xdt. (5.1)When solving partial differential equations on such time dependent domains, we needto adjust the equations for regions in the domain that are contracting or expanding,where solutions are augmented and reduced respectively, compared to the same equationresolved on a constant domain. It can be shown using the Reynolds Transport Theoremthat for a reaction-diffusion system of equations there is an extra divergence term toaccount for this∂u1∂t= D1∇2Ω(t)u1 + f(u1, u2)−∇Ω(t) · (u1~s),∂u2∂t= D2∇2Ω(t)u2 + g(u1, u2)−∇Ω(t) · (u2~s),(5.2)64where uj ∈ R, with j = 1, 2 being the number of equations, Dj are diffusion coefficients,∇2Ω(t) is the Laplace–Beltrami operator for the domain Ω at time t, ∇Ω(t) is the tangentialgradient on the domain Ω(t).In this chapter, we mostly use simple domains such as circles and spheres or sectionsof circle and spheres in our simulations. These domains work very well with the closestpoint method and their divergence terms are readily computable by hand. There aredifferent possibilities on how to evolve the domains in time. We list here a few of themand how they affect the system of partial differential equations.Note that in this chapter we will use the arrow notation for vectors in real space andbold typesetting for variables, vectors and operators defined on a discrete grid.5.1.1 Curves in Two-Dimensional SpaceLinear Domain GrowthOur first example is of a isotropically linearly growing circle. Isotropic growth meansdomains retain the same relative distance between any two points. The radius of thecircle follows the following functionR(t) = R0 · (1 + kt), (5.3)where R0 is the initial radius and k is the constant rate of growth. The velocity fieldwill be~s(t) = kxˆ, (5.4)where xˆ is the unit vector in the direction of ~x from the circle’s centre (usually chosento be the origin without loss of generality). The divergence term will be∇Ω(t) · (~su) = (∇Ω(t) · ~s)u = ku|~x| =kuR(t). (5.5)The first equality is due to the surface being always normal to ~s, so by the product rule∇Ω(t) · (~su) = (∇Ω(t) · ~s)u+∇Ω(t)u · ~s, (5.6)and the gradient of u has to be tangent to the surface, thus normal to ~s and the secondterm is equal to zero.Exponential Domain GrowthThe second example is of a circle growing isotropically and exponentially in time. Theradius function isR(t) = R0eρt, (5.7)65where R0 is the initial radius and ρ is a positive constant relative rate of change. Thevelocity field is~s = ρR(t)xˆ = ρ~x. (5.8)The divergence term is∇Ω(t) · (~su) = (∇Ω(t) · ~s)u = ρu. (5.9)Because the divergence term is a linear term it becomes much simpler to compute exactsolutions for some partial differential equations.Logistic Domain GrowthThe next example is another isotropic growth normal to the surface. The radius of thecircle follows a logistic functionR(t) = R0eρt1 + 1ζ(eρt − 1) , (5.10)where R0 is the initial radius, ρ is the initial relative rate of change and ζ is the relativesize of the limiting domain. The velocity field is~s(t) = ρR(t)(1− R(t)R0ζ)xˆ = ρ(1− R(t)R0ζ)~x. (5.11)The divergence term is thus∇Ω(t) · (~su) = (∇Ω(t) · ~s)u = ρ(1− R(t)R0ζ)u. (5.12)Domain Growth Along Bipolar TrajectoriesToroidal and bipolar trajectories (2.5)–(2.4) provide an example of a non-isotropic grow-ing domain. We introduce the bipolar coordinatesx =R sinh ηcosh η − cos ξ , y =R sin ξcosh η − cos ξ , (5.13)for η ∈ [0,∞) and ξ ∈ [pi/2, pi]. The curves of constant ξ will create a circle arc withendpoints at (−R, 0) and (R, 0) and centre at (0, R cot(ξ)). We may extend the arc intoa full circle by adding the arc of ξ+pi or equivalently extending the η values with η+pii.To find our divergence term we will use scaling factors, as in [48]∇Ω · ~s = ∂∂t(ln(h)), (5.14)66where h is the scale factor for the bipolar coordinates η used to describe our domainh =∣∣∣∣ ∂∂η~x(ξ, η)∣∣∣∣ . (5.15)We may use the η partial derivative on each coordinate and find the norm to findh =Rcosh(η)− cos(ξ) , (5.16)and then find our dissipation factor∇ · ~s = ∂∂tln(Rcosh(η)− cos(ξ))=−ξ˙ sin(ξ)cosh(η)− cos(ξ) =−ξ˙zR. (5.17)We may then use a function for ξ(t) according to our needs. For a given ξ value theradius of the circle isr =Rsin(ξ). (5.18)We can then solve for ξξ = arcsinRr. (5.19)If we know a function for r(t), then we can find ξ˙ξ˙ =−Rr˙r√r2 −R2 . (5.20)For example we can build an exponential growth circle such that at all times each pointmoves along bipolar trajectories. We haver(t) = r0eρt (5.21)ξ˙ =−Rρr0eρtr0eρt√r20e2ρt −R2 =−Rρ√r20e2ρt −R2 (5.22)[Note] This creates a singularity at ξ = pi/2, look into the substitution u = cos ξ andderive directly instead.In the case of logistic growth, we haveR(t) = R0eρt1 + 1ζ(eρt − 1) , (5.23)675.1.2 Surfaces in Three DimensionsWe can adapt each of the two-dimensional domain transformations into a three-dimen-sional domain by making a surface of revolution over a longitudinal angle ϕ ∈ [0, 2pi).As an example, we consider toroidal coordinates. In this case surfaces of constantξ values will create a spherical cap with circular base embedded in the xy-plane andradius R. We may append the ξ + pi surface to create a full sphere by extending the ηvalues with η + pii.The divergence term is now expressed by∇Ω · ~s = ∂∂t(ln(hηhϕ)), (5.24)where the hκ represent the scale factors for the variable κ = η, ϕhη =∣∣∣∣ ∂∂η~x(ξ, η, ϕ)∣∣∣∣ , hϕ = ∣∣∣∣ ∂∂ϕ~x(ξ, η, ϕ)∣∣∣∣ . (5.25)The derivatives and norms yieldhη =Rcosh(η)− cos(ξ) , hϕ =R sinh ηcosh(η)− cos(ξ) , (5.26)and the dissipation factor becomes∇ · ~s = ∂∂tln(R2 sinh η(cosh(η)− cos(ξ))2)=−2ξ˙ sin(ξ)cosh(η)− cos(ξ) =−2ξ˙zR. (5.27)Generally, for a three-dimensional domain that is a surface of revolution of a two-dimensional curve around a longitudinal angle and normal evolution, the dissipationfactor will be multiplied by two.5.2 Pattern Formation for Reaction-DiffusionSystemsIn this section we briefly review the pattern formation concepts that will be testednumerically in later sections. These patterns emerge from systems of reaction-diffusionpartial differential equations as in equation (5.2) for specific reaction functions f andg. Pattern formation models may have more than two reactants, hence more than twoequations, but here we demonstrate on systems of two equations.The two systems we use are the Shnakenberg model(f(u1, u2)g(u1, u2))=(a− u1 + u21u2b− u21u2), (5.28)68for a, b positive parameters and the Brusselator model(f(u1, u2)g(u1, u2))=(a− b1u1 + cu21u2b2u1 − cu21u2), (5.29)for a, b1, b2, c positive parametersThe classic model for Turing pattern formation [56] has a short range activatorreactant and a long range inhibitor, so D1 6= D2.Reaction-diffusion systems often feature a uniform steady-state solution, which wecall the patternless solution. For example, the Schnakenberg model has the followingpatternless solutionu01 = a+ b, u02 =b(a+ b)2. (5.30)For several configurations of model and system parameters, such as for small domains,only this patternless solution is stable. However, changing the parameters, for exampleincreasing the domain size, the patternless solution undergoes a Turing instability anda patterned state will emerge.Another bifurcation event on curves or surfaces are localized spike patterns, if thedifference between diffusion coefficients is much larger. In this case, the emerging pat-terns are sharp spikes surrounded by relatively large regions where solutions are almostconstant. In an evolving domain more spikes may coexist as the domain increases insize, which can lead in emergence or replication events [52].5.3 Adaptation of the Closest Point Method on anEvolving DomainIn this section we describe how one may use the closest point method to solve systemsof differential equations on evolving domains. We gave a brief description of the closestpoint method in Section 2.7.Let Ω = Ω(t) now be a domain evolving in time, with velocity vector field ~v. Thethree principles (2.71), which relate intrinsic surface differential operators to their Car-tesian counterparts, are preserved with the inclusion of the time dependent extensionoperator E(t):∇(E(t)u)(y) = ∇Ω(t)u(y), y ∈ Ω∇ · (E(t)g)(y) = ∇Ω(t) · g(y), y ∈ Ω (5.31)∇2(E(t)u)(y) = ∇2Ω(t)u(y), y ∈ Ω.Here the differential operators are defined as being equal to the constant domain ope-rators at any fixed time t.695.3.1 Time Derivatives of Closest Point ExtensionsWhen applying the closest point method on static surfaces, the extension operatorcommutes with the time derivative. With moving surfaces there are additional terms tobe considered because the closest point function cp(~x, t) is itself time dependent.Let v(~x, t) = E(t)u(~x, t), where E(t) is the extension operator at time t. This isequivalent to the expression v(~x, t) = u(cp(~x, t), t). Consequently, v is always constantin the normal direction to the moving surface and E(t) is idempotent, that is v = E(t)v.To relate the time derivative of v to that of u, we apply the chain rulevt = (Eu)t = (u(cp(~x, t), t)t= ut(cp(~x, t), t) +∇Ωu(cp(~x, t), t) · dcp(~x, t)dt= ut(cp(~x, t), t) + (∇v) (cp(~x, t), t) · dcp(~x, t)dt= E(t)ut + (E(t)∇v) · dcp(~x, t)dt,(5.32)for ~x in some embedding space around Ω(t). Here the first principle of (5.31) was used toreplace the surface gradient ∇Ω(t)u with derivatives of v. Note that while v is constantin the normal direction to the surface, ∇v is not generally constant in that direction.There are two possible approaches when computing the time derivative of the closestpoint function which are analogous to Eulerian and Lagrangian descriptions of velocityfields in fluid mechanics.5.3.2 Eulerian ApproachIn the Eulerian Case, points in the embedded space (which includes the band pointsintroduced in section 2.7.1) are considered fixed and thus the derivative of the closestpoint function is only the partial derivative of the functiondcp(~x, t)dt=∂cp(~x, t)∂t. (5.33)This derivative is the change in position of the closest point of a fixed band point as thesurface moves in time. This change is driven by two different types of surface motion.The first one is translation normal to the surface and changes the closest point functionin the same direction at the same rate and thus will be canceled by the inner productwith the surface gradient in equation (5.32). The other is local rotation of the surfaceat the closest point. A rotation will change the position of the closest point along thesurface in the direction of the sign of the rotation rate, as pictured in Figure 5.1.Demonstration: Aligning the x-axis on the surface tangent makes its slope m(t) =tan(θ(t)). The angle θ(t) is initially zero at time t. If the initial closest point is at the70d(t+∆t)d(t)θ˙(0, d(t))m(t) = 0m(t) = tan(θ(t))Figure 5.1: Diagram showing the evolution of the closest point function on a surfacerotating with respect to the closest point.origin, the new closest point after a small angle change may be computed using theclosest point to a line formulaX2 =m(t)d(t)m(t)2 + 1, Y2 =m(t)2d(t)m(t)2 + 1, (5.34)where d(t) is the distance of the band point to its closest point. Taking a time derivative,we get Y˙2 = 0 and thus our only change is in the x directionX˙2 = m˙(t)d(t) = θ˙d(t). (5.35)Finally, using small angles approximations the derivative of the angle or of the slope isthe directional derivative of the velocity field in the x directionθ˙ = ∂x(nˆ · ~s). (5.36)In surfaces embedded in three dimensions, the rate of change of the angle will be thesurface gradient of the normal component to the vector field.∂cp(~x)∂t= d(t) · ∇Ω(nˆ · ~s). (5.37)We note that in the numerical methods presented in the following sections the translationdue to this term is performed by “backtracking operators”, where instead of taking bandpoint values for previous time steps, we take the values at previously held positions ofthe current closest point and this translation is no longer needed.5.3.3 Lagrangian ApproachIn some of our numerical computations this drift of the closest point function will behandled automatically by the backtracking interpolation step, see Section 5.3.4, whichdoes the exact inverse process, in a Lagrangian fashion.71In the Lagrangian Case, we assume each band point moves along with its closestpoint,~x = ~x(t), ~˙x = ~s(cp(~x)). (5.38)The time derivative of the closest point function is thendcp(~x, t)dt=∂cp(~x, t)∂t+ Jcp~˙x. (5.39)The first term is the same as in the Eulerian Case, the second uses the Jacobian of theclosest point function, Jcp, which is equal to the projection on the surface [39]vt = ut(cp(~x), t) +∇Ωu(cp(~x), t) · ( ~Jcp~˙x), (5.40)where ~˙x is the closest point extension of the domain velocity field, thus ~Jcp~˙x = ~sΩ, thetangential velocity of the domain andvt = ut(cp(~x), t) +∇Ωu(cp(~x), t) · ~s. (5.41)This is the tangential part of the divergence term that is generated from the Reynoldstransport theorem. This arises because of the choice of tracking each position back intime along its domain evolution trajectory, which generates advection in the tangentialdirection to the surface.5.3.4 Time-Stepping in an Evolving DomainThe main challenge of adapting the closest point method in an evolving domain is toperform addition and substraction of solutions at two different time steps. In general, apoint in the embedding space at a specific time step no longer has the same closest pointon the surface at the next step. Moreover, when applying the method in simulations,solutions using this method are stored in a band of points and this band may increaseor decrease in size from one time step to another.To overcome this, we need to express solutions associated to previous time stepson a grid that has the same size as the current step and compensates for the movingclosest point positions. This is achieved by backtracking each closest point of the currenttime step to the previous surface, then interpolating past solutions on these backtrackedpoints using past bands. The way to perform this backtracking depends on whether weuse the Eulerian or the Lagrangian version of closest point function time derivative.In the Eulerian Case, the backtracking has to be normal to the past surfaces, whichcan be done by using the closest point function of a past surface evaluated at the currentclosest points and then interpolating using the values from the past band. The resultis thus an interpolation of the surface solution at the previous time that maps to thecurrent band points. In this case we numerically approximate the divergence term using72the velocity vector field and finite difference first derivative operators. Therefore thiscase is usually preferred in the cases where the velocity flow is known at each time step.Figure 5.2a illustrates the different positions and closest point functions.In the Lagrangian Case we need to interpolate the solutions on the past bandson the positions of the current closest points on the appropriate past times. Thiscan be done in a straightforward way if the domain evolution follows a parametrizedgrowth function. If no such function is available, we reverse the surface velocity field toadequately approximate the past surface positions. For example with the circle usingbipolar coordinates we may solve for each η coordinate using only the y coordinate ofeach closest points and the sign of their x coordinate. We then just need to generatenew points on past surfaces by changing ξ accordingly. Using this method we mustdiscard the tangential component of the divergence term as any tangential movementof the domain has already been considered by the backtracking. Figure 5.2b illustratesthis process for one point.The number of additional interpolation operators to be computed is equal to thenumber of past points required to use a specific time stepping algorithm. For examplethe explicit or implicit Euler methods use only one past point, so we only need one extraoperator. The semi-implicit backward differential formula of order n (SBDFn) [1] usesn past time steps, so we need n additional operators on each time step.5.3.5 Time-Stepping on Moving SurfacesIn this section we provide a few demonstrations of time-stepping schemes on the heatequation and a reaction-diffusion system. We work our way up from an explicit for-ward Euler scheme on the heat equation to an implicit-explicit SBDF scheme on areaction-diffusion system. We remind that in this section, variables in bold typesetting(for example u) represent the spatially discretized approximations of the solution andoperators or vectors restricted to the grid.Heat Equation with Forward Euler Time-SteppingWe start by showing the explicit forward Euler time discretization of the heat equationon a moving surface Ω(t) embedded in R3. Suppose u(x, t) is a function on Ω(t) thatsolvesut = ∇2Ω(t)u−∇Ω(t) · (~su). (5.42)If we extend both sides of this equationE(t)ut = E(t)(∇2Ω(t)u−∇Ω(t) · (~su)) ,and suppose that v = E(t)u, then applying the derivative of the extension operator, asin (5.41), we havevt = E(t)(∇2Ω(t)u−∇Ω(t) · (~su))− (E(t)∇Ω(t)u) · dcp(~x, t)dt . (5.43)73Ω(t) Ω(t−∆t)xcp(x, t)cp(cp(x, t), t−∆t)(a) Closest point backtracking, Eulerian version.Ω(t) Ω(t−∆t)xcp(x, t)Φ−∆t(cp(x, t))(b) Closest point backtracking, Lagrangian version.Figure 5.2: Comparison of the Eulerian and Lagrangian approaches to backtracking.74Notice there are still some surface functions u in this expression: we would like tocompute these using v, functions in the embedding space around Ω(t).Following the original approach of Ruuth and Merriman [53] we consider discretizingin time first, while staying continuous in space for the moment. Suppose at some timetn we have a function v(~x, tn), defined in a narrow band surrounding the surface whichis constant in the normal direction to the surface. For example, v(~x, tn) could be theextension of the exact solution E(tn)u(~x, tn). Suppose we also have ~se = E(t)~s which isthe closest point extension of the surface vector field s.Now consider the following expressionF (~x, tn) := ∇2v(~x, tn)−∇ · (~sev(~x, tn))− (∇v)(~x, t) · dcp(~x, t)dt, (5.44a)which is the right-hand side of (5.43), but with Cartesian differential operators insteadof the surface intrinsic ones. Because v(~x, tn) is constant in the normal direction, wecan apply the principles (5.31) to see thatF (~x, tn) = ∇2Ω(tn)u−∇Ω(tn) · (~su)−∇Ω(tn)u ·dcp(~x, t)dt, for ~x ∈ Ω(tn) (5.44b)Because these expressions match momentarily at time tn for points on the surface Ω(tn),we may apply a forward Euler step to compute a provisional approximate solution attime tn+1 = tn + ∆t by usingv˜(~x, tn+1) := v(~x, tn) + ∆t F (~x, tn). (5.44c)However, because F (~x, tn) is not constant in the normal direction, the provisional solu-tion will not be constant on normal directions to the surface. So we follow the forwardEuler step with an extension:v(~x, tn+1) := E(tn+1)v˜(~x, tn+1). (5.44d)This time-stepping scheme (5.44) is not yet a full discretization; we must next discretize(5.44a) in the Cartesian embedding space (R3 in this case)Fn = Lun −D1(s1un)−D2(s2un)−D3(s3un)−R(~x, tn), (5.45)where un represents the solution to w at the time step tn, L represents the discretizedLaplace operator, which could be computed from equation (2.74), Dj corresponds tothe first spatial derivative discretization in the xj direction, which could be computedfrom equation (2.75), sj corresponds to the jth coordinate of the velocity field and R isthe discretization of the last term in equation (5.44a)R(~x, tn) =D1u(cp(~x, tn), tn)D2u(cp(~x, tn), tn)D3u(cp(~x, tn), tn) · dcp(~x, tn)dt. (5.46)75In order to use a compact band around our time dependent surface that can poten-tially change between time steps, we apply the backtracking operator E(n)n+1 on un,1∆t(un+1 − E(n)n+1un) = E(n)n+1(Lun −D1(s1un)−D2(s2un)−D3(s3un)). (5.47)This matrix operator interpolates the solution at time n at the positions where closestpoints of step n + 1 previously were. As we previously mentioned, the use of thebacktracking operator removes the need to include the term R(~x, tn).We then re-extend following (5.44), we get the following fully discrete method toobtain un+1,u˜n+1 = E(n)n+1un + ∆tE(n)n+1(Lun −D1(s1un)−D2(s2un)−D3(s3un)), (5.48a)un+1 = En+1u˜n+1 (5.48b)to get the next time step on the band. By applying the interpolation operator matrixEn+1 for the surface at the time step n+ 1 we force the solution to remain constant inthe normal direction to the surface.Heat Equation with Backward Euler Time-SteppingConsider the extended system (5.43) from before. Now, assuming v = E(t)u for alltime, we may use the principles (5.31) to obtain the systemvt = ∇2v(~x, tn)−∇ · (~sev(~x, tn))− (∇v)(~x, t) · dcp(~x, t)dt. (5.49a)We apply the implicit backward Euler time discretization and spatial discretizationswhich yields1∆t(un+1 − un) = En+1Lun+1 − En+1D1(s1un+1)− En+1D2(s2un+1)− En+1D3(s3un+1)−R(~x, tn+1) + 2 · dimh2(I − En+1)un+1,(5.50)where the penalty term from [59], [36] has been added to stabilize the computation. Wecan collect this penalty term and lump it together with the extended Laplacian,Lˆ = En+1L +2 · dimh2(I− En+1), (5.51)where in the case of a surface embedded in R3, dim = 3. This Lˆ operator becomes ourdiscrete Laplace–Beltrami matrix operator.We can combine the terms in the right hand side of equation (5.50), excluding R, ina single linear discretized differential operator A, then apply backtracking interpolationon un,1∆t(un+1 − E(n)n+1un) = Aun+1. (5.52)76We can then obtain un+1 by solving the linear system,un+1 −∆tAun+1 = E(n)n+1un(I−∆tA)un+1 = E(n)n+1unAˆun+1 = E(n)n+1un,(5.53)where Aˆ = I−∆tA. We can then use the numerical solver GMRES to solve a solutionat each time step.This is a fully implicit treatment of the heat equation; we could instead treat theadvective terms explicitly. This would result in the SBDF1 or implicit-explicit Eulerscheme [1]; we discuss this more generally in the next section.Reaction-Diffusion System with SBDF2 Time-SteppingAs an example we show the time discretization of a reaction-diffusion of two agents ona moving domain Ω(t) using SBDF2 time stepping. The system of equations isut = d∇2Ω(t)u−∇ · (~s u) + f(u, v) (5.54)vt = ∇2Ω(t)v −∇ · (~s v) + g(u, v). (5.55)To apply the semi-implicit SBDF2 scheme on the embedded space, we treat the Laplace–Beltrami operator implicitly and the other terms explicitly, which yields12∆t(3ue(~x, tn+1)− 4ue(~x, tn) + ue(~x, tn−1))= d∇2ue(~x, tn+1) + 2fˆ(ue(~x, tn), ve(~x, tn))− fˆ(ue(~x, tn−1), ve(~x, tn−1))(5.56)12∆t(3ve(~x, tn+1)− 4ve(~x, tn) + ve(~x, tn−1))= ∇2ve(~x, tn+1) + 2gˆ(ue(~x, tn), ve(~x, tn))− gˆ(ue(~x, tn−1), ve(~x, tn−1)),(5.57)where ue, ve represent the closest point extension to the surface functions u and vrespectively. This explicit treatment of the diveregence terms will incur a Courant–Friedrichs–Lewy condition [14] for the time step size ∆t ≤ hmax(‖~s‖1) , which is not veryrestrictive for slow moving surfaces compared to ∆t ≤ O(h2) for the explicit treatmentof the Laplace–Beltrami operator.This, paired with the penalty term in the Laplacian operator as in (5.51) and paststeps interpolations yields12∆t(3un+1 − 4E(n)n+1un + E(n−1)n+1 un−1)= dLˆun+1 + 2fˆ(E(n)n+1un, E(n)n+1vn)− fˆ(E(n−1)n+1 un−1, E(n−1)n+1 vn−1)(5.58)12∆t(3vn+1 − 4E(n)n+1vn + E(n−1)n+1 vn−1)= Lˆvn+1 + 2gˆ(E(n)n+1un, E(n)n+1vn)− gˆ(E(n−1)n+1 un−1, E(n−1)n+1 vn−1),(5.59)77where un, vn represent the solution to u, v at the time step tn on the discretized gridspace, L represents the discretization of the Laplace–Beltrami operator and the hattedfunctions include the extra divergence operator terms, which we are treating explicitly,and E(i)j is the past step interpolation of time step i onto coordinates backtracked fromstep j. We then use the numerical solver GMRES to solve a solution at each time step.5.4 Convergence of the Heat EquationIn this section we have the result of a convergence analysis for the heat equationut = ∇2Ω(t)u−∇Ω(t) · (~su), (5.60)on a circle or a sphere undergoing different time dependent isotropic transformations.The different domain evolution functions are chosen such that exact solutions can becomputed, for example a circle undergoing an exponential scaling over time with relativerate of change ρ and pure mode initial condition U(0) = U0 cos θ, the heat equation hasthis exact solutionu = u0e−ρt+ e−2ρt−12R2ρ cos θ. (5.61)We can then compute the error by taking the difference between the simulation andexact solutions at each closest point and taking the maximum norm. We observe thatthe error converges to zero essentially quadratically in all tests. More details on thecomputation of the exact solution may be found in Appendix D.In this section we used a SBDF2 time-stepping scheme with a five-point (seven-pointfor the sphere) stencil (2.74) for the Laplacian operator discretization.We show here a series of tables describing a convergence analysis on the heat equationfor a few different evolutions.Table 5.1: ∆t = 0.2h, T = 10, Ri = 1, Rf = 2h Relative error Absolute error0.01 2.6950× 10−5 6.0252× 10−80.02 1.1396× 10−4 2.5480× 10−70.04 5.1260× 10−4 1.1465× 10−60.08 3.2686× 10−3 7.3311× 10−678Table 5.2: ∆t = 0.1h, T = 10, Ri = 1, Rf = 2h Relative error Absolute error0.01 2.9972× 10−5 6.7008× 10−80.02 1.3308× 10−4 2.9754× 10−70.04 6.4505× 10−4 1.4430× 10−6Table 5.3: ∆t = 0.05h, T = 10, Ri = 1, Rf = 2h Relative error Absolute error0.01 3.3827× 10−5 7.5625× 10−80.02 1.6239× 10−4 3.6310× 10−70.04 8.7219× 10−4 1.9516× 10−6Table 5.4: Uniform rotational motion x = eiωtx0, ∆t = 0.2h, T = 10, ω = 0.1,h Relative error Absolute error0.01 7.1022× 10−5 3.2242× 10−90.02 3.1193× 10−4 1.4164× 10−80.04 1.4148× 10−3 6.4313× 10−8Table 5.5: Exponential growth and uniform rotational motion x = ekt+iωtx0, ∆t = 0.2h,T = 10, ω = 0.1, Ri = 1, Rf = 2, Lagrangian approach.h Relative error Absolute error0.01 2.6152× 10−5 5.8460× 10−80.02 1.1064× 10−4 2.4734× 10−70.04 4.9922× 10−4 1.1165× 10−6Table 5.6: Exponential growth and uniform rotational motion x = ekt+iωtx0, ∆t = 0.2h,T = 10, ω = 0.1, Ri = 1, Rf = 2, Eulerian approach.h Relative error Absolute error0.01 2.6107× 10−5 5.8359× 10−80.02 1.1047× 10−4 2.4696× 10−70.04 4.9852× 10−4 1.1149× 10−6Table 5.7: Semicircle with homogeneous Dirichlet boundary conditions. ∆t = 0.2h,T = 10, Ri = 1, Rf = 2h Relative error Absolute error0.01 6.5217× 10−5 5.3564× 10−90.02 2.7851× 10−4 2.2880× 10−80.04 1.2877× 10−4 1.0589× 10−779Table 5.8: Exponentially growing sphere. ∆t = 0.2h, T = 10, Ri = 1, Rf = 2.h Relative error Absolute error0.04 2.3023× 10−3 1.1507× 10−80.08 9.6929× 10−3 4.8443× 10−80.16 4.1059× 10−2 2.0520× 10−7Table 5.9: Exponentially growing sphere. ∆t = 0.2h, T = 1, Ri = 1, Rf = 2.h Relative error Absolute error0.04 5.5087× 10−4 4.6674× 10−50.08 2.2395× 10−3 1.8975× 10−40.16 9.0355× 10−3 7.6555× 10−45.5 Schnakenberg RDE on Different Types of 2DDomain EvolutionIn this section we study a Schnakenberg reaction-diffusion equation on a circle usingexponential growth, logistic growth and then bipolar trajectories on a circle with expo-nential or logistical growth. Under the parameters we choose, there should be perioddoubling bifurcations occuring as the circle gets larger.In this section we used a SBDF2 time-stepping scheme with a five-point stencil (2.74)for the Laplacian operator discretization.The reaction functions for the Schnakenberg model aref(u, v) = a− u+ u2v, g(u, v) = b− u2v, (5.62a)and we chose the following parameter values for all simulationsa = 0.1, b = 0.9, d = 0.01. (5.62b)5.5.1 Exponential GrowthFor this trial we set time and space grid size h = 0.01, ∆t = 0.002, simulation timeT = 320 and the following parameter valuesR0 =12pi, ρ = 0.01. (5.63)We observe the bifurcations just like other similar simulations (such as [15], where thedomain growth is ten times slower) at similar times. Figures 5.4 and 5.5 show thesolutions on the band and as a function of the angle on the circle at different stages ofthe simulation. The vertical axis represents the arclength along the circle, starting atangle θ = −pi/2. Period doubling corresponds to the number of peaks doubling every80time the domain doubles in size. The first bifurcation takes twice as long because thesystem with homogeneous Neumann boundary conditions has a “half spike” pattern forsmall domains, where at one of the boundary points the solution is at the maximum ofa spike, and hence the zero-derivative condition is satisfied. These “half spike” patternscannot occur with periodic boundary conditions because we also need the solutions tomatch along with their derivative at the boundary points.Figure 5.3: Schnakenberg simulation on an exponentially growing circle. h = 0.01,∆t = 0.2 · h, ρ = 0.01. The arclength is measured counter-clockwise from the angleθ = −pi/2.5.5.2 Logistic GrowthFor the logistically growing domain we tried many configurations of maximal domainsize and time and space grid sizes, while keeping the initial radius and initial relativegrowth rate atR0 =12pi, ρ = 0.01. (5.64)We also set the total simulation time at T = 1400 to see if the final solution shapesare robust. To get results within reasonable times we had to increase the time step sizeconsiderably, which the implicit scheme should allow us to do. There is, however a limitto this increase, as we don’t want the band to change too much from one time step tothe next and effectively capture the action of the movement within the domain.Figure 5.6 shows a time-angle plot of the solution of a simulation with grid sizesh = 0.01, ∆t = 0.05 and relative maximum domain size ζ = 26. We again observe81(a) t = 16 (b) t = 144(c) t = 208 (d) t = 320Figure 5.4: Band solution of the Schnakenberg system on an exponentially growingcircle for four different time. Parameter values are the same used as in Figure 5.3.82(a) t = 16 (b) t = 144(c) t = 208 (d) t = 320Figure 5.5: Circle solutions of each snapshot in Figure 5.4, obtained by interpolatingthe band solution onto regularly angled points on the circle.83period doublings as the domain increases in size and as the domain stops increasing themaximal values of each peak is slightly reduced to a stable level. The times of transitioncorrespond very closely to results of Crampin et al. [15] on a constant domain computedusing time dependent diffusion coefficients, chosen to be equivalent to an exponentiallygrowing sphere.Figure 5.6: Schnakenberg Simulation on a logistically growing circle. h = 0.01, ∆t =0.05, ρ = 0.01, ζ = 26One should note that for isotropic evolution on a circle with periodic boundaryconditions any rotation of a solution along the circle also yields a solution. For randominitial conditions this will result in patterns that do not necessarily align with our choiceof reference angle (θ = −pi/2), whereas the equivalent solution on a growing line withhomogeneous Neumann boundary conditions, which does not have this symmetry, willnaturally align itself to the same solution each time.5.5.3 Exponential Growth Along Bipolar TrajectoriesIn this section we use an evolution function for the ξ parameter described in equations(5.20)–(5.21) such that each point on our domain goes through a bipolar trajectory (bluecurves in Figure 2.3), but the overall shape of the domain is the same as an exponentiallygrowing circle (red curves in Figure 2.3). We use the same initial radius and relativegrowth rate (5.63), but we start at time t0 = 0.1 to avoid a singularity at the ξ = pi/2value.We computed the divergence term using two different methods. Initially we usedan analytically precomputed form analoguous to (4.9) that yielded a term proportionalto the coordinate y and the second method used finite difference directional derivatives84applied to a product of the computed velocity field with the solution. Both methodsyield very similar results.We initially observe a clear anchoring behaviour with the solution aligning with theaxis of symmetry x = 0. This is due to the nonautonomous divergence term breakingrotational symmetry on our domain, but preserving reflexion symmetry on that axis.Then we see the transition of peaks with a higher y coordinate value happening earlierthan the lower peaks, which creates a six peak pattern around time 270. Figures 5.7and 5.8 show time plots of the solutions with the divergence term computed using thetwo different methods.Figure 5.7: Time plot of the solution on the circle as a function of the angle around itscentre with divergence term precomputed.5.5.4 Logistic Growth Along Bipolar TrajectoriesPairing bipolar trajectories with logistic total growth seems to be an efficient way tocontrol the final patterns. As in the previous exponential growth case, spikes closer to thetop of the sphere – where the domain is expanding more quickly – will split earlier, which,in this case, means that most splittings should occur there. In addition, the stabilizingevolution should create stable final patterns for pattern numbers that we cannot observewith uniform growth, provided we adjust the final domain size appropriately.Figures 5.9, 5.10 and 5.11 show numerical solutions for a decreasing final domainsize (from 26 to 20 to 16) and naturally, we see that the splittings happen later for thesmaller circle, until it vanishes for the domain sized 16, leading to a steady-state sixspiked pattern as the domain stabilizes.85Figure 5.8: Time plot of the solution on the circle as a function of the angle around itscentre with divergence term computed numerically.5.5.5 Stretching EllipseFor this simulation we started from a unit circle, then increased the major axis (x-axis)exponentially in time such that the final size is twice the initial size and decrease theminor axis (y-axis) also exponentially to three quarters of the initial size. For eachintermediate step the domain is an ellipse with known dimensions and the velocity fieldcan be computed directly component wise~s =(ρxxρyy), (5.65)and the divergence term is computed numerically, using the Eulerian backtracking met-hod.We first started by performing heat equation simulations on the stretching ellipsewith initial condition U(0) = cos(θ+pi/4), where the phase shift is to avoid symmetriesof the growth function. For comparison, we generated a reference solution using anon-embedded simulation using a finite difference method with an irregular grid∂ui∂x≈ ai+ui+1 + ai0ui + ai−ui−1, (5.66)86Figure 5.9: Schnakenberg Simulation on a logistically growing circle along bipolar tra-jectories. h = 0.01, ∆t = 0.05, ρ = 0.01, ζ = 26Figure 5.10: Schnakenberg Simulation on a logistically growing circle along bipolartrajectories. h = 0.01, ∆t = 0.05, ρ = 0.01, ζ = 2087Figure 5.11: Schnakenberg Simulation on a logistically growing circle along bipolartrajectories. h = 0.01, ∆t = 0.05, ρ = 0.01, ζ = 16withai+ =hi−hi+(hi+ + hi−)(5.67)ai0 =hi+ − hi−hi+hi−(5.68)ai− =−hi+hi−(hi+ + hi−), (5.69)where hi−, hi+ are the space step sizes to the left (clockwise) and right (counter clockwise)of the point xi on the unwound ellipse. The finite difference method uses 500 pointsthat are initially equidistant on the initial circle, the step size is ∆t = 2 · 10−6 and usesexplicit forward Euler. We then sampled 100 points for both solutions at the exact samepositions and compared their values by taking the relative error in infinity norm. Table5.10 shows convergence of both solutions at a superlinear rate, which is expected sincethe finite difference solution also contains a small error.We then performed simulations of the Schnakenberg system on the ellipse. Figure5.12 shows snapshots of the Schnakenberg reaction-diffusion system on an ellipse that isstretching in both the x and y directions, and Figure 5.13 shows snapshots of the samesystem for an ellipse that is growing similarly in the x direction, but slightly shrinks inthe y direction. The parameters of the Schnakenberg system are the same as equation(5.62b) and the total simulation time is Tf = 320 for both figures. We see that, especiallytowards the last peak splitting event (around t = 233 for the first, thicker ellipse and88∆x Relative error0.01 1.3748× 10−40.02 3.9788× 10−40.04 1.1640× 10−3Table 5.10: Comparison of heat equation solutions on a stretching ellipse with an initialcircle of radius one shape and final semi-major and semi-minor axes of lengths 2 and 0.75respectively. Other parameter values are ∆t = 0.2∆x, Tf = 3 and the initial conditionis u(0) = cos(θ + pi/4).t = 257 for the second, thinner one) peaks situated in areas that are stretched the most,near the y-axis split slightly before other peaks, which is to be expected, but overallsplittings still seem roughly simultaneous in both cases, with peaks doubling from 3 to6 to 12 peaks.(a) t = 75 (b) t = 175(c) t = 233 (d) t = 320Figure 5.12: Three snapshots of Schnakenberg solutions on a stretching ellipse. Initialshape is a circle of radius R = 1 and the axes are stretched by a factor of 8 for thesemi-major axis and a factor of 3 for the semi-minor axis.5.6 Numerical Examples in 3DIn order to demonstrate that our evolving domain closest point method works for higherdimension problems, we demonstrate here some experiments performed on a sphericalcap surface. The spherical cap is a section of a sphere when sliced by a plane. It is89(a) t = 75 (b) t = 175(c) t = 257 (d) t = 320Figure 5.13: Three snapshots of Schnakenberg solutions on a stretching ellipse. Initialshape is a circle of radius R = 1 and the axes are stretched by a factor of 8 for thesemi-major axis and a factor of 0.75 for the semi-minor axis.characterised by a circular boundary with radius R and a curvature index γ, given bythe sine of the angle the surface makes with the plane containing the circular boundary.Our main goal is to show that the method can handle the reaction-diffusion partialdifferential equation system that we previously used on a flattening spherical cap domainin Chapter 4. From this chapter we will use the predictions obtained through a centremanifold reduction for a spherical cap with constant radius R, and slowly decreasingcurvature γ(εt) for ε 1. In order to get numerical results in a reasonable amount oftime we need to apply the approximations discussed in the previous section.We used the SBDF2 time-stepping scheme with a seven-point stencil analogous toequation (2.74) for the Laplacian operator discretization.5.6.1 Schnakenberg SystemIn this section we repeat the Schnakenberg system (5.62a) computations on an evolvingsphere. We use the same parameter values (5.62b) as before.Exponential Evolution on a SphereFor this simulation we set space and time grid sizes h = 0.16, ∆t = 0.032, total simula-tion time T = 160, initial sphere radius R0 = 1 and relative growth speed of ρ = 0.01for the growing sphere and ρ = −0.01 for the shrinking sphere.90In the growing case, we observe the same type of peak splittings that was observedin the two-dimensional case as well as some three dimensional simulations. In theshrinking case we observe the selective decay of unstable peaks. This is to be expectedin the case of localized peak solutions, as was demonstrated by Rozada et al. [52] byincreasing or decreasing activator fuel (parameter a in the Schnakenberg model we use)on a Brusselator system, which is equivalent to increasing or decreasing the domain izerespectively. They also observed splittings on an effectively growing sphere computedusinf time-dependent diffusion coefficients. Overall the splittings and decays are not assymmetric as in the circle case, which might be due to the peaks needing more time tostabilize.(a) t = 121.6 (b) t = 133.76 (c) t = 150.4Figure 5.14: Three snapshots of Schnakenberg solutions on an exponentially growingsphere. We clearly see a splitting peak towards the centre of the figures. Initial radiusis R = 1 and the relative growing rate is ρ = 0.01.5.6.2 Brusselator System on a Linearly Growing SphereHere we study spike splitting events in the case of a sphere that is growing linearly. Weuse the following Brusselator reaction-diffusion equationsut = ε2∇2Ω(t)u−∇Ω(t) · (~su) + ε2E − u+ cu2v, (5.70a)κvt = d∇2Ω(t)v − κ∇Ω(t) · (~sv) +1ε2(u− u2v) (5.70b)We use the same methodology used by Rozada et al. [52] of using uniformly randominitial conditions, letting the stable pattern settle for 75 time units on a sphere of radiusR = 1. Then we grow the domain linearly over the next 50 time units, then keep thedomain radius at R = 2 from time t = 125 onwards. The same grid size ∆x = 0.025and initial conditions were used. The following parameter values were usedc = 0.8, ε = 0.12, d = 1.0, τ = 1/(0.8)2, E = 4.0. (5.71)91(a) t = 35.2 (b) t = 41.6 (c) t = 48.0Figure 5.15: Three snapshots of Schnakenberg solutions on an exponentially shrinkingsphere. We can see a decaying peak towards the centre of the figures. Initial radius isR = 5 and the relative growing rate is ρ = −0.01.(a) t = 105.0 (b) t = 115.5 (c) t = 126.0(d) t = 136.5 (e) t = 147.0 (f) t = 157.5Figure 5.16: Six snapshots of Schnakenberg solutions on an exponentially growing spherealong toroidal coordinates. We see more splittings at the top of the sphere, where thedomain is expanding the most, than at the bottom of the sphere, where the domain isrelatively constant. Initial radius is R = 1 and the relative growing rate is ρ = 0.01.92The main difference between the two experiments is that here we use an evolving domaininstead of using time-dependent parameters on a constant domain. We also included thedivergence term coming from the Reynolds transport theorem in the system. We alsoneeded to scale ε to 1/L in every term it is associated with, except for the Laplacian,in order to observe “peanut splittings”, as seen in [52].(a) t = 75.20 (b) t = 112.00 (c) t = 125.12(d) t = 150.40 (e) t = 184.96 (f) t = 236.16Figure 5.17: Six snapshots of Brusselator solutions on a sphere growing linearly betweent = 75 and t = 125, doubling its size and remaining constant outside the interval.Figure 5.17 shows the numerical solution of u with the divergence term present in thenumerical simulations. We observed the same two-fold “peanut” peak splittings from astable configuration of 4 peaks to another stable state of 8 peaks.5.6.3 Unevenly Stretching EllipsoidIn this section we extend the results on the ellipse obtained in Section 5.5.5 to an ellipsoidin three dimensions with the added feature that the domain evolution in the positivex-direction is different than in the opposite direction. We achieve this by computingclosest points on two ellipsoids corresponding to the two different evolution functions andonly keeping the valid coordinate values on either side of the yz-plane. Furthermore, we93kept the semi-major axes along y and z axes equal throughout the simulations. Figure5.18 shows the results.(a) t = 40 (b) t = 48(c) t = 56 (d) t = 64Figure 5.18: Snapshots of numerical solutions for the Schnakenberg system on an une-venly growing shape. We used a grid size of h = 0.04 at four time values. The initialconditions used are the homogeneous patternless solution with added uniformly randomnoise of amplitude 0.16 for u. We observe that peak splittings occur more often onthe right side (along the x-axis) of the domain, where the domain stretches faster thanon the left side. This is consistent with the previously obtained results on ellipses andspheres.By doing this we show that the closest point method on evolving domain may beadapted to domains with asymmetric evolution and shape and that qualitative charac-teristic of solutions are preserved locally. In this case we observe that peak splittingsoccur more often on the faster stretching positive x side of the domain than on theslower opposite side. This also shows that the method may be used on compositionsof connected subdomains, and could be optimized with better closest point functionsdefined on each subdomain separately.5.7 Brusselator System on a Flattening SphericalCapFinally we once again consider the computations performed in Section 4.3.5 on soluti-ons of the Brusselator system on a slowly flattening spherical cap. These simulationsconfirmed that we can skip some aspects of the growing domain (such as updating our94grids and matrix operators) on many time steps to improve calculation time and stillretain quadratic convergence, as long as the domain velocity is very small. In the speci-fic example shown in Figures 4.2, 4.3, 4.4 and 4.5 we ran the simulation over time stepsizes of ∆t = 0.1, but updated the domain every fifty steps. We kept the Reynold’sdivergence term in all time steps.Figure 5.19: Convergence analysis plot for the closest point method simulations. Thenumerical solutions are compared with a prediction generated from the centre manifoldequation. A dashed line corresponding to h2 times a constant is added to compare withquadratic convergence.For these numerical solutions, we developed a numerical surface integration methodin order to extract the (5, 1) mode part using an inner product projection method. Wefirst align the solution with the (5, 1) function on the bands, initially by using a rotationinterpolation operator to match each others maximum value position, then fine tunethis rotation to remove any residual phase between the functions using a binary searchalgorithm. We then performed the numerical integration on the spherical cap by using200 equidistant concentric circles each containing 200 equidistant sampling points. Wecomputed each circle integral with a midpoint rule sum then using another midpointrule sum on each circle integral weighted by the scale factor sin θ.By measuring the (5, 1) part of the numerical solution, we can then compare thenumerical solutions to a normal form reduction prediction calculated from a slowlyevolving centre manifold reduction in Chapter 4. Figure 4.4 shows the normal formcoefficient as a function of the curvature index γ as well as the equivalent coefficientscomputed from the numerical simulations for different grid sizes. We can see that thecurves corresponding to the numerical simulations seem to converge to the normal form95curve. For the final curvature value γ = 0.4515, we also computed a difference betweenthe numerical simulations and a centre manifold approximation solution. The L2 andmaximum norms of that difference are shown in Figure 5.19 along with a h2 line andconverge quadratically in both norms.96Chapter 6ConclusionsIn this thesis we presented several results relating to pattern formation on curved opensurfaces, especially on spherical caps, which can serve as a model for a plant tip. Wevalidated our findings with several numerical simulations, using a new adaptation of theclosest point method to evolving domains, which is then tested thoroughly and appliedto various other domain evolutions.We have first studied, in Chapter 3, bifurcating spatial patterns in a spherical capsurface, in a reaction-diffusion system that models morphogen concentrations in a de-veloping plant tip, such as a conifer embryo. In particular, we studied two cases ofmultiple bifurcations, where there are interactions between different patterning modes.In both cases there exist bifurcating nonlinear pure mode patterns, which are close tospherical cap harmonic functions, and nonlinear mixed mode patterns, which resemblespecific linear combinations of spherical cap harmonics. In what we call Case 1, wheretwo pure mode patterns have discrete rotational symmetries about the central axis ofthe spherical cap, the interaction is relatively simple to analyze, and for the parametervalues we considered the mixed modes are unstable and there is bistability of the puremodes. In Case 2, one pure mode pattern is annular and the interaction between themodes is more complicated to analyze. For the parameter values we considered, themixed mode patterns are stable. We have also computed finite element simulations thatsupport our analytical predictions.We began with a linear stability analysis of the unpatterned state, in terms of sp-herical cap harmonic functions. This was followed by reductions to center manifoldsand normal forms with symmetry. Finally the problems were reduced to the study ofamplitude equations. In Case 1, the amplitude equations were those of the classicaldouble pitchfork bifurcation with Z2⊕Z2-symmetry, which coincide with the amplitudeequations at a simple case of a Hopf-Hopf (or double Hopf) bifurcation. In Case 2 theamplitude equations can be thought of as those for a double pitchfork bifurcation withZ2 ⊕ Z2-symmetry, but perturbed by the addition of small symmetry-breaking termswhich retain only Z2-symmetry and a constraint that the origin remains an equilibrium.These amplitude equations coincide with those for a transcritical-Hopf bifurcation in aspecial case where the transcriticality is weak.As a model of a developing plant tip, our reaction-diffusion model shows that mor-phogen patterns that resemble observed initial growth patterns are stable. Pure modepatterns based on the (5, 1) eigenfunction and mixed mode patterns based on the (5, 1)and (3, 0) eigenfunctions both resemble the five-cotyledon configuration commonly seenin a conifer embryos [10]. Our bifurcation analysis shows that these patterns can be97stable in a somewhat wider parameter region than the one predicted from a single-modebifurcation analysis alone.In the study just mentioned, we have modeled pattern formation on a plant tipusing parameters such as tip curvature that are constant in time. In conifer embryos,patterning is observed to be preceded by tip flattening (corresponding to decreasing thetip curvature parameter), usually keeping the same tip radius R. Since the time scaleof pattern formation (on the order of hours or less) is much faster than the time scaleof tip flattening (over several days), as a first approximation a quasi-static approachto modeling seems reasonable, where one assumes the tip curvature parameter γ to beconstant while the pattern reaches a steady state.Next, in Chapter 4 we reduced the Brusselator reaction-diffusion system (2.8) ona time dependent non-isotropically slowly evolving spherical cap domain undergoing adelayed bifurcation. We did this using an asymptotic expansion of the normal formobtained from the constant (time independent) domain case. We use concepts frombifurcation theory such as centre manifold expansions and equivariant normal forms.This allows for the reaction-diffusion system to be approximated by a nonautonomousordinary differential equation. We obtained normal forms for m0 6= 0 modes and per-formed calculations using the m0 = 5, n0 = 1 mode. These reproduce a transition fromthe quasi-patternless state to the patterned state using significantly less computing re-sources compared to a numerical simulation of the partial differential equation system.This gives us a technique to better understand emergence of patterned solutions aroundcritical stability parameter values on an evolving domain.The existence and properties of exponentially attracting centre manifolds for systemsof nonautonomous partial differential equations are well established (e.g. [13]). Thecentre manifold expansion we use is from [49], adapted to slowly changing coefficients.It remains to prove, in future work, that our formal results using the WKB method(4.34) indeed imply the exponential dichotomy properties required for the existence ofa centre manifold. Instead of a proof we provide numerical results in Section 4.3.5 thatgive evidence toward its existence.There are, however some limitations to this method. First for the asymptotic ex-pansion to be valid we need slow domain evolution. Also, we need to start from thequasi-patternless state and go through a m0 6= 0 stability curve. More work would berequired to study the m0 = 0 case or transitions between patterned states. The trunca-tion of terms in the centre manifold expansion also limits the validity of predictions tosolutions close to the quasi-patternless solution. Finally, obtaining the order-ε term forthe cubic term C(εt, ε) of the normal form (4.77) could produce a better normal formapproximation and should be achievable through similar methods, but with considerablymore work.Other formulations for the curvature function may be used, for example γ(τ) =γ0 − ∆γ2 arctan τ , for −∞ < τ < ∞, that has a similar, almost linear, evolution nearτ = 0 and asymptotically approaches constant quantities when t approaches eitherinfinities. The latter formulation could be useful in future simulations where we wish to98give time for the patterns to develop in order to observe them, for example when usingquicker domain changes, and for investigating exponential dichotomies, which requireat least semi-infinite time intervals.Other future work should include more numerical simulations in order to study moretransitions between different modes or the quasi-patternless solution to other modes,such as the codimension two bifurcations studied in Chapter 3. In some cases, for regi-ons in the parametric portraits especially close to the origin, patterns evolve very slowlyand are only weakly stable, and could be overwhelmed by the effects of a slowly chan-ging parameter, typically including the phenomenon of delayed bifurcations, dependingon how slowly the parameter changes. However, we expect the more robust features ofthe codimension two bifurcations will persist when a slowly changing parameter is con-sidered. We could also account for different ways to change the domain. An even morecomplete and realistic model would have changes in the domain appearing as dynamicalconsequences of growth catalysts. These effects have been studied in simulations (e.g.[4], [16], [26], [28], [45], [57]) but much work remains to analytically characterize thepatterning mechanisms in this more realistic setting.In Chapter 5, we show that the Closest Point Method can be adapted to evolvingdomains, provided we have knowledge of the domain velocity vector and the appropriateclosest point function at every time step. We successfully showed that the method pre-serves pointwise quadratic convergence for the heat equation, as well as for the sphericalcap Brusselator model. We also demonstrated its use on many different evolving domainreaction-diffusion patterning systems to further validate the method and show resultson further non-isotropically evolving domains.Because it can work using only the velocity field this method could potentially beused for solution-dependent domain evolution, which are often used in morphogenesismodels. What we need to develop is a fast algorithm that could compute the closestpoint function of arbitrary curves and surfaces. The nonlocality of the closest pointfunctions do make these computations quite challenging, and the best alternative in theforeseeable future could be to compute the functions from scratch each time step.Another logical extension to this project would be to use it to simulate systems byadding a small thickness to the surface or working with domains with bulk in general.Simulations pairing bulk and surface reaction-diffusion systems have succeeded in re-laxing some of the activator-inhibitor restrictions for Turing patterns to emerge [37],as well as providing biological models for intracellular protein oscillations [25] and cellpolarization [17]. The Closest Point Method should be well suited for these types ofproblems as the algorithms can easily separate band points that are in the bulk andthose that have a closest point on the surface of the domain.Moreover, there are morphogenesis models that use a domain whose evolution isdictated by the solution of a similar system, but incorporating biomechanical concepts.For example Brinkmann et al. [8] have presented a tissue surface model in three di-mensions with small thickness and can show multiple morphogenesis patterns seen inmany organisms. In the example of the conifer embryo, for example, we could add some99thickness to the exterior tissue, or expect domain growth to be promoted in the regionsof maximal concentration of one of the chemical agents. Any such features added tothe model would make the asymptotic analysis much more complicated, but interestingnumerical experiments could be made, though the closest point function would have toadapt to the irregularly evolving domain.100Bibliography[1] Ascher, U., Ruuth, S., and Wetton, B. Implicit-explicit methods for time-dependent partial differential equations. SIAM Journal on Numerical Analysis 32,3 (1995), 797–823. 58, 73, 77[2] Auer, S., and Westermann, R. A semi-lagrangian closest point method fordeforming surfaces. Computer Graphics Forum 32, 7 (2013), 207–214. 2[3] Axelsson, O., and Barker, V. Finite Element Solution of Boundary ValueProblems. Classics in Applied Mathematics. Society for Industrial and AppliedMathematics, 2001. 33[4] Baker, R. E., and Maini, P. K. A mechanism for morphogen-controlled domaingrowth. Journal of Mathematical Biology 54, 5 (May 2007), 597–622. 99[5] Barreira, R., Elliott, C. M., and Madzvamuse, A. The surface finiteelement method for pattern formation on evolving biological surfaces. Journal ofMathematical Biology 63, 6 (Dec 2011), 1095–1119. 46[6] Berrut, J.-P., and Trefethen, L. N. Barycentric lagrange interpolation.SIAM Review 46, 3 (2004), 501–517. 20[7] Bilsborough, G. D., Runions, A., Barkoulas, M., Jenkins, H. W., Has-son, A., Galinha, C., Laufs, P., Hay, A., Prusinkiewicz, P., and Tsian-tis, M. Model for the regulation of Arabidopsis thaliana leaf margin development.Proceedings of the National Academy of Sciences 108, 8 (2011), 3424–3429. 1, 3[8] Brinkmann, F., Mercker, M., Richter, T., and Marciniak-Czochra,A. Post-Turing tissue pattern formation: Advent of mechanochemistry. PLOSComputational Biology 14, 7 (07 2018), 1–21. 99[9] Budrene, E. O., and Berg, H. C. Dynamics of formation of symmetricalpatterns by chemotactic bacteria. Nature 376, 6535 (Jul 06 1995), 49–53. 1[10] Butts, D., and Buchholz, J. T. Cotyledon numbers in conifers. Transactionsof the Illinois State Academy of Science 33, 2 (1940), 58–62. 34, 97[11] Carr, J. Applications of Centre Manifold Theory. Applied Mathematical Sciences.Springer New York, 1981. 24101[12] Chen, Y., Kolokolnikov, T., and Zhirov, D. Collective behaviour of largenumber of vortices in the plane. Proceedings of the Royal Society A 469, 2156(2013), 20130085. 1[13] Chicone, C., and Latushkin, Y. Center manifolds for infinite dimensionalnonautonomous differential equations. Journal of Differential Equations 141, 2(1997), 356 – 399. 98[14] Courant, R., Friedrichs, K., and Lewy, H. On the partial difference equa-tions of mathematical physics. IBM Journal of Research and Development 11, 2(March 1967), 215–234. 77[15] Crampin, E. J., Gaffney, E. A., and Maini, P. K. Reaction and diffusion ongrowing domains: Scenarios for robust pattern formation. Bulletin of MathematicalBiology 61, 6 (Nov 1999), 1093–1120. 80, 84[16] Crampin, E. J., Hackborn, W. W., and Maini, P. K. Pattern formation inreaction-diffusion models with nonuniform domain growth. Bulletin of Mathemati-cal Biology 64, 4 (Jul 2002), 747–769. 99[17] Cusseddu, D., Edelstein-Keshet, L., Mackenzie, J., Portet, S., andMadzvamuse, A. A coupled bulk-surface model for cell polarisation. Journal ofTheoretical Biology (2018). 99[18] Eftimie, R., Perez, M., and Buono, P.-L. Pattern formation in a nonlocalmathematical model for the multiple roles of the tgf-b pathway in tumour dynamics.Mathematical Biosciences 289 (2017), 96 – 115. 1[19] Erneux, T., Reiss, E. L., Holden, L. J., and Georgiou, M. Slow passagethrough bifurcation and limit points. asymptotic theory and applications. In Dyn-amic Bifurcations (Berlin, Heidelberg, 1991), E. Benoˆıt, Ed., vol. 1493 of LectureNotes in Mathematics, Springer Berlin Heidelberg, pp. 14–28. 51[20] Fujita, H., Toyokura, K., Okada, K., and Kawaguchi, M. Reaction-diffusion pattern in shoot apical meristem of plants. PLoS ONE 6, 3 (03 2011),1–13. 1, 3[21] Garzo´n-Alvarado, D. A., Martinez, A. M. R., and Segrera, D. L. L.Appearance and formation of seed and pericarp may be explained by a re-action–diffusion mechanism? a mathematical modeling. Mathematical and Com-puter Modelling 55, 3 (2012), 853 – 860. 46[22] Golubitsky, M., and Schaeffer, D. M. Singularities and Groups in Bifurca-tion Theory. Applied Mathematical Sciences. Springer-Verlag, 1985. 27102[23] Guckenheimer, J., and Holmes, P. Nonlinear Oscillations, Dynamical Sys-tems, and Bifurcations of Vector Fields. Applied Mathematical Sciences. SpringerNew York, 2002. 13, 14[24] Haberman, R. Slowly varying jump and transition phenomena associated withalgebraic bifurcation problems. SIAM Journal on Applied Mathematics 37, 1 (1979),69–106. 18[25] Halatek, J., and Frey, E. Highly canalized mind transfer and mine sequestra-tion explain the origin of robust mincde-protein dynamics. Cell Reports 1, 6 (2012),741 – 752. 99[26] Harrison, L. G., Adams, R. J., and Holloway, D. M. Dynamic regulation ofgrowing domains for elongating and branching morphogenesis in plants. Biosystems109, 3 (2012), 488 – 497. Biological Morphogenesis: Theory and Computation. 99[27] Harrison, L. G., and von Aderkas, P. Spatially quantitative control of thenumber of cotyledons in a clonal population of somatic embryos of hybrid larchLarix x leptoeuropaea. Annals of Botany 93, 4 (2004), 423–434. 3[28] Holloway, D. M., and Harrison, L. Pattern selection in plants: Couplingchemical dynamics to surface growth in three dimensions. Annals of Botany 102,1 (2008), 361 – 374. 1, 33, 99[29] Holloway, D. M., Rozada, I., and Bray, J. J. H. Two-stage patterning dy-namics in conifer cotyledon whorl morphogenesis. Annals of Botany 121, 3 (2018),525–534. 1[30] Hong, Y., Zhu, D., Qiu, X., and Wang, Z. Geometry-based control of firesimulation. The Visual Computer 26, 9 (Sep 2010), 1217–1228. 2[31] Kim, T., Tessendorf, J., and Thu¨rey, N. Closest point turbulence for liquidsurfaces. ACM Transactions on Graphics 32, 2 (Apr. 2013), 15:1–15:13. 2[32] Kuznetsov, Y. A. Elements of Applied Bifurcation Theory. Applied Mathema-tical Sciences. Springer, 2010. 14, 27[33] Langford, W. F. Periodic and steady-state mode interactions lead to tori. SIAMJournal on Applied Mathematics 37, 1 (1979), 22–48. 30[34] Lobry, C. Dynamic bifurcations. In Dynamic Bifurcations, E. Benoˆıt, Ed.,vol. 1493 of Lecture Notes in Mathematics. Springer Berlin Heidelberg, 1991, pp. 1–13. 16[35] Macdonald, C. B., Brandman, J., and Ruuth, S. J. Solving eigenvalue pro-blems on curved surfaces using the closest point method. Journal of ComputationalPhysics 230 (2011), 7944–7956. 20, 60103[36] Macdonald, C. B., and Ruuth, S. J. The implicit Closest Point Method forthe numerical solution of partial differential equations on surfaces. SIAM Journalon Scientific Computing 31, 6 (2009), 4330–4350. 20, 76[37] Madzvamuse, A., and Chung, A. H. W. The bulk-surface finite elementmethod for reaction-diffusion systems on stationary volumes. Finite Elements inAnalysis and Design 108 (2016), 9 – 21. 99[38] Madzvamuse, A., Gaffney, E. A., and Maini, P. K. Stability analysis of non-autonomous reaction-diffusion systems: The effects of growing domains. Journalof Mathematical Biology 61, 1 (2010), 133–164. 46, 64[39] Ma¨rz, T., and Macdonald, C. Calculus on surfaces with general closest pointfunctions. SIAM Journal on Numerical Analysis 50, 6 (2012), 3303–3328. 19, 72[40] Meinhardt, H., Prusinkiewicz, P., and Fowler, D. The Algorithmic Beautyof Sea Shells. Virtual laboratory. Springer, 2003. 1[41] Miura, T., Shiota, K., Morriss-Kay, G., and Maini, P. K. Mixed-modepattern in doublefoot mutant mouse limb—turing reaction–diffusion model on agrowing domain during limb development. Journal of Theoretical Biology 240, 4(2006), 562 – 573. 46[42] Murray, J. D. Mathematical Biology. II, third ed., vol. 18 of InterdisciplinaryApplied Mathematics. Springer-Verlag, New York, 2003. Spatial models and bio-medical applications. 7[43] Nagata, W., Zangeneh, H. R. Z., and Holloway, D. M. Reaction-diffusionpatterns in plant tip morphogenesis: Bifurcations on spherical caps. Bulletin ofMathematical Biology 75, 12 (2013), 2346–2371. v, 1, 3, 4, 33[44] Narisawa, S. An implicit closest point method for solving convection-diffusionequations on a moving surface. PhD thesis, National Chiao Tung University, Hsin-chu, Taiwan, 2013. 2[45] Neville, A. A., Matthews, P. C., and Byrne, H. M. Interactions betweenpattern formation and domain growth. Bulletin of Mathematical Biology 68, 8 (Nov2006), 1975–2003. 99[46] Palin, R. J. A comparison of cell wall properties of Arabidopsis thaliana. PhDthesis, University of Birmingham, July 2011. 3[47] Petras, A. The closest point method for the numerical solution of partial differen-tial equations on moving surfaces. PhD thesis, Simon Fraser University, Burnaby,Canada, 2016. 64104[48] Plaza, R., Sa´nchez-Gardun˜o, F., Padilla, P., Barrio, R., and Maini,P. The effect of growth and curvature on pattern formation. Journal of Dynamicsand Differential Equations 16, 4 (2004), 1093–1121. 46, 47, 64, 66[49] Po¨tzsche, C., and Rasmussen, M. Taylor approximation of integral manifolds.Journal of Dynamics and Differential Equations 18, 2 (Apr 2006), 427–460. 51, 55,98[50] Prigogine, I., and Lefever, R. Symmetry breaking instabilities in dissipativesystems. ii. The Journal of Chemical Physics 48, 4 (1968), 1695–1700. 1[51] Reynolds, O., Brightmore, A. W., and Moorby, W. H. Papers on mecha-nical and physical subjects, vol. 3: The Sub-Mechanics of the Universe. Cambridge:The University Press, 1903. 46[52] Rozada, I., Ruuth, S., and Ward, M. The stability of localized spot patternsfor the Brusselator on the sphere. SIAM Journal on Applied Dynamical Systems13, 1 (2014), 564–627. 2, 64, 69, 91, 93[53] Ruuth, S., and Merriman, B. A simple embedding method for solving partialdifferential equations on surfaces. Journal of Computational Physics 227 (01 2008),1943–1961. 2, 19, 58, 61, 75[54] Schnakenberg, J. Simple chemical reaction systems with limit cycle behaviour.Journal of Theoretical Biology 81, 3 (1979), 389 – 400. 1[55] Stannard, A. Dewetting-mediated pattern formation in nanoparticle assemblies.Journal of Physics: Condensed Matter 23, 8 (2011), 083001. 1[56] Turing, A. M. The chemical basis of morphogenesis. Philosophical Transactionsof the Royal Society of London. Series B, Biological Sciences 237, 641 (1952), 37–72.1, 7, 69[57] van Mourik, S., Kaufmann, K., van Dijk, A. D. J., Angenent, G. C.,Merks, R. M. H., and Molenaar, J. Simulation of organ patterning on thefloral meristem using a polar auxin transport model. PLoS ONE 7, 1 (01 2012),1–9. 1, 3, 99[58] von Aderkas, P. In vitro phenotypic variation in larch cotyledon number. In-ternational Journal of Plant Sciences 163, 2 (2002), 301–307. 23[59] von Glehn, I., Ma¨rz, T., and Macdonald, C. B. An embedded method-of-lines approach to solving partial differential equations on surfaces. arXiv preprintarXiv:1307.5657 (07 2013). 19, 20, 76[60] Wen, A. Pattern Formation on Growing Domains. Master’s thesis, University ofOxford, Oxford, UK, 2011. 2, 64105Appendix AComputation of Normal FormCoefficients for Codimension 2Bifurcations in a Constant SphericalCapIn this appendix we show some details of the computation of the first several Taylorseries coefficients for the centre manifolds and normal forms used in Chapter 3. Todetermine stability and bifurcation of equilibria of the amplitude equations, it sufficesto evaluate these coefficients at the critical parameter values that give the double zeroeigenvalues σ1 = σ2 = 0 for the linearization of the Brusselator about the patternlesssteady state solution.For any two, possibly complex, vector functions Uj = (Uj(θ, ϕ), Uj(θ, ϕ))ᵀ, j = 1, 2,we define their inner product as〈U1,U2〉 =∫ 2pi0∫ θmax0[U1(θ, ϕ)U2(θ, ϕ) + V1(θ, ϕ)V2(θ, ϕ)]sin θ dθ dϕ, (A.1)where an overline denotes complex conjugation. Corresponding to the critical eigen-functions (3.1), we find adjoint eigenfunctions in the formU(∗)j =(u(∗)mj ,njv(∗)mj ,nj)eimjϕP|mj |λmj,nj (γ)(cos θ), j = 1, 2, (A.2)where the coefficients u(∗)mj ,nj , v(∗)mj ,nj solve the adjoint algebraic eigenvalue problem forσ = 0, (−DXµmj ,nj + k1 k3k2 −DY µmj ,nj + k4)(u(∗)mj ,njv(∗)mj ,nj)=(00). (A.3)The coefficients are normalized asu(∗)mj = k3/N(∗)j , v(∗)mj= k5,j/N(∗)j ,N(∗)j = spi√Mj(k2k3 + k25,j)/N(0)j ,(A.4)106so that〈U(0)j ,U(∗)k 〉 = δjk, j, k = 1, 2, (A.5)where δjk is the Kronecker delta. This allows us to define projections onto the centerand stable subspaces.A.1 Case 1: m1 6= 0 and m2 6= 0In Case 1, we define the projection Pc onto the centre subspace Ec, of any real vectorfunction U = (U(θ, ϕ), V (θ, ϕ))ᵀ that satisfies the homogeneous Dirichlet boundaryconditions (2.37), byPcU = z1U(0)1 + z2U(0)2 + z¯1U(0)1 + z¯2U(0)2 , (A.6)wherez1 = 〈U,U(∗)1 〉 ∈ C, z2 = 〈U,U(∗)2 〉 ∈ C. (A.7)Since the stable subspace Es is complementary to Ec, the projection onto Es is I−Pc,where I is the identity operator.Now for U in the centre manifold W c we have the Taylor series expansionU =z1U(0)1 + z2U(0)2 + z¯1U¯(0)1 + z¯2U¯(0)2 + z21U2000 + +z1z2U1100,+ z1z¯1U1010 + z1z¯2U1001 + z22U0002 + z2z¯1U0110 + z2z¯2U0011+ z¯21U0020 + z¯1z¯2U0011 + z¯22U0002 +O(‖(z1, z2, z¯1, z¯2)‖3),(A.8)with the Ujklm belonging to the stable subspace Es. We use projections Pc and I−Pc tosplit the system (2.36) onto its components in the centre and stable subspaces, substitutethe Taylor series expansion (A.8), and use the fact that W cloc is locally invariant. Thencollecting powers of (z1, z2, z¯1, z¯2) in the projection of the system (2.36) onto the stablesubspace Es, at second order we obtainA0U2000 = −(I−Pc)B0(U(0)1 ,U(0)1 ),A0U1100 = −2(I−Pc)B0(U(0)1 ,U(0)2 ),(A.9)and eight more similar equations for the remaining Ujklm. These are linear nonhomo-geneous systems of partial differential equations, with homogeneous Dirichlet boundaryconditions, that can be solved for the vector functions Ujklm by expanding the righthand sides in infinite series of the spherical cap harmonics. In practice, we use trunca-ted finite series expansions (eight, twelve and sixteen terms and comparing results toassess numerical convergence) to obtain approximations for the vector functions Ujklm.We note that not all ten vector functions at second order are needed to find the requirednormal form coefficients Cjk.107The normal form coefficients are found by substituting (A.8) into the projection ofthe system (2.36) onto the centre subspace Ec. We obtain the formulasC11 =〈2B0(U2000,U(0)1 ) + 2B0(U1010,U(0)1 ) + 3C0(U(0)1 ,U(0)1 ,U(0)1 ),U(∗)1 〉,C12 =〈2B0(U1100,U(0)2 ) + 2B0(U1001,U(0)2 ) + 2B0(U0101,U(0)1 ) + 3C0(U(0)1 ,U(0)2 ,U(0)2 ),U(∗)1 〉,C21 =〈2B0(U1100,U(0)1 ) + 2B0(U0110,U(0)1 ) + 2B0(U1010,U(0)2 ) + 3C0(U(0)1 ,U(0)2 ,U(0)1 ),U(∗)2 〉,C22 =〈2B0(U0200,U(0)2 ) + 2B0(U0101,U(0)2 ) + 3C0(U(0)2 ,U(0)2 ,U(0)2 ),U(∗)2 〉.(A.10)For the critical parameter values (2.52), we use the mathematical software package Mapleto evaluate the formulas for the normal form coefficients, obtaining the numerical values(3.13) given in Section 3A.2 Case 2: m1 6= 0 and m2 = 0In Case 2, the eigenfunction U(0)1 is complex while U(0)2 is real, and we define theprojection Pc onto the centre subspace Ec byPcU = zU(0)1 + z¯U(0)1 + xU(0)2 , (A.11)wherez = 〈U,U(∗)1 〉 ∈ C, x = 〈U,U(∗)2 〉 ∈ R. (A.12)As in Case 1, we expand U in W cloc in a Taylor series, but now with the formU =zU(0)1 + z¯U(0)1 + xU(0)2 + z2U200 + zz¯U110 + zxU101 (A.13)+ z¯2U020 + z¯xU011 + x2U002 +O(‖z, z¯, x‖3). (A.14)Following the same procedure as in Case 1, we obtain a linear nonhomogeneous systemfor each Ujkl, solved the same way as in Case 1.We have formulas for the quadratic normal form coefficientsB11 = 〈2B(U(0)1 ,U(0)2 ),U(∗)1 〉,B21 = 〈2B(U(0)1 ,U(0)1 ),U(∗)2 〉,B22 = 〈2B(U(0)2 ,U(0)2 ),U(∗)2 〉.(A.15)These do not involve any of the second order coefficients Ujkl from the centre manifoldexpression. For critical parameter values (2.52) for the intersection of the (m1, n1) =108(5, 1) and (m2, n2) = (0, 3) marginal stability curves, we numerically evaluate the for-mulas for the quadratic coefficients Bjk to obtain the numerical values (3.19). For thecubic normal form coefficient, we obtainC11 = 〈2B(U200,U(0)1 ) + 2B(U110,U(0)1 ) + 3C(U(0)1 ,U(0)1 ,U(0)1 ),U(∗)1 〉,C12 = 〈2B(U002,U(0)1 ) + 2B(U101,U(0)2 ) + 3C(U(0)1 ,U(0)2 ,U(0)2 ),U(∗)1 〉,C21 = 〈2B(U101,U(0)1 ) + 2B(U011,U(0)1 ) + 2B(U110,U(0)2 ) + 6C(U(0)1 ,U(0)1 ,U(0)2 ),U(∗)2 〉,C22 = 〈2B(U002,U(0)2 ) + C(U(0)2 ,U(0)2 ,U(0)2 ),U(∗)2 〉,(A.16)and evaluating these formulas at the critical parameter values (2.52) gives the numericalvalues (3.20).109Appendix BComputation of Bifurcation Curves,for Codimension 2 Bifurcations in aConstant Spherical Cap, Case 2In this appendix we give some details of the analysis of the amplitude equations (3.17),in Case 2r˙ = σ1r +B11xr + C11r3 + C12rx2x˙ = σ2x+B21r2 +B22x2 + C21r2x+ C22x3.(B.1)The first equation is odd in r, and the x-axis r = 0 is invariant under the flow. Weobserve that the values of the quadratic coefficients are much smaller than their cubiccounterparts, thus the amplitude equations can be viewed as a small perturbation ofthe double pitchfork amplitude equations (3.11) seen in Case 1. Sufficiently far from theorigin, the (σ1, σ2)-parametric portrait in Case 2 resembles one for a double pitchforkbifurcation with Z2⊕Z2-symmetry. Closer to the origin in the (σ1, σ2)-plane, the effectsof the small symmetry-breaking terms become more apparent.B.1 The Quadratic Normal FormIf we truncate the Taylor series of the reduced system (3.14) at quadratic order, theresulting amplitude equations arer˙ = σ1r +B11xr,x˙ = σ2x+B21r2 +B22x2,(B.2)with truncated remainders O(‖(r, x)‖3). There is a pitchfork bifurcation of the origin(r, x) = (0, 0) and mixed modes for σ1 = 0, and a transcritical bifurcation of the originand non-trivial x-modes for σ2 = 0. In addition, there are pitchfork bifurcations ofnontrivial x-modes and mixed modes, for parameters (σ1, σ2) on the line B22σ1 = B11σ2.The parametric portrait and phase portraits for (B.2) are valid for the reduced system(3.14) in neighbourhoods of the origin in the (σ1, σ2)-plane and the (r, x)-half-plane,respectively. However, these neighbourhoods turn out to be very small (see Figure3.6). Because the cubic coefficients are much larger than the quadratic coefficients, theyinfluence the dynamics, even moderately near the origin.110B.1.1 Bifurcations of the OriginThe origin (r, x) = (0, 0) is an equilibrium for the amplitude equations (3.17) for all(σ1, σ2). Investigating the linearized stability and bifurcations of the origin, we findthat the origin is asymptotically stable for parameters in the quadrant σ1 < 0, σ2 < 0,that the σ1-axis (σ2 = 0) corresponds to transcritical bifurcations of x-modes, and thatthe σ2-axis (σ1 = 0) corresponds to pitchfork bifurcations of mixed modes.B.1.2 Bifurcations of x-modesSince the x-axis is invariant for (3.17), setting r = 0 we obtainx˙ = σ2x+B22x2 + C22x3 = x(σ2 +B22x+ C22x2), (B.3)on the x-axis, thus we see again the transcritical bifurcation of x-modes at r = 0, x = 0when σ2 = 0. Nontrivial x-modes exist for σ2 +B22x+C22x2 = 0, and a fold bifurcationof x-modes occurs at r = 0, x = −B22/(2C22) when σ2 = B222/(4C22) = −1.0685× 10−5.These two bifurcations are arranged as parts of a perturbed pitchfork bifurcation, suchas shown in Figure 3.10a. The x-mode perturbation line σ2 = −1.0685 × 10−5 can beseen in the parametric portrait in Figure 3.5b, where it appears as a dash-dotted cyanhorizontal line. In Figure 3.5a, this fold bifurcation line is too close to the transcriticalbifurcation line σ2 = 0 for any difference between them to be visible at this larger scale.When we study the linearized stability of nontrivial x-modes (0, x), x 6= 0, we findthat the linearization of (3.17) is diagonal, with eigenvalues corresponding to the r- andz-directions. The eigenvalue in the x-direction is the standard one corresponding toequilibria of the one-dimensional equation (B.3). The eigenvalue in the r-direction isσ1+B11x+C12x2, and when this eigenvalue changes sign it signals a pitchfork bifurcationof a mixed mode (r, x), r > 0, from a nontrivial x-mode. Thus to detect these pitchforkbifurcations we simultaneously require that nontrivial x-modes exist and their eigenvaluein the r-direction is zero:σ2 +B22x+ C22x2 = 0,σ1 +B11x+ C12x2 = 0.(B.4)This system of equations (B.4) can be viewed as the parametric equations for a curvein the (σ1, σ2)-plane with parameter x. This dashed red line appears in Figures 3.5a-3.5c. The curve passes through the origin, and intersects the x-mode fold bifurcationcurve with a quadratic tangency. We verify that pitchfork bifurcations indeed occur forparameter values along this curve, and check the linearized stability of the bifurcatingmixed modes.111B.1.3 Bifurcations of Mixed ModesMixed mode equilibria of the amplitude equations (3.17) bifurcate from the origin inpitchfork bifurcations, and in general can be found as solutions (r, x) of the systemσ1 +B11x+ C11r2 + C12x2 = 0,σ2x+B21r2 +B22x2 + C21r2x+ C22x3 = 0,(B.5)with r > 0. We can solve the first equation for r2 and substitute into the second equationto obtainf0 + f1x+ f2x2 + f3x3 = 0, (B.6)withf0 =−B21C11σ1,f1 = σ2 − C21C11σ1 − B21B11C11,f2 = B22 − C21B11 +B21C12C11,f3 = C22 − C12C21C11.(B.7)This can be viewed as an unfolding of a cusp singularity. To locate the two curves offold singularities that meet at the cusp, we take the derivative with respect to x andsolve the systemf0 + f1x+ f2x2 + f3x3 = 0,f1 + 2f2x+ 3f3x2 = 0.(B.8)Now we can express σ1 and σ2 as functions of x and use these parametric equations totrace curves in the (σ1, σ2)-plane. However, we need to recall that r2 > 0, so with thisconstraint only portions of the fold bifurcation curves survive, seen as the dotted bluecurves in Figures 3.5a-3.5c. Solving for r2 = 0 on these curves shows that the endpointsof the fold bifurcation curves fall on one of the pitchfork bifurcation curves.112Appendix CTime Derivative of theLaplace–Beltrami EigenfunctionIn the expansion of the WKB solution (4.35) we are required to use the time derivativeof the harmonic function Φmn(τ) in (4.45). As the domain evolves with time so do theLegendre functions associated to the surfaces. This means that any point on the surfaceevolving in with velocity vector v will have a slightly different value for the functionΦmn. This is the time derivative we need to use.Using domain evolution along toroidal coordinates, this corresponds to an evolvingξ while keeping η and φ constant. We can express Φmn asΦmn(τ) = eimφP|m|λmn(τ)(ζ(τ)), (C.1)which will be useful when using the symmetries of the equation. The Legendre functionargument ζ(τ) is expressed asζ(τ) =1− cosh η cos ξ(τ)cosh η − cos ξ(τ) , (C.2)in toroidal coordinates. We can then find the time derivative using the chain ruleΦ′mn(τ) = eimφ[∂∂λP|m|λmn(τ)(ζ(τ))∂λmn(τ)∂τ+∂∂ζP|m|λmn(τ)(ζ)∂ζ(τ)∂τ]. (C.3)The derivative is thus split in two terms: the λ derivative and the ζ derivative. The firstterm does not have a practical expression that can be easily implemented in Maple, sowe used finite difference approximations to find the derivative∂∂λP|m|λ (ζ(τ)) ≈P|m|λ+h1(ζ(τ))− P |m|λ−h1(ζ(τ))2h1(C.4)∂λmn(τ)∂τ≈ λmn(τ + h2)− λmn(τ − h2)2h2, (C.5)for h1, h2 fixed, positive and very small. The derivatives were then tested by graphicallycomputing linear approximations to validate their values. Note that we were alreadyusing a numerical solver to find the λ values.113The ζ derivative term does have a useful identity∂∂ζP|m|λ =λζP|m|λ (ζ)− (λ+m)P |m|λ−1(ζ)ζ2 − 1 (C.6)∂ζ∂τ=sin ξ(τ) sinh2 η(cosh η − cos ξ(τ))2 · ξ′(τ) (C.7)=−γ′(τ)γ(τ)√1− γ(τ)2 ·γ(τ)2 sinh2 η(cosh η +√1− γ(τ)2)2 , (C.8)where we used the identitiesξ′ =γ′cos ξ, andcos ξ = −√1− γ2.(C.9)We can then use inner products of the new functions in order to find the series expression(4.45) of Φ′mn.114Appendix DExact Solutions of the HeatEquation on an Evolving DomainIn this appendix we demonstrate the analytical exact solutions of the heat equation onexponentially growing domains used in Section 5.4.D.1 Exponentially Growing CircleLet Ω(t) be an exponentially growing circle centered at the origin with radius functionR(t) = R0eρt. (D.1)As we have seen before the heat equation for this evolving domain isut = ∇2Ω(t)u− ρu, (D.2)where the Laplacian operator is∇2Ω(t)u =uθθR(t)2. (D.3)The Laplacian operator will have, up to arbitrary rotations, an orthogonal basis ofeigenfunctionsΦm(θ) = cosmθ, (D.4)for m any non-negative integer. Taking our initial condition to be any one of theseeigenfunctionsu(0, θ) = cosmθ, (D.5)the effect of the Laplacian is given by∇Ω(0)u(0, θ) = −m2R20u(0, θ). (D.6)Separating the variablesu(t, θ) = u˜(t) cos(mθ), (D.7)115we obtainu˜t cos(mθ) =−m2R(t)2u˜ cos(mθ)− ρu˜ cos(mθ), (D.8)which leads to a non-autonomous linear ODE for u˜u˜t =(−m2R(t)2− ρ)u˜ =(−m2e−2ρtR20− ρ)u˜, u˜(0) = 1, (D.9)which has the solutionu˜ = ee−2kt−12R20k−kt. (D.10)D.2 Exponentially Growing SphereFor the exponentially growing circle of radius R(t) the heat equation isut = ∇2Ω(t)u− 2ρu, (D.11)with Laplacian operator∇2Ω(t)u =1R(t)2 sin θ∂∂θ(sin θ∂u∂θ)+1R2 sin2 θ∂2u∂ϕ2, (D.12)where θ and ϕ are the polar and longitudinal angles respectively.Spherical harmonic functions Y ml (θ, ϕ) for m integer and l > |m| are eigenfunctionsof the Laplacian operator and create an orthogonal basis, the simplest of which wouldbe m = 0, l = 1Y 01 (θ, φ) = cos θ. (D.13)If we choose this function as our initial condition u0(θ), then its Laplacian is∇Ω(0)u0(θ) = −2R20u0(θ). (D.14)Using separation of variables as before we obtain the following ODE for u˜u˜t = −2(e−2ρtR20+ ρ)u˜. (D.15)Its solution beingu˜ = ee−2ρt−1ρR20−2ρt. (D.16)116
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Pattern formation on curved surfaces
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Pattern formation on curved surfaces Charette, Laurent 2019
pdf
Page Metadata
Item Metadata
Title | Pattern formation on curved surfaces |
Creator |
Charette, Laurent |
Publisher | University of British Columbia |
Date Issued | 2019 |
Description | Patterns emerge in various growing biological organisms, like conifer embryos, often from an homogeneous preceding state. The embryo tip, initially hemispherical shaped gradually flattens as cotyledons arranged in a roughly regular pattern emerge. A common way to model these patterns is through reaction-diffusion systems of partial differential equations. This thesis relays results obtained studying such systems and provides various new results about pattern formation on curved surfaces via diffusion driven bifurcations. We first describe situations where, for certain critical values of domain or differential equation parameters, the unpatterned solution is unstable to two different linear normal modes. We use centre manifold and normal form reductions to analyze the existence and stability of pure and mixed modes of nonlinear patterned solutions of the reaction-diffusion system, for parameters near two cases of critical values. In one case, the system reduces to a well known example of mode interaction. In the other case, the mode interaction is new, due to very small quadratic terms in the normal form. We then perform a reduction of a nonautonomous Brusselator reaction-diffusion system of partial differential equations on a spherical cap with time dependent curvature using an asymptotic series expansion on the centre manifold reduction. Parameter values are chosen such that the change in curvature would cross critical values which would change the stability of the patternless solution in the constant domain case. The non-isotropic nature of the domain evolution insert a small patterned component to the previously patternless state, which we call the `quasi-patternless solution'. The evolving domain functions and quasi-patternless solutions are derived as well as a method to obtain this nonautonomous normal form. The obtained reduction solutions are then compared to numerical solutions. Finally we provide an adaptation of the closet point method to evolving domains. We perform several convergence analysis experiments of the heat equation on test surfaces and obtain quadratic convergence. Then, a few reaction diffusion simulations are shown on evolving surfaces, some featuring non-isotropic evolution. The previously shown application of the Brusselator system on an evolving spherical cap are compared to a centre manifold reduction and also show quadratic convergence. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2019-03-28 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0377645 |
URI | http://hdl.handle.net/2429/69326 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2019-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2019_may_charette_laurent.pdf [ 7.94MB ]
- Metadata
- JSON: 24-1.0377645.json
- JSON-LD: 24-1.0377645-ld.json
- RDF/XML (Pretty): 24-1.0377645-rdf.xml
- RDF/JSON: 24-1.0377645-rdf.json
- Turtle: 24-1.0377645-turtle.txt
- N-Triples: 24-1.0377645-rdf-ntriples.txt
- Original Record: 24-1.0377645-source.json
- Full Text
- 24-1.0377645-fulltext.txt
- Citation
- 24-1.0377645.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
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"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0377645/manifest