Finite to infinite steady state solutions, bifurcations of an integro-differential equation

We consider a bistable integral equation which governs the stationary solutions of a convolution model of solid--solid phase transitions on a circle. We study the bifurcations of the set of the stationary solutions as the diffusion coefficient is varied to examine the transition from an infinite number of steady states to three for the continuum limit of the semi--discretised system. We show how the symmetry of the problem is responsible for the generation and stabilisation of equilibria and comment on the puzzling connection between continuity and stability that exists in this problem.

For an overview of the use of (1) in materials science, see [10]. There are many papers dealing with the mathematical analysis of this equation, which examine existence and stability of travelling waves [3], the structure of the stationary solutions set [2], propagation of discontinuities [11], coarsening [9] and long time behaviour [15,19,24,20]. Note, in particular, that in [15] it is shown that if the diffusion coefficient ε is sufficiently large, a "Conway-Hopf-Smoller" type result holds: the only stable steady state solutions, say, in L ∞ (R), are the constant stable steady states of the kinetic equation u t = f (u). Thus, if we choose f (u) = u(1 − u 2 ), the stable states are u = 1 and u = −1. On the other hand, if ε = 0, (1) admits an uncountable set of equilibria: let X, Y and Z be any disjoint sets such that A ∪ B ∪ C = R, then a function u(x) that is equal to 1 on X, −1 on Y and 0 on Z is a steady state solution. Note that if Z = ∅, all the resulting equilibria are stable in L ∞ (I). Furthermore, it is shown in [9] that there exists an ε 0 > 0 which depends on the kernel J ∞ , such that for all 0 < ε < ε 0 the set of steady state solutions of (1) is in one-to-one correspondence with the set of equilibria of u t = f (u). Hence, in view of the above, it is of interest to perform a bifurcation analysis of the set of steady states of (1), as we decrease ε from some initially large value to zero, and investigate the transition from a finite to infinite set of solutions.
To the best of our knowledge, such a study has not been performed before. The object of this paper is precisely such a study of the spatially discretised version of (1). For simplicity, here we restrict ourselves to 1-periodic patterns.
If we choose spatially one-periodic initial data u(x, 0), then from (1) it is clear that for all x ∈ R and t ∈ R + u(x, t) = u(x + 1, t).
Then from (1) we have where and x ∈ [0, 1]. Thus, for 1-periodic initial data we only need to solve the problem (1) on the interval Ω = [0, 1] with the kernel J(x). For the kernel given by (2), J ∞ (x) and J(x) are plotted in Figure 1.
For J defined by (6) the following two properties hold.
dx. Property 1. above has an important influence on the spectrum of the matrix governing the semi-discretised version of (1) as we explain in the next section. From now on we work on [0, 1] and use the kernel J given in (6)

