Characterization of Traveling Waves Solutions to an Heterogeneous Diffusion Coupled System with Weak Advection

: The aim of this work is to characterize Traveling Waves (TW) solutions for a coupled system with KPP-Fisher nonlinearity and weak advection. The heterogeneous diffusion introduces certain instabilities in the TW heteroclinic connections that are explored. In addition, a weak advection reﬂects the existence of a critical combined TW speed for which solutions are purely monotone. This study follows purely analytical techniques together with numerical exercises used to validate or extent the contents of the analytical principles. The main concepts treated are related to positivity conditions, TW propagation speed and homotopy representations to characterize the TW asymptotic


Problem Description and Objectives
Typically, the models involving spatial diffusion have been considered as coming from simple physical principles. This is the case of the Fick law that establishes a relation between the variable flux in a media and the gradient of the variable concentration. The application of such law leads to the classical gaussian order two diffusion. Nonetheless, in applied areas such as biology, optics, structures or materials, the gaussian diffusion has been extended to account for new ways of modeling introducing high order diffusion operators. Such operators are currently subjected to intensive research, as an example, Bonheure and Hamel have shown the De Giorgi's conjecture for a fourth order Allen-Cahn equation together with classical solutions bounds [1]. The Allen-Cahn elliptic equation is used to model stationary bi-stable systems in physics, chemistry or biology.
In some practical cases, the fourth-order operators emerge from already known order two diffusion. As an example, the classical Fisher-Kolmogorov equation was proposed to study the interaction of different populations in a biological environment: It has been observed an onset of instabilities near degenerate points given by expression (1) ( [2] and references listed there), which led to propose the Extended Fisher-Kolmogorov equation to model the behaviour of bi-stable systems. These systems can be defined as those with only two uniform states and a solutions "traveling" from one stable solutions to the other, either forming a heteroclinic or homoclininc orbit [3]. In [4,5], Peletier and Troy, on one hand, and Bonheure, on the other hand, showed the existence of oscillatory spatial patterns for the Extended Fisher-Kolmogorov equations. Additionally, they exhibited examples of oscillating heteroclinic (the author also called it kinks) and homoclinic orbits (pulses) in the spatial domain. The instabilities were found to be permanent oscillations, leading to think that there shall be evolution flows hidden by the regularity of the second order diffusion. Therefore, the original Fisher-Kolmogorov equation was perturbed with a fourth order spatial derivative, leading to the Extended Fisher-Kolmogorov equation: where ∆ 2 = ∂ 4 ∂x 4 . In the classical sense, the Extended Fisher-Kolmogorov equation requires solutions to have continuous derivatives up to the fourth order. One can think, preliminary, that oscillating functions (such as sine, cosine or a combination of both) may be appropriate candidates to constitute solutions. Peletier and Troy showed the existence of oscillating solutions [5]. In addition and making use of a development in the exponential bundle of solutions, Rottschäfer and Doelman showed the nature of such oscillations [2].
A previous work [6] developed a set of analysis about the existence of minimal heteroclinic orbits for a class of fourth-order ODE systems (not necessarily cooperative) with variational structure. In the present analysis, we develop hetereclinic orbits for a cooperative system within the PDE theory and making use of analytical and numerical evidences to account for solution profiles in the Traveling Waves domain. In addition, the analysis focuses on the construction of exponential bundles of solutions that are represented through homotopy graphs. To illustrate the relevance of the exponential bundles of solutions, recently, Dang [7] has provided a general method on the complex plane to analyze exponential solutions (and others rational and elliptic) for the (2 + 1)-dimensional, the (3 + 1)-dimensional Boiti-Leon-Manna-Pempinelli equations and the (2 + 1)-dimension Kundu-Mukherjee-Naskar equation. In addition, in [8] the author studies Traveling solutions to the non-local Fisher-KPP equation considering only the kernels for which the spatially uniform steady state u = 1 is stable. Further, the search of waves propagation for a damped wave equation has been explored in [9] for a fractional Laplacian with nonlinear source. The existence of solutions is assessed based on Galerkin approximations combined with the potential well theory to show the decay behaviour of solutions. Along the presented analysis, such decay behaviour is presented within the Traveling Waves solutions and exponential bundles. In addition, it is the intention to search for the most appropriate TW-solution with positivity and homogeneous convergence towards stationary solutions. In this sense, some previous works shall be mentioned. In [10], the authors develop multiscale methods and asymptotic analysis to understand the homogenization in domains with heterogeneous strips. With the same intention, an analysis in [11] aims to characterize homogeneous processes in heterogeneous reaction-diffusion environments. In [12], compactness criteria are employed to characterize the homogenization of a diffusionconvection equation with divergence-free velocities. Further, some applications to other disciplinary subjects of homogenizing techniques in heterogenous porous media can be consulted in [13,14] that illustrate the relevance of the topic.
Fourth-order operator equations have been assessed by Lyapunov stability approaches [15], in which the existence of bifurcation branches for even, periodic solutions in both the Swift-Hohenberg and extended Fisher-Kolmogorov equations have been considered. During the present analysis, we use the homotopy analysis in stead of pure Lyapunov methods.
The high-order operator induces a set of instabilities in the proximity of the stationary solutions. Recently, Díaz and Naranjo [16] have shown the oscillating behaviour of selfsimilar solutions and have characterized regions of positivity for a class of high order cooperative system with no advection. This work permits to introduce insights on the instabilities shown by spatially inhomogeneous structures when they are modeled by diffusion (see in [17] and references there in for further details). Particularly, it is the intention to characterize the propagation features of the heteroclinic orbits connecting two spatially homogeneous solutions anticipated by the cooperative system formulation. Note that cooperativeness shall be understood as the synergistic collaboration between species to prosper and grow in a territory. Our observations will be made in the traveling wave domain, and mainly, in the traveling waves fronts and tips where the transition involving exponential bundles of solutions happen. The set of equations can be summarized as where ∈ R and sufficiently small for our purposes, with initial conditions The minus sign before the fourth order spatial derivative is set to account for a regular asymptotic stable system. As discussed, such degenerate diffusion aims to introduce oscillatory patterns close the equilibrium so as to model a center manifold (in the sense of oscillatory) behaviour. The terms v and u in the u-equation and the v-equation forcing terms, respectively, account for the coupled cooperation between species. The advection is introduced to account for a certain preferred direction in the space, for instance, a direction of food and/or shelter.
It is to be noted that terms u(a − u) and v(b − v) have been dealt previously, for a single-order two-diffusion equation, by Kolmogorov, Pretrovskii and Piskunov, leading to classical and so-called KPP problem of order two [18]. The KPP term is typical in biological systems (also called Allee effect) to model the birth, growth and death in species.
This analysis pretends to go beyond existence and uniqueness, providing some solutions profiles obtained by a combination of numerical and analytical methods. Furthermore, this study provides a stability analysis in the Traveling wave domain by homotopy representations.

