Self-replicating spots in reaction-diffusion systems

In this article we discuss the self-replicating spot, a p ticlelike phenomenon that occurs in reaction-diffusion ~RD! systems@1#. The spots consist of localized regions in whi the concentrations of the reactants differ from the surrou ing concentration fields. They grow, reaching a critical s at which time they divide in two. The two resulting spo again grow and divide. This process, which is visually sim lar to cell division, continues indefinitely. The long-time b havior depends on the precise values of the external con parameters, but typically consists of a chaotic ‘‘soup’’ which many spots compete for resources as illustrated in 1. Those spots that find adequate resources continue to and divide. Those that are unable to find adequate resou decay into the background. The spots observed in @1# were found during an attempt to model labyrinthine patterns in ferrocyanide-iodate-sulfate~FIS! reaction @2#. Since then, replicating spot patterns have been observed both num cally and experimentally in the FIS reaction @3,4#. There are obvious differences in the Gaspar-Showa model of the FIS reaction @5# and of the models discussed b us and others. This fact suggests that replication is a gen feature characterizing a broad class of reaction-diffusion s tems. In@6#, we presented some arguments in support of proposition. These arguments included both a heuristic scription of the process of replication and demonstrations analytic features common to several related model RD s tems. It turns out that replication is more general than analysis accounts for. Nevertheless, we think it worthwh to spell out the details of the theory presented in @6#. We remark here that various aspects of the replica phenomenon have been discussed by other authors @7,8#. Kerner and Osipov have a large body of work on larg amplitude dissipative structures including an analysis of static division of one-dimensional pulses as the system is changed. Gurevich and Mints have a body of work replication of thermal hot spots in composite supercondu ors. In the article by Petrov, Scott, and Showalter replicat


I. INTRODUCTION
In this article we discuss the self-replicating spot, a particlelike phenomenon that occurs in reaction-diffusion ͑RD͒ systems ͓1͔.The spots consist of localized regions in which the concentrations of the reactants differ from the surrounding concentration fields.They grow, reaching a critical size at which time they divide in two.The two resulting spots again grow and divide.This process, which is visually similar to cell division, continues indefinitely.The long-time behavior depends on the precise values of the external control parameters, but typically consists of a chaotic ''soup'' in which many spots compete for resources as illustrated in Fig. 1.Those spots that find adequate resources continue to grow and divide.Those that are unable to find adequate resources decay into the background.The spots observed in ͓1͔ were found during an attempt to model labyrinthine patterns in the ferrocyanide-iodate-sulfate ͑FIS͒ reaction ͓2͔.Since then, replicating spot patterns have been observed both numerically and experimentally in the FIS reaction ͓3,4͔.
There are obvious differences in the Gaspar-Showalter model of the FIS reaction ͓5͔ and of the models discussed by us and others.This fact suggests that replication is a generic feature characterizing a broad class of reaction-diffusion systems.In ͓6͔, we presented some arguments in support of this proposition.These arguments included both a heuristic description of the process of replication and demonstrations of analytic features common to several related model RD systems.It turns out that replication is more general than our analysis accounts for.Nevertheless, we think it worthwhile to spell out the details of the theory presented in ͓6͔.
We remark here that various aspects of the replication phenomenon have been discussed by other authors ͓7,8͔.Kerner and Osipov have a large body of work on largeamplitude dissipative structures including an analysis of the static division of one-dimensional pulses as the system size is changed.Gurevich and Mints have a body of work on replication of thermal hot spots in composite superconductors.In the article by Petrov, Scott, and Showalter replication is described in a manner qualitatively consistent with our analysis.The article by Doelman, Kaper, and Zegeling is complementary to and is of the greatest relevance to the present analysis.They use geometric singular perturbation methods to prove the existence of ''a plethora of periodic stationary solutions'' to the Gray-Scott model.They also prove the nonexistence of any rigidly traveling spot solutions: traveling spot solutions deform as they travel.This is consistent with our analysis in which the traveling spot solutions ultimately undergo an instability that leads to replication.
Our main goal here is to present in detail and to generalize the asymptotic analysis that was sketched in ͓6͔.Section II introduces the Gray-Scott model and expands upon the qualitative picture of spot replication presented in ͓6͔.Sec-FIG. 1. Snapshot of the replicating spot phenomenon in two space dimensions ͑from two-dimensional simulations described in ͓1͔͒.
tion III is devoted to developing in detail the onedimensional asymptotic analysis of the Gray-Scott model that we sketched in ͓6͔.Section IV will generalize our scaling results to kinetics other than the ones considered in detail in Sec.III.In Sec.V we summarize the main points of this paper.

