TURING BIFURCATION IN THE SWIFT–HOHENBERG EQUATION ON DETERMINISTIC AND RANDOM GRAPHS

. The Swift–Hohenberg equation (SHE) is a partial differential equation that explains how patterns emerge from a spatially homogeneous state. It has been widely used in the theory of pattern formation. Following a recent study by Bramburger and Holzer [2], we consider discrete SHE on deterministic and random graphs. The two families of the discrete models share the same continuum limit in the form of a nonlocal SHE on a circle. The analysis of the continuous system, parallel to the analysis of the classical SHE, shows bifurcations of spatially periodic solutions at critical values of the control parameters. However, the proximity of the discrete models to the continuum limit does not guarantee that the same bifurcations take place in the discrete setting in general, because some of the symmetries of the continuous model do not survive discretization. We use the center manifold reduction and normal forms to obtain precise information about the number and stability of solutions bifurcating from the homogeneous state in the discrete models on deterministic and sparse random graphs. Moreover, we present detailed numerical results for the discrete SHE on the nearest-neighbor and small-world graphs.


Introduction
Networks of interacting dynamical systems (a.k.a.interacting particle systems) form an important class of models of natural and technological systems.Examples include neuronal networks, swarm of fireflies, coupled lasers, and power grids, to name a few [10,33,34,36].
When the number of particles is large, taking a limit as the number of particles goes to infinity is an effective tool for analyzing network dynamics.The continuum limit in the form of a nonlocal PDE has been very successful for studying synchronization and pattern formation in large systems of coupled oscillators on a variety of graphs [25,28,29,40].Therefore, it is important to understand the accuracy of approximation of large dynamical networks by a continuum equation.For solutions of initial value problems on convergent families of graphs, this was done in [23,24,26].
In many theoretical studies as well as in practical applications, valuable information about system dynamics is gained by studying regimes bifurcating from simpler solutions under the variation of the control parameter.Therefore, it is of interest to understand how well and under what conditions, the bifurcation structure of large networks can be obtained from its continuum limit.Specifically, suppose the continuum model undergoes a bifurcation at a certain value of the control parameter.What can be said about the discrete system?Will it undergo a bifurcation at a close parameter value?Will the bifurcating solutions resemble those obtained in the continuum case?These are nontrivial questions in general, because certain features (such as symmetries) that are present in the continuum limit may not survive discretization.The questions become even more challenging if the large network system is random.
In this paper, following the work of Bramburger and Holzer [2], we address these questions in the context of the Turing bifurcation in the discrete Swift-Hohenberg equation (SHE) on deterministic (Cayley) and random graphs.We postpone the discussion of how our approach and results differ from [2] and turn to the formulation of the discrete and continuous models next.
The classical SHE plays a prominent role in the theory of pattern formation (see [8] and references therein).It has the following form where u(t, x) ∈ R and γ is a control parameter.The normal form reduction near the bifurcation, which makes use of the symmetries present in the system, shows existence of a two-parameter family of stable stationary nontrivial solutions bifurcating from the trivial solution u ≡ 0 (see Section 2.4.3 in [15]): for small positive γ and every δ ∈ T.
The nonlocal SHE is obtained by replacing the second derivative ∂ 2 x with a nonlocal operator L W : 3) where (1.4) Here, W : T 2 → [0, 1] is a measurable function that is symmetric a.e. on T 2 , and f ∈ L 1 (T).
In this paper, we assume that W has the following form for a given even function S ∈ L 1 (T).In this case, d W is independent of x ∈ T. Graphons of the form (1.5) arise as limits of convergent sequences of Cayley graphs.For this reason, they are referred to as Cayley graphons.It is instructive to study this case first, because the nonlocal SHE (1.3) with (1.4) and (1.5) is the closest nonlocal analog of the classical model (1.1) on a periodic domain.Since K W has a discrete (real) spectrum with the only accumulation point at zero, the bifurcation of stable stationary nontrivial solutions occurs if κ is set to κ = −d W + λ k 0 , k 0 ∈ N, where {λ k } k∈Z are eigenvalues of K W satisfying λ −k = λ k , k ∈ N. Similar to (1.2), we are interested in solutions bifurcating from u = 0 for small γ.The normal form analysis of (1.3) closely resembles the normal form analysis of (1.1) and illustrates the role of translation invariance in the continuum model, see Theorem 3.8 below.
Along with the nonlocal model (1.3) we consider a discrete model on a graph Γ given by u = − (L Γ − κ) 2 u + γu − u 3 , ( where u(t) ∈ R n , the cubic nonlinearity u 3 is understood in the componentwise sense, and the graph Laplacian L Γ is defined as follows which is associated with the adjacency matrix A Γ = (a ij ) 1≤i,j≤n .
In this paper, we study (1.6) with (1.7) on two families of graphs, one is deterministic and the other is random.Both families are constructed using a given W and converge to W almost surely.We refer the reader to [23,26] for the overview of the ideas from the graph limit theory that are relevant here and to the references in the literature on this subject.For either graph model, one can show that the solution of the IVP for (1.3) approximates the solutions of the IVPs of the discrete problems on finite time intervals (cf.[23,24,26]).Thus, (1.3) provides a common continuum limit for the discrete models on both deterministic and random graphs.Note that this however does not guarantee that the solutions bifurcating from the trivial solution will resemble (1.2).In fact, in general, it is not clear that the discrete models will undergo a bifurcation at all.The reason for this is that the translation symmetry present in the continuum model does not survive discretization.The lack of the translation symmetry affects computations of the normal forms and thus the bifurcations.Our analysis predicts the exact number of solution families bifurcating from the trivial solution in the deterministic model, see Theorem 4.2, based on the discrete group of symmetries.On the other hand, the symmetries are broken in the random model and we are only able to find the lower and upper bounds on the number of solution families bifurcating from the trivial solution, see Theorem 5.1.
The organization of the paper is as follows.In Section 2, we complete the description of the discrete model (1.6) with (1.7) by providing the details on the deterministic and random graphs.The former are weighted Cayley graphs and the latter are W-random graphs (cf.[21]).
After that we turn to the analysis of the continuous model in Section 3. Specifically, we demonstrate the well-posedness of the continuous model and analyze the bifurcation of the spatially homogeneous solutions.The bifurcation analysis is based on the normal form of the system in the Fourier coordinates.The derivation of the normal form uses the symmetries present in the system and resembles the analysis of the classical SHE in [15].The existence of a family of spatially periodic solutions is invariant with respect to the continuous spatial translations.
In Section 4, we study discrete SHE on weighted Cayley graphs.The bifurcation analysis follows the same lines as in the continuous case, where in place of Fourier transform we now use the discrete Fourier transform.An important distinction of the discrete model is that the invariance under any translations in the continuous system is replaced by the invariance with respect to a discrete group of translations.This results in finitely many solutions bifurcating from the trivial solution instead of the continuous family (1.2).
In Section 5, we move on to study the discrete model on random graphs, which is the main problem addressed in this paper.In the random setting, the discrete model does not possess the discrete translational invariance and the normal form approach, which worked for Cayley graphs, is no longer applicable.To overcome this obstacle, we use the proximity of the system on a sufficiently large random graph to that analyzed in Section 4 to derive a leading order approximation for the normal form of the system at hand.This allows us to study the bifurcation for the system on random graphs and to obtain precise bounds on the number of solution families bifurcating in each random realization of the discrete graph.This is where our approach is different from the approach taken in [2].The normal form for the bifurcation on a random graph is compared with the one on the associated discrete deterministic graph instead of the one on the associated continuous nonlocal model.We determine the translational parameter precisely from the normal form, whereas the translational parameter is defined implicitly in [2, Theorem 4.1] from a linear transformation of eigenvectors of the linearized equations.Additional comments on the differences between our work and [2] can be found in Remark 5.6.
In Section 6, we present numerical experiments with SHE on small-world graphs.This example illustrates the selection of random stationary patterns in SHE on random graphs.We also explain the implications of the lack of continuous translational symmetry in the corresponding averaged SHE on deterministic Cayley graph.In Section 7, we discuss the utility of graphons in the analysis of dynamical systems on graphs and potential applications of our techniques to related network models.
The approximation result needed for the analysis of SHE on random graphs is given in the Appendix.It is derived in using the method of [12] based on the concentration inequality for adjacency matrices of W -random graphs.

Discretization
The goal of this section is to complete the formulation of the discrete model (1.6) by supplying the details on the families of deterministic and random graphs.In what follows, we assume that n is even in (1.6) and denote N = n/2.This assumption is used to simplify computations of the normal forms.We will denote the deterministic graph by Γ N W and the random graph by ΓN W .
2.1.The discrete SHE on deterministic graphs.To define the family of deterministic graphs (Γ N W ), we fix N ∈ N and discretize T as follows.Let h .= 1 2N , x i = ih, and Γ N W is a weighted graph on 2N nodes indexed by integers from [−N +1, N ].An edge between nodes i and j is supplied with a weight In addition, we assume Since W (x, y) = S(x − y) according to (1.5), we have a Toeplitz matrix A = (a ij ) with If S k := a k0 , then a ij = S i−j and the graph Laplacian on Γ N W can be rewritten in the form: where d N W is a constant independent of i.The SHE on the deterministic graph Γ N W has the following form: where 2.2.The discrete SHE on random graphs.The second family of graphs ( ΓN W ) is random and corresponds to Γ N W , N ∈ N. Denote the adjacency matrix of ΓN W by ÃN = (ã ij ) −N +1≤i,j≤N .We postulate that two distinct nodes of ΓN W i and j are connected with probability a ij , i.e., In addition, ãii = 0 and ãji = ãij .The SHE on the random graph ΓN W is given in the form: where where DN W is a diagonal matrix.
Graphs defined by (2.4) are dense almost surely.In dense graphs the number of edges scales quadratically with the number of vertices, i.e., under this model Γ N W has O(N 2 ) edges.In this paper, we extend the random graph model to allow for sparse graphs.To this end, we introduce a nonincreasing sequence α N ∈ (0, 1] and modify (2.4) as follows As before, ãii = 0 and ãji = ãij .
If α N ≡ 1 then (2.7) reduces to (2.4) and we obtain the sequence of dense random graphs as above.On the other hand, if α N ↘ 0 then the expected number of edges in ΓN W is N (N − 1)α N ≪ N 2 , which implies that ΓN W is a sparse graph with probability 1.By varying the rate of convergence of α N to zero, one can control the degree of sparseness of ΓN W .We need to impose the following technical condition on (α N ): for some M > 0 dependent of N .
We illustrate the families of graphs, which we just defined, with the following example.
Remark 2.3.The deterministic SHE model (2.3) is a Galerkin approximation of the continuous SHE model (1.3).On the other hand, it is also related to the random model (2.5) via averaging, because E ãij = a ij by construction.Therefore, one can view the continuous SHE model (1.3) as a continuum limit of either of the discrete models (2.3) or (2.5).For a related class of nonlocal models, it is known that the initial-value problems for the deterministic and random discrete models approximate that for the continuum one [23,26].The same techniques apply to the models at hand and the approximation results continue to hold for the discrete and continuum SHEs.However, the analysis in the remainder of this paper does not depend on the validity of these results.

The continuum SHE
In this section, we study the nonlocal SHE model (1.3), which serves as a continuum limit for the discrete SHE models (2.3) and (2.5) on deterministic and random graphs Γ N W and ΓN W respectively.Our objective is to obtain a spatially dependent steady state via a Turing bifurcation of the trivial solution.For W in the form (1.5), the nonlocal SHE on Cayley graphon can be written in the form where κ and γ are real parameters and Throughout this section, we assume that S is a given integrable function on T, i.e., S ∈ L 1 (T).Furthermore, S is assumed to be even almost everywhere.By Young's convolution inequality, K S is a bounded operator from L p (T) to L p (T) for any 1 ≤ p ≤ ∞ with the bound Remark 3.1.Since K S is a bounded operator from L 2 (T) to L 2 (T) and T is compact, the spectrum σ(L S ) is purely discrete and consists of eigenvalues {λ k } k∈Z obtained from the Fourier modes: Since S is an even function, consists of real eigenvalues.However, unless S is even on T, K S is not a self-adjoint operator in L 2 (T) and the eigenvalues {λ k } k∈Z only satisfy the reduction λ k = λ−k for k ∈ N which does not exclude the possibility of complex eigenvalues.
The following lemma gives the well-posedness result of the time-evolution problem (3.1) in the phase space X := L ∞ (T), which is continuously embedded into L 2 (T) due to the bound ∥f ∥ L 2 (T) ≤ ∥f ∥ L ∞ (T) .Lemma 3.2.Assume that S ∈ L 1 (T) and its periodic extension to R is an even function.For every u 0 ∈ X , there exists a unique global solution u ∈ C 1 ([0, ∞), X ) such that u(0, •) = u 0 .The unique solution u ∈ C 1 ([0, τ 0 ], X ) is Lipschitz continuous with respect to u 0 ∈ X for every finite τ 0 > 0.
Proof.By (3.4) with p = ∞, K S is a bounded operator from X to X .In addition, the nonlinear term of (3.1) is closed in X due to the bound is a C 1 map from X to X .By the standard results of the semi-group theory [4, Chapter 3], for every u 0 ∈ X , there exists a unique local solution u ∈ C 1 ([0, τ 0 ], X ) for some τ 0 > 0.
Thanks to the repulsive cubic nonlinearity, we have the following bound: ). Lipschitz continuity of the local solution u ∈ C 1 ([0, τ 0 ], X ) with respect to u 0 ∈ X for every finite τ 0 > 0 follows from Gronwall's inequality.□ Remark 3.3.The nonlinear term of (3.1) is not closed in L 2 (T).However, it is closed in H 1 per (T) given by since H 1 per (T) is a Banach algebra with respect to pointwise multiplication.Moreover, H 1 per (T) is continuously embedded into a space of bounded and continuous functions satisfying the periodic boundary conditions.Compared to H 1 per (T), functions in the phase space X = L ∞ (T) do not have to be continuous or to satisfy the periodic boundary conditions.This is more suitable in the context of solutions of the discrete SHE on deterministic and random graphs.
Proof.The continuous SHE (3.1) admits the following symmetries: • the spatial translation x → x + h, ∀h ∈ R due to periodic conditions, • the spatial reflection x → −x due to even S, • the sign reflection u → −u due to odd nonlinearity, which can be easily confirmed.The new solutions are generated by the symmetries.□ Remark 3.5.The nonlocal SHE (3.1) has two real parameters γ and κ.Parameter γ is small and is used to characterize Turing bifurcation of a spatially dependent steady state from the zero solution.On the other hand, parameter κ is the tuning parameter defined by the bifurcation condition according to the following definition.
Assumption 3.6.We fix k 0 ∈ N, assume that λ k ̸ = λ k 0 for every k ∈ N\{k 0 }, and choose Remark 3.7.Since K S is a compact operator on L 2 (T), eigenvalues {λ n } n∈Z satisfy λ n → 0 as |n| → ∞.Hence, Assumption 3.6 implies that there is C 0 > 0 such that The following theorem presents the main result of this section.Although it is an exercise from [15, Section 2.4.3],we write the computational details explicitly, since they are useful for analysis of Turing bifurcation in the discrete SHE.
Theorem 3.8.Under Assumption 3.6, there exists γ 0 > 0 and C 0 > 0 such that for every γ ∈ (0, γ 0 ) there exists a non-trivial time-independent solution u γ (• + δ) in X of the nonlocal SHE model (3.1),where u γ is an even function satisfying and δ ∈ T is an arbitrary translational parameter.The orbit of time-independent solutions {u γ (• + δ)} δ∈T is asymptotically stable in the time evolution of the nonlocal SHE in X .
Proof.In order to use the center manifold theorem and to derive slow dynamics along the center manifold, we use the Fourier series and obtain the evolution problem in the form where we have used that κ = λ k 0 − d S .By standard results from Fourier analysis, u(t, The nonlinear term of (3.7) is closed since ℓ 1 (Z) is a Banach algebra with respect to the convolution sum.
Since u is real, we have a −k = āk for all k ∈ Z. Due to the three symmetries identified in Lemma 3.4, the vector field in (3.7) is equivariant under the transformation and under the transformations: where k 0 ∈ N is defined in Assumption 3.6.
By Assumption 3.6, the rest of eigenvalues {λ k } k∈Z\{±k 0 } are bounded away from λ k 0 = λ −k 0 with the bound (3.5).By the center manifold theorem [15, Theorem 2.9], there exists a center manifold in Xbif spanned by A := a k 0 ∈ C and Ā := a −k 0 .Since the system is closed on (3.10), the center manifold can be expressed as graphs of functions: The dynamics on the center manifold can be expressed by the amplitude equations where F 1 is a C ∞ function in A and Ā with γ-dependent coefficients which commutes with the symmetries of (3.7).Due to the equivariance (3.8), the amplitude equations can be transformed to the normal form: ) where P 1 is an analytic function in |A| 2 with γ-dependent coefficients.Moreover, due to the symmetry with respect to the transformation a k → a −k , P 1 has real-valued coefficients.Similarly, we express functions Ψ m in the form: where P m is a C ∞ function in |A| 2 with γ-dependent real-valued coefficients.
Due to the cubic nonlinearity in (3.7), we obtain where the remainder terms of the order of O(|A| 4 ) is defined by Ψ 3 since Ψ m with m ≥ 5 give a higher-order contribution of O(|A| 6 ) to the normal form (3.13).It follows from (3.13) that there exists a time-independent solution of the form where δ is an arbitrary parameter.This time-independent solution (3.14) yields a non-trivial time-independent solution of the nonlocal SHE (3.1) in X satisfying the expansion in (3.6).
To determine stability of the time-independent solution u γ (• + δ) for every given δ ∈ R, we note that all eigenvalues in the spectrum of −(K S − d S − κ) 2 are located in the lefthalf plane of the complex plane with the exception of the double zero eigenvalue.Hence, there is no unstable manifold of the system (3.7).The time-independent solution (3.14) is orbitally asymptotically stable in the time evolution of the reduced equation (3.11).By the standard decomposition near the orbit {u γ (• + δ)} δ∈T , a perturbation of the initial solution u γ (• + δ 0 ) defines a time-dependent solution which approaches as t → +∞ exponentially fast to the final solution u γ (• + δ ∞ ) with δ ∞ being close to δ 0 .Hence, the orbit {u γ (• + δ)} δ∈T is asymptotically stable in the time evolution of the nonlocal SHE in X .□ Remark 3.9.The leading-order approximation for Ψ 3 in (3.12) can be obtained from (3.7) and (3.13).Substitution yields which is solved by expanding P 3 (|A| 2 ) in |A| 2 with the leading order: where λ 3k 0 ̸ = λ k 0 by the assumption.Similarly, one can find the leading-order expansions of Ψ m in (3.12) for m ≥ 5.
Related to the time-independent solution u γ (• + δ) of the nonlocal SHE in Theorem 3.8, we can introduce the linearized operator in the form where we have used that κ = λ k 0 − d S according to Assumption 3.6.Since u γ is a bounded function on T, L γ is a bounded operator from L 2 (T) to L 2 (T) for all sufficiently small γ for which u γ is defined.Hence, the spectrum of L γ is purely discrete and consists of real eigenvalues, see Remark 3.1.The following lemma uses the smallness of γ and gives a precise information on the location of these eigenvalues.
Proof.The existence of Λ 1 (γ) = 0 follows from the translational invariance of the nonlocal SHE given by (3.1) since δ ∈ R is a free parameter of the steady state u γ (• + γ) ∈ X in Theorem 3.8.Since u ′ γ ∈ L 2 (T) is assumed, we obtain by direct differentiation that The rest of the spectrum of L γ in L 2 (T) follows from the perturbation theory for selfadjoint operators with purely discrete spectrum which implies continuity of eigenvalues with respect to small parameter γ ∈ (0, γ 0 ).The eigenvalue Λ 2 (γ) coincides with the linearization of the slow motion at the center manifold given by (3.13): which are strictly negative and bounded away from 0 by Assumption 3.6.Hence Due to the presence of zero eigenvalue Λ 1 (γ) = 0, the operator L γ is not invertible.This creates difficulties in the persistence argument when the limiting nonlocal model (3.1) is replaced by the discrete SHE models on the deterministic or random graphs, which destroy the continuous translational symmetry.

