Exact Propagating Wave Solutions in Reaction Cross-Diffusion System

Reaction-diffusion systems with cross-diffusion terms in addition to, or instead of, the usual self-diffusion demonstrate interesting features which motivate their further study. The present work is aimed at designing a toy reaction-cross-diffusion model with exact solutions in the form of propagating fronts. We propose a minimal model of this kind which involves two species linked by cross-diffusion, one of which governed by a linear equation and the other having a polynomial kinetic term. We classify the resulting exact propagating front solutions. Some of them have some features of the Fisher-KPP fronts and some features of the ZFK-Nagumo fronts.

1. Introduction. Reaction-diffusion systems are models that are used widely to model physical, chemical, biological, ecological processes. Realistic models of such processes are typically quite complicated and can only be dealt with numerically. However qualitative understanding of the most important features benefits from analytical approaches, even if that requires sacrfices in quantitative accuracy. This may be achieved by using asymptotic methods and/or considering "toy models".
One of the first and famous "toy" reaction diffusion systems is Fisher's [4] model of propagation of an advantageous gene, which is a special case of the equation considered by Kolmogorov, Petrovsky and Piskunov [7]. We refer to it as Fisher-KPP model. Another early archetypal reaction-diffusion equation was a model of flame propagation considered by Zeldovich and Frank-Kamenetsky [17], which later became known also as Schlögl model [10] and Nagumo equation [8]. We refer to it as ZFK-Nagumo equation. Both models have monotonic propagating wavefront solutions of similar appearances, but each has its own distinct mechanism. The Fisher-KPP model shows the transition from an unstable resting state to a stable resting state, while the ZFK-Nagumo model shows the transition from one stable resting state to another stable resting state. Another qualitative difference between them is that ZFK-Nagumo model exhibits a unique, up to a constant shift in time or space, propagating front solution with a fixed speed and shape, whereas Fisher-KPP model has a family of solutions with a continuous range of possible speeds. The importance of these toy models goes well beyond providing simplest examples. For instance, the ZFK-Nagumo equation can be considered as the fast subsystem in describing pulse waves in the FitzHugh-Nagumo and similar systems using singular perturbation techniques [15].
In the last three decays, several interesting phenomena have been discovered in reactiondiffusion systems which have cross-diffusion of the reaction species in addition or instead of their self-diffusion [14,2,13,12]. For example, propagating waves in reaction-cross-diffusion systems (RXD) with excitable reaction kinetics could penetrate each other on collision, a behaviour that is unusual for excitable systems with self-diffusion.
The properties of solutions in RXD systems have been mostly studied numerically. An analytical approach has been attempted e.g. in [2]. In that work, fast-slow separation between reaction kinetics of two reacting species is assumed. The fast subsystem has piece-wise linear kinetics and linear cross-diffusion, and admits exact analytical solutions in the form of propagating fronts. Unlike the Fisher-KPP and ZFK-Nagumo fronts, these front solutions are oscillatory. They can be matched asymptotically with slow pieces to obtain complete asymptotic description of propagating pulses. The fast subsystem in this approach is different from the Fisher-KPP and ZFK-Nagumo equations in two aspects: that it is two component and it is piece-wise linear, as opposed to the "classical" toy models which are both one-component and with polynomial nonlinearity of the kinetics. At least two components are of course required to have cross-diffusion. However, the necessity of piecewise defined kinetics is less obvious.
In the present work, we aim to clarify the situation, and investigate the possibility of having exact front solutions in a cross-diffusion system with polynomial kinetics. The idea is to postulate the solutions and deduce the governing equations from there. For simplicity and as the first step, we only consider here monotonic fronts, similar to those found in the ZFK-Nagumo equation.
The paper is organized as follows. The history and the problem formulation are in section 2. In section 3, we consider the possibilities of chosing polynomial nonlinearity for the reaction term. In section 4, we discuss the simpest aspects of stability of possible solutions. Then we show the correspondent polynomial function suitable for solutions in of the wavefront type. These are presented in section 5. We proved the possibility to use the wavefront solution of the system as generalisation for Fisher-KPP also we analyse the right choices of the parameters to imitate Fisher-KPP model. These analyses are shown in section 6 and section 7. We return to the question of stability, now for the selected wavefront solution, in section 8. Results of numerical simulation are presented in section 9. These simulations show that the wavefront are unstable. The type of this instability is investigated in section 10 and the paper is concluded by discussion in section 11.
The system (2.1) is a reaction self-diffusion system when D uu > 0 and D uv = D vu = 0. Otherwise, D uv = 0 and/or D vu = 0, we have reaction-cross-diffusion system. If = 0, v ≡ 0, the self-diffusion system degenerates to the ZFK-Nagumo equation [17,10,8] for u(x, t), with an exact propagating front solution. A piece-wise linear N-shaped variant of f (u) also leads to exact propagating front solution [8]. Qualitative properties of this equation, including existence of propagating front solutions, persist for a generic N-shape, and for 0 < 1, these solutions can form a basis of asymptotic description, see for instance the review [15].
A similar asymptotic approach for 0 < 1 was considered for the cross-diffusion case of (2.1) in [2]. To make the problem analytically tractable, the consideration there was restricted to a piecewise linear N-shaped function f (u) and pure cross-diffusion, with self-diffion totally absent, D uu = D vv = 0.
In this paper we consider the same system as was dealt with in in [2], namely and intend to extend the methodology of [15] and [2] for a polynomial function f (u). In absence of self-diffusion terms, we abbreviate D vu = D u and D uv = D v . We start by recapitulation of the approach of [2] to set the scene and introduce notation and terminology. Direct numerical simulations of (2.2) with no-flux boundary conditions and cubic f (u) produces, in particular, solutions in the form of propagating pulses of a fixed shape, as illustrated in Figure 1. For small , the width of the pulse grows as O −1 . This means that in the limit → 0, the wave front and the wave back of the pulse go apart. Both the front and back represent transitions between two distinct equilibrium points, say (u 1 , v 1 ) and (u 2 , v 2 ).  Figure 1. (a) The direct numerical simulation of (2.2) the reaction cross-diffusion system with a no-flux boundary exhibits a propagating pulse with f (u) = u(u − 0.3)(1 − u) and the values of parameters are = 0.001, Du = 5 and Dv = 0.5. (b) At small distance and time, the front of the pulse of reaction cross-diffusion system with a cubic nonlinearity approach to two asymptotic states (u1, v1) and (u2, v2).
In the limit → 0 the system (2.2) turns into The two equilibria (u 1 , v 1 ) and (u 2 , v 2 ), the asymptotic states of the wave front and the wave back, satisfy f (u j ) = v j . Letû(ξ) = u(x, t) andv(ξ) = v(x, t) be a propagating wave solution of (2.3), where ξ = x − ct and c > 0. Substituting this into the system (2.3) yields As the front is an asymptotically approaching distinct steady states, we havê Integrating (2.5) with respect to ξ gives When ξ → ±∞, we obtain from (2.8) that v =v 1 =v 2 and then equation (2.4) turns into We have from equation (2.6) that cv = D uû , hencev = D uû /c. Substituting this into (2.5) yields whereû is a wave solution for the reaction cross-diffusion system (2.3). This differential equation is deduced by applying the wave variable on the reaction crossdiffusion system (2.3). Biktashev and Tsyganov [2] have replaced f (û) by a piecewise linear function. The fronts that are obtained from the piecewise linear function are oscillatory fronts that are similar to the fronts of the wave solution seen in numerical simulations with cubic f (û). We seek to consider a polynomial function for f (û) instead of piecewise linear function, which would still yield explicit analytical solutions for propagating fronts.
We aim that functionf (û) is a polynomial. This can be assured if y(û) is a polynomial. Let us find the possible degree of the polynomials y(û) andf (û). Let P n be the set of polynomials of degree n. If y ∈ P n , then y A y 2 + yy + B =f (û) ∈ P 3n−2 .
If n = 1 thenf (û) is linear, which is not of interest for us, as this cannot produce two distinct solutions for (2.9). If n = 2 thenf (û) is quartic. This quartic polynomial is comparable to cubic, in that it can describe bistability, if it has at least three simple roots. Therefore, y ∈ P 2 ,f ∈ P 4 is the simplest suitable choice.
The travelling wave differential equation for ZFK-Nagumo can be solved analytically by a reduction of order [8]. Incidentally, in that solution y(û) is also quadratic. It is convenient to factorise the quadratic polynomial y(û), for some constants k = 0, g and h. Note that due to (2.6), (2.7) and (3.2), we have {u 1 , u 2 } = {g, h}. From (3.2) and (3.4), we obtain where C is an arbitrary constant. The front wave described by (3.5) is illustrated in Figure 2. Onceû(ξ) is known, we can findv(ξ) using (2.8), as Obviously, the profile of componentv represents not a wave front but a pulse. In accordance with (2.7), we havev (±∞) = v .
4. On the stability of the front solutions. The stability of any front solution we seek shall depend, in particular, on the stability of its asymptotic steady states; the latter is easily done using standard methods. The system (2.3) can be written in the matrix form Suppose w * = [u * , v * ] T is an equilibrium, i.e. F(w * ) = 0. We perturb this point, and in the linear approximation we have where F = ∂F/∂w is the Jacobian matrix. By separation of variables, particular solutions of (4.1) bounded in space can be written as linear combinations of Substituting (4.2) in (4.1), gives and eigenvalue problem where f = ∂f /∂u, and the eigenvalues are see Figure 3. Therefore, if f (u * ) is positive, then Re (λ 1,2 ) ≥ 0 and the steady state (u * , v * ) is unstable, and if f (u * ) is negative then Re (λ 1,2 ) < 0 for all µ = 0, and the state is stable in linear approximation. Of course, even if both asymptotic states are stable, the stability of the whole front solution will still depend on the discrete spectrum; this is outside our scope.