II. THE GRAY-SCOTT MODEL
Both the asymptotics and the qualitative arguments will be developed in relation to the Gray-Scott model.The model ͓9͔ is given by ‫ץ‬u ‫ץ‬t ϭٌ 2 uϪuv 2 ϩA͑1Ϫu ͒, ͑1͒ ‫ץ‬v ‫ץ‬t ϭ␦ 2 ٌ 2 vϩuv 2 ϪBv.
Here u(x,t) and v(x,t) are scalar fields representing the concentrations of two chemical species.The ratio of their diffusion coefficients is ␦ 2 and A and B are parameters describing a feed from an external reservoir with the fixed concentrations uϭ1 and vϭ0.We now present the heuristic picture of spot replication that describes the behavior observed in Eq. ͑1͒.
We begin by considering a region of ground where some flammable liquid fuel (u) is continuously seeping in from a reservoir maintained at unit concentration.This effect is modeled by the term A(1Ϫu) in Eq. ͑1͒.This diffuses rapidly relative to the temperature v.If the fuel is depleted locally it will diffuse in from the sides to bring the level back up.We refer to this effect as the ''lateral'' or ''diffusive'' feed.
We now consider the effect of increasing the fuel's temperature v, which at equilibrium is a constant, say, vϭ0.The fuel seeps up so slowly ͑i.e., A is small͒ that if the entire domain is ignited the fuel will flare up and burn out.͑A mathematical model would have only one fixed point: fuel concentration equal to 1 and temperature equal to 0.͒ The kinetics of the fuel and fire are excitable, i.e., if the fuel is warmed by the sun, it will not do anything special, relaxing back to its equilibrium value when the sun goes down, but if it is perturbed by a blowtorch, we expect it to start burning and the temperature to increase ͑which is to say that the system exhibits autocatalytic behavior͒.Such behavior is, for example, given by the local kinetics u ˙ϭϪuv 2 ϩA(1Ϫu), v ˙ϭuv 2 ϪBv.
Finally, any fire that starts will spread diffusively across the field.We will require that the diffusion constant for this process is much smaller than that of the fuel spreading; for this reason, we include the small parameter ␦ 2 in front of the Consider what happens when a localized region of the fuel field is ignited.The fire will slowly begin to spread, until it is consuming all the fuel it can get.As the fuel is depleted in the burning region, it creates a gradient in the fuel concentration, which results in the fuel diffusing into the burning spot from the sides.Note that the lateral feed is essential since the external feed is not sufficient by itself to keep the fire burning.In essence, the nonburning regions around the spot act as integrators of the feed, channeling the collected fuel into the spot.Thus every spot has a quiescent zone around it, which collects the fuel required to keep it burning.
Consider what occurs as the strength of the external feed is varied.For low feed strengths, the spot remains localized, consuming all the fuel that diffuses into it.A finite region can support several spots provided the distance between them is great enough that they can collect sufficient fuel.As the feed strength is increased, an isolated spot increases in size in order to consume all of the fuel it receives.As the spot widens, it will eventually reach a point where the lateral feed is inadequate to maintain its center in the hightemperature state.In this situation, the temperature at the center of the spot will drop, leading to a double-bumped spot.This structure may be thought of as two nascent spots, each with an asymmetric flux distribution.The spot either moves towards regions of higher flux until the flux distribution is equalized or increases to the point where an additional replication occurred.
This model gives, within a reaction diffusion framework, the classical mechanism of size selection through competition between the rates of growth of the surface area and the volume ͓10͔.This mechanism is related to the ''activatorinhibitor'' model ͓11͔, with the fuel playing the role of inhibitor.It is important to note, however, that this is not a Turing instability, but rather a different type of structure in an excitable system.We note that this description is not married to any particular set of kinetics; any kinetics that provide these qualitative features should do.We will in fact demonstrate replicating structures in other models.

III. ANALYTIC SOLUTIONS IN THE LIMIT ␦Ӷ1
We now turn to analytic solutions of the system ͑1͒ in the limit ␦Ӷ1 that correspond to the evolution observed in Figs.
2 and 3. We can see from these figures that the spatial domain is divided into ''inner'' or spot regions where v is large and u is small and ''outer'' regions where u and v are closer to their fixed point values ͑uϭ1, vϭ0͒.In the outer regions, all the extrema of u and v are maxima and minima, respectively.Thus all the maxima of v and minima of u occur in the inner regions.In this section we determine how the different quantities scale with ␦ and obtain the equations of motion that the rescaled quantities satisfy in the different regions.We also show that the solutions in the different regions can be ''matched'' so that they go smoothly into one another.We discuss the main features of the solutions that we obtain in this way and compare them with those found via numerical simulations of the original set of equations ͑1͒ for very small ␦.This comparison shows very good agreement.

A. Scalings
In order to determine how u, v, x, and the velocity of the spots C scale with ␦, we performed a series of numerical simulations for different small values of ␦.Then, by plotting the characteristic length scales x and the maximum and minimum values of u and v in the different regions and the velocity C as functions of ␦ we extracted the scaling relations.͑The scaling can also be determined analytically with a few assumptions about the nature of the solutions.͒We find that v in ϳ␦ Ϫ1 , Cϳ␦, u in ϳx in ϳ␦, u out ϳx out ϳ1, and v out is transcendentally small.The subscripts ''in'' and ''out'' refer, respectively, to the inner and outer regions.
The simulations also show that there are two characteristic time scales.The spots evolve slowly over long-time intervals and then divide and replicate on a fast time scale.
Since the spots travel an O(1) distance at an O(␦) velocity, we conclude that the slow time scale is O(␦ Ϫ1 ).The solutions that we construct are valid only during these long slow intervals, but they do predict when and why replication occurs.As we show later, these slowly varying solutions are such that u has a unique extremum in each outer and inner region ͑maxima in the outer and minima in the inner re-gions͒.

B. Rescaled quantities and equations of motion
Taking into account all these facts, we derive now a simplified set of equations that determine the evolution of a single spot in the interval ͓x 0 M ,x 1 M ͔.Here x 0 M and x 1 M are the locations of two successive maxima of u that are, in general, functions of the slow ͓O(␦ Ϫ1 )͔ time variable.Any point in the spot can identify its position.We will call this point x 1 m .We will show that it differs from the location of the unique minimum of u in this interval by an O(␦) quantity.The point x 1 m is also a function of the slow time variable.The interval ͓x 0 M ,x 1 M ͔ is divided into three regions: the outer regions to the right and left of the spot and the inner region centered at x 1 m .As we show later, the single-spot solution can be used to construct an array of moving spots on a larger ͑possibly infinite͒ domain.
Given the observed scalings, we define in ϭBϵ␦Bt, ͑2͒ in the inner region and out ϭϵ␦t ͑4͒ in the outer ones.We also introduce the following expansions of u and v in powers of ␦ : and assume v out ϭ0. ͑8͒ We introduce the expansions ͑5͒-͑8͒ and the rescaled coordinates ͑2͒-͑4͒ in the evolution equations ͑1͒.Equating terms with equal powers in ␦ yields the following hierarchy of equations.In the outer regions, for O(␦ 0 ), and for O(␦ n ), nу1, In the inner regions, for O(␦ Ϫ1 ), where we have defined the linear operator

