Front propagation in spatially modulated media

J. Armero, 1 A. M. Lacasta, 2 L. Ramı́rez-Piscina, 2 J. Casademunt, 1 J. M. Sancho, 1,3 and F. Sague ́s Departament d’Estructura i Constituents de la Mate `ria, Universitat de Barcelona, Avenida Diagonal 647, E-08028 Barcelona, Spa Departament de Fı ́sica Aplicada, Universitat Polite `cnica de Catalunya, Avenida Dr. Gregorio Maran ̃ón 50, E-08028 Barcelona, Spain Institute for Nonlinear Science, Department 0407, University of California, San Diego, 9500 Gilman Drive, La Jolla, California 92093-0407 Departament de Quı ́mica Fı́sica, Universitat de Barcelona, Avenida Diagonal 647, E-08028 Barcelona, Spain ~Received 31 July 1996 !


I. INTRODUCTION
The study of propagating fronts has been a problem of great interest in a rich variety of very different situations ͓1-3͔.The steady-state velocity and shape of the front are problems not yet satisfactorily solved in many situations.Although different kinds of fronts can be defined in many contexts, it is commonly accepted that they can be classified and then studied in a general framework.In this paper we will consider stable fronts which propagate with a fixed velocity and flat shape if the medium is isotropic and homogeneous.This type of situation can be modeled by a reaction-diffusion equation for the order parameter.The phenomenology of this situation is well known ͓4-10͔.Here we will consider the case of a front moving in a nonhomogeneous medium.Within this general framework we will assume that some parameter controlling the local front velocity presents a transversal spatial modulation.This has been the case for several experiments of chemical waves propagating into modulated excitable media ͓11-13͔, in which stationary patterns with a fixed velocity were obtained by maintaining spatial modulations of the chemical excitability of the medium.
Our present study shows that the selected stationary pattern and the corresponding propagation velocity result from a nontrivial global competition process between different maxima of the local velocity, which are coupled through local curvature effects.An example of such competition is shown in Fig. 1 ͑top͒, where there exists a modulation ͓Fig. 1 ͑bottom͔͒ in the x direction of the local front velocity u(x).Starting from a planar front, the system initially mimics the modulating function by developing as many front maxima, which hereinafter we will call fingers, as local maxima presented by that function.The evolution of each one of these fingers, moving with different velocities, turns out to depend on the local details of u(x) around the maxima.That gives rise to a slow competition process where some of the competing fingers are eliminated and some others survive.Finally, the resulting stationary pattern may be quite different from the initial one.The main objective of this paper is to characterize and explain these facts within an analytical framework based on singular perturbation techniques, and provide simple analytical criteria to predict the velocity and the qualitative shape of the stationary front shape for different kinds of spatial modulations.An experimental text of this selection problem was presented, for a photosensitive version of the Belousov-Zhabotinsky reaction, in Ref. ͓13͔.
The paper is organized as follows.In Sec.II, the interface dynamics in generic reaction-diffusion systems is briefly reviewed, and the extension of the theory to smooth modulated media is performed.Section III is devoted to the main problem, namely the description of the competition process.The final stationary state is characterized in Sec.IV, where the number of surviving fingers and the stationary velocity of the whole front are determined for the case of smooth modulations.A particular case in which the spatial variation is not smooth, and consists of a homogeneous stripe of a larger velocity surrounded by a homogeneous area with a smaller one, is analyzed in detail in Sec.V.In Sec.VI, we summarize our main results and make some final remarks.Numerical results are presented along the paper showing an excellent agreement with the analytical predictions.