5.
Fixing the polynomial reaction term. In this section we will find the particular form of the polynomial functionf (û), as well as the parameters A and B that satisfy (3.3). To achieve this, we substitute (3.4) into (3.3), which gives We take our quartic polynomialf (û) in the following form: where {u 1 , u 2 } = {g, h} and without loss of generality σ = ±1; a different scaling of f would just result in a change of the spatial and temporal scale of the solutions. By substituting (5.2) into (5.1), we obtain By equating like terms we obtain This imposes five constraints onto a set of 11 parameters k, g, h, σ, D u , D v , u 1 , u 2 , u 3 , u 4 and c; hence we can describe all solutions of this system by assigning six of these parameters as free, and then finding the remaining five parameters as dependent on these six free parameters. We restrict consideration to real values of parameters in both groups, except possibly the roots u 3,4 . Moreover, as parameters g and h fix the positions of the pre-and post-front resting states of the solution (3.5), it convenient to have these two among the free parameters; note also that we have already constrained σ to ±1.
6. Possible types of solutions. As discussed in the Introduction, this study is not motivated by any real-world applications leading to specific examples of reaction cross-diffusion systems. Rather, we are interested in theoretical possibilities achievable within a certain class of models. With this in mind, we want to see if we can make the reaction cross-diffusion system with quartic polynomial to look like generalizations, in one sense or another, of other well-known models, from the much better studied class of systems with self-diffusion. We shall say that we "imitate" those models. The models that we want to imitate are Fisher-KPP and ZFK-Nagumo. Those models exhibit propagating front solutions with asymptotics If we identify the scalar field u here with the namesake first dynamic variable in our system, then this property can be achieved by letting g = 0 and h = 1 in (3.5).
We found in the previous section that the stability of a spatially uniform steady state depends on the sign of the derivative of the quartic polynomial at that state. In terms of stability, to imitate the ZFK-Nagumo wave, we would need a stable pre-front state and a stable post-front state, and consequently an unstable equilibrium in between. To imitate the Fisher-KPP wave front we would need to have an unstable pre-front state and a stable postfront state, with zero or two equilibria in between. In this respect, the possibilities for front waves from the reaction cross-diffusion system with quartic polynomial are constraint by the following proposition. Proof. From (5.3), among the roots off (û) we have {u 1 , u 2 } = {g, h}, and the other two roots, u 3,4 , are the roots of the quadratic in the square brackets, which is equivalent tô , this implies that either g and h are two inner roots while u 3 and u 4 are the two outer roots, or vice versa. If u 3 = u 4 the g and h are the two outer roots out of the three, and if u 3,4 ∈ R, then g and h the only two, therefore automatically the outer, roots.
From Proposition 6.1, we conclude that of the resting states of the front wave solution, only one can be stable but not both. That means, in the considered reaction cross-diffusion system with the quartic polynomial, it is impossible to imitate ZFK-Nagumo wave in terms of the stability of the resting states, but there is a chance to imitate Fisher-KPP wave. We note, however, that for any given set of parameters of the model, the speed of the front solution is in any case uniquely fixed by (7.4), and this feature is characteristic of ZFK-Nagumo fronts rather than Fisher-KPP fronts.
7. Choice of Signs to Imitate Fisher-KPP. We have found that there is a possibility to imitate Fisher-KPP front wave, in terms of the stability of the pre-front and post-front equilibria, by reaction cross-diffusion system with quartic polynomial nonlinearity. In this section, we will turn this possibility into reality, by identifying appropriate parameter choices.
Firstly, let us make sure that solution given by (3.5) satisfies the asymptotic boundary conditions of Fisher-KPP front wave, In section 5 we found that six parameters in (5.4) can be arbitrary assigned. We choose k, g and h as three of such free parameters, in order to satisfy (7.1). We have already committed ourselves to the choice {g, h} = {0, 1}, and we require k = 0. Table 1 lists the resulting four a priori possibilities. Table 1 Examining possible choices to imitate Fisher-KPP front. The symbols ( ) and ( ) mean that χ(ξ) is an increasing or decreasing function, respectively.