͑15͒
the rescaled velocity c( in )ϵͱB(‫ץ‬x 1 m ‫ץ/‬ in )ϭC/␦ͱB, and where ͓u in v in 2 ͔ (n) is the sum of the terms of O(1) in the expansion of u in v in 2 /B 3/2 ␦ n .The boundary conditions, required to complete the definition of L, are specified in the following subsection and Appendix B.

C. Matching and boundary conditions
Now we need to solve each equation of the hierarchy in the appropriate region and then ''match'' the inner and outer solutions where the regions meet.The slow time variables out and in enter only parametrically.Therefore, each equation is an ordinary differential equation ͑ODE͒ with coefficients ͑e.g., c͒ that are actually functions of out or in .In order to solve the ODE's we need to supplement them with boundary conditions, which will also be functions of out or in .The matching of the solutions in the different regions determines these boundary conditions.For this reason, we will first analyze how the solutions behave in the ''matching'' region and then deduce the boundary conditions that allow the construction of a solution uniformly valid in the domain ͓x 0 M ,x 1 M ͔.For the clarity of the presentation we concentrate in this section on the leading-order equations.All higher-order calculations are discussed in Appendix A.
From both the left and right outer regions, the matching regions correspond to the limit xϷx 1 m .In the inner region the left and right matching regions correspond to the limits x in →Ϯϱ.The condition v out ϭ0 implies that v in →0 as x in →Ϯϱ.Thus, if in Eqs.͑11͒ we neglect the terms that are nonlinear in v in , we find Equations ͑16͒ and ͑17͒ imply that where L Ϯ , M Ϯ , and v Ϯϱ are functions of in .Using Eqs.
͑18͒ and ͑5͒, we find that in the matching regions ͓12͔.
In order to determine the behavior of the outer solutions in the matching regions xϷx 1 m , we Taylor expand them around xϭx 1 m .We define and where the subscripts Ϫ and ϩ correspond, respectively, to the outer regions to the left and right of the spot.Using these definitions and Eq.͑3͒, we then write the outer solutions in the matching regions as from which we obtain

͑24͒
Thus, in order to have u out ϭu in in the matching regions, we find from Eqs. ͑20͒ and ͑24͒ that Thus, this discussion shows that the appropriate boundary conditions for Eqs.͑11͒ are while those for Eqs.͑9͒ and ͑10͒ are ‫ץ‬u outϪ given the fact that we are looking for solutions such that u out has maxima at x 0 M and x 1 M .Equations ͑25͒-͑27͒ determine how the boundary conditions of the inner and outer equations are related.

D. Construction of the single-spot solution
We now show how to obtain the single spot solution up to O(␦) in u and O(␦ Ϫ1 ) in v for a case in which x 0 M and x 1 M do not change with time.This holds, in particular, when x 0 M ϭ0 and x 1 M →ϱ ͑i.e., a spot moving into an infinite me-dium͒.In order to simplify the notation, since we focus on the leading-order terms, from now on we drop the superscripts ͑1͒ from the inner solutions, unless otherwise noted.The construction of the single-spot solution to all orders is discussed in Appendix B.
Equations ͑9͒ and ͑10͒ subject to the boundary conditions ͑29͒ can be solved explicitly to all orders and for all out 's.
To leading order we find where we have used Eq.͑25͒.Once we have the solutions u outϮ (0) ( out ), we can calculate the fluxes of u into the spot to leading order in ␦, L Ϯ ( 1) , by means of Eq. ͑22͒.As the reader will notice, they provide a link between the analytic solution that we discuss in this section and the heuristic explanation of spot replication presented in Sec.II.If we further use Eq.͑26͒, we find Notice that, since x 0 M рx 1 m рx 1 M , then 0рL Ϯ рA 1/2 /B.As one would expect intuitively, the fluxes L Ϯ are monotonic bounded functions of the distance from the spot to the adjacent maxima.
Given the fluxes L Ϯ into the spot at a particular time ˜, the boundary conditions ͑28͒ are completely specified.Therefore, in principle, we can solve Eqs.͑11͒.We do this numerically by a shooting method.Solving Eqs.͑11͒ determines the values of c and M Ϯ ͑see Appendix B͒.
Using Eqs.͑32͒ and the definition of c, c( in ) ϭͱB(‫ץ‬x m ‫ץ/‬ in ), we find

͑33͒
Since c is uniquely determined in terms of L Ϯ by solving Eqs.͑11͒, then Eq. ͑33͒ can be integrated to find L Ϯ at any time in which in turn can be used to obtain x 1 m ( in ) and c( in ).Therefore, by this procedure, we know L Ϯ , c, and all their derivatives at any time.Knowing this and M Ϯ , which we also get by solving Eq. ͑11͒, we can completely specify u outϮ (1) ( ˜).The derivatives L Ϯ (2) of u outϮ (1) can then be used to integrate Eq. ͑12͒ as we explain in Appendix B. In this way, the solution up to order ␦ in u and ␦ Ϫ1 in v is obtained as a function of .

