Lecture 7 Notes: Linear Stability Analysis and Pattern Formation

18.354J Nonlinear Dynamics II: Continuum Systems 7 7.1 Lecture 7 Linear stability analysis and pattern formation Linear stability analysis of fixed points for ODEs Consider a particle (e.g., bacterium) moving in one-dimension with velocity v(t), governed by the nonlinear ODE d v(t) = −(α + βv 2 )v =: f (v). dt (148) We assume that the parameter β is strictly positive, but allow α to be either positive or negative. The fixed points of Eq. (148) are, by definition, velocity values v∗ that satisfy the condition f (v∗ ) = 0. For α > 0, there exists only p one fixed points v0 = 0. For α < 0, we find the three fixed points v0 = 0 and v± = ± −α/β. That is, the system undergoes pitchfork bifurcation at the critical parameter value α = 0. To evaluate the stability of a fixed points v∗ , we can linearize the nonlinear equation (148) in the vicinity of the fixed points by considering small perturbations v(t) = v∗ + δv(t). (149) By inserting this perturbation ansatz into (148) and noting that, to leading order, f (v∗ + δv) ' f (v∗ ) + f 0 (v∗ ) δv = f 0 (v∗ ) δv, (150) we find that the growth of the perturbation δv(t) is governed by the linear ODE d δv(t) = f 0 (v∗ ) δv(t), dt (151) which has the solution δv(t) = δv(0) ef 0 (v ∗ )t . (152) If f 0 (v∗ ) > 0, then the perturbation will grow and the fixed point is said to be linearly unstable. whereas for f 0 (v∗ ) < 0 the perturbation will decay implying that the fixed point is stable. For our specific example, we find f 0 (v0 ) = −α , 2 f 0 (v± ) = −(α + 3βv± ) = 2α (153) This means that for α > 0, the fixed point v0 = 0 is stable, indicating that the particle will be damped to rest in this case. By pcontrast, for α < 0, the fixed point v0 becomes unstable and the new fixed points v± = ± −α/β become stable; that is, for α < 0 the particle will be driven to a non-vanishing stationary speed. Equation (148) with α < 0 defines one of the simplest models of active particle motion. 7.2 Stability analysis for PDEs The above ideas can be readily extended to PDEs. To illustrate this, consider a scalar density n(x, t) on the interval [0, L], governed by the diffusion equation ∂n ∂2n =D 2 ∂t ∂x (154a) ∂n ∂n (0, t) = (L, t) = 0. ∂x ∂t (154b) with reflecting boundary conditions, This dynamics defined by Eqs. (154) conserves the total ‘mass’ Z L N (t) = dx n(x, t) ≡ N0 , (155) 0 and a spatially homogeneous stationary solution is given by n0 = N0 /L. (156) To evaluate its stability, we can consider wave-like perturbations n(x, t) = n0 + δn(x, t) , δn =  eσt−ikx . (157) Inserting this perturbation ansatz into (154) gives the dispersion relation σ(k) = −Dk 2 ≥ 0, (158) signaling that n0 is a stable solution, because all modes with |k| > 0 become exponentially damped. 7.3 Swift-Hohenberg theory of pattern formation As a simple generalization of (154), we consider the simplest isotropic fourth-order PDE model for a non-conserved real-valued order-parameter ψ(x, t) in two space dimensions x = (x, y), given by (159) ∂t ψ = F (ψ) + γ0 ∇2 ψ − γ2 (∇2 )2 ψ, where ∂t = ∂/∂t denotes the time derivative, and ∇ = (∂/∂x, ∂/∂y) is the two-dimensional Laplacian. The force F is derived from a Landau-potental U (ψ) F =− ∂U , ∂ψ U (ψ) = a 2 b 3 c 4 ψ + ψ + ψ , 2 3 4 (160) where c > 0 to ensure stability. The appearance of higher-order spatial derivatives means that this model accounts for longer-range effects than the diffusion equation. This becomes immediately clear when one writes a (159) in a discretized form as necessary, for example, when trying to solve this equation numerically on a space-time grid: second-order spatial derivatives require information about field values at nearest neighbors, whereas fourth-order derivatives involves field values at next-to-nearest neighbors. In this sense, higher-thansecond-order PDE models, such as the Swift-Hohenberg model (159), are more ‘nonlocal’ than the diffusion equation (154). The field ψ could, for example, quantify local energy fluctuations, local alignment, phase differences, or vorticity. In general, it is very challenging to derive the exact functional dependence between macroscopic transport coefficients (a, b, c, γ1 , γ2 ) and microscopic interaction parameters. With regard to practical applications, however, it is often sufficient to view transport coefficients as purely phenomenological parameters that can be determined by matching the solutions of continuum models, such as the one defined by Eqs. (159) and (160), to experimental data. This is analogous to treating the viscosity in the classical Navier-Stokes equations as a phenomenological fit parameter. The actual predictive strength of a continuum model lies in the fact that, once the parameter values have been determined for a given set-up, the theory can be used to obtain predictions for how the system should behave in different geometries or under changes of the boundary conditions (externally imposed shear, etc.). In some cases, it may also be possible to deduce qualitative parameter dependencies from physical or biological considerations. For instance, if ψ describes vorticity or local angular momentum in an isolated ‘active’ fluid, say a bacterial suspension, then transitions from a > 0 to a < 0 or γ0 > 0 to γ0 < 0, which both lead to non-zero flow patterns, must be connected to the microscopic self-swimming speed v0 of the bacteria. Assuming a linear relation, this suggests that, to leading order, a0 = δ −αv0 where δ > 0 is a passive damping contribution and αv0 > 0 the active part, and similarly for γ0 . It may be worthwhile to stress at this point that higher-than-second-order spatial derivatives can also be present in passive systems, but their effects on the dynamics will usually be small as long as γ0 > 0. If, however, physical or biological mechanisms can cause γ0 to become negative, then higher-order damping terms, such as the γ2 -term in (159), cannot be neglected any longer as they are essential for ensuring stability at large wave-numbers, as we shall see next. Linear stability analysis The fixed points of (159) are determined by the zeros of the force F (ψ), corresponding to the minima of the potential U , yielding ψ0 = 0 (161a) and b ψ± = − ± 2c r b2 a − , 4c2 c if b2 > 4ac. (161b) Linearization of (159) near ψ0 for small perturbations δψ = 0 exp(σ0 t − ik · x) (162) σ0 (k) = −(a + γ0 |k|2 + γ2 |k|4 ). (163) ψ = ψ± + ± exp(σ± t − ik · x) (164)   σ± (k) = − −(2a + bψ± ) + γ0 |k|2 + γ2 |k|4 . (165) gives Similarly, one finds for the dispersion relation In both cases, k-modes with σ > 0 are unstable. From Eqs. (163) and (165), we see immediately that γ2 > 0 is required to ensure small-wavelength stability of the theory and, furthermore, that non-trivial dynamics can be expected if a and/or γ0 take negative values. In particular, all three fixed points can become simultaneously unstable if γ0 < 0. Symmetry breaking In the context biological systems, the minimal model (159) is useful for illustrating how microscopic symmetry-breaking mechanisms that affect the motion of individual microorganisms or cells can be implemented into macroscopic field equations that describe large collections of such cells. To demonstrate this, we interpret ψ as a vorticitylike 2D pseudo-scalar field that quantifies local angular momentum in a dense microbial suspension, assumed to be confined to a thin quasi-2D layer of fluid. If the confinement mechanism is top-bottom symmetric, as for example in a thin free-standing bacterial film, then one would expect that vortices of either handedness are equally likely. In this case, Figure 1: Numerical illustration of structural transitions in the order-parameter ψ for symmetric (a) mono-stable and (b) bi-stable potentials U (ψ) with b = 0. (c) Snapshots of the order-parameter field ψ at t = 500, scaled by the maximum value ψm , for a mono-stable potential U (ψ) and homogeneous random initial conditions. (b) Snapshots of the orderparameter at t = 500 for a bi-stable potential. For γ0  −(2π)2 γ2 /L2 , increasingly more complex quasi-stationary structures arise; qualitatively similar patterns have been observed in excited granular media and chemical reaction systems. (159) must be invariant under ψ → −ψ, implying that U (ψ) = U (−ψ) and, therefore, b = 0 in (160). Intuitively, the transformation ψ → −ψ corresponds to a reflection of the observer position at the midplane of the film (watching the 2D layer from above vs. watching it from below). The situation can be rather different, however, if we consider the dynamics of microorganisms close to a liquid-solid interface, such as the motion of bacteria or sperms cells in the vicinity of a glass slide (Fig. 2). In this case, it is known that the trajectory of a swimming cell can exhibit a preferred handedness. For example, the bacteria Escherichia coli and Caulobacter have been observed to swim in circles when confined near to a solid surface. More precisely, due to an intrinsic chirality in their swimming apparatus, these organisms move on circular orbits in clockwise (anticlockwise) direction when viewed from inside the bulk fluid (glass surface). Qualitatively similar behavior has also been reported for sea urchin sperm swimming close to solid surfaces. Hence, for various types of swimming microorganisms, the presence of the near-by noslip boundary breaks the reflection symmetry, ψ 6→ −ψ. The simplest way of accounting for this in a macroscopic continuum model is to adapt the potential U (ψ) by permitting values b 6= 0 in (160). The result of a simulation with b > 0 is shown in Fig. 2a. In contrast to the symmetric case b = 0 (compare Fig. 1c), an asymmetric potential favors the formation of stable hexagonal patterns (Fig. 2a) – such self-assembled hexagonal vortex lattices have indeed been observed experimentally for highly concentrated spermatozoa of sea urchins (Strongylocentrotus droebachiensis) near a glass surface (Fig. 2b). Figure 2: Effect of symmetry-breaking in the Swift-Hohenberg model. (a) Stationary hexagonal lattice of the pseudo-scalar angular momentum order-parameter ψ, scaled by the maximum value ψm , as obtained in simulationsof Eqs. (159) and (160) with b > 0, corresponding to a broken reflection symmetry ψ 6→ −ψ. Blue regions correspond to clockwise motions. (b) Hexagonal vortex lattice formed spermatozoa of sea urchins (Strongylocentrotus droebachiensis) near a glass surface. At high densities, the spermatozoa assemble into vortices that rotate in clockwise direction (inset) when viewed from the bulk fluid. 7.4 Reaction-diffusion (RD) systems RD systems provide another generic way of modeling structure formation in chemical and biological systems. The idea that RD processes could be responsible for morphogenesis goes back to a 1952 paper by Alan Turing (see class slides), and it seems fair to say that this paper is the most important one ever written in mathematical biology. RD system can represented in the form ∂t q(t, x) = D∇2 q + R(q), (166) where • q(t, x) as an n-dimensional vector field describing the concentrations of n chemical substances, species etc. • D is a diagonal n × n-diffusion matrix, and • the n-dimensional vector R(q) accounts for all local reactions. 7.4.1 Two species in one space dimension As a specific example, let us consider q(t, x) = (u(t, x), v(t, x)), D = diag(Du , Dv ) and R = (F (u, v), G(u, v)), then ut = Du uxx + F (u, v) (167a) vt = Dv vxx + G(u, v) (167b) In general, (F, G) can be derived from the reaction/reproduction kinetics, and conservation laws may impose restrictions on permissible functions (F, G). The fixed points (u∗ , v∗ ) of (167) are determined by the condition   F (u∗ , v∗ ) R(u∗ , v∗ ) = = 0. (168) G(u∗ , v∗ ) Expanding (167) for small plane-wave perturbations     u(t, x) u∗ = + (t, x) v(t, x) v∗ (169a) with  = ˆ e σt−ikx   ˆ σt−ikx = e , η̂ (169b) we find the linear equation   ∗  k 2 Du 0 Fu Fv∗ ˆ ≡ M ˆ, ˆ + σˆ = − 0 k 2 Dv G∗u G∗v  (170) where Fu∗ = ∂u F (u∗ , v∗ ) , Fv∗ = ∂v F (u∗ , v∗ ) , G∗u = ∂u G(u∗ , v∗ ) , G∗v = ∂v G(u∗ , v∗ ). Solving this eigenvalue equation for σ, we obtain   q 1 2 2 ∗ ∗ ∗ ∗ ∗ ∗ 2 σ± = −(Du + Dv )k + (Fu + Gv ) ± 4Fv Gu + [Fu − Gv + (Dv − Du )k ] . (171) 2 In order to have an instability for some finite value k, at least one of the two eigenvalues must have a positive real part. This criterion can be easily tested for a given reaction kinetics (F, G). We still consider a popular example. 7.4.2 Lotka-Volterra model This model describes a simple predator-prey dynamics, defined by F (u, v) = Au − Buv, (172a) G(u, v) = −Cv + Euv (172b) with positive rate parameters A, B, C, E > 0. The field u(t, x) measures the concentration of prey and v(t, x) that of the predators. The model has two fixed points (u0 , v0 ) = (0, 0) , (u∗ , v∗ ) = (C/E, A/B), (173) with Jacobians  Fu (u0 , v0 ) Fv (u0 , v0 ) Gu (u0 , v0 ) Gv (u0 , v0 )   = A 0 0 −C  (174a) and    Fu (u∗ , v∗ ) Fv (u∗ , v∗ ) A − BEC = Gu (u∗ , v∗ ) Gv (u∗ , v∗ ) C −A −C + AE B  . (174b) It is straightforward to verify that, for suitable choices of A, B, C, D, the model exhibits a range of unstable k-modes.