The Speed of Traveling Waves in a FKPP-Burgers System

Abstract

We consider a coupled reaction–advection–diffusion system based on the Fisher-KPP and Burgers equations. These equations serve as a one-dimensional version of a model for a reacting fluid in which the arising density differences induce a buoyancy force advecting the fluid. We study front propagation in this system through the lens of traveling waves solutions. We are able to show two quite different behaviors depending on whether the coupling constant $\rho$ is large or small. First, it is proved that there is a threshold $\rho _0$ under which the advection has no effect on the speed of traveling waves (although the advection is nonzero). Second, when $\rho$ is large, wave speeds must be at least $\mathcal {O}(\rho ^{1/3})$. These results together give that there is a transition from pulled to pushed waves as $\rho$ increases. Because of the complex dynamics involved in this and similar models, this is one of the first precise results in the literature on the effect of the coupling on the traveling wave solution. We use a mix of ordinary and partial differential equation methods in our analytical treatment, and we supplement this with a numerical treatment featuring newly created methods to understand the behavior of the wave speeds. Finally, various conjectures and open problems are formulated.

A Bounding the Minimum Wave Speeds Numerically

Here we briefly discuss the numerical bounding algorithm that is used to produce the values in Table 1. In complete generality the method relies on obtaining auxiliary functions defined by convex differential inequalities whose existence provides trapping boundaries defined by the level sets of the auxiliary function. To be explicit, suppose we are given an ODE

$$\begin{aligned} {\dot{y}} = F(y), \end{aligned}$$

(A.1)

with $F:\mathbb {R}^n \rightarrow \mathbb {R}^n$ and some subset of phase-space $\Omega \subseteq \mathbb {R}^n$. We aim to find a continuously differentiable function $H:\Omega \rightarrow \mathbb {R}$ whose zero level set forms a trapping boundary of (A.1) in $\Omega$. Following [4], one way to obtain such a function is to determine the existence of a constant $\lambda > 0$ such that

$$\begin{aligned} \lambda F(y)\cdot \nabla H(y) \leqq -H(y), \quad \forall y\in \Omega . \end{aligned}$$

(A.2)

Indeed, the term $F(y)\cdot \nabla H(y)$ is exactly the derivative of H along trajectories of (A.1), and so we see that the region $\{y\in \Omega :\ H(y) \leqq 0\}$ is forward invariant with respect to the dynamics of (A.1). Crucially, searching for H over some convex class of functions gives that the set of functions satisfying (A.2) is convex as well.

In the above analysis we have used the fact that traveling waves of (1.2) correspond to heteroclinic orbits of a related spatial ODE. To use auxiliary functions to either confirm existence or non-existence of heteroclinic orbits requires coupling the existence of a function H satisfying (A.2) with other constraints that are specific to the ODE in question. These conditions are detailed explicitly below as they apply to the ODEs (3.6) and (4.5). In general, providing existence of a heteroclinic orbit using (A.2) would require that the target and source equilibria of the desired heteroclinic orbit both satisfy $H(y) \leqq 0$, while the boundaries of $\Omega$ for which the unstable manifold of the source equilibrium could escape $\Omega$ have $V(y) > 0$. Proving non-existence of a heteroclinic orbit is slightly simpler since we could have the source equilibrium lying in the interior of the forward invariant region defined by $H(y) \leqq 0$, while the target lies in the complement of this set.

In the present investigation the wave speed $c > 0$ functions as an external parameter in the ODEs in which we wish to confirm the existence or non-existence of heteroclinic orbits. Given the existence of a minimum wave speed, one approach would be to find the extremal value of c for which an auxiliary function H can be obtained to confirm the existence or non-existence a desired connection in phase-space. Unfortunately, this leads to a non-convex optimization problem (see [4, Section 2c]) and so, to implement the process numerically, we perform computations at multiple fixed values of c. For example, to find the smallest value of c for which a heteroclinic connection exists, we perform the following iterative procedure. We begin with a sufficiently large value of c for which a heteroclinic orbit exists and another sufficiently small value for which it does not. We can then repeatedly bisect in c, attempting to find a (potentially) different auxiliary function H that confirms the existence of a heteroclinic orbit in $\Omega$ at each new value of c. The smallest such c for which this can be performed then becomes an upper bound on the minimum wave speed. Confirming non-existence is similarly performed through such a bisection method in the wave speed c. All other parameters in the system (i.e. $\nu ,\rho > 0$) are held fixed.