E. The inner equations and their solutions
We now describe the properties of Eqs.͑11͒ and of the solutions that satisfy the boundary conditions ͑28͒.We first note that the feed term of u, A(1Ϫu), is completely absent from Eqs. ͑11͒.This agrees with our heuristic picture of spot formation: a spot's structure is determined by the lateral fluxes of fuel into it, which to this order are L Ϫ (1) and L ϩ (1) .The absence of the feed term also allows us to scale the rate B out of the equations, such that the solutions for different values of B are related by a simple rescaling of the variables.Another feature of the equations is that ‫ץ‬ 2 u in /‫ץ‬x in 2 у0 for all x in due to the fact that u in and v in are never negative.This means that u in can only have one minimum and no other extremum.Thus we cannot describe pulse division of the u field using the inner equations alone.However, we can obtain solutions with different number of extrema in v in .In particular, we are able to predict spot replication since the splitting is always preceded by v in going from having one to having three extrema.
We mentioned before that solving Eqs.͑11͒ and ͑28͒ gives the velocity c as a function of L Ϯ .In general, there is not a unique c for each pair ͑L Ϫ , L ϩ ͒, but a discrete set of values.Of all the possible solutions, we will describe only those that we think are relevant for spot replication, unless otherwise noted.Even if we restrict ourselves to these cases, the surface c(L Ϫ , L ϩ ) has various sheets that are intertwined in a complicated fashion.We describe its structure using cuts of constant L ϩ .Note that due to the symmetry of the equations and of the boundary conditions, if there is a solution with L Ϫ ϭL 1 , L ϩ ϭL 2 and cϭc 0 then there is also a solution with L Ϫ ϭL 2 , L ϩ ϭL 1 and cϭϪc 0 .Thus plots with constant L Ϫ can be obtained from those with constant L ϩ by simply changing the sign of c.We show in Fig. 4͑a͒ a cut of constant L ϩ ϭ1 and in Fig. 4͑b͒ another one with L ϩ ϭ1.7901 ͑these values were the ones used for the simulations shown in Figs. 3 and 2, respectively͒.In both cases we observe two curves, but while in Fig. 4͑a͒ they intersect at L Ϫ ϭL ϩ , cϭ0, in Fig. 4͑b͒ they do not intersect at all.Moreover, the value of c for the solid line curve in Fig. 4 is never zero.We show how this intersection disappears in Fig. 5, where we have again plotted c vs L Ϫ for various values of L ϩ between 1 and 1.7901.A smooth change occurs as L ϩ is increased between L ϩ ϭ1 and L ϩ ϭ1.334: the dashed curve goes from crossing the line cϭ0 only once ͑at L Ϫ ϭL ϩ ͒ to crossing it three times ͑preserving the crossing at L Ϫ ϭL ϩ ͒.We do not show this change in Fig. 5, but it is not hard to imagine how it occurs.We do show how both curves behave at L ϩ ϭ1.334 and 1.3рL Ϫ р1.38 in Fig. 5͑a͒.There we can see that for some range of values of L Ϫ there are four different solutions: three on the dashed curve and one on the solid one.The two curves still cross at L Ϫ ϭL ϩ , cϭ0, which means that there are two stationary solutions ͓13͔.These solutions are both symmetric, one single bumped and the other double bumped.Figure 5͑b͒ corresponds to L ϩ ϭ1.3385.There we see that there are two intervals of L Ϫ with four coexisting solutions.In one of these intervals three solutions lie on the dashed curve and one on the solid one, while in the other interval the situation is reversed.
The two curves approach one another as L ϩ is increased and so do the intervals with four solutions.At L ϩ Ϸ1.3387 ͓Fig.5͑c͔͒ the intervals merge, the two curves come together, and the solutions on both of them become equal at the new point of intersection.As L ϩ is increased further, two ''new'' solid and dashed curves can be distinguished, which do not intersect.They separate as shown in Fig. 5͑d͒, which corresponds to L ϩ ϭ1.34.Note that half of each of the new curves shown in this figure ''comes from'' the solid curve of Fig. 5͑a͒ while the other half comes from the dashed one.As L ϩ is increased further, the ''loop'' of the upper curve shrinks and finally disappears.This occurs at L ϩ Ϸ1.347, which is shown in Fig. 5͑e͒.Also at this point the two symmetric solutions with L Ϫ ϭL ϩ and cϭ0 become equal and disappear as L ϩ is further increased.The absence of symmetric solutions with cϭ0 may be observed in Fig. 5͑f͒; which corresponds to L ϩ ϭ1.38.We also observe that the singularities in the derivative of the upper curve disappear as L ϩ is increased.

F. Putting the pieces together
We now match the inner solutions to the outer ones to determine the dynamics of the patterns.As before, we restrict ourselves to cases in which x 0 M and x 1 M are time independent.For simplicity, we further assume that x 1 M →ϱ, the extension to finite x 1 M being straightforward.Under these assumptions, Eq. ͑32͒ implies that L ϩ is also time independent.Thus we can readily use Figs. 4 and 5 to study the evolution of the whole solution.Equations ͑32͒ also determine that 0 ϽL Ϫ рL ϩ ϭA 1/2 /B at all times.Therefore, the inner solutions with L Ϫ ϾL ϩ are not relevant for this case.
Let us first consider a case with A 1/2 /Bϭ1 ͓Fig.4͑a͔͒ and inner solution lying on the solid curve.Assuming that initially x 1 m is finite, then L Ϫ (0)ϽL ϩ and c(0)Ͼ0.Thus, ac- cording to Eq. ͑33͒, L Ϫ will increase steadily, L Ϫ →L ϩ , and c will decrease, c→0, as →ϱ.Therefore, the system approaches a symmetric stationary solution asymptotically in time.A similar situation holds if the inner solution initially lies on the dashed curve of Fig. 4͑a͒, but in this case the system would go to a double-bumped solution, if it were stable ͑which it turns out not to be͒.
Let us now consider a case with A 1/2 /Bϭ1.7901͓Fig.4͑b͔͒.If initially the inner solution lies on the solid curve and x 1 m is finite, then L Ϫ (0)ϽL ϩ and c(0)Ͼ0 as before.However, in this case, while L Ϫ also increases, approaching L ϩ asymptotically in time, c goes to a nonzero value.This means that the spot, were it stable, would continue to move forever at an almost constant speed.On the other hand, if the inner solution lies initially on the dashed curve and x 1 m is finite, then v(x) approaches an asymmetric stationary solution for which cϭ0 and L Ϫ ϭL*Ϸ1.256.If initially L Ϫ ϽL*, then L Ϫ will increase and c will decrease with increasing time, while the situation will be reversed if initially L Ϫ ϾL*.
How do these behaviors compare with the numerical simulations of Eqs.͑1͒?We can compare the L ϩ ϭ1 case with the simulation of Fig. 3 and the L ϩ ϭ1.7901 case with the simulation of Fig. 2 before the first splitting.͑The comparison is done on half the spatial range spanned by the simulation͒.In both cases we observe that after the initial transient dies out, the system approaches a situation in which the inner solution lies on the solid curves of Figs.4͑a͒ and  4͑b͒.Indeed, the analytical solutions are in very good agreement with the numerical ones, as we will show.Now, there is a major difference between the cases L ϩ ϭ1 and L ϩ ϭ1.7901.While in the first case the analytical solution describes the evolution observed in the simulation for all times, in the second case, the agreement lasts only between splittings.As mentioned before, the dynamics during splittings cannot be described within the assumptions of our analytical solutions since the splitting transition occurs on a fast time scale.
Given these results, we might then ask the following two questions.First, why does the numerical simulation always ''choose'' an inner solution that lies on the upper curves?Second, why does spot splitting occur for L ϩ ϭ1.7901?In addition, related to this, does the analytical solution provide any hint on when a splitting is about to occur?The answer to these question depends on the stability of the solutions that we discuss in the next subsection.