The discrete SHE on Cayley graphs
Here we study the discrete SHE (2.3) on the deterministic Cayley graph Γ N W .By using (2.2), we rewrite the discrete model in the form: where N is integer, κ and γ are real parameters, and the linear operator A N W : Z N → Z N is given by the convolution sum where {S j } j∈Z N satisfies S −j = S j for all j ∈ Z N .System of differential equations (4.1) can be viewed as an evolution equation on Z N .The classical solutions are interpreted as elements of C 1 (R, R Z N ), the space of continuously differentiable vector functions of t ∈ R.
Lemma 4.1.If {u j (t)} j∈Z N is a solution of the discrete SHE (4.1), so are Proof.Similarly to the continuous SHE (3.1), the discrete SHE (4.1) admits the following three symmetries: • the discrete spatial translation j → j + m, ∀m ∈ Z N due to periodic conditions, • the spatial reflection j → −j due to even {S j } j∈Z N , • the sign reflection u → −u due to odd nonlinearity, which can be easily confirmed.The new solutions (4.3) are generated from {u j (t)} j∈Z N by symmetries.□ Our objective is to obtain a spatially dependent steady state of the discrete SHE (4.1) via a Turing bifurcation of the zero solution.By using the discrete Fourier transform, one can obtain eigenvalues of A N W in the form {λ N k } k∈Z N with We again use parameter γ in (4.1) to characterize Turing bifurcation and parameter κ to satisfy the bifurcation condition.
We will locate a spatially dependent steady state bifurcating from the zero solution by adapting the proof of Theorem 3.8 with the discrete Fourier transform replacing Fourier series.One key new feature of the center manifold analysis is that the translational parameter δ ∈ T which was arbitrary in Theorem 3.8 takes exactly 4N admissible values under some non-degeneracy conditions.In comparison with Assumption 3.6, we set k 0 = 1 for the bifurcating mode in order to simplify computations of the normal form.The following theorem presents the main result of this section.15), then there exists γ 0 > 0 and C 0 > 0 such that for every γ ∈ (0, γ 0 ) and every integer N ≥ 3 there exist two non-trivial time-independent solutions u N γ , v N γ ∈ R Z N of the discrete SHE (4.1), where u N γ is symmetric about j = 0 and satisfies and v G γ is symmetric about the mid-point between j = 0 and j = 1 and satisfies One of the two solutions is asymptotically stable in the time evolution of the discrete SHE in C 1 (R, R Z N ) and the other one is unstable.These solutions generate (2N ) asymptotically stable and (2N ) unstable solutions on Z N via the discrete group of spatial translations.
Proof.We use the discrete Fourier transform with real a 0 , a N , and possibly complex a −k = āk for 1 ≤ k ≤ N − 1.The Fourier amplitudes are extended periodically with the period 2N as the sequence {a k } k∈Z N .The discrete SHE (4.1) transforms to the evolution problem in the form where the summation terms are adjusted due to the (2N )-periodicity of the Fourier modes in the Fourier space and we have used that κ = λ N 1 − d N W . Symmetries in Lemma 4.1 imply that the system (4.8) is equivariant under the transformation and under the transformations: a k → a −k and a k → −a k .In particular, the system (4.9) is closed on the subspaces We apply the center manifold theorem [15, Theorem 2.9] under the assumption that λ N k ̸ = λ N 1 for k ̸ = ±1.Similar to the proof of Theorem 3.8, there exists a center manifold of the system (4.8)spanned by a 1 ≡ A ∈ C and a −1 = Ā.Since the system is closed on (4.11), the center manifold can be expressed as graphs of functions: (4.12) The slow dynamics on the center manifold can be expressed by the amplitude equations where F 1 is a C ∞ function in A and Ā with γ-dependent coefficients which commutes with the symmetries of (4.8).Compared to the normal form for the amplitude equations in (3.11) and (3.12) the discrete Fourier modes are (2N )-periodic so that and Āe Following the classification of normal forms under the symmetry (4.9) in [7], we transform the amplitude equations to the normal form: where , and Ā2N with γ-dependent coefficients.Due to the symmetry with respect to the transformation a k → a −k , Q 1 and R 1 have real-valued coefficients.Similarly, we express functions Ψ m in the form: where Q m and R m are C ∞ functions in |A| 2 , A 2N , and Ā2N with γ-dependent real-valued coefficients.
We assume the following non-degeneracy condition: where r N is a γ-dependent real-valued coefficient.Similarly to (3.13), we obtain from the cubic nonlinearity in (4.8) that if N ≥ 3.By using the polar form A = ρe iθ , we write If r N ̸ = 0, there exist two distinct time-independent solutions for θ = 0 and θ = π 2N on interval [0, π N ).If N ≥ 3, both solutions can still be expressed at the leading order in the form: with either δ = 0 (which corresponds to θ = 0) or δ = 1 2 (which corresponds to θ = π 2N ).For δ = 0, we get the real solution {a k } k∈Z N on the subspace (4.10) that corresponds to the solution u N γ of the discrete SHE (4.1) which is symmetric about j = 0 and satisfies (4.5).For δ = 1 2 , we get a complex solution {a k } k∈Z N that corresponds to the real solution v N γ of the discrete SHE (4.1) which satisfies (4.6).The leading order and hence the solution v N γ is symmetric about the mid-point between j = 0 and j = 1 due to the symmetry with respect to the transformation u j → u 1−j in Lemma 4.1.
To obtain the stability conclusion for the time-independent solutions, we observe that all eigenvalues in the spectrum of −(A N W − λ N 1 ) 2 are located in the left-half plane of the complex plane with the exception of the double zero eigenvalue.If r N ̸ = 0, then linearization of the time-dependent equation (4.13) at the leading-order solution (4.16) for N ≥ 3 yields for the perturbations (ρ ′ , θ ′ ): where the upper sign corresponds to the solution u N γ and the lower sign corresponds to the solution v N γ .Dynamics of (4.17) in ρ ′ is asymptotically stable, whereas dynamics of (4.17) in θ ′ is either asymptotically stable or unstable.This yields the conclusion that one of the two solutions u N γ and v N γ is asymptotically stable in the time evolution of the discrete SHE in R 2N  per and the other one is unstable.
The two solutions u N γ and v N γ generate (4N ) solutions by using the discrete group of translations u j → u j±m for every m ∈ Z N by Lemma 4.1.The stability of the translated solutions coincide with the stability of u N γ and v N γ .□ Remark 4.4.The graphons defined by the graphs Γ N W used the formulation of the discrete SHE (4.1) converge to W , the kernel used in the continuous SHE (3.1), in the cut-norm as N → ∞.Thus, λ N k → λ k for every fixed k ∈ Z [37].From this we conclude that Assumption 3.6 with k 0 = 1 implies the corresponding assumption of Theorem 4.2, i.e., λ N k ̸ = λ N 1 , k ̸ = ±1 for sufficiently large N .Remark 4.5.If the non-degeneracy condition (4.15) is not satisfied, one needs to expand functions Q 1 and R 1 in (4.13) to the higher orders and to obtain the higher-order nondegeneracy conditions.If it happens that R 1 ≡ 0 and Q 1 (|A| 2 , A 2N , Ā2N ) = P 1 (|A| 2 ), then the time-independent solution (4.16) exists with arbitrary δ ∈ R/Z.In this degenerate case, there exists a family of non-trivial time-independent solutions u N γ,δ ∈ R Z N of the discrete SHE (4.1) with arbitrary parameter δ ∈ [0, 1], where u γ,δ satisfies sup The orbit of time-independent solutions {u N γ,δ } δ∈R/Z is asymptotically stable in the time evolution of the discrete SHE in C 1 (R, R Z N ).Although such degenerate cases may exist in other discrete models, see, e.g., [16], the explicit computations for the particular discrete SHE model (4.1) show that r N ̸ = 0 for every N ≥ 3.
Remark 4.6.Cases N = 1 and N = 2 are exceptional.In both cases, the discrete Fourier transform in the subspace (4.11) is given by the sum of two terms If N = 1, then A in (4.19) is real and satisfies from which the stable nontrivial solutions at A = ± √ γ/2 exist in addition to the unstable solution A = 0 for γ > 0. If N = 2, then A in (4.19) is complex and satisfies Using the polar form A = ρe iθ , this equation is reduced to the system ρ = γρ − 3ρ 3 − ρ 3 cos(4θ), θ = ρ 2 sin(4θ), from which the two nontrivial time-independent solutions are given by .
The linearization shows that the first solution is linearly unstable and the second solution is asymptotically stable.
The following two examples give details of computations of the center manifold reductions for N = 3 and N = 4, for which the discrete Fourier transform in the subspace (4.11) is given by the sum These details show that the non-degeneracy condition (4.15) is satisfied for N = 3, 4.
Example 4.7.If N = 3, then A is complex and B is real.The system (4.1) in the decomposition (4.20) is reduced to the system of two equations We compute the center manifold reduction B = Ψ 3 (A, Ā) in powers of A according to (4.14) with real B by writing where c 0 is a real coefficient that depends on γ.From the system of differential equations for A and B, we find at the cubic order that , which agrees with the expansion (3.15).Substituting B = Ψ 3 (A, Ā) into the first equation of the system, we obtain consistently with (4.13) that the slow dynamics of A is given by Since r N =3 = −3c 0 > 0, the non-degeneracy condition (4.15) is satisfied.
Example 4.8.If N = 4, then both A and B are complex.The system (4.1) in the decomposition (4.20) is reduced to the system of two complex-valued equations We compute the center manifold reduction B = Ψ 3 (A, Ā) in powers of A according to (4.14) by writing where c 0 , c 1 , and b 0 are real coefficients that depend on γ.From the system of differential equations for A and B, we find recursively at the cubic and quintic powers of A that .
Substituting B = Ψ 3 (A, Ā) into the first equation of the system, we obtain consistently with (4.13) that the slow dynamics of A is given by ) < 0, the non-degeneracy condition (4.15) is satisfied.Remark 4.9.One can show with the method of mathematical induction that the nondegeneracy condition (4.15) is satisfied for every N ≥ 3.
The following lemma is important for the persistence argument, when the discrete SHE is perturbed by a small correction term.Lemma 4.10.Let u N γ , v N γ ∈ R Z N be defined by Theorem 4.2 for γ ∈ (0, γ 0 ) under the nondegeneracy condition (4.15).Then, the matrix operator , where (u γ ) 2 is a diagonal matrix computed on the squared entries of u N γ , v N γ , is invertible.
Proof.This follows from the available information about the eigenvalues of the linearized system (4.8) at the time-independent solutions {a k } k∈Z N constructed from (4.13) with (4.12) and (4.14).Eigenvalues of the linearized system (4.17) are bounded away from zero and so are eigenvalues −(λ N k − λ N 1 ) 2 for k ̸ = ±1.□ Remark 4.11.The inverse matrix A −1 γ in Lemma 4.10 behaves poorly as γ → 0 because the linearized system (4.17) has one eigenvalue −2γ +O(γ 2 ) and the other eigenvalue of O(γ N −1 ) under the non-degeneracy condition (4.15).As a result,