The semi-discretised system
We discretise in space using piecewise-constant functions [9] and collocating at the uniformly spaced element mid-points, x = x j+ 1 2 , j = 0, 1, 2, · · · , N − 1. Setting u j = u(x j+ 1 2 , t), we have the semi-discrete approximation of (1) given by where now u(t) ∈ R N , supplemented with some initial condition u(0) = u 0 ∈ R N . The nonlinearity F : R N → R N is given by F j (u) = f (u j ). It remains to specify the N × N matrix A N . If we put h = 1 N , its elements are given by From Lemma 1 it follows that A N is a symmetric circulant matrix generated by the elements a 1,1 , . . . , a 1,N . Hence the theory of circulant matrices can be used to characterise its spectrum precisely. Let W k be the N distinct roots of z N − 1 = 0, so W k = exp i2πk N , for k = 0, 1, 2, ..., N − 1. Then the following theorem holds: . Let A N be the circulant matrix defined by a 1,1 , a 1,2 , ..., a 1,N . Then −A N is diagonalisable with eigenvalues with corresponding eigenvectors v k = 1, Let us see what this implies in our case for the spectrum of the discretisation.
Lemma 2. The following three properties hold for the spectrum of A N 1. λ 0 = 0; 2. λ k = λ N −k ; 3. Let I N be the convex hull of the set of non-zero eigenvalues of A N . As N → ∞, I N converges in the Hausdorff metric to the set Before we prove this lemma, let us explain what it means. First of all, we must have a zero eigenvalue with a constant eigenvector, because, like in the case of the Neumann Laplacian, the equation conserves mass. Secondly, the pairing of the eigenvalues is simply the consequence of the symmetry J(x) = J(1 − x) inherited from the evenness of the kernel J ∞ . Finally, the third part of the lemma implies that as N → ∞, the spectrum accumulates at the point 1 0 J(x) dx. Note that in the case of J ∞ (x) = 100/π exp(−100x 2 ), we explicitly have Proof. 1. From (8), putting W 0 = 1, we immediately obtain from (9) that λ 0 = 0. 2. From part 1. of Lemma 1 it follows that for all j = 2, . . . , N , so that the matrix A N is symmetric. Hence its eigenvalues λ k are real. But then taking complex conjugates of 3. Finally, by taking the limit as N → ∞ in (9) we immediately obtain that Our aim is to examine bifurcations in this system and, below, we perform a numerical pathfollowing of solution branches. Some of these will, by symmetry, arise in pitchfork bifurcations from the trivial solution u = 0. Here we examine analytically the values of ε where such bifurcations may occur in the semi-discrete system and later we can compare to the numerically found values. Linearising around the zero solution, we have the eigenvalue problem and hence bifurcations from the zero solution will only occur if µ = 0, or in other words, if Thus, for the semi-discrete system (7) we can fully characterize the values of ε where bifurcations of the zero solution occur, namely For example, for N = 32, J ∞ (x) = 100 π exp(−100x 2 ) and f (u) = u(1 − u 2 ), we have using (9), the results of Lemma 2 and the formula (11) that bifurcations from the zero solution are expected at the values of ε as in Table 2. Note that for this case of N = 32, the value of ε 1 agrees to 12 decimal points with the limiting value of ε 1 , 1/(1 − exp(−π 2 /100)) as N → ∞, (see part 3 of Lemma 2). ε 1 = ε 31 10.6403416996149 ε 9 = ε 23 1.00033746722800 ε 2 = ε 30 3.06584313146254 ε 10 = ε 22 1.00005172586163 ε 3 = ε 29 1.69885748860222 ε 11 = ε 21 1.00000650971543 ε 4 = ε 28 1.25968856776757 ε 12 = ε 20 1.00000067252291 ε 5 = ε 27 1.09266327932320 ε 13 = ε 19 1.00000005703325 ε 6 = ε 26 1.02948119722480 ε 14 = ε 18 1.00000000397031 ε 7 = ε 25 1.00800141737791 ε 15 = ε 17 1.00000000022729 ε 8 = ε 24 1.00180943793728 ε 16 1.00000000002128 Table 1: For N = 32 Bifurcation values in terms of ε of the zero solution.
Let us examine the eigenvectors of −A N in some more detail. Since both v k and v N −k are eigenvectors, we immediately have that Re (v k ) and Im (v k ) are eigenvectors. Define the cyclic This follows since if v is an eigenvector, then so is e 2iπ/N v.
Remark. In the above argument, we can pass to the limit as N → ∞ and arrive at the somewhat startling conclusion that cos(2πkx) and all their translates are eigenfunctions of is as long as it has the right symmetry property. Of course, cosines are also the eigenfunctions of the Neumann Laplacian. It is very pleasing to obtain such a result via a semi-discretisation. Finally we note that fixed points of the semi-discrete problem satisfy Thus at ε = 0 stable solutions are given by 1]. Unstable solutions at ε = 0 are given by where X ∪ Y ∪ Z = [0, 1] with some nonempty Z.
We use this to define solutions with different numbers of interfaces. When X = [0, α) and That is loosely speaking we have at ε = 0 one jump in the solution upto cyclic shift. Two-interface, three-interface solutions, etc., are defined similarly. Thus, for example, the branch of solutions corresponding to orbit A in Table 2 are of one-interface and those corresponding to E are of three-interface.

