Global stability of stationary patterns in bistable reaction-diffusion systems.

We study a piecewise linear version of a one-component, one-dimensional reaction-diffusion bistable model, with the aim of analyzing the effect of boundary conditions on the formation and stability of stationary patterns. The analysis proceeds through the study of the behavior of the Lyapunov functional in terms of a control parameter: the reflectivity at the boundary. We show that, in this example, this functional has a very simple and direct geometrical interpretation.


I. INTRODUCTION
The subject of self-organization or nonequilibrium pattern formation has captured the attention of researchers during more than a decade, and is by now one of the most active fields in the physics of complex systems [l]. The extremely rich variety of nonequilibrium systems that one can consider calls for different descriptions.
Among them, the reaction-diffusion (RD) approach has shown to be a very fertile source of models for interesting phenomena in natural and social sciences, where structures can arise and last for longer or shorter periods of time according to their degree of stability [2]. In particular, the possibility of a solitary pattern without propagation in an infinite medium was established [3].
For oneand two-component systems, boundary conditions (BC's) have been recently shown to play a relevant role in the formation and stability of patterns [4,5]. We are concerned here with the role of BC's in pattern selection, and particularly with the global stability of the resulting structures. As it was indicated in Ref. [4], the choice of Neumann or Dirichlet BC corresponds physically to complete reAection or absorption at the boundaries, respectively. The more realistic case of a partially reflecting (or absorbing) boundary is adequately taken into account by the albedo BC [4,5].
The specific model we shall focus on belongs to a family of one-component models with a broad range of applications [7]. These are so-called bistable reaction-diffusion models, whose general formulation is =DV T+F(T), Bt where D is the difFusion coeKcient, T(z, t) is a real field that represents the magnitude of interest, and the nonlinear term F(T) is the source (or dissipation) term. In the bistable case, F(T) is generally a cubiclike function (with three roots, two stable and one unstable fixed points).
As in Ref. [4,5], we shall concentrate on a bounded one-dimensional system and impose the most general linear homogeneous BC, namely, the albedo BC. These relate the normal derivative of the field with its value at the boundary I, through an albedo parameter k: Their physical meaning is that the boundary acts as a partially absorbing or re6ecting medium. They have Neumann's BC (totally refiecting, i.e. , k~0) and Dirichlet s BC (totally absorbing, i.e. , k~00) as limits, and can describe most real situations in a more accurate way.
In this paper we want to present an analysis of the effect of the BC on the formation and stability of stationary patterns through an overview of the global stability of the structures. In order to reach such a goal we shall profit by using the concept of Lyapunov functional and through its study get information on the form that the albedo parameter controls the shape of the patterns and their global stability. This kind of approach, based on the use of the concept of the Lyapunov functional, has not been exploited in the realm of reaction-diffusion systems because it is usually not possible, insofar as some potential conditions are not fulfilled, to obtain a Lyapunov function in a general problem. However, Graham and collaborators [6], who have been pioneers in introducing the concept of a nonequilibrium potential, have investigated the form of such Lyapunov-like functionals in several problems of systems far from equilibrium, such as spatially extended systems described by the complex 1063-651X/95/52(1)/129 (8) and about global stability even though the system does not fu1611 the above potential conditions. In this work, as we study a one-component system, the form of the functional is well known [7]. Nevertheless, we perform such an analysis in order to gain insight into the way that the variation of the albedo parameter changes the value of the Lyapunov functional, altering the relative stability by including changes between the stable and metastable states. The stable states will correspond to attractors (minima) of the functional while the unstable ones will be saddle points, defining the barrier height between attractors.
Those forms, as well as their first derivatives, must be matched at the spatial coordinates of the boundaries between the activated and unactivated regions, that we denote generally by z". Through that matching procedure we get the general solution for the stationary case.
In order to identify the matching points z", we have to solve the transcendental equations given by the conditions T(z") = T,h. These result in self-consistent implicit equations for z" in terms of k.
The local stability can be analyzed by means of the linearized equation for the perturbation, in a standard way. The eigenvalues of this equation allow us to distinguish the unstable patterns from the stable ones. Writing the solution as II. THE BALLAST RESISTOR T(z, t) = T,"(z)+P(z)e (4) We particularize our analysis to the ballast resistor model [8,9]. The ballast resistor has been used as a current stabilizer for a long time. It consists of a wire (or a thin ribbon) through which an electric current flows, immersed in a thermal bath whose temperature is taken as reference, i.e. , Tz =0. For the sake of concreteness we restrict ourselves here to its superconducting version, described by the so-called "hot-spot" model and used as a model for thermal transfer in superconducting microbridges [10]. Below the critical temperature T, )0 there is no Joule dissipation; above T" the latter occurs at a basically constant rate (independent of T). The piecewise linearization results from the above-mentioned characteristics, as a natural approximation.
After a suitable variable rescaling, the energy-balance equation can be written in the form aT= 2 Bt =7' T -T+8(T -T ), ch where 6 is the Heaviside function, and T, i, = T, /Th, T& )0 being the ratio between the dissipation rate and the transfer of heat to the bath. T(r, t ) is the scaled temperature field, which determines the (nonequilibrium) state of the system. It turns out [9] that, when T, h & 1, the system has only two uniform stationary states: T =0 and T =1. When albedo BC's are imposed [4], only T =0 survives as a stationary and homogeneous solution. There are however inhomogeneous stationary solutions to Eq. (3) where two kinds of regions coexist: those where T )T, t, (the system is said to be activated), and those where T & T,h (the system is called unactivated). In this work we restrict ourselves to the one-dimensional case (wire) and consider the system to be bounded in the interval -zL (z (zL. In the near future we shall consider other geometries. We keep the left end either isolated (Neumann BC) or in contact with the heat bath (Dirichlet BC) and assume that the right end is partially isolated, so that the heat Aux at this end is partially refiected (albedo BC).
Equation (3) can be solved, in the stationary case, following the technique already discussed in Ref. [4]. DifFerent  does not distinguish between the stable and metastable states. This is an important question when there are two or more locally stable structures for a given BC, as in our case. In order to perform a global stability analysis, it becomes necessary to write Eq. (3) in a potential form: where 7 is the Lyapunov functional [7], defined as (7) and satisfying We see that V vanishes for stationary distributions and is negative in any other case. Every stable stationary distribution corresponds to a minimum of this functional.
During its temporal evolution, V decreases until it reaches one of its minima. Hence, the evolution of this system consists in the approach to one of its stable stationary distributions.
9' is of great usefulness in discussing the global stability of the patterns. The unstable structures are related to extrema which are not minima of x In the present case, 9 is defined by and we need to impose the Albedo BC at +zL. Integrating by parts the kinetic term in Eq. (7), we find that where the terms upper and lower refer, respectively, to the areas of the activated region above and below T, h (see Fig. 1). In two-region patterns, that region is located between z, and zz, and in three-region patterns, it is located between z, &and z,2. Despite the fact that its validity seems to be restricted to the present model, this is a remarkable result: the value of V is uniquely determined by this area difference, which is easily calculable.
As it is to be expected, linearly stable solutions are relative minima of X The stable stationary solution corresponds to the absolute minimum of V. The other linearly stable solutions correspond to metastable states. The dependence of V upon k will give us information about the way in which the boundary conditions affect the global stability of the structures.
The Lyapunov functional so obtained can be explored in the neighborhood of a reference state. We particularize the analysis to two-region patterns (to be discussed later). In order to analyze the nature of the resulting solution we make a functional expansion of Eq. (8) around a stationary solution T"(z): leading finally to + Zl 7[T"(z)+P(z, t) ] = 2[ T"(z)] -f P dz, (14) which shows explicitly that the unstable patterns are associated with saddle points in the functional space. It is the eigenfunction in the unstable manifold that produces the increase of the functional, showing the saddle-point character. The value of V at the unstable patterns gives the barrier height between the uniform and nonuniform attractors (the separator between the respective attraction basins). This information is extremely relevant in order to study the effects of Auctuations on this type of systems [6,11].

IV. STATIONARY STRUCTURES AND THEIR STABILITY
Here we discuss the formation and global stability of several types of stationary patterns. We focus our attention on simple patterns. Five cases will be considered: patterns of two regions with an activated zone on the right and (a) a Dirichlet BC or (b) a Neumann BC at -zL, patterns of three regions with a central activated zone and (c) a Neumann BC or (d) a Dirichlet BC at zt ', and (e) periodic structures. We impose the albedo BC at the right end in all the cases except (e), where the BC's are naturally Neumann BC's in both boundaries.
After some algebra, the general solution in the two-region case can be written in terms of the values of the temperature profile and its derivative at z =z"T(z,), and T'(z, ), as P'/T, ', (z, ) .
The albedo BC at zL gives the relationship between k and T, '. The method can be readily generalized for structures with many zones.
A. Two-region patterns In the case of two-region patterns, there is a maximum allowed value of k, beyond which three-zone structures appear. In Fig. 2 the behavior of z, vs k can be seen for the patterns arising in case (a) (the curves are parameterized with the value of T,I, ). For small k values there are two branches of z, (Td"k) and therefore two nonuniform structures exist, associated with this boundary condition.
The patterns corresponding to the lower branch have a large activated zone. T,I, appears as a control parameter rendering qualitative changes in the shape of the branches. For small values of T,h, they are "open;" the situation changes when T, h increases, and the branches are "closed. " It is stability which selects between these branches. As we shall see later, the lower branch turns out to be the stable one. Figure 3 shows the forms of typical patterns for this case. Curves 1 and 2 correspond, respectively, to a lower and upper branch of z, ( T,h, k). The behavior of z, vs k for the patterns in case (b) is depicted in Fig. 4. We see that there exists only one branch for each value of T,I, and k. Figure 5 [4,5], the branches lying above the marginal-stability line (determined by the vertices of the closed curves) correspond to unstable stationary solutions, whereas those 2.00 Here also a maximum allowed k exists, but in contrast with case (a) the curves are single valued. The values of T,z for each curve are (1) 0.045, (2) 0.225, (3) 0.315, (4) 0.405, (5) 0.45, (6) 0.5, (7) 0.55, and (8) (2) 0.34, (3) 0.382, (4) 0.404, (5) 0.42, (6) 0.43, (7) 0.435, (8) 0.44, (9) [4,5], we can infer that the stable patterns are associated with the existence of an important activated region, which is equivalent to large dissipation. Comparing the areas above and below the T, & level in the activated region, we see that stable patterns are associated with a larger upper area.
If we remember that T =0 is a stable stationary homogeneous solution, it is necessary to make a global stability analysis in order to resolve stable from metastable states. value of z, . The same feature for the "closed" branches can be seen in Fig. 7, where the cusp corresponds to the marginal stability point. The global analysis confirms the results obtained by linearization and, in addition, allows us to distinguish the stable from the metastable patterns.
In fact, the inhomogeneous stationary structure is the stable solution for the whole physical k range, for values of T,z lower or of the order of 0.32 (see Fig. 6). This situation changes when T,z increases, and there are values of T,h and k such that T=O is the stable state and the nonuniform patterns are metastable (see Fig. 7). We see that k appears as a control parameter whose variation originates a change of the relative stability between locally stable structures. In Fig. 8 we show the behavior of the Lyapunov functional vs k for patterns in case (b). The curves are monotonically increasing with k, so the unstable patterns will become more unstable as k grows.

B. Three-region patterns
In the case of three-region patterns, there are no limitations on the value of k. As already stated, in the limit k -+Do, we recover the Dirichlet BC. Three-region pat- FIG. 7. The value that the Lyapunov functional takes on the stationary patterns of case (a), plotted against k, for the "closed" branches in Fig. (2). The values of T,h for each curve are (1) 0.44, (2) 0.45, (3) 0.46, and (4) 0.47. terns have two matching points indicating the limits between the activated and the unactivated regions. In Figs. 9 to 12 we show the behavior of z, vs k, together with typical patterns for cases (c) and (d), respectively. Figure   9 shows that there exists only one pair of matching points for each value of k and T,&, that is, only one nonuniform stationary solution associated with the external conditions (the same result was obtained for two-region patterns under an identical BC). These structures can only exist in a finite interval of k. Figure 10 shows typical patterns for this case. In Fig. 11, the behavior of z, vs k is very different: we see that there exists a region where for each value of k and T, & there are two pairs of values (z, &, z,z) (that is, two nonuniform stationary patterns).
(ii) Although twoand three-region patterns cannot coexist for the same values of T, I, and k, it can be seen that as k decreases (for fixed T, s ), the associated patterns become more stable. As a consequence, two-region patterns associated with the Dirichlet-Neumann BC result to be the most stable. Also, the "barrier" between these attractors and the homogeneous one ( T =0) is the largest.

C. Periodic patterns
As was indicated in [9], the system has also analytic periodic solutions as long as the applied BC's (the same at both ends) are not Dirichlet BC's. Their local stability was analyzed in Ref. [9], finding that all these structures are unstable. With the aim of exploring whether the number of nodes of the proposed solution has any effect on the stability of the periodic patterns, we have analyzed the Lyapunov functional. Figure 15 shows the behavior of V as a function of the number of nodes. A saturation effect appears as the number of nodes increases, for a fixed length. In all the cases, the value of the functional remains positive, showing the unstable character of this type of solution.

V. FINAL REMARKS
In this paper we have analyzed the effect of the boundary conditions on the formation and stability of stationary patterns and, in particular, the global stability of the resulting structures. We want to remark that the albedo BC has a very important inQuence on the formation and global stability of dissipative structures, which has been overlooked until recently [4,5,12]. The albedo parameter controls the shape of the patterns and their global stability. Its variation changes the value of the Lyapunov functional, and then induces changes between the stable and metastable states (it alters the relative stability). The unstable structures correspond to saddle points of the functional, and define the barrier between attractors. The numerical exploration of the extrema allows us to identify the nature of the solutions as an alternative to the linear stability analysis. This line of thought has been pioneered by Graham and collaborators [6], who have introduced the concept of nonequilibrium potentials. They have investigated the form of such Lyapunov-like functionals in several problems. The present work is a first step towards bringing these concepts into the realm of reaction-diffusion systems.