Selection , shape , and relaxation of fronts : A numerical study of the effects of inertia

We study the problem of front propagation in the presence of inertia. We extend the analytical approach for the overdamped problem to this case, and present numerical results to support our theoretical predictions. Specifically, we conclude that the velocity and shape selection problem can still be described in terms of the metastable, nonlinear, and linear overdamped regimes. We study the characteristic relaxation dynamics of these three regimes, and the existence of degenerate ~‘‘qu nched’’! solutions.


I. INTRODUCTION
Front propagation is a subject of intense study due to its relevance in a variety of contexts ͓1͔.In almost all cases, this problem was modeled by a parabolic partial differential equation.This research was devoted to the case of a globally stable state that invades an unstable or metastable state, focusing on the velocity and profile selection problems ͓2-11͔.A relevant result is that the phenomenology can be classified into three regimes: metastable, nonlinear, and linear ͑see Sec.II below for further details͒.
An important question regarding this problem is the relevance of inertia ͑i.e., when can a parabolic equation be used͒ and what are its main effects on the overdamped scenario.To this end, we will consider a damped, hyperbolic partial differential equation.We emphasize that models of this type arise when one considers a more realistic jump process for individuals whose probability density is described by the partial differential equation.Moreover, hyperbolic equations of this kind describe many actual physical phenomena, such as, e.g., population dynamics, nonlinear transmission lines, cell motion, branching random walks, dynamics of ferroelectric domains, and others ͓12-18͔.Therefore, it is physically crucial in all those contexts to know whether or not one can use the simpler, parabolic equation.
To our knowledge, the above issue has not been systematically addressed in the past.Important results were obtained by Hadeler ͓15͔, who found a mathematical proof of the relationship of both equations, although without an explicit discussion of the three main regimes and the shape of the selected front.More recently, in Ref. ͓17͔ a damped hyperbolic equation was studied, but only in the linear regime.Finally, there was a brief mention of the hyperbolic problem in Ref. ͓19͔ in connection with the algebraic relaxation of the front velocity, although the corresponding asymptotic analytical prediction was not checked numerically.
In this paper, we address the question of how the selection problem is affected, concerning both the velocity and shape of the front edges, the existence of degenerate solutions, and the relaxation to a steady state in the different regimes.We study the problem from a numerical point of view, trying to obtain accurate results in every regime.In addition, we show how theoretical predictions can be obtained in a simple way within a scheme similar to that in Ref. ͓8͔.The results so obtained are relevant in a more general context, with the same phenomenology but without exact solutions, providing a way to identify the different regimes from numerical experiments.The paper is organized as follows: In Sec.II, we introduce the model and briefly review the results for the overdamped case.Section III is devoted to a theoretical and numerical study of the inertial dynamics of fronts.Conclusions are summarized in Sec.IV.