The discrete SHE on W -random graphs
We now turn to the discrete SHE model (2.5) on the W -random graph ΓN W .Using (2.6), we rewrite it as follows We are interested in the Turing bifurcation of the trivial solution of (5.1).Appendix A shows that the matrix ÃN W is close to the matrix A N W in the cut-norm.Consequently, eigenvalues of ÃN W are close to the eigenvalues {λ N k } k∈Z N of A N W given by (4.4) [37].The following theorem presents the main result on the steady-state solutions in (5.1).
Fix γ ∈ (0, γ 0 ) with γ 0 given in Theorem 4.2 and with r N ̸ = 0 in (4.15).There exist N 0 ≥ 3 such that for every N ≥ N 0 there exist at least 4 and at most 4N values of δ such that the discrete SHE (5.1) admits time-independent solutions Half of solutions are asymptotically stable in the time evolution of the discrete SHE (5.1) in C 1 (R, R Z N ) and the other half of solutions are unstable.
Remark 5.2.The system (5.1)does not have any symmetries except of the symmetry with respect to the sign reflection u → −u.Theorem 4.2 gives existence and stability of two distinct state u N γ and v N γ for sufficiently small γ > 0, which are translated to every point of the lattice chain by the discrete translational symmetry of Lemma 4.1.Since the lattice chain in Z N has 2N sites, we can count 4N distinct steady-state solutions in the discrete SHE model (4.1), of which 2N are stable and 2N are unstable.Compared to this conclusion, we do not have an exact count of the number of steady-state solutions on the random graph because of the broken symmetries.The number of steady solution is a random number divisible by 4 between 4 and 4N.
Remark 5.3.Recall that the matrix A γ in Lemma 4.10 has a very small eigenvalue of the size O(γ N −1 ) for large N ≥ 3. Due to this small eigenvalue, we cannot not prove that the steady-state solutions of (4.1) persist as the steady-state solutions of (5.1) because the perturbation is not sufficiently small, see Appendix A, and the implicit function theorem cannot be used for the persistence argument.To overcome this problem, we develop again the approach based on the center manifold reduction, where the main difference is that the linear part of (5.1) is no longer diagonalizable by the discrete Fourier transform.The cubic nonlinear term still enjoys the same transformation under the discrete Fourier transform as in the proof of Theorem 4.2.The other distinction from the deterministic setting is that we can no longer use γ > 0 as a small continuation parameter.Instead, we consider the small parameter γ as fixed in (0, γ 0 ) with γ 0 given in Theorem 4.2 and continue the solution with respect to an additional small parameter µ induced by randomness which is only small for sufficiently large N ≥ N 0 .
Proof of Theorem 5.1.As in the proof of Theorem 4.2, we use discrete Fourier transform To apply the result in Appendix A, we write u = F a, F = ω jk −N ≤j,k≤N −1 , ω = e iπ/N , with F −1 = (2N ) −1 F * , where F * stands for the conjugate transpose of F .By applying the inverse Fourier transform to both sides of (5.1), we have where Recall that similarity transformation A N W → F −1 A N W F diagonalizes the linear part of the deterministic model (4.1).Thus, where and we have used that κ = λ N 1 − d N W . Ã is a symmetric matrix, whose entries above the main diagonal are independent random variables.Further, E Ã = A. From these facts, by Bernstein's inequality, it follows that with probability at least 1 − O 25 −N , we have max (cf.Lemma A.1). Pick µ > 0 small and fixed.It follows from (5.6) that with high probability for sufficiently large N , we have max Proceeding along the lines of the proof of Theorem 4.2, we express the slow manifold of the system (5.3) as the graph of functions a 0 = Ψ 0 (A, Ā), a 1 = A, a −1 = Ā, and Note that the difference from (4.12) that the general L couples the odd-numbered and evennumbered Fourier amplitudes.
The functions Ψ k (A, Ā), k ̸ = ±1 can be obtained from the condition that γ and µ are small, whereas (λ N k − λ N 1 ) 2 in Lkk are strictly positive and bounded away from zero for every k ̸ = ±1.The dynamics on the slow manifold can be expressed by the amplitude equation: where F 1 is a C ∞ function in A and Ā with (γ,µ)-dependent coefficients.The Taylor series expansion includes only odd powers of A and Ā due to the only symmetry of the random model (5.1) with respect to the sign reflection u → −u.
The power expansion of F 1 (A, Ā) is different from those in the proof of Theorem 4.2 due to the presence of the perturbation terms in Lκ which are of the order of O(µ).In addition, it is no longer true that F 1 (A, Ā) can be written in the form (4.13).Nevertheless, the nonzero terms in the expansion of (4.13) remain dominant terms in the expansion of (5.8) if µ is sufficiently small.Thus, we have with some (γ, µ)-dependent coefficients which are bounded as |γ| + |µ| → 0. Since γ ∈ (0, γ 0 ) is fixed and µ in (5.7) can be chosen sufficiently small, we have γ + µα 1 > 0, −3 + µβ 1 < 0, and r N + O(µ) ̸ = 0.By using the polar form A = ρe iθ , we write ( Since µ is selected to be much smaller than γ, there exists only one positive root of equation given by independently of the value of θ.The value of θ is defined from the second equation of system (5.10) with ρ 2 = O(γ).Since there is no discrete translational symmetry of the system (5.1),compared to the system (4.1),we have to consider θ defined on [0, 2π).Roots of θ are defined from equation where ρ is expressed from (5.11).Since the left-hand-side is a trigonometric polynomial in 2θ, there exist at least four roots of θ in [0, 2π) and at most 4N roots since r N + O(µ) ̸ = 0.The total number of roots is divisible by 4. For each root of θ, the root of ρ in (5.11) is uniquely defined and the bound (5.2) with δ := N π θ follows.The stability conclusion of Theorem 5.1 follows from the linearization of system (5.10)near each root and from the fact that all other eigenvalues of − Lκ are strictly negative of the order of O(1) for small γ and µ.Using notations (ρ ′ , θ ′ ) for perturbation terms to the root (ρ, θ), we write the linearized equations in the form: Linearized evolution in ρ ′ is asymptotically stable, whereas linearized evolution in θ ′ is asymptotically stable for half of solutions and is unstable for the other half of solutions.□ Remark 5.4.Coefficients α 1 and α 2 in (5.9) can be easily computed from the linear part of system (5.3).Since Ψ k (A, Ā) = O(µ)|A| for k ̸ = ±1, we have at the leading order: Remark 5.5.If α 2 ̸ = 0, we have exactly four time-independent solutions, from which two are asymptotically stable and two are unstable.The two stable (or unstable) solutions are related to each other by the sign reflection symmetry u → −u of the discrete SHE (5.1).This follows from the fact that the center manifold reduction relies on a trigonometric polynomial in 2θ for which θ 0 and π + θ 0 are equivalent points.If ∆ = 0 in (5.4), then we have 4N time-independent solutions, identically to the outcome of Theorem 4.2.
Remark 5.6.In comparison with Theorem 4.1 in [2], we do not include the quadratic terms in the discrete SHE models to avoid technical computations of the near-identity transformations.We also specify the particular case k 0 = 1 in Theorems 4.2 and 5.1 to simplify computations of the normal forms.On the other hand, we give a precise statement of how the translational parameter δ of Theorem 3.8 is determined in the case of the discrete graphs (both in the deterministic and random cases) and how many time-independent solutions exist for the discrete graph models.In addition, the proof of Theorem 4.1 in [2] is incomplete.The analysis of the Fourier mode w 2 corresponding to the small eigenvalue l 2 = −δ 2 N 2 ρ 2 is not included in the proof.In our setting, the dynamics of w 2 is captured by the equation for θ in (5.10).This equation is important, because it determines the number of branches bifurcating from the spatially homogeneous equilibrium.

The discrete SHE on small-world graphs
In this section, we illustrate the bifurcation analysis of the discrete SHE models on the deterministic and random graphs with numerical results.To this end, we use the family of small-world graphs from Example 2.1.This is a representative example, for which Assumption 3.6 can be verified analytically.
Recall the definition of the small-world graphon W (x, y) = S(x − y) with S ∈ L 1 (T) given by where p ∈ [0, 1] and r ∈ 0, 1 2 .The eigenvalues of the Hilbert-Schmidt operator K S in (3.2) are known explicitly (cf.[5]): By Theorem 3.8, for small γ > 0, the continuous SHE model (3.1) has a continuous family of asymptotically orbitally stable solutions {u γ (• + δ)} δ∈T , where u γ is approximated by We shall now consider how the stable solutions with the expansion (6.3) persist in the discrete SHE models on the deterministic and random graphs.For the discrete SHE model (4.1) with (4.2) and (6.1), we compute the eigenvalues of A N W in the form For small γ > 0, Theorem 4.2 yields existence of the two discrete families of solutions {σ m u N γ } m∈Z N and {σ m v N γ } m∈Z N , where σ m is the shift operator defined by (σ m u) j = u j+m , j ∈ Z N and the profiles of u N γ and v N γ are approximated respectively by ) for j ∈ Z N .One of the two solutions is asymptotically stable and the other solution is unstable.
The random SHE model (5.1), we use the same κ and define the symmetric matrix ÃN W with zeros on the main diagonal and with the entries above the main diagonal ãN ij being independent random variables such that , where a N ij is a coefficient of the adjacency matrix of the weighted Cayley graph.Theorem 5.1 implies the existence of asymptotically stable solutions with the approximation for small γ > 0 and appropriate fixed δ.Both the deterministic and random models were initialized with the leading order term on the right-hand side of (6.3).The shift δ was set to −π/2 in plot a and to 0 in plot b.The two models were integrated for 100 units of time.The solutions of the deterministic models are plotted in red and the solutions of the random model are plotted in blue.
On finite time intervals, the solutions of the discrete SHE on deterministic and random graphs are expected to remain close, provided the initial data are sufficiently close and N is large (cf.[26]).This is clearly seen in the simulations shown in Figure 3.The snapshots of the trajectories of the deterministic and the random models, shown in red and blue respectively, lie very close to each other.For the deterministic model, we know that the trajectory relaxes to one of the 2N stable states lying in the vicinity of the initial condition.As to the snapshots of the SHE on random graph, we see that they relax to the form predicted by Theorem 5.1.We do not know whether the solutions of the random SHE shown in Figure 3 are close to the actual equilibria, because the evolution in the translational direction is extremely slow due to the O(γ N −1 ) eigenvalue, and cannot be resolved by numerics for large N .To illustrate the effects of the slow drift in numerical simulations, we have performed simulations of the deterministic and random SHE models for N = 5. Figure 4 shows outcomes of the numerical simulations for the discrete SHE model (4.1) with two different initial conditions given by the red and cyan lines.The two initial data are given by (6.3) with a half of the amplitude and shifted to the right and to the left relative to the grid points (shown by magenta dots).The two solutions quickly converge to the near-equilibrium solutions with the correct amplitude of 0.2 (blue and black curves) which then slowly shift further from the grid points towards the mid-point between the two grid points.The snapshot on the left panel was stopped at time t = 100, but the slow drift of the zeros of the two solutions on the much longer time interval is shown on the right panel.Zeros are computed from the linear interpolation between the grid points and the color scheme is the same as that for the final states on the left panel.This numerical experiment shows that the asymptotically stable solution is given by a discrete translation of the solution v N γ given by the expansion (6.5). Figure 5 shows outcomes of the numerical simulations for the random realizations of the discrete SHE model (5.1) with five different initial conditions shown on the top left panel.The drift generally occurs faster in the random model because the relevant small eigenvalue is of the size O(µ) if α 2 ̸ = 0, where µ is defined by (5.7).Only two time-independent solutions are stable if α 2 ̸ = 0 and this generic case is shown on the top right panel.The two stable solutions are given by (6.6) with some specific δ and they are related by the sign reflection of each other.The color scheme of the final states corresponds to that of the initial conditions and the missing colors correspond to the final states identical to those shown in red and yellow.
If α 2 is zero or close to zero, then four, six, eight, or ten time-independent solutions could be stable according to Theorem 5.1 and since half states are related by the sign reflection of the other half states, five numerical conditions could be used to detect all cases.However, the cases with more stable states become rare events.Our numerical simulations showed random configuration of the discrete SHE model with four and six stable states (bottom left and right panels, respectively) but we did not find random configurations with eight and ten stable states.The four configurations on the bottom left panel shows that two states are given by the sign reflection of the other two states (red versus yellow and blue versus cyan).The same is true for the bottom right panel (red versus cyan and black versus yellow) but there is one more stable state (blue line) which is not matched by the sign reflected state (since our simulations only involve five initial conditions).