Choices
Results Clearly, choices that comply with (7.1) are (I) and (III). In both cases, equation (3.4) gives The quartic polynomialf (û) posited in (5.2) allows σ = 1 or σ = −1. Remember that the equation for the coefficients atû 4 in (5.4) states If σ = 1 then the solution (3.5) will not satisfy the condition (7.1): since D u , D v and c are positive, equation (7.3) implies k < 0, which is inconsistent with (7.2). So, we must choose σ = −1, which together with {g, h} = {u 1 , u 2 } = {0, 1} turns the system (5.4) to Previously, we let variables g, h and k be free parameters. We now add to that list D u and D v . The rest of the variables will be dependent on those as follows: 3 + 36ρ, (7.5) The quartic polynomial now has the form where u 3,4 is given by (7.5) and (7.6). We expect that, in principle, if the quartic polynomial is substituted into the system (2.3) , i.e. Remember that by virtue of (7.8) and (7.4), this means that the location of the roots u 3,4 is determined by the three parameters k, D u and D v .
8. Stability of the Resting States. Previously, we have linearised the system (2.3) for general function f (u) about the equilibrium and derived the formula of the eigenvalues (4.4). Substituting the quartic polynomial function (7.9) into the function of the eigenvalue yields that, the eigenvalues of the equilibrium u 1 = 0 are given by while the eigenvalues of the equilibrium u 2 = 1 are given by In the "inner roots" case I, the two roots u 3 and u 4 have different signs, and are to opposite sides of 1. Thus, from (8.1) and (8.2) we deduce that the pre-front u 1 = 0 is stable and the post-front u 2 = 1 is unstable.
The similarity between Fisher-KPP and inner roots case is that both systems have two consecutive roots off (u) that coincide with the resting states of a wave front. The difference between them is that the pre-front in Fisher-KPP is unstable and the post-front is stable, while in inner roots case the pre-front is stable and the post-front is unstable.
In the "outer roots" cases III and IV as well as "the only two roots" case V, wee see from (8.1) and (8.2) that the pre-front u 1 = 0 is unstable and the post-front u 2 = 1 is stable. This matches the stability of the equilibria in Fisher-KPP model.
The marginal case II gives Re (λ 1,2 ) = 0 so the stability of the resting states cannot be established in linear approximation, and requires separate consideration. We leave this outside the scope of this paper. Table 2 sums up the results of above analysis. Table 2 The stability of the resting states in the front wave depends on the choice of the roots of the quartic polynomial. In the next section we will show the result of the numerical simulation for each case.