Results
We take for our computations the kernel function with f (u) = u(1 − u 2 ) and vary the parameter ε. For small values of N it is possible to enumerate all possible solutions of the semi-discrete system (12) with ε = 0 and to analyse their continuation to ε > 0 using the theory of bifurcation with symmetry. This we do below for N = 4 and these analytic results were used to check the validity of our numerics. We implemented in Matlab a standard pseudo arc-length continuation algorithm with step size control as described in [14,22,13] for the discrete problem (12). Since A N is a circulant matrix, we take advantage of reducing storage costs as the full information of A N can be obtained storing one row or column only, see [4] and references therein. Furthermore the use of the FFT for each matrix vector multiplication reduces the computational cost. We detect bifurcation points by observing where eigenvalues of the Jacobian J = D u F of the nonlinear system F (u, β) = 0 cross the imaginary axis and perform branch switching at those points by perturbing in the direction of the associated eigenvector. The arc-length ℓ of u(x) ∈ C 1 (R) is defined in the standard way and we approximate the arc-length of u(x) with the mid-point rule and using the standard forward difference approximation for the derivative. With a uniform discretization we get where h = x j+1 − x j and u j ≈ (x j ). Note that although ℓ only makes sense for u ∈ C 1 however we can evaluate ℓ h even when u is discontinuous at grid points. Then, for N = 32 we compute the bifurcation diagram numerically and gain insight into the structure of the bifurcation diagram of the original continuous problem. Finally, we examine the large N limit and formulate the results of the numerics as two conjectures concerning the interplay of continuity and stability and the behaviour of saddlenode bifurcations as α → 1/2. For the continuous system, the symmetry group is O(2) × Z 2 , and so for a finite number of nodes N , we use Γ N = D N × Z 2 equivariance structure [12,18].

The N = 4 case
Inverses of the nonzero eigenvalues of the 4 × 4 matrix A 4 are {91.82, 183.63, 183.63}, so we expect primary branches to bifurcate from the zero solution at those values of ε. Note that all primary branches have zero mean, but the converse is not true. Since here we know all the solutions at ε = 0 and their stability, and since symmetry properties are conserved on primary branches, we can cut down the work considerably by looking only at orbits of solutions under Γ 4 . In the table 2, we collect all the orbits, their lengths and the corresponding isotropy subgroups Σ x . There, f stands for the group generated by f ∈ D 4 × Z 2 . Now we can immediately draw the bifurcation diagram using the following three rules [12,18]. First a bifurcating branch must have the isotropy subgroup which is a subgroup of the isotropy subgroup of the primary branch; secondly dimensions of unstable manifolds have to match at a bifurcation point to satisfy the principle of exchange of stability, and thirdly at ε = 0; the number of nodal domains must increase from one bifurcation point to the next. With these rules there is only one way to construct the bifurcation diagram; see Figure 2 (a) and (b), where the y-axis is not to any scale, and is only intended to make clear the end-points of various branches at ε = 0. These figures show the bifurcation structure arising from bifurcations of the zero solution.
We would like to make the following observations. The stable non-zero-mean branches corresponding to the orbit K have to arise through a saddle-node bifurcation. Numerically, this happens at a value of epsilon ≈ 49.294 that is smaller than the value ε 4 = 122.432 at which the branches of the orbit A become stable, see Figure 3 which shows the numerically     Figure 3: Numerically computed bifurcation diagram from the zero solution from the first (a) and second (b) bifurcation points. Solid lines represent stable and broken lines unstable unstable solutions respectively. We indicate above each branch the dimension of the unstable manifold. We show the homogeneous solutions or the orbits J and K connected by a saddlenode bifurcation on (b). Compare to the theoretical prediction in Figure 2. computed diagram. We will see the equivalents of these statements in higher dimensional discretisations. Finally, we did not perform a Liapunov-Schmidt calculation to determine the order of bifurcations at the double eigenvalue point ε = 183.63, but the opposite assignment of stabilities cannot be reconciled with the above rules of bifurcation.