II. LOCAL EQUATION OF MOTION
Our starting point is a field model of a two-dimensional system with scalar order parameter (r ជ ,t) governed by the generic reaction-diffusion equation ͓2-9͔ ‫ץ‬ ‫ץ‬t ϭٌ 2 ϩF͑͒.͑1͒ The reaction kinetic term F() is assumed to be a nonlinear continuous function that allows the existence of two homogeneous stationary states or phases, 1 and 2 , i.e., F( 1)ϭF( 2 )ϭ0.The interface between these two phases is supposed to be thin as compared to its typical scale of curvature, and it is identified with our local front.We are interested here in those situations where the propagating front describes the invasion of the 2 state, either metastable or unstable, by the globally stable 1 state.
As it is well known, Eq. ͑1͒ has a planar front solution propagating at a well-defined velocity u 0 ͓8͔.When the 2 state is metastable, the velocity u 0 is uniquely determined ͓14͔.If 2 is unstable, the steady-state version of Eq. ͑1͒ does not uniquely determine the propagation velocity.In this case, the asymptotic front speed is selected dynamically and, for sufficiently localized initial conditions, it corresponds to the velocity of the front propagating with the steepest decay to 2 .When the linear-marginal-stability criterion holds ͑linear regime͒, the front velocity approaches asymptotically the value u 0 ϭ2ͱFЈ(0) ͓4,5͔.However, for some parameter regime near the metastable region ͑nonlinear regime͒, the linear-marginal-stability criterion fails ͓8͔.
When the front is not planar, curvature effects correct the actual front velocity, tending to restore the planar shape if the medium is homogeneous, and for a modulated medium providing the spatial coupling that determines the transient dynamics and the final steady state.By projecting Eq. ͑1͒ on the interface and using a quasistatic approximation ͓15,16͔, the front dynamics can be well approximated by the local equation where v n is the local normal front velocity and denotes the local curvature of the front.The normal velocity is taken as positive if the globally stable 1 state is invading the 2 state, and the curvature is taken negative at the tip of a finger of phase 1 .
Others types of models could have been considered instead of the simple model represented by Eq. ͑1͒, with two coupled equations, such as in excitable media ͓17͔ or in solidification systems ͓18͔, presenting a richer phenomenology.
In order to study the influence of a spatial modulation of the medium on the front dynamics we consider an explicit spatial dependence in the reaction term of Eq. ͑1͒.This dependence is introduced through a modulation of an external control parameter ␣.By considering an orthogonal fixed frame (x,y) with the y axis along the direction of propagation of the front, we define our model equation as This particular way to introduce the modulation has two restricting features.First, the modulation appears as a multiplicative factor in the reaction term.This simplifies the model by preserving the homogeneous stationary states, and modulates only the strength of the driving force acting on the front.Second, and most importantly, the modulation has only a spatial dependence in the transverse direction to the front propagation.This situation is interesting by itself, as shown by existing experiments designed under these conditions ͓11-13͔.
When such a modulation of the reaction term is present, the planar front is no longer a solution, and Eq.͑2͒ should be modified.If the spatial modulation is sufficiently smooth, within the quasistatic approximation, and following the usual projection methods ͓15͔, Eq. ͑2͒ can be generalized to The explicit relation of u(x) with the external modulation ␣(x) given by u(x)ϭu 0 ͱ␣(x) can be derived in the follow- ing way.In the reference frame of the moving onedimensional stationary front, zϭyϪut, Eq. ͑3͒, reads

ЉϩuЈϩ␣F͑͒ϭ0, ͑5͒
where a prime denotes differentiation with respect to z.By rescaling z in Eq. ͑5͒ by zϭ/ͱ␣, we obtain This is the equation of a planar front arranged to the spatially homogeneous reaction term F() propagating at a velocity given by u/ͱ␣.This velocity is nothing but u 0 , so uϭu 0 ͱ␣.
Since this relation is verified at each point of the front, we have

͑7͒
For numerical integration, it is convenient to write Eq. ͑4͒ as an equation for the front position yϭy(x,t), In Fig. 2, we compare numerical integrations of the starting field model Eq.͑3͒ and of the effective equation of motion Eq. ͑8͒ for the same temporal evolution of an initially planar front.In our simulations we have used a standard finitedifference Euler algorithm with ⌬xϭ0.5 and ⌬tϭ0.1.The reaction kinetic term considered is F()ϭ 2 Ϫ 3 , and the external modulation ␣(x) is given by the two-maxima function shown in the inset of Fig. 2. Different choices of the function F() may slightly affect the accuracy of the projected equation ͑4͒.We will comment on this in Sec.VI.In Fig. 2, fronts are shown at different times in a frame moving at the propagation speed of the fastest finger.We see that in both models the slow finger on the left is eliminated, and the front reaches a stationary shape with just one finger ͑local maximum͒, that is, a qualitatively different pattern from that of the external modulation ␣(x).This simple case illustrates how the effective local equation for the front Eq. ͑8͒ captures the competition mechanisms present in the original field model, with the advantage of a much simpler analytic and numeric treatment ͑see Sec.III͒. Figure 1 presents also the same phenomena but with a more complicated modulation.