G. Stability: How spot splitting arises
We interpret the splitting of the spots as an instability of the spot solutions just described that occurs on the fast time scale t.To demonstrate this, we calculate the stability of the solutions of Eq. ͑1͒.We introduce a perturbation to the base solution: u(x,)→u(x,)ϩ(x,)exp͓(/B)t͔, v(x,)→v(x,)ϩ(x,)exp͓(/B)t͔, where ϳO(␦) and ϳO(␦ Ϫ1 ) in the inner regions and ϳO(1) and ϭ0 in the outer ones.Note that the perturbations evolve on the fast time scale, whereas the base solutions u and v vary only on the slow time scale .Linearizing Eq. ͑1͒, transforming to the inner coordinates, and keeping only terms to leading or- In the outer region, the linearized equations allow only the null solution to leading order in ␦.Matching this to the inner region perturbations yields the boundary conditions ‫ץ‬ ‫ץ‬x in →0 as x in →Ϯϱ, ͑35͒ ͑x in ͒→0 as x in →Ϯϱ.
For a given value of L ϩ and L Ϫ , we determine the stability of the solution by calculating u in and v in numerically and then computing the eigenvalues of a spatial discretization of the operators ͑34͒ and ͑35͒.A positive implies instability.We know that there is always a zero eigenvalue due to the translational symmetry of the equations.This provides a useful check on our numerics since the zero eigenmode is proportional to the derivatives ‫ץ‬u in /‫ץ‬x in and ‫ץ‬v in /‫ץ‬x in .
In particular, we have looked at the stability of the solution branches plotted in Figs. 4 and 5.We show some of the results in Fig. 6, where we have plotted the maximum eigenvalue ␣ of the operators ͑34͒ and ͑35͒ as a function of L Ϫ for various values of L ϩ , with solid ͑dashed͒ curves corresponding to solid ͑dashed͒ curves of Figs. 4 and 5. Figure 6͑a͒ corresponds to the solutions of Fig. 4͑a͒, for which L ϩ ϭ1.There we can see that the solid curve solution is stable while the other one is unstable.This explains why only this solution is observed in numerical simulations of the set ͑1͒.This situation ͑i.e., a stable solid curve and an unstable dashed one͒ continues to hold if we further increase L ϩ until we reach the point at which there are three different solutions on the solid curve for the same values of L Ϫ and L ϩ ͓see, e.g., Figs.5͑b͒ and 5͑c͔͒.At this point, some of the solutions on the solid curve are unstable, as may be observed in Fig. 6͑b͒.We also see that all the solutions on the lower curve of Fig. 5͑c͒ are unstable.When the solid and dashed curves of the c vs L Ϫ plot are rearranged ͑see Fig. 5 and its explana-tion͒, the stability also changes.The solutions on the new solid curves are stable for low enough values of L Ϫ , but become unstable for large values of L Ϫ .This may be observed in Fig. 6͑c͒, where we have plotted the values of ␣ for the solutions of Fig. 5͑e͒.There we see that although the solutions on the dashed curve of Fig. 5͑e͒ remain unstable for L Ϫ ͓1.3, 1.38͔, the values of ␣ decrease continuously from ␣ϭ0.15 at L Ϫ ϭ1.3 to ␣Ϸ0.01 at L Ϫ ϭ1.38.On the other hand, there are stable solutions on the solid curve of Fig. 5͑e͒ for all L Ϫ Ͻ1.3476, but all the solutions on this curve are unstable for L Ϫ Ͼ1.3476.This explains why at L ϩ ϭ1.7901, when there is a one to one correspondence between solutions on the solid curve and values of L Ϫ ͓see Fig. 5͑f͔͒, the solutions on this curve go from being stable at small values of L Ϫ to being unstable at large values of L Ϫ ͓see Fig. 6͑e͔͒.As we have already mentioned, the case L ϩ ϭ1.7901 corresponds to the numerical simulation of Fig. 2 before the first splitting.Thus we conclude that the splitting occurs because the solution becomes unstable.
The stability analysis in this section is somewhat different from the standard linear analyses that are the foundation of bifurcation theory.Here the analysis is of the model that arises when the evolution occurs on a slow time scale.We obtain the base solution assuming that time scales as O(␦ Ϫ1 ) and we thus approximate it by a solution of the first set of equations in the hierarchy ͑9͒-͑14͒.We then look for instabilities that grow on an O(1) time scale, and for this reason we modify the equations accordingly, as described before.An instability in this case implies that it is impossible to split the original set of equations ͑1͒ into the hierarchy ͑9͒-͑14͒.It also implies that our approximate base solution is no longer ''close'' to the actual solution of the full set ͑1͒.This breakdown of the model precludes the type of analyses that are usually done in bifurcation theory ͑indices, etc.͒.