II. MODEL AND BACKGROUND
The front dynamics subject to inertia is generically modeled by the dimensionless hyperbolic equation ͓15,17,19͔ ⑀ tt ϩ t ϭ xx ϩ f ͑ ,a͒.͑1͒ The parabolic or overdamped limit is obtained by letting ⑀ →0 ͓15,17͔, which leads to ͑note that it is a singular pertur-bation͒ which, being well known, will be taken as the reference scenario.As a representative example of a nonlinear reaction term we take The parameter a acts as an external control parameter.The allowed interval for a, (Ϫ1/2,1͔, ensures the global stability of the 2 state.We then have the following regimes. ͑i͒ For a(Ϫ1/2,0), 1 is a metastable state.It can be proven in this case that the kinklike solution which fulfills the boundary conditions, is unique, and that it invades the 1 state with velocity v nl and spatial asymptotic decay k nl : This is called the metastable regime.͑ii͒ For a͓0,1͔, 1 is an unstable state.In this case, the situation is somewhat more rich.
͑a͒ a͓1/2,1͔.The linear-marginal-stability criterion applies ͓8͔: Initial conditions with compact support evolve toward a solution with minimum velocity, whose spatial decaying profile is given by This is the linear regime, in which the front velocity relaxes to a stationary value following an algebraic law, instead of the exponential decay found in other instances ͑see Ref. ͓19͔ and references therein͒.͑b͒ a͓0,1/2͔.The linear-marginal-stability criterion now fails, and a full nonlinear solution of Eq. ͑2͒ is needed, so we are in the nonlinear regime with solutions ͑5͒ and ͑6͒.It can be shown that initial conditions with kϾk*ϭͱ2a propagate, after a short transient, with the nonlinear velocity v nl of the unique solution of the metastable regime ͑i͒ and k nl .
In the last two regimes higher propagation velocities can exist by choosing suitable initial conditions with an asymptotic spatial decay given by k i Ͻk l or k i Ͻk*, respectively.In fact, there is a continuous degeneracy of solutions for steadily propagating fronts as well as, correspondingly, a continuum of possible velocities.All those solutions behave as e Ϫkx as x→ϱ ͑the front travels to the right in view of the boundary conditions͒, with a k-dependent velocity.Such fronts maintain their shape as far as the asymptotic decay is concerned, and propagate with a velocity vϭ(k i 2 ϩa 2 )/k i .These degenerate solutions are usually referred to as quenched fronts.

III. RESULTS
We begin with a study of the effects of inertia corresponding to Eq. ͑1͒.We expect that the system will exhibit frontlike solutions (xϪvt).In accordance to this, let us assume that Eq. ͑1͒ has a generic frontlike solution of the form ͓17͔ ͑x,t͒ϭ͑BxϪAt͒ϭ͑ y ͒. ͑8͒ Substitution of this expression into Eq.͑1͒ leads to an equation for (y), given by yy ϩA y ϩ f ͑ ,a͒ϭ0, ͑9͒ which is the same equation one would obtain beginning with Eq. ͑2͒.The parameters A and B have to be taken as where v(a) is the front velocity of the parabolic case ͓Eq.͑2͔͒.
It is important to note that, at this level of theoretical analysis, both models ͓Eqs.͑1͒ and ͑2͔͒ are esentially the same.What makes them different is the presence of the parameter B and, very likely, the dynamics, because it is assumed that the selection problem is a dynamical one and then the necessary information cannot be obtained from Eq. ͑9͒ alone.
From Eqs. ͑8͒-͑10͒, we find that the velocity is

͑11͒
in agreement with the rigorous result of Ref. ͓15͔.
Let us now focus on the shape of the frontlike solution, an aspect which has not been studied before for the hyperbolic case, to our knowledge, but which is immediately understood within the present analysis.The front behaves as e Ϫkx for x→ϱ: As k(a) is the spatial decay of the leading front, described by the solution of Eq. ͑9͒, the corresponding k(⑀,a) has to be k͑⑀,a ͒ϭBk͑ a ͒ϭͱ1ϩ⑀v 2 ͑ a ͒k͑ a ͒.͑12͒ Of course, for ⑀ϭ0 we recover the parabolic or overdamped values ͓Eqs.͑6͒ and ͑7͔͒, whereas in the opposite limit ⑀ →ϱ and v(⑀,a)ϳ⑀ Ϫ1/2 .In principle, these two equations ͑11͒ and ͑12͒ give us the information we were looking for, provided v(a) and k(a) are known, which implies that the parabolic problem has to be solved for the three regimes.A remarkable conclusion of our theoretical analysis is that the boundary between linear and nonlinear regimes (aϭ1/2) coincides for both models, independently of the inertia parameter ⑀.However, these being nonrigorously proved analytical results, we have to check whether these solutions are observable through numerical simulations of the full partial differential equation ͓Eq.͑1͔͒.
We begin the discussion of the numerical results with discussing Fig. 1, which contains data regarding the verification of the theoretical prediction for v(⑀,a) ͓Eq.͑11͔͒.In all cases, the initial profile (x,0) is a steplike function, which is allowed to evolve numerically under a fourth-order Runge-Kutta scheme ͓20͔.The plotted velocities are measured after transients have died out.As can be seen, the agreement between theory and simulations is excellent, even for values of ⑀ϭ10, which can hardly be regarded as a small perturbation.As for the front shape, we observed that during the numerical evolution, the initial step transformed into a tanh-like front, similar to expression ͑5͒.Again, there is a nice agreement between theory and the simulations, a little worse for very large values of ⑀ and aϷ1 because of the larger time needed for transients to die out.
We have to note that in both the linear and the nonlinear regimes there is the possibility of two front solutions with different values of k but the same velocity.This can be obtained from the linear approximation of Eq. ͑1͒ by substituting the leading contribution of the front, ϳe Ϫk(xϪvt) , xӷvt, ͑13͒ into Eq.͑1͒.In this way we find ⑀k 2 v 2 ϩkvϭk 2 ϩa, ͑14͒ an equation which has two solutions in k for every v value.
In the same way, we can analyze the spatial tail of the back-end front, which is characterized by a spatial decay k s .In the metastable and nonlinear regimes, Eq. ͑5͒ is a solution, and then k s ϭk nl (⑀,a).In the linear regime, we can proceed as in Eq. ͑14͒ by performing a linear analysis around the stable state ϭ1.Substituting the relaxing or back-end contribution of the front, ϳ1Ϫe k s (xϪvt) , xӶvt, ͑15͒ into Eq.͑1͒, we find another expresion for v(k s ), different from Eq. ͑14͒: It is worth mentioning here that the nonlinear solution ͑6͒ is a solution of both Eqs.͑14͒ and ͑16͒.Expression ͑16͒ gives us either v(k s ) or k s (v) as a function of the model parameters: In Table I we present the results of numerical simulations which confirm these predictions with only minor round-off errors.We have also verified that the same analysis works for the model of Ref. ͓21͔.However, we note that this is not a general approach as, for instance, the model with f () ϭϩd 3 Ϫ 5 proposed in Ref. ͓8͔ does not fulfill the condition k s ϭk nl ͓8͔.
As for the degenerate or quenched solutions, we have found that, much in the same way as in the parabolic model, we also have the possibility of generating metastable solutions which move faster than the selected one in the hyperbolic problem.This phenomenon appears for aϾ0.The idea is to prepare an initial front profile with a small value of k i in Eq. ͑5͒.In the nonlinear regime the condition is kϽk*, where k* is evaluated in Eq. ͑14͒ using v nl (⑀,a).For the linear regime, the condition is k i Ͻk l (⑀,a).We have found numerically that the velocities of these solutions agree with Eq. ͑14͒ up to an accuracy better than 0.1%, and that their spatial decay remains the same as that of the initial condition.For these particular solutions, it can be analytically shown that the back-end spatial relaxation k s is different from k i .Furthermore, for ⑀Ͼ1/2, the possibility arises that the relaxation of the back end to the stable state is oscillatory instead of simply exponential.Both predictions have been verified satisfactorily by numerical simulations.
To complete our study of the selection problem, we have also considered the way the velocity approaches its asymptotic value.In the overdamped case, the steplike profile accelerates and changes its shape until the stationary shape and velocity are reached.In the underdamped problem, we have found that, in general, the scenario is the same, although the approach to the asymptotic velocity is nonmonotonic, i.e., the velocity increases initially to values larger than the asymptotic one, and then relaxes to it in a damped, oscillatory manner, the number and amplitude of the oscillations depending on the model parameters.Leaving aside these initial oscillations, we have checked that v(t) relaxes exponentially in the metastable regime.This behavior, predicted by theory ͓8͔, can be seen in the inset of Fig. 2, which shows that it is very clear in the parabolic model, whereas in the hyperbolic model it is masked by the above mentioned oscillations ͑see the inset͒.However, the relaxation time is of the same order as that of the parabolic problem, and therefore the behavior is roughly the same.Com-FIG. 1.Comparison of analytical ͑lines͒ and numerical results ͑symbols͒ The solid line corresponds to the overdamped case (⑀ ϭ0), whereas dashed lines represent the analytical results ͓Eq.͑11͒ for the velocities ͑upper panel͒ and Eq.͑12͒ for k ͑lower panel͔͒ for the inertial case.M, NL, and L indicate the metastable, linear, and nonlinear regimes, respectively.Symbols: ⑀ϭ0 ͑squares͒, ⑀ϭ1 ͑circles͒, and ⑀ϭ10 ͑triangles͒.versely, in the linear regime, we have found a power law relaxation with an exponent ϳ1, for the parabolic model, and some evidence of smaller exponents ͑decreasing with increasing ⑀) in the hyperbolic model.Examples of this are shown in Fig. 2 ͑upper curves͒.We note that the ⑀ 0 curves bend slightly at long times; this could imply a deviation from the power law behavior, but we cannot exclude that it might come from the manipulation of the numerical data.Hence, we are more confident of the initial part of the curve, and thus we conclude that the decay is of power law form in the linear regime.Finally, for the nonlinear regime ͑the lower curves in Fig. 2͒, both models give a slow relaxation ͑not exponential͒, whose functional form is not really clear, although it might be fitted with some accuracy by a stretched exponential.In any event, the data indicate that the decay is neither power law nor exponential in this case.
Regarding this problem of the approach to the asymptotic velocity, it is interesting to check the theoretical prediction for the case aϭ1 in Ref. ͓19͔, namely, where v*ϭv l (⑀,1), kϭk l (⑀,1), and Dϭ(1ϩ4⑀) Ϫ2 .In view of the fact that Ref. ͓19͔ did not contain any numerical check of this claim, we have verified this prediction in simulations.
Our result is shown in Fig. 3 for ⑀ϭ1, and the agreement is very good.In fact, as can be seen from the plot, if one substitutes v(⑀,1) with the numerically obtained asymptotic value ͑as discussed in Ref. ͓19͔͒, the theoretical prediction overlaps the numerical results.We have verified that this agreement holds for several values of ⑀ and, in particular, is much better for ⑀ϭ0, for which there is no need to modify the theoretical value v* in Eq. ͑18͒.In addition, in this overdamped limit, it is clearly seen the need for the t Ϫ3/2 term in Eq. ͑18͒ for the prediction to be accurate.This is less dra-matic when ⑀ 0 because D is very small, and then the asymptotic is much slower and practically irrelevant, as shown in Fig. 3.