III. ANALYTICAL DESCRIPTION OF THE COMPETITION PROCESS
The main question we address in this section is the analytical description, in the context of a singular perturbation scheme, of the competition process ͓19͔.To this aim, we write Eq. ͑4͒ in terms of the angle variable (x,t), defined by tanϭ‫ץ‬y/‫ץ‬x ͑see Fig. 2͒.The following geometrical relations hold: This can be substituted into Eq.͑4͒ to obtain Given that y(x,t)ϭ͐ x dxЈtan, we differentiate Eq. ͑9͒ with respect to x to obtain

͑10͒
In view of Eq. ͑10͒ a simple picture of the competition process between fingers arises.The value of the expression in parentheses in that equation is the local velocity in the y direction.After a short transient, this quantity reaches a roughly constant value for each finger, but different from that of other fingers, in such a way that the spatial derivative in Eq. ͑10͒ is very small, and the shape of each finger does not change significantly.Only the contact points between fingers, which move with different velocities, will present values for the spatial derivative very different from zero, and hence will do so for the time derivative.Therefore, the competition dynamics is basically governed by these contact regions.
Our approach here will be to build a perturbative scheme valid for sufficiently smooth modulations.That will be equivalent to perturbing on the curvature term in Eq. ͑10͒, which, being the highest order derivative, defines a singular perturbation.It is expected that the curvature in most of the front will scale with the typical length scale of the modulation, while the contact regions between fingers referred to above will behave as boundary layers in the perturbative scheme.The matching order by order of the corresponding inner and outer expansions will define the actual solution of the problem.
Without loss of generality we take the modulation of the system as periodic with period L, and consider 1/L as our perturbative parameter.From now on the modulation will be described by the periodic fixed function u(x/L), in such a way that L becomes a parameter that controls the smoothness of the modulation.We start by obtaining the equation for the outer solution by rescaling variables as