H. Comparison with the numerical simulations
We have compared the analytical solutions with those obtained from numerical simulations of the original equations ͑1͒.We have found very good agreement both in the ''form'' of the spots u in (1) (x in ) and v in (1) (v in ) and in their velocity c( in ).We have also observed that spot splitting occurs after the onset of the instabilities.This means that the stability calculation allows the prediction of the occurrence of splittings.All these properties are an indication that the analysis is essentially correct.
We illustrate some of these results with Figs. 7 and 8. Figure 7 is a plot of the maximum growth rate of the instability ␣ as a function of the location of the spot in rescaled coordinates x m for the solution on the solid curve of Fig. 4͑b͒, which has L ϩ ϭ1.7901.This analytical solution is to be compared against simulations with Aϭ0.02 and Bϭ0.079. Figure 8 is a plot of the spot location in rescaled coordinates x m vs in for this analytical solution; a simulation with A ϭ0.02, Bϭ0.079, and ␦ 2 ϭ0.01; and a simulation with A ϭ0.02, Bϭ0.079, and ␦ 2 ϭ0.0001.We can immediately observe the good agreement among all curves before the occurrence of splittings.This means that the analytical solution gives a good approximation of the numerical one between splittings.On the other hand, if we compare Figs.7 and 8, we observe that the onset of the instability of the analytical solution occurs before the splitting and that the time delay between both events decreases with ␦.This behavior highlights the predictive value of the stability calculation.The stability calculation is done under the assumption that the instabilities grow on a fast time scale.It is clear that this assumption is violated immediately after the onset, where the growth rate is as small as we wish.Therefore, there is a time interval on which our simplified model with two separate time scales does not hold.Since the slow time scale is of O(␦ Ϫ1 ), we expect the simplified model to be valid after the growth rate gets bigger than ␦.As is clear from Fig. 8, the rescaled time at which this occurs decreases with ␦.This means that the transient with no separate time scales gets smaller with ␦.Thus we expect the stability analysis to be more accurate as ␦ decreases.In fact, this may be observed in Fig. 8, where it is clear that the onset of the instability gives a better estimate for the time of splitting as ␦ gets smaller.

I. Construction of the N-spot solution
So far we have considered a single spot moving into an infinite medium.However, our solutions describe correctly the evolution observed between splittings in simulations with any number of spots.Suppose that we want to describe the evolution of N spots in the interval ͓0,ϱ͔ that are located at x i m ( in ), 1рiрN.We are interested in the case in which u has a maximum at xϭx 0 M ϭ0 and another one at x FIG. 7. Plot of the maximum growth rate of the instability ␣ as a function of the location of the spot in rescaled coordinates x m for the solution on the solid curve of Fig. 4͑b͒, which has L ϩ ϭ1.7901.We observe that the onset of the instability occurs at x m Ϸ7.8.

ϭx N
M →ϱ for all in .The coordinates of all the other maxima of u,x i M , 0ϽiϽN, are functions of in .Assuming the ordering x iϪ1 M Ͻx i m Ͻx i M , 1рiрN, the various x i M 's and x i m 's are related by Thus, using an extension of Eqs.͑32͒ to the case of N spots and Eq.͑36͒ we can write the rescaled fluxes L Ϯi ϵ ϮL Ϯi (1) /B into the ith spot as where, as before, c i ϭͱB(‫ץ‬x i m ‫ץ/‬ in ) and c 0 ϭ0.The N-spot solution can be constructed in almost the same way as we did with the single spot one.Given the distances between spots or the fluxes of u into each of them, we can solve Eqs.͑11͒ for each spot separately.This determines the values of c i and of the other parameters of the solution as a function of the fluxes.This information can then be used to integrate Eqs.͑38͒.Here is the main difference with respect to the single spot case: now the time evolution of the various spots is coupled by Eqs.͑38͒.Once the fluxes and the velocities are known as functions of in , we can go to higher orders as in the single-spot case.

IV. OTHER MODELS
An analysis similar to the one presented in the preceding section can be carried out for a whole class of models of which Eqs. ͑1͒ are a particular example.The evolution equations for the models in this class can be written as where ␣у0, ␤у0, ␥Ͼ0, ⑀Ͼ1, and ␣(⑀Ϫ1)р␤␥.Also in this more general case we can divide the space into inner and outer regions of width O(␦) and O(1), respectively.In the outer regions u is of O(1) and v is transcendentally small, while in the inner ones, uϳ␦ 2(pϪ1)/͓␤␥Ϫ(␣Ϫ1)(⑀Ϫ1)͔ and v ϳ␦ 2␥/͓␤␥Ϫ(␣Ϫ1)(⑀Ϫ1)͔ .We can obtain a hierarchy of equations equivalent to Eqs. ͑9͒-͑14͒ if we introduce the same rescaled coordinates as before and power series expansions for u and v consistent with the new scalings.In this way, the equations to leading order in ␦ are Eq.͑9͒ in the outer and a simulation with Aϭ0.02, B ϭ0.079, and ␦ 2 ϭ0.0001, and the analytic solution that corresponds to these values of A and B ͑i.e., L ϩ ϭ1.7901͒.A comparison with Fig. 7 shows that the splitting occurs after the onset of the instability.
in the inner one where the velocity c defined as before.
The outer equations can be solved to all orders.However, we cannot ensure the existence of a solution to Eqs. ͑40͒ in the general case.If such a solution exists, then the matching and the construction of the whole solution can be carried out as before.On the other hand, we cannot ensure that this solution would be ''attracting'' in the same way it is for Eqs.͑1͒.However, in all the cases we tried we observed spot replication in the simulation of the corresponding Eqs.͑39͒.We show an example in Fig. 9, where we have plotted v(x) at different times for a model with ␣ϭ0, ␤ϭ1, ␥ϭ1, ⑀ϭ2, Aϭ0.3, Bϭ0.08, and ␦ 2 ϭ0.005.