Choice of roots Pre-front Post-front Matching with Fisher-KPP
9. Numerical Simulations.
9.1. General settings. We simulate numerically the reaction cross-diffusion system and u 3 and u 4 are dependent parameters defined in (7.5) and (7.6). We apply no-flux boundary conditions, and the initial condition taken from the analytical solution, that is whereû andv are defined in (3.5) and (3.6).  We will show the results of the simulation for cases I, III, IV and V identified above. For each case, we pick an appropriate set of values of the free parameters to satisfy the correspoinding conditions. Table 3 lists the parameter values used and the corresponding equilibria. Note that the value of D u for Case IV in the table is given to three decimal places; in fact it was determined from the exact condition that ρ = −1/12, which implies The numerical simulations are done using finite differences, fully explicit first order for time and second order central for space. The space discretization interval is [−a, b] = [−37. 5,150] and the discretisation steps are ∆x = 0.15 and ∆t = 4 × 10 −6 unless otherwise stated. The choice of the discretization steps is motivated by the numerical stability and accuracy analysis of the scheme, which will be presented later.
9.2. The inner roots case. As shown above, in this case the pre-front equilibrium u 1 = 0 is stable, while the post-front equilibrium u 2 = 1 is unstable. Hence we expect in simulations that the post-front state evolves to another, stable equilibrium. This is indeed what happens in simulations, see Figure 4.
For the parameters used in this simulations, the unstable equilibrium u 2 = 1 is surrounded by the pre-front equilibrium u 1 = 0 and the upper stable equilibrium u 4 = 1.5. Thus in this case we expect the post-front state attracted to either of these two stable equilibria.
In fact, the solution curiously does both, i.e. attracted first to the upper stable equilibrium, u 4 = 1.5, but does not stay there for long and departs for the lower stable equilibrium, u 1 = 0. As a result, a pulse-shaped solution develops, with the pre-front and post-front states at u 1 = 0, and the plateau state near u 4 = 1.5. This phenomenology is similar to that observed in [2] for excitable (i.e. one stable equilibrium) cross-diffusion systems, incluiding oscillatory front and oscillatory back, both trigger waves from one stable equilibrium to anotherand is of course very far from the initial condition which is a monotonic front from a stable equilibrium to an unstable one.