͑12͒
As stated before, starting from a planar front, each local maximum of the function u(z) forms a local maximum or finger of the front.The stationary shape of each competing finger i moving at a velocity v i is then given by Substituting the following expansions in the above Eq.͑13͒: we obtain at the lowest order, FIG. 2. Same temporal sequence of the propagation of an initially planar front under the two-maxima modulation ␣(x)ϭ1Ϫ0.125sin͓(2/L)x͔Ϫ0.375 cos͓(4/L)x͔ ͑see the inset͒.Solid lines correspond to numerical integration of the field model Eq.͑3͒ with F()ϭ 2 Ϫ 3 and dashed lines to numerical integration of the effective equation of motion Eq. ͑8͒.Fronts are plotted every ⌬tϭ100 in the frame moving at the propagation velocity of the fastest finger.
and, for the velocities, where u m i and u m i Љ are the value and the second derivative of the modulating function u(z) at its maximum in the finger i.
In view of Eq. ͑19͒ the velocity is corrected by the length scale given by the spatial variations of the modulating function u near its local maximum.This result can be directly translated to the original field model Eq.͑3͒.In particular, given the relation u(x)ϭu 0 ͱ␣(x), we obtain, up to the low- est orders in the inverse system size, where ␣ m i is the value of the spatial modulation ␣(x) at the maximum of the ith finger, and ␣ m i Љ its second derivative at the same maximum.
The nature of the competition process among fingers is clearly illustrated in Fig. 3, where the local velocity of the front in the y direction, vϭ‫ץ‬y/‫ץ‬t, is plotted for the same spatial modulation and the same planar initial condition as in Fig. 2. Figure 3 shows that the slower finger is actually being invaded by the faster one, as if a one-dimensional front were propagating laterally, in the x direction, with a certain velocity c.This effective front is well defined within our perturbative scheme, and its dynamics can be obtained from the inner solution of the equations in the boundary layers placed at the contact regions between fingers.We write then Eq. ͑10͒ in the original variables x, t and, up to the lowest order in 1/L, consider u as a constant.We obtain Inner solutions coming from this equation, valid inside the boundary layer, have to be matched with the outer solutions found before, which are valid in the regions outside the boundary layer.That gives, as boundary conditions for a layer placed between two fingers moving at velocities v Ϫ and v ϩ , the following We next look for a stationary solution of Eq. ͑21͒ of the form (x,t)ϭ(xϪct)ϭ(), which satisfies where a prime denotes differentiation with respect to ϭxϪct.An integration of this equation gives where K is a constant.Its value is obtained by imposing the boundary condition Eq. ͑22͒, with the result So we obtain, for the velocity of the lateral invading front, This invading front that describes the competition between the two fingers does not propagate uniformly since, according to Eqs. ͑16͒ and ͑17͒, Ϯ ride on the local details of the modulated velocity u(x).This invasion process is clearly seen in recent experiments ͓13͔.Making use of Eq. ͑23͒, the spatial dependence of the invasion velocity c may be explicitly written as The propagation velocities of the two fingers, v Ϯ , are given by Eqs.͑18͒ and ͑19͒ when evaluated at the two local maxima of u(x).
The analytical prediction for the invasion velocity c, Eq. ͑25͒, is not valid in the initial short transient when fingers are still not formed.Once the front forms its initial fingers, as dictated by the form of u(x), the competition process that arises is well described by Eq. ͑25͒ as long as v Ϯ 2 Ϫu 2 (x) is non-negative.
In Fig. 4 we compare the analytical predictions of Eq. ͑25͒ for the case F()ϭ 2 Ϫ 3 , to numerical integrations of the effective front motion Eq. ͑8͒ for systems of two different sizes L but the same two-maxima spatial modulation ␣(z) and the same planar initial front.

IV. STATIONARY STATE
A result from the perturbative analysis above is that, for large enough L, only one finger survives, even if there exist some other local maxima of the modulating function.From the lowest-order approximation, Eq. ͑4͒ gives negative curvature only for the absolute maximum u M of the modulating function.In that case, only the fastest finger survives, and its velocity is adopted by the whole front.This velocity is in the original, fully dimensional, variables The shape of the front is again given by Eqs.͑16͒ and ͑17͒, but applied to the values u M and u M Љ of the absolute maxi- mum of u(x).
In terms of the external modulation ␣(x), this velocity reads This perturbative result gives the more accurate values for the whole front velocity the larger is L or, equivalently, for a fixed L, the smoother is the external modulation.In Fig. 5 we show the stationary front velocity for F()ϭ 2 Ϫ 3 and for ␣(x) given by a Gaussian centered in a Lϭ500 system.
The Gaussian modulation has a fixed maximum ␣ M ϭ1.5, so the curvature at the tip, ͉ M ͉ϭ͉␣ M Љ ͉ provides a measure of the length scale of the spatial variation of the modulation, (͉␣ M Љ ͉/␣ M ) Ϫ1/2 .The broken line in Fig. 5 corresponds to the analytical result Eq. ͑27͒ and, for sufficiently smooth modulation, agrees accurately with results from numerical integration of both the local ͑effective͒ and the two field models.
In general, the single-finger stationary state will only occur if the velocity v selected by the system and given, up to the first order, by Eq. ͑26͒, is greater than the value of the other local maxima of the u function.This condition is al-ways satisfied for sufficiently large L. Conversely, any secondary local maximum with u greater than v will present a negative curvature ͓see Eq. ͑4͔͒, and thus will form a secondary finger.In this way, the comparison between the value of v and that of the different local maxima of u(x) does provide a criterion to identify the surviving fingers in the final stationary state, such as those for which u m i Ͼv.On the other hand, to obtain from the perturbative analysis a good prediction for the selected velocity v one has to take the largest value taken by Eq. ͑14͒ at the different local maxima u m i of u(x).Notice that the selected velocity calculated this way and, therefore, the final stationary front shape, depend not only on the values of u(x) at their different local maxima but also depend on the second derivative of the modulation at those maxima.
With this picture in mind, a more detailed discussion of Fig. 1 is now in order.This figure shows the case of an eight-maxima modulating function.The steady state presents only five maxima or fingers, instead.Finger numbers 2, 5, and 7 have been eliminated during the transient.In the upper part of the figure we see how the initial planar front develops the eight maxima of the modulation ͑bottom of the same figure͒.The competition process then sets in on a much slower time scale.In the bottom of Fig. 1 we also see how maximum 4 has the largest predicted velocity according to Eq. ͑26͒ ͑denoted by crosses͒, which determines the final velocity ͑dashed line͒.Maxima numbers 1, 3, 6, and 8 also have a larger value of u(x) at their tip than the final selected velocity, and therefore they survive all the way to the steady state.This example also shows how the final shape of the front can be very different from that of the modulation.