Discussion
Dynamical principles underlying formation of patterns in complex systems are important for understanding a range of phenomena arising across multiple disciplines from morphogenesis to autocatalytic reactions, to firing patterns in neuronal networks [30].Traditionally, pattern formation has been studied in the continuous setting through the framework of reaction-diffusion partial differential equations.Recently, the interest has shifted towards understanding patterns in discrete systems, driven by the ubiquity of networks in contemporary science.Interestingly, already in his seminal paper on morphogenesis, Turing considered a reaction-diffusion model on a lattice, i.e., a discrete system [38].There is an extensive literature on lattice dynamical systems, which covers pattern formation and propagation phenomena [14].In the past decade, there have been many studies exploring pattern formation in complex networks including random networks.The key analytical challenge in dealing with this class of models, which was not present in studies of partial differential equations or lattice dynamical systems, is handling network topology, which can be random.Until recently most studies of the dynamics of complex and random networks had to resort to heuristic and numerical arguments (see [31,41,17,1,18] for a representative albeit limited sample of studies of Turing patterns in networks).
The situation has changed with the development of the theory of graphons [20].The use of graphons allows a rigorous derivation of the continuum limit for interacting dynamical systems on a large class of graphs including random graphs [23,24,26], which can be used effectively for studying dynamics on large networks.For example, the use of graphons led to the breakthrough in the analysis of synchronization and pattern formation in systems of coupled phase oscillators on networks [6,27], interacting diffusions on graphs [32,22], mean-field games [3], and graph signal processing [35,11].
Graphons provide multiple analytical benefits.Oftentimes, graph limits possess additional symmetries that are not present in the individual realizations of random graphs.For example, the graph limit of the small-world family of graphs used in the present work is isotropic, in contrast to the graphons corresponding to the realizations of small-world graphs.This symmetry enables the effective use of the Fourier transform for analyzing the limiting model (cf.Section 3).Likewise, the deterministic (averaged) discrete model (4.1) is shift invariant and can be studied using discrete Fourier transform (cf.Section 4).This in turn can be used to understand the dynamics of the random model (5.1).In general, the proximity in cut norm of a network to a symmetric network in can facilitate the analysis of the original nonsymmetric model.
In the context of the bifurcation problems considered in this paper, the relationship between the spectra and eigenspaces of the deterministic and random discrete models, as well as their continuum counterpart, is crucial.The theory of graph limits offers very efficient tools for tracking this relation.The convergence of kernel operators in cut-norm automatically implies the proximity of the corresponding eigenvalues and eigenspaces [37].This has significant implications for the analysis of dynamical systems.For instance, once the proximity in cut norm of the linear operators corresponding to the interaction terms of the random model and its deterministic counterpart was verified (cf.Appendix A), the proximity of the eigenvalues and the corresponding eigenspaces followed automatically.In contrast, the analysis in [2] based on operator norm topology and Davis-Kahane estimates requires substantial efforts.Importantly, the analysis in the present paper extends to sparse networks.
Our techniques apply naturally to other pattern-forming systems on random graphs, including Gierer-Meinhardt model [17], Mimura-Murray model of interacting prey-predator populations [31] and many other activator-inhibitor systems.An interesting area of potential applications are neural fields [9].Fourier methods played an important role in the analysis of pattern in nonlocal neural field models, which are very similar to the continuum limit analyzed in Section 3 [19].We expect that the methods developed this paper may lead to interesting results for discrete neural fields with random connectivity.Finally, the concentration estimate for the random linearized operator in Fourier space (cf.Appendix A) may be useful in graph signal processing.This completes the analysis of S 1 , the first term on the right hand side of (A.7).The remaining terms are analyzed similarly and result in either O (α n n) −1/2 bound as for S 2 above or in O (α 3 n n) 1/2 bound for S 1 .Thus, ∥∆∥ ∞ = O (α 3 n n) 1/2 as claimed.□ Remark A.3.Note that O (α 3 n n) 1/2 come from quadratic terms like S 2 .In the second order nonlocal spatial operator in (3.1) is replaced with the first order operator (3.2), as one encounters in the neural field type models, then the bound on ∥∆∥ ∞ can be improved to O (α n n) 1/2 .This means that the results of this paper would hold for α n = O(n −1 ), i.e., for graphs of bounded degree.