We refer the reader to [4] for a more complete discussion of the above method of bounding the wave speed and how to implement this procedure numerically. Briefly, after determining the inequalities that a desired auxiliary function H must satisfy to confirm existence or non-existence of a heteroclinic orbit, we search for H numerically as a degree $d \geqq 1$ polynomial in y with tunable real-valued coefficients. These convex inequalities for H in the space $\mathbb {R}[y]_d$, the set of all polynomials with real-valued coefficients and degree $\leqq d$, defines a semidefinite program that can be relaxed to a series of sum-of-squares constraints that are numerically tractable. The numerical implementation may either find admissible values for the coefficients of H or return that no such values exist. The results in Table 1 were obtained using the MATLAB software YALMIP (version R20190425) [27, 28] to translate the sum-of-squares constraints into semidefinite programs which are then solved using Mosek (version 9.0) [34]. The code to reproduce the values in Table 1 is available at the repository GitHub/jbramburger/FKPP-Burgers. Computations are performed by optimizing over $\lambda > 0$ at each fixed degree $d \geqq 1$ of H. The degree is successively increased until convergence in the bound is observed, typically around $d = 10$ to $d = 14$.

We conclude this section with a brief discussion of the specific conditions put on H to bound the wave speed from above and below in the cases $\nu = 0$ and $\nu \ne 0$, respectively.

$\mathbf{Upper Bounds on }$ $c_*(\rho )$: From the work of Section 3, our ODE of interest is the planar system (3.6) and the region of interest is $\Omega = \mathcal {B}$, defined in (3.8), for each fixed $\rho > 0$ and $\nu = 0$. From Proposition 1 we have that the unstable manifold of $(1,u_0)$ can only leave $\mathcal {B}$ by crossing $U = 0$, and so we impose the conditions

$$\begin{aligned} \begin{aligned}&-\lambda \bigg [(c-U)\bigg (-cT + UT + \frac{U}{2\rho }(2c - U)\bigg )H_T(T,U) \\&+ \rho T(T-1)H_U(T,U)\bigg ] \\&\geqq (c- U)H(T,U), \\&H(T,0) \geqq \varepsilon T(1-T), \\&-\varepsilon \geqq H(1,u_0) \\ H(0,0)&= 0, \end{aligned} \end{aligned}$$

(A.3)

for all $(T,U)\in \mathcal {B}$, where subscripts denote partial differentiation. Indeed, the first condition is (A.2) for the specific ODE (3.6), multiplied through by the positive term $c-U$ to maintain that the inequality is entirely stated in terms of polynomial functions when V is polynomial. For any $\varepsilon > 0$, the second condition guarantees that the set $\{(T,0):\ 0< T < 1\}\subset \mathcal {B}$ is such $H(T,0) > 0$ for all $0< T < 1$, thus confirming that the unstable manifold of $(1,u_0)$ cannot cross $U=0$. The results in Table 1 always use $\varepsilon = 10^{-4}$. We comment that the presence of $\varepsilon > 0$ comes from the fact that strict inequalities cannot be handled numerically. The third condition gives that the equilibrium $(1,u_0)$ lies in the interior of the forward invariant region, thus confirming that its unstable manifold in $\mathcal {B}$ also lies in this set. Finally, the fourth condition gives that the origin lies on the boundary of the forward invariant region. Finally, we note that we cannot have $H(0,0) \leqq -\varepsilon$ since this would be inconsistent with the second condition, thus leading to the requirement that $H(0,0) = 0$, as stated above.

$\mathbf{Lower Bounds on }$ $c_*(\rho )$: Again we take $\Omega = \mathcal {B}$ and consider the ODE (3.6). The conditions to confirm non-existence of a heteroclinic orbit for each fixed $\rho > 0$ are then

$$\begin{aligned} \begin{aligned}&-\lambda \bigg [(c-U)\bigg (-cT + UT + \frac{U}{2\rho }(2c - U)\bigg )H_T(T,U)+ \rho T(T-1)H_U(T,U)\bigg ]\\&\quad \geqq -(c- U)H(T,U), \\&\quad -H(T,0) \geqq \varepsilon U(1-U), \\&H(1,u_0) = 0, \\&H(0,0) \geqq \varepsilon , \end{aligned} \end{aligned}$$