V. STRIPED CASE
In this section we consider the particular case in which u(x) takes a constant value u M in a central stripe of width W along the y direction, and a smaller constant value u m outside.This is clearly an exception to the perturbative treatment of Sec.IV because the modulating function is not smooth.We address this case here because it has been previously studied both theoretically and experimentally in Ref.
͓11͔, which directly motivated the present study.In that reference the theoretical predictions were obtained numerically.Here we will provide some analytical results with the hope of gaining some insight into a situation where the perturbative approach is not appropriate, particularly from the dependence of the steady state on the parameters of the problem.
According to Eq. ͑8͒, the steady-state shape of the front must satisfy where we have used the fact that vϭ‫ץ‬y/‫ץ‬t is the constant velocity of the whole front in the y direction.This equation can be solved piecewise for constant u(x), so an explicit solution for this case can be found by appropriate matching at the boundaries between regions.
The fact that the perturbative prediction fails here is obvious from the fact that all orders except the lowest one vanish, since all derivatives of u(x) are zero.However, notice that the zeroth-order solution does describe the correct large-L limit, which corresponds to three straight pieces, one horizontal in the central stripe, and two inclined ones on the sides of the stripe, with angles given by Eq. ͑16͒, and a selected velocity equal to u M .
In order to obtain explicit analytic expressions for the selected velocity we will find a rigorous lower bound using the explicit solutions for constant u(x).This bound will indeed turn out to be a very good estimate of the actual selected velocity.
To solve Eq. ͑28͒ by reduction of order, we seek a solution of the form yЈϭsinh. Its substitution into Eq.͑28͒ gives a first-order equation for , which is easy to solve by direct integration.Without loss of generality let us assume that the central stripe is placed at the center of the system of size L with periodic boundary conditions.Then the symmetry of the problem enables us to solve Eq. ͑29͒ for the central and lateral regions separately.The solution y c for the central band 0ϽxϽW/2 with yЈ(0)ϭ0 is then

͑30͒
with a c ϭv/u M Ͻ1.The solution y l for the lateral region W/2ϽxϽL/2 with yЈ(L/2)ϭ0 is given by

͑31͒
with a c ϭv/u m Ͼ1.For the lateral region there exists also the trivial family of planar solutions given by cos p ϭu/v.
The solution at the central band, Eq. ͑30͒, is periodic.The most relevant feature is that it intersects itself, giving rise to periodic loops, separated by mostly planar regions.These loops are unphysical, so the matching points with the external solutions Eq. ͑31͒ must be such that there are no loops between them.This physical requirement provides a lower bound to the selected velocity, by imposing that the distance between the turning points with infinite slope be greater or equal than the size of the central band W. The bound for a c is then given by the equation The lower bound for the steady front velocity obtained from the previous equation turns out to be very close to the actual velocity selected from the complete matching of the solutions, since the solution near the turning point has a large curvature, so the point with the right slope should be very close to it.Given that Notice that the dependence on W displayed by Eq. ͑33͒ is quite different from that the system size dependence obtained for the smooth case in the framework of the singular perturbative scheme.The excellent accuracy of these approximated results can be seen in Fig. 6.An exhaustive theoretical study of a model system with two bands of different excitability ͑local velocity͒ was presented very recently in Ref. ͓20͔.Their results are in agreement with the ones obtained in this section.