Existence and Traveling Waves Structure
The TW formulation for problem Q (3) consists on operating in the TW variable (y) and profiles ( f , g), given by the following expressions [19]: where λ 0 is the propagation speed.
The system (Equation (3)), then, reads as Call λ = λ 0 + . The TW formulation is subjected to the pseudo-boundary conditions: First, and aiming to understand the basic TW profiles features, we enunciate the Lemma 1: Lemma 1. The traveling waves moves from y → −∞ to y → ∞. In other words, the wave speed λ is positive. (6) by f , the second one by g and integrating by parts (note the calculations are presented only for the first equation):

Proof. Multiplying the first equation in Equation
such that Each of the involved integrals needs to be assessed with certain conditions at −∞ and ∞ to be specified: The integral is assessed between −∞ and +∞, where we admit the following approximations: Therefore, We continue assessing the next involved integral: such that: The process is followed for the next integral involved: In the assessment of the previous integral, we are interested on determining the sign character rather than a precise value. The cooperative state makes the solutions to evolve closely upon selection of an appropriate value in the TW-speed (to this end, refer to the analysis shown afterwards in Figures 1 to 5). Then, assume that close the equilibrium conditions at −∞ and ∞ the following holds f g ∼ f g . Returning to Equation (15): Finally, we repeat the same integration by parts for the last of the integrals involved: Note that the integral on right has been already assessed in Equation (14). Then, The compilation of the assessed integrals yields Considering a > 0, immediately λ > 0. The same process can be followed for the second equation of Equation (6) to obtain Moreover, considering b > 0, λ > 0 as well.
To further characterize the TW existence, the system in Equation (6) is converted into a linear system by the standard change of variables: This can be expressed as Note that the partial derivatives with respect to f i and g i for i = 1, 2, 3, 4 of each component on the right hand vectorial function in Equation (23) are continuous which is needed to ensure the Lipschitz condition. In addition, and to support the global existence of TW profiles f , g, the cooperative system in Equation (6) is solved with a numerical algorithm.
One of the key questions to answer is whether it exists a minimal TW speed, λ, such that solutions are stable in the proximity of the stationary solutions. The minimal speed existence, for which non-oscillating behaviour is given, is common to the KPP order two problems [18]. Nonetheless, for higher order systems, it is not possible to ensure the existence of a minimal speed with a monotone decay at infinity, due to the existence of oscillating behaviour. Following the KPP order two philosophy [18], this minimal speed is the critical speed at which any oscillation in the solution vanishes. Nonetheless, in this case, as we are dealing with fourth-order operators, this is not guaranteed, being a notable difference compared to the KPP models involving a second order parabolic operator. Thus, if a minimal speed is not guaranteed, we can search for other common property of the TW profiles in the KPP order 2 problems. In particular, when the TW moves at the minimal-critical speed, the solution is positive everywhere for all y ∈ R. Our target can be translated into finding a suitable value of the TW-speed for which the first minimum in y > 0 is positive. Nonetheless, we cannot ensure the positivity of the solution for y ∈ R as the natural instabilities in the high order operator impedes the possibility of a maximum principle.
The estimation of a TW-speed, at which the first minimum in y > 0 is positive in Equation (6), is done via a numerical algorithm. The numerical analysis has been done over a sufficiently large y interval [−1000, 1000] to avoid the influence of the pseudo-boundary conditions (Equation (7)) in the integration domain. The relative and absolute errors for each iteration has been considered as 10 −6 and the number of nodes varies from 10 4 to 10 5 .
Additionally, the numerical results provide the evidence related with the global structure of the both f , g for different values of the TW speed λ.
First, the TW-speed λ is assumed to be equal for both solutions f , g in which the condition a = b is considered. Previous to any formal proof, Figures 1 and 2 suggest that the oscillatory character of the TW decreases for increasing values of λ. This property is common with the KPP-2 problems [18].