V. CONCLUSIONS
We have presented the calculations contained in ͓6͔ in detail and shown how they generalize to cover a class of systems that can be described by the fuel and fire picture.
The key element for spot replication is the multivaluedness of c(L Ϫ ,L ϩ ) and the disappearance of the cϭ0 branch of solutions when the fuel flux exceeds a critical value.
Ours is not the only theory of spot replication.The existence of nontrivial stationary solutions that approach the fixed value uϭ1, vϭ0 as ͉x͉→ϱ in the Gray-Scott model has also been analyzed by Doelman, Kaper, and Zegeling ͓7͔.They prove the existence of single-and multiple-pulse solutions in the infinite line.They also prove that traveling pulses of the same type cannot exist.The authors conclude that this nonexistence is somehow related to splitting: ''while the numerically observed moving pulses begin to resemble the non-existing ͓sic͔ traveling solitary pulses more and more, they must undergo some transformation, such as pulse-splitting ͓sic͔.''This is similar to what we have found in our stability analysis, since in our case a fast developing instability means the nonexistence of a solution with the time scale separation we had assumed.
Hagberg and Meron ͓8͔ also see spot replication in a modified Fitzhugh-Nagumo model.This model has also been analyzed in this context by Muratov and Osipov ͓14͔.The analysis in ͓8͔ is centered around the nonequilibrium Ising-Bloch transition: a bifurcation from a static front to a pair of counterpropagating Bloch waves.This transition is in essence a pitchfork bifurcation in which the cϭ0 solution loses stability to solutions cϭϮͱ, where is the distance from the bifurcation.In their work spot replication occurs when the velocity of a propagating curved front changes sign as a function of curvature along the front.Thus the variation of curvature along the front is a key element in the analysis of Hagberg and Meron.However, one of the experimental observations is of a growing disk of high pH.Eventually the center of the disk collapses, leaving a propagating annulus.This effect has not been observed in the modified Fitzhugh-Nagumo model, which forms the basis of the work by Hagberg and co-workers.It seems unlikely that it will be observed since no replication has been observed in the onedimensional modified Fitzhugh-Nagumo equations.The fact that the annuli have curved fronts should not make much difference since the curvature is constant along the front.
In a key respect our analyses are similar in that c is a multivalued function of either extrinsic or intrinsic parameters ͑such as flux and curvature͒ and the stability of the c ϭ0 solution changes at a bifurcation point.However, in our case the c surface has many different sheets that are intertwined in a complicated fashion as can be seen by the cuts shown in Fig. 5, whereas in the case considered by Hagberg and Meron the static cϭ0 solution undergoes a pitchfork bifurcation.As far as qualitatively matching the experiments, we do not do a particularly good job either since our spots correspond to spots of high concentration of the autocatalytic species and the experimentally observed spots are spots of high pH.Since H ϩ is presumably the autocatalytic species our picture does not compare well.
Lee and Swinney ͓4͔ have performed numerical simulations of the four-variable Gaspar-Showalter ͓5͔ model in one space dimension and find spot replication also.There are qualitative differences between Lee and Swinney's simulations and their experiments.In his experiments, colliding fronts repel and in their simulations they annihilate.Another difficulty in reconciling the experiments and simulations is that in the simulations replication was found only when the diffusion coefficient of H ϩ was set lower than the other species.The absence of a mechanism that slows down the diffusion of H ϩ relative to the other species makes it doubtful that the simulations provide an accurate mirror of the experiments.We also remark that Lee and Swinney's simulations are not well described by either our work or the work of Hagberg and Meron.The null clines of the two-variable Gaspar-Showalter model bear little in common with either the Gray-Scott null clines or the Fitzhugh-Nagumo null clines.
Currently there are no results, either analytical or numerical, in full qualitative agreement with the experimental results.Since there are qualitative differences between the analyses, the experiments, and Lee and Swinney's numerical simulations we conclude that the replication phenomenon is more general than demonstrated by the current analyses.and ЈϭwЈϩx 1 M Ϫx 1 m .We now focus on the solutions to the inner equations, for which we need to specify the boundary conditions.According to the discussion in Appendix A, the appropriate boundary conditions are where nу1 and the various w Ϯ (n) must be calculated as in Eq. ͑A10͒ with u outϮ (k,n) , respectively.It follows from Eqs. ͑B1͒-͑B3͒ and ͑A10͒ that the derivatives ‫ץ‬w Ϯ (n) /‫ץ‬x in , nу1, do not depend on M n .In fact, they can be completely specified in terms of the parameters A and B, the functions ͕M Ϯ (l) , lϽn͖, the distances x 1 m Ϫx 0 M and x 1 m Ϫx 1 M ͓or the fluxes L Ϯ through Eq. ͑32͔͒, and their derivatives.Thus the boundary conditions for u in (n) can be written completely in terms of the same quantities.In particular, for nϭ1, they are given by Eq. ͑28͒.Then, since solving the inner equations for (u in (n) ,v in (n) ) determines the values of M Ϯ (n) ͑see the discussion below͒, we need to solve the inner equations at each order for all times before we can move to the next order.
Although we cannot solve the inner equations analytically, we can still describe how to carry on the construction of the whole solution up to any desired order as follows.Let us choose, at one particular time ˜, a value for the distances x 1 m Ϫx 0 M and x 1 m Ϫx 1 M .Given these distances, we can immediately calculate the outer solutions u outϮ (0) , using Eqs.͑30͒ and ͑31͒, and the fluxes L Ϯ , using Eq.͑32͒, at that same time.Given the fluxes, we can integrate numerically Eqs.͑11͒ in order to find a particular solution (u*,v*) that satisfies the boundary conditions ͑28͒ at that time.Let us write this solution in a way such that in the matching regions x in →Ϯϱ.Suppose, without loss of generality, that the origin of coordinates is such that ‫ץ‬u in */‫ץ‬x in ϭ0 at x in ϭ0.Then, any pair of functions of the form is also a solution of Eqs.͑11͒ with asymptotic behavior where Thus any translation of the original solution of the form ͑B6͒ is a solution of Eqs.͑11͒ that satisfies the boundary conditions ͑28͒ and can be matched to the left and right outer solutions up to this order.So a general solution of Eqs.͑11͒ is a function of of the form ͑B6͒-͑B8͒, where is the point at which ‫ץ‬u in (1) /‫ץ‬x in ϭ0.It is clear that while M ˜Ϯ and v ˜Ϯϱ change with , c and L Ϯ remain invariant.Therefore, as we mentioned before, we can integrate Eq. ͑33͒ without knowing the exact value of .We proceed by noting that a numerical integration of Eq. ͑12͒ requires that we specify seven constants: the boundary parameters L Ϯ , M ˜Ϯ , v ˜Ϯϱ and the velocity c.In addition, due to the translational invariance of Eq. ͑12͒, there is a degenerate family of solutions differing only in their shift from the origin .The numerical procedure consists of shooting the solutions into the origin from their asymptotic behaviors at x in ϳϮϱ ͓Eq.͑B7͔͒ ͑using a suitably large value of x in in place of x in ϳϮϱ͒.At the origin, there will be four matching conditions: the continuity of u ˜,v ˜and their first derivatives.We introduce a fifth matching condition to break the translational degeneracy by requiring ‫ץ‬u ˜/‫ץ‬x in ϭ0 at x in ϭ0 ͑this has the effect of setting ϭ0 and M ˜ϮϭM Ϯ * ͒.If we now fix the values of L Ϯ ͑by specifying the locations of the spots͒ we see that we have five free parameters to vary such that the five boundary conditions are satisfied, which implies that there will be a discrete spectrum of acceptable values for M Ϯ * , v ˜Ϯϱ , and c.Numerically we accomplish this by using an Adams-Bashforth solver to integrate the equations, while varying the parameters using a Powell method ͑see ͓15͔ for details on these methods and the subroutines used͒.
This integration gives c and L Ϯ ͓and L Ϯ (1) via Eq.͑26͔͒ as functions of in .Also, we can formally write M Ϯ and v Ϯϱ as in Eq. ͑B8͒ for all times with an unknown function of in .Thus, using Eqs.͑26͒, ͑27͒, ͑32͒, ͑B1͒, and ͑B3͒, we can write the solutions u outϮ (1) and its derivatives L Ϯ (2) , at all times, in terms of L Ϯ (1) ( in ), c( in ), and the unknown function ( in ).In particular, the L Ϯ (2) read Thus, using Eqs.͑B4͒, ͑A1͒, and ͑B9͒, we can write the boundary conditions for Eq.͑12͒: v in ͑ 2 ͒ →0 as x in →Ϯϱ at all times in terms of the same quantities ͓15͔.Now, Eq. ͑12͒ is linear and we know that Lwϭ0 for The existence of this zero mode of L is due to the translational invariance of Eqs.͑11͒.Due to this zero mode, we can also see that if there is a solution (u**,v**) T of Eq. ͑12͒, then it is not unique since any function of the form ͓u** ϩ␣(‫ץ‬u in (0) /‫ץ‬x in ), v**ϩ␣(‫ץ‬v in (0) /‫ץ‬x in )͔ T with ␣ an arbitrary real function of in is also a solution.Now, the boundary conditions are not ''translationally invariant'' in the following sense.In the hypothetical case that we could solve the equations at a given time ˜for an interval of values of L Ϯ ( 2) , then the different solutions would not be related by a spatial translation.Therefore, we expect Eq. ͑12͒ to be solvable only for certain values of L Ϯ ( 2) or, equivalently, certain values of ( in ).
This discussion shows that we first need to ensure that Eq. ͑12͒ is indeed solvable.In order to do this, we transform our problem so as to have homogeneous boundary conditions.For this purpose we define