IV. CONCLUSIONS
To summarize this work about inertial effects on front dynamics, we have numerically checked that fronts governed by the hyperbolic equation ͑1͒ behave in a fashion qualitatively very similar to the parabolic case, the inertial term yielding mostly quantitative corrections.In particular, the separation into metastable, linear, and nonlinear regimes holds for the underdamped case, even for very large values of the inertial parameter.This includes the selected velocity for the propagating front, the spatial decay constants ͑both in the back-end and the forward parts͒, the possibility of quenched fronts, and the approach to the asymptotic velocity ͑with special emphasis on the aϭ1 case͒.
Moreover, it is possible to understand this separation with the same approach and the ideas of the linear-and nonlinearmarginal-stability criteria, the boundaries between them being the same as in the overdamped case.This is an important result for two reasons: First, it shows that, for actual physical systems where the overdamped approximation may or may not be accurate, the picture of the phenomenology that one should expect is not changed qualitatively by neglecting inertial terms.On the other hand, this result provides a starting point to study the effect of any other perturbation that may affect the front dynamics, allowing one to stick to the well known framework worked out for the overdamped front dynamics to interpret the effects of possible additional perturbing terms.
Finally, it is worth commenting on the applicability of the results reported here to related problems.The theoretical predictions and their excellent agreement with the numerical results have permitted us to establish very firmly the exis-FIG.2. Approach to the asympotic velocity in the nonlinear ͑lower curves, aϭ0.25) and linear regimes ͑upper curves, a ϭ0.75).Solid lines are the overdamped cases, and dashed lines correspond to ⑀ϭ1.v asymp is the asymptotic velocity reached in the simulation.The inset is a plot of the metaestable regime (aϭ Ϫ0.25), exhibiting its exponential decay.FIG. 3. Approach to the asympotic velocity for aϭ1 ͑linear case͒ and ⑀ϭ1.Circles: numerical solution.Dot-dashed line: Eq. ͑18͒ up to first order.Dotted line: Eq. ͑18͒.Full line: Full Eq. ͑18͒ with the theoretical velocity replaced by the asymptotic velocity obtained from the numerical simulation.Note that the difference between full and dotted lines is less than 0.2%.tence, range of validity, and main features of the three regimes for the underdamped case.We believe that this is a general result.In this case, the characterization of the different temporal decays of the front velocities can then be useful to identify the border between the metastable, nonlinear, and linear regimes in other, non tractable models directly from simulations.

TABLE I .
Back-end (k s ) and forward ͑k͒ spatial decay coefficients obtained from numerical simulations and compared to theoretical predictions (k s th , k th ).Numerical values for k are computed from exponential fits to the tail.