The case of N = 32
Though an analysis similar to that in the case of N = 4 can be attempted here, the numbers of orbits are astronomical, and we rely on our numerical continuation method, the results of which match exactly the predictions of the analysis in the case N = 4. In Figure 4 we plot in (a)-(d) sample solution branches of the bifurcation diagram with N = 32. If we start with a large value of ε we see in (a) and (b) the first bifurcation arises at ε ≈ 10.64 as predicted by the theory in Table 2. In (a) we show the continuation of the zero mean one-interface which undergoes a pitchfork bifurcation at ǫ ≈ 2.012. In (b) we have plotted the one-, three-, five-and seven interfaces and their stabilization. In (c) we show details of the bifurcation structure close to the pitchfork at ε ≈ 2.012 (note for clarity one branch of the pitchfork seen in (a) is not plotted). As α → 0.5, the saddle-node bifurcation points converge to ε 3,2 . This structure is repeated for the other n-interface solutions and is illustrated in (d) for the three-interfaces solutions. Here we see that the zero-mean one-interface solution branches stabilize at ε 3,2 ≈ 2.012. Below we will formulate a conjecture concerning the limiting value which we call ε 0 1 at which the one-interface branch with zero-mean stabilizes as N → ∞.

The limiting problem and conclusions
It is not hard to prove (see for example [3]) that if ε > 1, steady state solutions of (4) are continuous, since the function −εu + f (u) is monotone. Hence it is interesting to understand when the solutions lose continuity (certainly, for ε = 0 there are no non-constant continuous solutions).
The non-trivial stable one-interface zero-mean solution branches (α = 0.5) that originate at ε = 0 were investigated in detail as we change N . If we define M by then for a C 1 function this converges to max x∈[0,1] |u x | and so we can identify where the solution is continuous. Figure 5 plots in (a) M against ε along a branch of one-interface zero-mean solutions for N = 2 p , p = 4, 5, 6,7,8,9,10. If we let ε d 1 be the value of ε at which this branch of solutions becomes discontinuous, then this figure suggests ε d 1 = 1. This is supported in (b) which shows for different ε convergence of the derivative M with N on a log log scale. Furthermore the loss of continuity appears to coincide with the loss of stability. In Figure 6 we show numerically that the bifurcation values converge to ε 0 1 = ε 0 2 = ε 0 3 = 1 as N → ∞. Now we can collect our observations and form two conjectures. First we consider the zeromean interface branches. Let ε 0 m be the value at which the zero-mean 2m−1-interface solution becomes stable. Let ε d m be the value of ε at which this branch becomes discontinuous. Then we have Conjecture 1: ε 0 m = ε d m . We can prove a very weak form of this conjecture for m = 1. From the results of [3] it follows that discontinuous stationary solutions will exist for any ǫ such that the function is non-monotone. On the other hand, from Theorem 2.1 of [6] it follows if g(u) is monotone, there are no nonconstant minimizers of the energy functional (3). Hence we have ε s 1 ≤ ε d 1 . However we do not have the inequality the other way. We now consider the saddle-node bifurcation of the non-zero mean interface solutions. Now, let u s be a branch of 2m − 1-interface stable solutions of (4) with mean s, and let ε b,s m be the value of ε at which the saddle-node bifurcation giving rise to the branch occurs. Then we have Conjecture 2: lim s→0 ε b,s m = ε 0 m . These two conjectures, if true, would lead to the bifurcation picture sketched in Figure 7. In (a) we plot the zero-mean one-interface branch and have indicated the continuum of saddlenode bifurcations ε s 1 that approach the bifurcation at ε s 1 = ε d 1 = 1. In (b) we indicate the first four branches of the infinite number that bifurcate from zero, the branches of associated saddle-node bifurcations and here we have that lim s→0 ε b,s m = ε 0 m = ε d m . In addition our numerical investigation seems to indicate that ε s 1 = ε d 1 = 1 = ε s 2 = ε d 2 . Finally let us consider the stable solutions -that is the solutions we expect to see from any simulation. Thus we have for ǫ > 1 two stable solutions, then a region of parameter space with an infinite number of stable solutions of one and three interface type, then a region of parameter space with one, three and five interfaces and so on. In conclusion the diffusion coefficient ε determines the number and type of stable solutions.     We show the primary branch and the continuum of saddle-node from the saddle-node bifurcations. In (b) we indicate that this structure is then repeated for the three-interface, five-interface,... solutions. The solid circles represent lim s→0 ε b,s m = ε 0 m = ε d m . The bifurcations for the one and three interface solutions both occur at ε = 1.