The Result of Simulation of Distinct Real Roots, Double Roots and Complex
Roots. The behaviour of the propagating wave front for the distinct real roots case and double roots case is quite similar. The simulation shows that the numerical propagating wave remains close to the analytical wave for a period of time. Then an oscillation appears in the onset of the front. After that the oscillation grows as the time evolves, which causes the numerical solution to break up. The results of the simulation of distinct real roots case is shown in Figure 5 while the results of double roots case is shown in Figure 6.
For complex roots case, we observe that the instability occurs earlier than all previous cases (inner roots case, outer roots case and double roots case). Moreover, the numerical front does not last as long as those front waves in the other cases, see Figure 7.
10. The instability of the solution. In the previous sections we have shown the results of direct numerical simulation on reaction cross-diffusion system (9.1) where the initial condition is an exact analytical wave solution. This analytical solution presents a monotonic wave front.
We have considered four cases, corresponding to different positions of the roots of the quartic polynomial. In all four cases considered, there are oscillations which appear in the onset of the wave front. These oscillations grows as time evolves, which obviously means the propagating wave front is not stable. We now would like to address the question whether this was due to dynamical instability in the underlying partial differential equations, or numerical instability, i.e. artefact of the numerical scheme used.
Our plan on how to distinguish numerical instability from the numerical is as follows. If the instability is numerical, then its features shall significantly depend on details of the numerical scheme. For instance, the oscillations could be reduced by changing the discretisation steps. Conversely, the dynamical instability the behaviour of the solution may be affected by refining the discretisation steps only slightly, if the simulation is "resolved".   A crude theoretical analysis of numerical stability of the scheme we use can be achieved by removing the kinetic terms from system (9.1). In this way, we obtain the following For the forward-time, central-space discretization on the grid x ∈ ∆x Z, t ∈ ∆t Z, using the standard von Neumann stability analysis, for the Fourier component (u, v) ∝ e iqx we find the amplification factor ν, such that which means that the numerical scheme is unstable as the condition |ν| ≤ 1 will not be satisfied, in principle, for any choice of discretization steps. However, let us look at the quantitative aspect of the numerical instability. Namely, let us estimate the time it takes for the numerical instability to grow to macroscopic value. Supposing, for a crude estimate, that the seed of the instability comes from round-off errors, so is of the order of machine epsilon ε, and it will become significant when it grows to an order of 1. Then, with the amplification factor ν(q), the number of time steps required for that will be at least ln |1/ε| / max q (ln |ν(q)|). Taking the leading order approximation for the ln |ν(q)| in (10.1), we get the time interval required for the instability to grow to macroscopic size as By substituting the values of parameters we used in our simulation, we see that in all cases T inst is much bigger than the time T break taken for the numerical waves to break up. Table 4 clarifies more by numbers. We took ε = 10 −15 . This comparison suggests that even though the numerical scheme is unstable, the numerical instability is highly unlikely to be the reason for the behaviour observed in our numerics, and we must consider the possibility of a dynamical instability.
So, according to our plan, we have verified the plausibility of a dynamical instability by repeating the simultions at different discretization steps. We have repeated each of the simulations, once with bigger discretization steps and once with smaller discretization steps. We have found that the behaviour of the solution does not significantly change even after we refine the discretisation. More precisely, once the oscillations appear, we have found the growth rate of the oscillation is the same in all different discretisation steps. Figure 8 illustrates that for the "outer roots" case: even though the moment of onset of the instability depends on the discretization, its growth rate is not affected by it.
The same thing happened in double roots case and complex roots case. Change of discretisation steps changes the time of the onset of the instability, but not the growth rate of the instability, as can be seen in see Figures 9 and 10. For the "inner roots" case, the initial condition is a front of invasion of an unstable state into a stable state, and the numerical simulation show behaviour different from other cases: now the instability appears, at first, as the elevation of the u-field right behind the front. So we observe how this instability changes with different discretization steps. The result is shown in Figure 11. We see, again, that the time of the onset of the instability does depend on the discretization steps, but the growth rate remains the same. The subsequent behaviour of the solution also remains qualitatively similar, involving formation of a propagating pulse with a plateau and a back -even though shifted in time and differing in detail, which is of course only expectable for a solution affected by a dynamical instability.
We can conclude that insofar as it may be established by numerical simulations, the analytical front solutions are dynamically unstable.    11. Discussion. The main purpose of the paper, which has been successfully achieved, was to demonstrate the feasibility, and provide an example, of constructing a PDE model of a certain class which has desirable analytical solutions. Regardless of the utility of the particular example we have considered, we hope that the technique we used may be helpful in other problems similarly formulated.
More specifically, our aim has been a reaction-cross-diffusion system with a polynomial nonlinearity, which would have solutions in the form of a propagating front. We have found that to achieve that, the nonlinearity must be at least quartic, in which case the system may indeed have solutions in the form of a monotonic propagating fronts. The situation is similar to ZFK-Nagumo model rather than Fisher-KPP model in that for given parameters of the system, the speed and shape of the front solution are uniquely defined.
We have further established that in terms of stability of pre-front and post-front equilibria, the proposed model may be likened to the Fisher-KPP system (one of the equilibria is stable and the other unstable) but not ZFK-Nagumo (with both equilibria stable).
The quartic nonlinearity can be of various diffierent classes depending on behaviour of its four roots: when the asymptotic equilibria are two inner roots, two outer roots out of four, two outer roots out of three, the only two simple roots (with the other two being complex) and two double roots.
We have made simulations of selected examples of the proposed model belonging to different algebraic classes, and in all of these examples it happened that the analytical solutions are dynamically unstable, with some of the instabilities distinct from those related to the unstable pre-front equilibrium. Since the conclusion about instability of the solutions is based only on direct numerical simulations of arbitrarily selected examples, it requires further investigation, perhaps both theoretically and numerically, and wider parametric searches. A good survey of the relevant theory can be found in [9], and examples of numerical tools suitable for this task are AUTO [3] and WAVETRAIN [11].
Returning to feasibility of proposed PDE system as a model of real processes, we recall that KPP-Fisher is a viable model despite the unstable pre-front state. As it is well known, there are two inter-related reasons for that. One reason is the positivity of the equation: nonnegative initial conditions guarantee that the solution will remain non-negative at all times. Since the linearly unstable pre-front state is 0, i.e. at the border of the domain invariant under the system, this motivates restriction on the class of perturbations considered to those that would respect the positivity. The other reason is also related to the fact that the prefront state is 0, but is of physical rather than mathematical nature: it motivates applications in which the dynamic field represent an essentially non-negative quantity with the meaning of a concentration of some kind; specifically, in the seminal papers [4,7] it was population density. With that physical sense of the dynamic field, the magnitude of physically feasible perturbations related to fluctuations must decay as the system gets closer to the pre-front state, and exactly vanish at that state. This motivates consideration of solutions in specially constructed functional spaces that take this issue into account, in which the solution may be stable -despite the formal instability of the pre-front state in the sense of generic dynamical systems theory. In this context, the possibility of, and, as numerics show, preference for, the non-monotonic fronts is only possible because the class of model we consider does not possess the positivity property. Here we note that the models with linear cross-diffusion cannot have that property in principle, see e.g. [5].
The above consideration motivates possible continuation of the present work: • ZFK-Nagumo type fronts, i.e. monotonic fronts with stable pre-front and stable postfront states, may be sought for in models with polynomial nonlinearity of degrees higher than four; • Reasonably stable monotonic fronts switching from a zero pre-front state may be observed in models with nonlinear cross-diffusion, e.g. "pursuit-evasion" type mutual taxis of the components; • As the fronts actually observed in numerical simulations of cross-diffusion models so far are typically oscillatory, search of exact solutions of that kind would involve "inventing" an ansatz more sophisticated than that given by (3.4) and (3.5). All that should be considered in the context that the problem addressed in this paper is about the "fast subsystem" in (2.2), and encompasses just the first step in the singular perturbation theory in the limit → 0.
Finally, we note that the oscillatory pulses observed in the numerical simulations with cross-diffusion models look similar to models with "snakes and ladders" bifurcation diagrams [6,1], including reaction-diffusion systems with excitable kinetics [16]. This suggests that a similar bifurcation diagram may be expected for reaction-cross-diffusion models, which is a yet another possible direction for future research.