Conjecture 1.
Let consider a = b in the set of Equation (6) and that the TW-speed is common to both solutions f , g. Let define the set M as representing the location of the TW front along the y axis and let define expressing a set whose elements locate the minimum points in solution f (y) beyond the TW-front (i.e., in the TW-tail). In these conditions there exist a value of λ for which f (min(m)) > 0. This value for λ represented by λ f (min(m))>0 has been sharply estimated to be λ f (min(m))>0 = 2.394.
Note that the conjecture is postulated based on the numerical evidences results presented in Figure 3.
Conjecture 1 is particularly relevant and expresses some analogy with the minimal TW-speed (λ min ) typical of the TW associated to the KPP problems of second order [18]. In the high-order case, we cannot directly consider the minimal TW-speed; nonetheless, we have shown the existence of an equivalent λ f (min(m))>0 in the cooperative system. Indeed, the most remarkable difference between the KPP-2 (KPP order two) TW and the cooperative system with a high order operator relays on the fact that we shall change the concept of finding a λ min by finding a λ f (min(m))>0 . For a KPP-2 TW moving at λ min , the profile does not oscillate as the second order profiles changes from the sub-critical solution to the critical one [18]. In the fourth-order cooperative system, we use the concept of λ f (min(m))>0 , for which we know that the oscillatory behavior is observed in the TW-tail when y >> 1 (see Figure 3). This is a common behavior to all high-order parabolic operators. Thus, we have found a λ f (min(m))>0 for which, in an appropriate inner region (inner compared to y >> 1), it is possible to express the high order TW profile in a similar manner as done in the KPP-2 problem. Note that the results obtained, up to now, consider that both TW profiles f , g move with the same speed λ. It is considered, now, the possibility of different TW-speeds resulting in the following cooperative system: We keep on the philosophy of finding suitable TW speeds (in this present case different for f , g) for which the conditions f (min(m)) > 0, g(min(m)) > 0 hold. The property of finding λ f (min(m))>0 and λ g(min(m))>0 for each TW profile ( f , g) is now more subtle. We start by assuming that λ 1 = 2.394 and we try to answer the following question: Is there an interval in λ 2 for which there exist λ 1, f (min(m))>0 and λ 2,g(min(m))>0 ?
To answer this question, the cooperative system has been numerically modeled. The following proposition in the form of conjecture compiles the numerical evidences.
Admit the following discussion to support the Conjecture enunciation. Consider that a = b in the set of Equation (6). The proof of this proposition is given in Figure 4 (left). It is convenient to highlight that there exists TW profiles for other λ 2 > λ 2,g(min(m))>0 = 2.437 and λ 1 > λ 1, f (min(m))>0 = 2.394 (see Figures 4 (right) and 5 as examples), nonetheless these profiles are oscillatory from the first minimum. Now, the intention is to explore the effect of the advection term over the positivity regions in the TW characterization. To this end, admit that the ∼ λ 0 , i.e., the advection is considered to interact significantly with the TW propagation. As a consequence, consider λ 1 = λ 0,1 + , λ 2 = λ 0,2 + , then while keeping positivity in the first minimum as described in Figure 4.