Figure 1 .
Figure 1.(a) W takes values 1 − p and p over the black and white regions respectively.(b),(c) Pixel plots of the adjacency matrix of Γ N W and its random counterpart ΓN W .

Remark 4 N
of (2N ) states obtained from either (4.5) or (4.6) with the discrete group of translations are the sign reflections of the other N states, in accordance with the sign reflection symmetry in Lemma 4.1.

. 1 )
where u ∈ C 1 (R, R Z N ) and ÃN W = (ã ij ) is the adjacency matrix of ΓN W .To make the setting for the model on a random graph consistent with the model on the deterministic Cayley graph, we have kept the periodic setting Z N in (5.1).The matrix ÃN W and the diagonal matrix DN W acts on components of u ∈ Z N inside [−N, N − 1] only.The same notations ÃN W and DN W are used to denote the linear operators on Z N and on [−N, N − 1].

Figure 3 .
Figure 3. Numerical solutions of the SHE on deterministic SHE with kernel (6.1) with parameters N = 400, p = 0.1, r = 0.2, and γ = 0.05 plotted in red and its random counterpart plotted in blue.Both models were initialized with the discretization of the leading order term in (6.3) and integrated for 100 units of time.The shift δ was set to −π/2 in plot a and to 0 in plot b.

Figure 3
Figure 3 presents results of numerical simulations of the discrete models derived from the nonlocal SHE with kernel (6.1) with parameters N = 400, p = 0.1, r = 0.2, and γ = 0.05.Both the deterministic and random models were initialized with the leading order term on the right-hand side of (6.3).The shift δ was set to −π/2 in plot a and to 0 in plot b.The two models were integrated for 100 units of time.The solutions of the deterministic models are plotted in red and the solutions of the random model are plotted in blue.

Figure 4 .
Figure 4. Numerical approximations of the discrete SHE model (4.1) for N = 5, p = 0.9, r = 0.2, and γ = 0.03.a Two initial conditions (red and cyan curves) and the outcomes of numerical simulations at t = 100 (blue and black curves).b Numerically obtained zeros of the two solutions versus time relative to the midpoint between two grid points (red dotted line).

Figure 5 .
Figure 5. Numerical approximations of random realizations of the discrete SHE model (5.1) for N = 5, p = 0.9, r = 0.2, and γ = 0.03.a Five initial conditions shifted relative to each other.b-d Outcomes of numerical simulations with two, four, and six stable solutions.