VI. CONCLUDING REMARKS
The results presented in this paper refer basically to the local model derived from standard projection techniques from a reaction-diffusion field model, where the nonlinear function defining the reaction term has not been specified.In the simulations reported here we used F()ϭ 2 Ϫ 3 , for which we have quantitatively checked the validity of the local approximation.The question of the degree of indepen- dence of the result on the nonlinear function F() naturally arises.It is well known that the presence of a linear term such as in F()ϭϪ 3 changes the nature of the velocity selection problem ͓8͔.For F()ϭϪ 3 , linear marginal stability theory applies, and the transient decay to the steady state is of power-law nature.For the case F()ϭ 2 Ϫ 3 , instead, nonlinear theory applies, and the transient regime decays exponentially.This difference is important since the quasistatic approximation involved in the projection of the problem into a local equation is not as well justified in the linear case due to its long transients.Although simulation of the linear model may give apparently stationary shapes which are qualitatively similar to the ones obtained with the nonlinear model, in the former case it is always much more difficult to conclude about stationary states because of contributions to the velocity vanishing as t Ϫ1 may produce ar-bitrary large deviations in the front position.In Fig. 7 we show the convergence to the steady value in the two cases of F().It is also shown the effect of the spatial discretization, which may produce errors comparable to the transient effects.
If a more complicated model than Eq.͑3͒ is chosen, but presenting a front structure obeying Eq. ͑2͒, then our conclusions apply in the same way.Our theoretical results are in agreement with recent experiments in front propagation in excitable media ͓11,13͔.
In conclusion, we have studied fronts propagating through two-dimensional media modulated in the transverse direction.An effective local equation for the motion of the front has been derived and it has been used to explain, in the context of a singular perturbation scheme, the dynamics of the competition process leading to the nontrivial stationary solution.Explicit criteria to determine the qualitative shape of the stationary front and quantitative estimations of the steady front velocity have been found for the generic case of smooth modulations and for a particular case of nonsmooth modulation with relevance to experiments.A generic picture of competition between fingers has naturally arisen in terms of lateral fronts propagating in the transverse direction describing the invasion of slower fingers by their neighboring faster ones.Finally, the reliability of the effective local equation for different reaction terms of the original reactiondiffusion model has been discussed.

FIG. 1 .
FIG. 1. Temporal evolution of an initially planar front.Fronts are shown in the frame comoving with the fastest finger ͑4͒.At early times, the front y(x,t) mimics the modulated local velocity u(x) with eight maxima ͑bottom͒.Three slow fingers are eliminated before the front reaches the stationary shape with only five fingers.See more details in the text.The front is plotted at times 0,10, 20, 200, 2500, 15 000, and 30 000 ͑solid line͒.Also see the discussion at the end of Sec.IV.

FIG. 3 .
FIG. 3. Competition process between two fingers.The vertical component of the local velocity of the front is plotted every ⌬tϭ50for the same evolution shown in Fig.2.The effective transversal velocity c ͓see Eq. ͑25͒ in the text͔ is also indicated.

FIG. 4 .
FIG. 4. Scaled velocity c of the transversal invading front for two systems of different sizes L ͑cf. figure legend͒ under the same external modulation shown in Fig. 2. Thick solid lines correspond to the analytical prediction, Eqs.͑25͒, for each size L.

FIG. 6 .
FIG. 6. Stationary velocity for fronts propagating in striped media with u M ϭ0.866 and u m ϭ0.707 as a function of the central stripe size W. Dashed line is the analytical prediction Eq. ͑33͒.

FIG. 7 .
FIG. 7. Transient effects of the nonlinearities in the numerical evolution of the front velocity.Full symbols correspond to the F()ϭϪ 3 model, and hollow symbols to the F()ϭ 2 Ϫ 3 model.Triangles correspond to ⌬xϭ1, and circles to ⌬xϭ0.5.The cross is the analytical value of the velocity.