(A.4)

for all $(T,U)\in \mathcal {B}$. The first condition is a variant of (A.2), in this case of the form

$$\begin{aligned} \lambda F(y)\cdot \nabla H(y) \leqq H(y), \end{aligned}$$

(A.5)

which again gives that $H \leqq 0$ is forward invariant with respect to the dynamics of F. The reason for this change in conditions is to maintain that the third and fourth conditions, which separate the target and source equilibria, are consistent with the condition guaranteeing forward invariance of $H \leqq 0$. Indeed, using (A.2) in the place of the first condition above would require that $H(0,0) = 0$, thus not necessarily separating the equilibria and in turn not necessarily prove non-existence of a heteroclinic orbit since both equilibria would belong to the closed set $H \leqq 0$. The second condition gives that the line $\{(1,U\, 0< U < 1\}\subset \mathcal {B}$ lies in the forward invariant region, and since $(1,u_0)$ is a saddle, a local analysis near this equilibrium reveals that its unstable manifold in $\mathcal {B}$ must enter into the forward invariant set $H \leqq 0$. As in the upper bounds on $c_0(\rho )$, the second condition guarantees that $H(1,u_0)$ cannot be taken to be negative, and therefore we can only have $H(1,u_0) = 0$, in turn necessitating the alternative first condition.

$\mathbf{Upper bounds on }$ ${\overline{c}}_*(\nu ,\rho )$: Here now we consider the ODE (4.5) with $\Omega = \mathcal {C}$ for $\nu ,\rho > 0$ fixed. From Proposition 2 we have that the unstable manifold of $(1,u_0,c-u_0)$ that enters into $\mathcal {C}$ can only leave by crossing $V = 0$. We then seek $H:\mathbb {R}^3\rightarrow \mathbb {R}$ satisfying

$$\begin{aligned} \begin{aligned}&-\lambda \bigg [(-cT + UT + V)H_T(T,U,V) + \frac{1}{\nu }\bigg (-cU + \frac{1}{2}U^2 \\&\quad + \rho V\bigg )H_U(T,U,V) + T(T-1)H_V(T,U,V)\bigg ] \geqq H(T,U,V), \\&H(T,U,0) \geqq \varepsilon T(1-T) + \varepsilon U(u_0 - U), \\&\quad -\varepsilon \geqq H(1,u_0,c-u_0), \\&H(0,0,0) = 0, \end{aligned} \end{aligned}$$

(A.6)

for all $(T,U,V)\in \mathcal {C}$. The first condition is exactly (A.2) with F specified by the right-hand-side of (4.5). The remaining conditions are analogous to those for the upper bounds on $c_0(\rho )$.

$\mathbf{Lower bounds on }$ ${\underline{c}}_*(\nu ,\rho )$: Again we consider the dynamics (4.5) with $\Omega = \mathcal {C}$ for fixed $\nu ,\rho > 0$. To determine non-existence of a heteroclinic connection from $(1,u_0,c-u_0)$ that remains in $\mathcal {C}$ we seek a function $H:\mathbb {R}^3\rightarrow \mathbb {R}$ satisfying

$$\begin{aligned} \begin{aligned}&-\lambda \bigg [(-cT + UT + V)H_T(T,U,V) + \frac{1}{\nu }\bigg (-cU + \frac{1}{2}U^2\\&+ \rho V\bigg )H_U(T,U,V) + T(T-1)H_V(T,U,V)\bigg ] \geqq H(T,U,V), \\&\quad -\varepsilon \geqq H(1,u_0,c-u_0), H(0,0,0) = 0, \end{aligned} \end{aligned}$$

(A.7)

for all $(T,U,V)\in \mathcal {C}$, along with the condition

$$\begin{aligned} H(T,U,V) \geqq 0, \quad \forall T,U,V\in [0,\delta ]. \end{aligned}$$

(A.8)

The first condition above represents (A.2), while the following two conditions put the source equilibrium $(1,u_0,c-u_0)$ in the interior of the forward invariant region and the origin on the boundary. For some sufficiently small $\delta > 0$, the condition (A.8) works to guarantee that a region around the origin lies outside of the forward invariant region $H \leqq 0$ when H is a polynomial in (T, U, V). In our numerical implementations we take $\delta = 0.05$; larger values of $\delta$ lead to less precise bounds, while smaller values show little change in the lower bound.