͑B12͒
where f Ϯ (x)→1 for x→Ϯϱ.For example, we may choose f Ϫ (x)ϭ(1Ϫtanh x)/2 and f ϩ (x)ϭ(1ϩtanh x)/2.Therefore, we now need to solve the equation subject to the boundary conditions ͓16͔ ͩ ‫ץ‬p ‫ץ‬x in q ͪ ϭ0 for x in →Ϯϱ.͑B14͒ The zero mode w satisfies these new boundary conditions, namely, ‫ץ‬w 1 /‫ץ‬x in ϭ0ϭw 2 as x i n→Ϯϱ.Then, the operator adjoint to L, L* also has a zero mode w* such that ‫ץ‬w 1 */‫ץ‬x in ϭ0ϭw 2 * as x i n→Ϯϱ.Then, for Eq.͑B13͒ to be solvable we need its right-hand side to be orthogonal to w* ͑see, e.g., Ref. ͓17͔ and Krischer and Mikhailov ͓7͔͒.The orthogonality condition will establish a relationship between L Ϫ (2) and L ϩ (2) at each value of in , from which, in principle, we should be able to determine the value of ( in ).
Assume now that we have determined ( in ) and found a solution to Eq. ͑B13͒ or, equivalently, to Eq. ͑12͒.Then, as we mentioned before, we can add an arbitrary multiple of the zero mode of L and still have a solution.This arbitrary multiple can be determined at the next order by exactly the same technique that we have just described, namely, by requiring solvability of Eq. ͑13͒, which is linear and involves the same operator L. Clearly, at each order we will have the same arbitrariness, in the sense that we can add any multiple of the zero mode of L and still have a solution.A solvability condition at the following order will then determine this multiple.