Traveling Waves Homotopy Discussion
In the asymptotic case with y → ∞, the cooperative system (Equation (6)) can be linearized, leading to the following expression: First, to simplify the operations involved and without loss of generality, we consider a = b = 1. The linearization exercise permits to represent the system (Equation (28)) by making use of the first order system autonomous matrix: Or equivalently, The parametric analysis in the matrix A is very complex as it involves four parameters over a eighth order matrix. Thus, we will consider some specific values for the parameters involved, so that we account for a sufficient set of evidences to determine the asymptotic behaviour as y → ∞. Additionally, the fact of having linearized, resulting in the matrix A, permits to proceed via homotopy graphs to further determine the behaviour of the exponential bundles associated to the matrix A eigenvalues.
For this case, the characteristic polynomial of A for different values of λ reads as For a generic λ the characteristic polynomial of the matrix A reads We represent now the polynominal solutions for different values in λ to check the homotopy evolution starting from λ = 0. The results are summarized in Lemma 2.

Lemma 2. (A)
For any λ > 0, the linearized system (Equation (28)) (obtained in the limit with y → ∞) with a = b = 1 and λ 1 = λ 2 = λ presents a 6-D stable family of solutions (corresponding to two pairs of complex conjugates: one constant and one real negative eigenvalues); and a 2-D unstable family of solutions (corresponding to two different real solutions). (B) The eigenvalues tend to accumulate into four different clusters with increasing distance among them when λ → ∞. One cluster is formed by the two eigenvalues with Re < 0, Im > 0, another cluster with Re < 0, Im < 0 and another one with Re > 0, Im = 0, (see Figure 6, and for increasing values of λ, see Figure 7 with λ = 60 and Figure 8 with λ = 100 and λ = 1000). (C) There exists the null eigenvalue for any value of λ.   Proof. The determination of the eigenvalues in the matrix (Equation (29)) provides the existence of complex eigenvalues which introduce oscillatory bundles in the proximity of the null critical point. This feature is common to all high-order operators (see in [20] and references therein) and expresses the difficulty to determine a pure monotone TW at a suitable speed λ = 0).
Under Lemma 2, it is possible to check that the oscillatory solutions cannot be avoided for a suitable value of the TW-speed λ. Even further, one can check by a careful look to Figures 6-8 that the imaginary parts of the complex eigenvalues increase when the TW-speed increases.
The characteristic polynomial for the matrix A adopts the following structure: In this case ab = 1, the null eigenvalue is not a solution of Equation (33), as in the case when ab = 1. Again, the asymptotic behaviour of solutions is obtained for particular values of the parameters a, b and λ due to the difficulty of performing a parametric root analysis in Equation (33).
Note that, with the aim of understanding the different solution bundles for different actual values of a, b and λ, we analyze first if Equation (33) has a pure imaginary root of the form µ = ki for λ = 0 and k = 0 ∈ R: rearranging terms: From the imaginary part, the two real solutions k = ( a+b 2 ) 1 4 and k = −( a+b 2 ) 1 4 are obtained. Substituting the first root into the real part we arrive to the following equation: In the limit with a → 0 the following equations holds which has not positive real roots (b > 0). One can check, after doing the first derivative, that Equation (38) is always positive for b > 0.
On the opposite side, if we consider values for both a → ∞ and b → ∞, the leading term in Equation (37) is the product ab, therefore the function J(a, b) evaluated under these circumstances is negative.
Given the positivity of the function J(a, b) for sufficiently small values of a and b, the negativity of the same function for sufficiently high values of a, b and the continuity properties, it can be concluded the existence of a combination of a, b satisfying the condition J(a, b) = 0. Given such values for (a, b), the imaginary solution is given by The existence of such pure imaginary solution in Equation (37) can be further characterized by representing the homotopy graph for a fixed random value of the TW-speed (assumed to be λ = 2). Figures 9-11 confirm, indeed, the existence of such pure imaginary root in Equation (37) for a → ∞ and b → ∞.
In addition, and for λ = 0 the homotopy graph presents pure imaginary solutions, see Figure 12.    The oscillatory character of the solutions f , g increases with λ → 0. Figure 13 represents the evolution of a TW for a sufficiently small value of λ together with the homotopy diagram. The eigenvalues satisfying Re < 0 has Re → 0. The consecutive figures 14 and 15 give the evolution of the TW and the homotopy graph for different values of the TW-speed.   For λ → ∞ and a = b in the same order a b, the homotopy graph tends to dense and converges to four eigenvalues with the following properties: • Two eigenvalues densing to a homotopy point with Re < 0 and Im > 0. • Two eigenvalues densing to a homotopy point with Re < 0 and Im < 0. • Two eigenvalues densing to a homotopy point with Re > 0 and Im = 0. • Two eigenvalues densing to a homotopy point with Re → 0 and Im → 0.
The analysis considers the case with a = 1, b = 2 and different values for λ 1 = λ 2 within the same order of magnitude and with substantially different order of magnitudes.
The case with a = 1 and b = 2 is analyzed considering two TW-speed cases: the case with λ 1 = 0.1 and λ 2 = 100 is represented in Figure 16 and the case with λ 1 = 100 and λ 2 = 0.1 in Figure 17.
The homotopy graphs, for this third case, are formed by two eigenvalues with Re > 0 and Im = 0 (Unstable) and six eigenvalues with Re < 0.

Conclusions
The objectives highlighted in relation with the searching, finding and characterization of TW profiles have been proposed based on analytical and numerical evidences. The main question related with the existence of TW speeds for which the natural instabilities of the four order operator are minimized has been answered. This fact has permitted to extend the classical concept of TW speed for order two diffusion as done by the KPP models [18]. The weak advection term has been shown to be bounded by the own TW propagation speeds determined so as to keep positivity in the first minimum. In addition, the homotopy representation has permitted to characterize the TW tail in the proximity of the stationary solutions in the cooperative system Q and to determine the asymptotic behaviour of TW solutions via homotopy analysis.
Funding: This research has be supported by the Universidad Francisco de Vitoria Research Directorate.
Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Data is made available upon request.

Conflicts of Interest:
The author declares no conflicts of interest.