Multidimensional Stability of Planar Traveling Waves for Competitive–Cooperative Lotka–Volterra System of Three Species

Abstract

We investigate the multidimensional stability of planar traveling waves in competitive–cooperative Lotka–Volterra system of three species in n-dimensional space. For planar traveling waves with speed $c>c^{*}$, we establish their exponential stability in $L^{\infty }\left({\mathbb{R}}^n\right)$, which is expressed as $t^{-\frac{n}{2}}{\mathrm{e}}^{-{\epsilon }_{\tau }\sigma t}$, where $\sigma >0$ is a constant and ${\epsilon }_{\tau }\in (0,1)$ depends on the time delay $\tau >0$ as a decreasing function ${\epsilon }_{\tau }=\epsilon \left(\tau \right)$. The time delay is shown to significantly reduce the decay rate of the solution. Additionally, for planar traveling waves with speed $c=c^{*}$, we demonstrate their algebraic stability in the form $t^{-\frac{n}{2}}$. Our analysis employs the Fourier transform and a weighted energy method with a carefully chosen weight function.

1. Introduction

There have been a lot of achievements in the traveling wave solutions of reaction–diffusion equations since the pioneering contributions of Fisher [1] and Kolmogorov [2]. The parabolic development equation has been studied widely because it can describe many reaction–diffusion phenomena in combustion theory, ecology, infectious diseases, and so on (see Buckmaster and Ludford [3], Cheng and Yuan [4], da Silveira and Fontanari [5], Gallay [6], Schaaf [7], Smith and Zhao [8], Volper et al. [9], Wu et al. [10], Yu and Yuan [11], and Zhao and Ruan [12]). One of the most fundamental problems within the field of traveling waves is the stability of traveling wave solutions. This is a pivotal issue in the studies of traveling waves and represents a highly intricate and challenging area of investigation. Over the past decade, the study of multidimensional stability in traveling wave solutions has attracted considerable interest. As to the following reaction–diffusion equation

$$\frac{\partial u(t,x)}{\partial t}=\Delta u(x,t)+f\left(u(x,t)\right),x\in {\mathbb{R}}^n,t>0,$$

(1)

for some $\theta \in (0,1/2)$, Xin [13] first considered the multidimensional stability of planar traveling waves in the case of a bistable reaction–diffusion equation. Specifically, he obtained that in the case where the perturbation of a planar traveling wave is small enough in $H^m\left({\mathbb{R}}^n\right)\cap L^1\left({\mathbb{R}}^n\right)(m\ge n+1,n\ge 4)$, the solution of the initial value problem will converge to the planar traveling wave in $H^m\left({\mathbb{R}}^n\right)$ at the rate of $O(t^{-\frac{n-1}{4}})$, based on the linear operator semigroup theory. Here, a perturbed state is a condition of a system that is stable under small perturbations but can transition to a more stable state if sufficiently disturbed. This state exists between stability and instability, allowing a system to remain unchanged over time under different initial perturbations. Aiming at the same problem, Levermore and Xin [14] further proved the stability of the planar traveling wave in $L_{loc}^2\left({\mathbb{R}}^n\right)$ with $n\ge 2$. The maximum principle and spectral theory, which makes a difference in the construction of a Lyapunov function, plays the key role in the proof. Smith and Zhao [8] established the global stability of traveling waves with regard to a reaction–diffusion equation, with time delay using squeezing methods, comparison principle, and upper and lower solutions. In addition, Matano et al. [15] showed that the planar traveling waves of (1) are asymptotically stable under almost periodic perturbations or sufficiently large initial perturbations that decay at spatial infinity. However, they also discovered a specific solution that oscillates indefinitely between two planar traveling waves, revealing that planar traveling waves lose asymptotic stability under more general perturbations. Later, Matano and Nara [16] expanded these results, proving that planar traveling waves are asymptotically stable in $L^{\infty }\left({\mathbb{R}}^n\right)$ under spatially ergodic perturbations, a broader class that includes quasi-periodic and almost periodic perturbations. Recently, Volpert and Petrovskii [17] made a complete review of the development trend and the latest contribution of reaction–diffusion wave in biology. More precisely, they provided a deeper insight into wave propagation, stability, and bifurcation. They also discuss the new model and its influence on understanding the process of biological invasion or disease spread. For a more comprehensive understanding, the following findings provide crucial insights, including Aronson and Weinberger [18], Bates and Chen [19], Kapitula [20], Matano and Nara [16], Matano et al. [15], Sheng [21], Volpert [22], Zeng [23], Vidyasagar [24] and the references therein.

Introducing time delay into the nonlinear reaction term of reaction–diffusion equations plays a critical role in capturing the dynamics of processes, where the response of the system depends on past states. Time delays are particularly relevant in modeling biological, chemical, and ecological systems where delays arise naturally, including in the gestation of biological species, chemical reaction processes, or resource replenishment in ecosystems. By incorporating delays, models account for the temporal gap between cause and effect, offering a more realistic depiction of the underlying phenomena. Mei and Wang [25] studied a Fisher-KPP type reaction–diffusion equation with time delay given by

$$\frac{\partial u(t,x)}{\partial t}=D\Delta u(t,x)-d\left(u(t,x)\right)+{\int }_{{\mathbb{R}}^n}{f}_{\alpha }\left(y\right)b\left(u(t-\tau ,x-y)\right)\mathrm{d}y,x\in {\mathbb{R}}^n,t>0,$$

(2)

where $D>0$ represents the diffusion rate and $d\left(u\right),b\left(u\right)$ are nonnegative nonlinear functions. By the method of the Fourier transform and weighted energy method, they proved that all noncritical planar traveling waves are exponentially stable, while critical planar traveling waves exhibit algebraic stability of the form $t^{-\frac{n}{2}}$. Huang et al. [26] further generalized these results to nonlocal diffusion equations, expanding the understanding of stability in systems with time delays and nonlocal effects. More related results can be found in Freedman and Gopalsamy [27], Gourley [28], Li et al. [29], Lin et al. [30], Yu and Mei [31], Faye [32] and so on. Considering both time delay and diffusion, we usually obtain time delay systems with nonlocal interaction. This occurs because populations at a specific spatial location at time $t-\tau$ become distributed across all spatial locations at time t as a consequence of diffusion. Chen and Shi [33] developed a reaction–diffusion equation incorporating nonlocal effects and logistic growth as follows, offering a more realistic framework for modeling complex spatiotemporal dynamics,

$$\left\{\begin{array}{cc}\frac{\partial u(x,t)}{\partial t}=d\Delta u(x,t)+\lambda u(x,t)\left(1-{\int }_{\Omega }K(x,y)u(y,t-\tau )\mathrm{d}y\right),\hfill & x\in \Omega ,t>0,\hfill \\ u(x,t)=0,\hfill & x\in \partial \Omega ,t>0,\hfill \end{array}\right.$$

(3)

where $\Omega$ represents a connected bounded open domain in ${\mathbb{R}}^n$, $\partial \Omega$ is a smooth boundary, $\lambda$ is the specific constant, and $K(x,y)$ denotes the possible diffusion of mature population. They obtained the stability and Hopf bifurcation of (3). For more from the literature, see [34,35,36,37,38,39].

On the Lotka–Volterra competitive reaction–diffusion system, Ruan and Zhao [40] conducted an in-depth study focusing on two-species models with time delay, including the predator–prey and competitive systems. They successfully established criteria for uniform persistence and global extinction. Building on this foundation, Wang and Lv [41] further investigated entire solutions of a diffusive Lotka–Volterra competitive model with nonlocal delays. Through the use of comparison principles and the upper–lower solution method, they demonstrated the existence of traveling waves. In terms of the cooperative system, Abdurahman and Teng [42] considered an n-species Lotka–Volterra cooperative system and proved the sufficient conditions of uniformly strong persistence, uniformly weak average persistence, and uniformly strong average persistence of population. In what follows, Li and Wang [43] studied the Lotka–Volterra model of diffusion cooperation with distributed delay and nonlocal spatial effects. By choosing different kernels, they use an iterative technique to establish sufficient conditions for the existence of traveling wave solutions, which connect the zero equilibrium and the positive equilibrium. In addition to these works, extensive research on the Lotka–Volterra competitive or cooperative system has been carried out, as detailed by Chermiha and Davydovych [44], Tian and Zhao [45], Du and Ni [46], Han et al. [47], Lin and Li [48], Ma et al. [49], Tang et al. [50], Wang et al. [51], who offer valuable insights into the broader dynamics and applications of these models. Although it is abstract, studying ecological problems in $n>3$ dimensions holds significant theoretical and practical value. Higher-dimensional models can incorporate complex factors such as environmental gradients, species traits, or temporal dynamics, providing a more comprehensive understanding of ecological processes like adaptive evolution or biodiversity maintenance. Furthermore, many ecosystems are best represented as complex networks, where high dimensional abstractions capture intricate interactions such as species dispersal or habitat fragmentation. The study of such systems in higher dimensions also reveals universal scaling laws and critical thresholds, which enhance the generality and predictive power of ecological theories. While the multidimensional stability of planar traveling waves for scalar reaction–diffusion equations has been extensively studied, there is limited research on the stability of such waves in systems, particularly those involving nonlocal nonlinearity in higher-dimensional spaces. In population biology, competition and cooperation are fundamental interactions between species. Thus, in this paper, we study the following competitive and cooperative Lotka–Volterra system with the diffusion and nonlocal nonlinearity

$$\left\{\begin{array}{c}\frac{\partial u_1(t,x)}{\partial t}=d_1\Delta u_1(t,x)+r_1{u}_1(t,x)[1-a_{11}{u}_1(t,x)+a_{12}g\ast u_2(t-\tau ,x)-a_{13}g\ast u_3(t-\tau ,x)],\hfill \\ \frac{\partial u_2(t,x)}{\partial t}=d_2\Delta u_2(t,x)+r_2{u}_2(t,x)[1+a_{21}g\ast u_1(t-\tau ,x)-a_{22}g\ast u_2(t,x)-a_{23}{u}_3(t-\tau ,x)],\hfill \\ \frac{\partial u_3(t,x)}{\partial t}=d_3\Delta u_3(t,x)+r_3{u}_3(t,x)[1-a_{31}g\ast u_1(t-\tau ,x)-a_{32}g\ast u_2(t-\tau ,x)-a_{33}{u}_3(t,x)],\hfill \end{array}\right.$$

(4)

for $x\in {\mathbb{R}}^n,t>0$, where $d_i>0,r_i>0,\tau \ge 0,a_{ij}\ge 0$ are constants and $u_i(t,x)$ denotes the density of species i at the location $x\in {\mathbb{R}}^n$ and time $t>0$, respectively $(i,j=1,2,3)$. In this system, the species $u_1$ and $u_2$ are cooperative with each other, and the species $u_3$ is competitive to both species $u_1$ and $u_2$. The nonlocal interaction term with time delay $g\ast u_i(t-\tau ,x)$ is defined as

$$g\ast u_i(t-\tau ,x)={\int }_{{\mathbb{R}}^n}g\left(y\right)u_i(t-\tau ,x-y)\mathrm{d}y,i=1,2,3.$$

The kernel function g is a continuous, even, and nonnegative function such that ${\int }_{{\mathbb{R}}^n}g\left(x\right)\mathrm{d}x=1$. The spatial convolution arises because species diffuse, meaning they were not located at the same spatial point x at the earlier time $t-\tau$. Consequently, interspecific competition or cooperation for resources depends on more than the population density at a single spatial point but a weighted average involving values at all points in space (see Gourley and Britton [35], Gourley et al. [36]). When the nonlocal interaction is very narrow, in the limiting case where g is the delta function centered at the origin (i.e., $g\left(x\right)=\delta \left(x\right)$), the system (4) can be reduced to the following time-delayed system

$$\left\{\begin{array}{c}\frac{\partial u_1(t,x)}{\partial t}=d_1\Delta u_1(t,x)+r_1{u}_1(t,x)[1-a_{11}{u}_1(t,x)+a_{12}{u}_2(t-\tau ,x)-a_{13}{u}_3(t-\tau ,x)],\hfill \\ \frac{\partial u_2(t,x)}{\partial t}=d_2\Delta u_2(t,x)+r_2{u}_2(t,x)[1+a_{21}{u}_1(t-\tau ,x)-a_{22}{u}_2(t,x)-a_{23}{u}_3(t-\tau ,x)],\hfill \\ \frac{\partial u_3(t,x)}{\partial t}=d_3\Delta u_3(t,x)+r_3{u}_3(t,x)[1-a_{31}{u}_1(t-\tau ,x)-a_{32}{u}_2(t-\tau ,x)-a_{33}{u}_3(t,x)].\hfill \end{array}\right.$$

(5)

We can refer to [49,52,53] and references in it to learn more details about the research of the competitive–cooperative system (5).

A planar traveling wave of system (5) takes the form of $u_i(t,x)={\varphi }_i(x·\nu +ct)$ ($i=1,2,3$, where $\nu \in {\mathbb{R}}^n$ is a fixed unit vector), which connects two equilibria of (5). In one-dimensional space, Hu et al. [54] investigated the stability of traveling waves for the Lotka–Volterra competition system with three species in $\mathbb{R}$. The methods they adopted were the weighted functional space, spectrum approach, and squeezing theorem. The existence and stability of traveling waves for system (5) have also been studied in [49,52,53]. However, to our knowledge, few results address the multidimensional stability of planar traveling waves in competitive–cooperative systems like (5) in higher-dimensional spaces. Multidimensional analysis allows for a deeper understanding of stability as those transverse perturbations perpendicular to the propagation direction can lead to phenomena like wave breaking, turbulence, or pattern formation, which one-dimensional models cannot capture. Furthermore, higher dimensions introduce unique geometric and topological properties, such as wavefront curvature and the effects of heterogeneous environments, which significantly influence wave propagation. This paper aims to investigate the multidimensional stability of planar traveling waves in system (5). We demonstrate that planar traveling waves with speed $c>c^{*}$ (where $c^{*}$ is defined in (15)) are exponentially stable in $L^{\infty }\left({\mathbb{R}}^n\right)$, with the decay rate expressed as $t^{-\frac{n}{2}}{\mathrm{e}}^{-{\epsilon }_{\tau }\sigma t}$ by the weighted energy method and Fourier transform, where $\sigma >0$ and $\epsilon \left(\tau \right)\in (0,1)$ is a decreasing function of the time delay $\tau >0$. Notably, time delay significantly slows the decay rate. For waves with $c=c^{*}$, we prove that planar traveling waves exhibit algebraic stability, with a decay rate of $t^{-\frac{n}{2}}$ in $L^{\infty }\left({\mathbb{R}}^n\right)$.

The structure of this paper is as follows. Section 2 introduces the necessary notations and discusses the existence and stability of planar traveling waves. Section 3 establishes the multidimensional stability of planar traveling waves, including the case where $c=c^{*}$. Finally, Section 4 provides numerical simulations to support the main findings on the basis of Appendix A, which presents exact planar traveling waves without time delay.

2. Preliminaries and Main Results

First, we elaborate on some necessary notations throughout this paper. $C>0$ denotes a generic constant and $C_i(i=0,1,2,\cdots )$ represents a specific constant. Let $∥·∥$ and ${∥·∥}_{\infty }$ denote the 1–norm and ∞–norm of the matrices (or vectors), respectively. Let $\Omega$ be a domain, typically $\Omega ={\mathbb{R}}^n$ and $L^p(\Omega )$ be the Lebesgue space of the integrable functions defined on $\Omega$. $\alpha =({\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_n)$ denotes a multi-index with nonnegative integers ${\alpha }_i\ge 0(i=1,2,\cdots ,n)$. The derivatives for a multi-dimensional function $f\left(x\right)$ are denoted as ${\partial }^{\alpha }f\left(x\right)={\partial }_{x_1}^{{\alpha }_1}\cdots {\partial }_{x_n}^{{\alpha }_n}f\left(x\right)$. $W^{k,p}(\Omega )(k\ge 1,p\ge 1)$ is the Sobolev space where the function $f\left(x\right)$ is defined on $\Omega$ and its weak derivatives ${\partial }^{\alpha }f\left(x\right)\left(\right|\alpha |\le k)$ also belong to $L^p(\Omega )$. We denote $W^{k,2}(\Omega )$ as $H^k(\Omega )$. Further, $L_w^p(\Omega )$ denotes the weighted $L^p$ space with a weighted function $w\left(x\right)>0$. Its norm is defined by

$${\parallel f\parallel }_{L_w^p(\Omega )}={\left({\int }_{\Omega }w\left(x\right){\left|f\left(x\right)\right|}^p\mathrm{d}x\right)}^{1/p}.$$

$W_w^{k,p}(\Omega )$ is the weighted Sobolev space with the norm

$${∥f∥}_{W_w^{k,p}(\Omega )}={\left(\sum _{\left|\alpha \right|\le k}{\int }_{\Omega }w\left(x\right){\left|{\partial }^{\alpha }f\left(x\right)\right|}^p\mathrm{d}x\right)}^{1/p}.$$

Fourier transform is defined as

$$\mathcal{F}\left[f\right]\left(\eta \right)=\widehat{f}\left(\eta \right):={\int }_{{\mathbb{R}}^n}{\mathrm{e}}^{-\mathrm{i}x·\eta }f\left(x\right)\mathrm{d}x.$$

The inverse Fourier transform is given by

$${\mathcal{F}}^{-1}\left[\widehat{f}\right]\left(x\right):=\frac{1}{{\left(2\pi \right)}^n}{\int }_{{\mathbb{R}}^n}{\mathrm{e}}^{\mathrm{i}x·\eta }\widehat{f}\left(\eta \right)\mathrm{d}\eta .$$

In what follows, we recall the solution formula and the decay rate of the linear differential equation with time delay, which makes a difference in proofs of decay rates in Section 3. Now let us consider the following delayed differential system,

$$\left\{\begin{array}{c}\frac{\mathrm{d}}{\mathrm{d}t}z\left(t\right)=Az\left(t\right)+Bz(t-\tau ),t>0,\hfill \\ z\left(s\right)=z_0\left(s\right),s\in [-\tau ,0],\hfill \end{array}\right.$$

(6)

where $A,B$ are constant $n\times n$ matrices, $z\left(t\right)\in {\mathbb{R}}^n$, and $\tau >0$ denotes a time delay. In [55], Khusainov and Ivanov obtained the solution formula of (6) in the case of $n=1$. In [56], Ma et al. presented the solution formula of system (6) and gave a sufficient condition of the global stability for the trivial solution of the linear delayed system (6) in the general case of $n\ge 2$ shown as follows.

Lemma 1 

([56]). If the initial data satisfy $z_0\left(s\right)\in C^1([-\tau ,0],{\mathbb{R}}^n)$, then the solution of system (6) can be shown as

$$z\left(t\right)={\mathrm{e}}^{A(t+\tau )}{\mathrm{e}}_{\tau }^{B_{1}t}{z}_0(-\tau )+{\int }_{-\tau }^0{\mathrm{e}}^{A(t-s)}{\mathrm{e}}_{\tau }^{B_1(t-s-\tau )}[z_0^{\prime }\left(s\right)-A{z}_0\left(s\right)]\mathrm{d}s,$$

where $B_1=B{\mathrm{e}}^{-A\tau }$ and ${\mathrm{e}}_{\tau }^{B_{1}t}$ represents the delayed exponential function, which takes the following form

$${\mathrm{e}}_{\tau }^{B_{1}t}=\left\{\begin{array}{cc}0,\hfill & -\infty <t<-\tau ,\\ 1,\hfill & -\tau \le t<0,\\ 1+B_{1}t,\hfill & 0\le t<\tau ,\\ 1+B_{1}t+\frac{B_1^2}{2!}{(t-\tau )}^2,\hfill & \tau \le t<2\tau ,\\ ⋮\hfill & ⋮\\ 1+B_{1}t+\frac{B_1^2}{2!}{(t-\tau )}^2+\cdots +\frac{B_1^m}{m!}{[t-(m-1)\tau ]}^m,\hfill & (m-1)\tau \le t<m\tau ,\\ ⋮\hfill & ⋮\end{array}\right.$$

(7)

According to Ma et al. [56], we can give a conclusion on the global stability for the trivial solution of the linear delayed system (6).

Lemma 2 

([56]). Suppose $\mu \left(A\right):=\frac{{\mu }_1\left(A\right)+{\mu }_{\infty }\left(A\right)}{2}<0$, where ${\mu }_1\left(A\right)$ and ${\mu }_{\infty }\left(A\right)$ denote the matrix measure of A induced by the matrix 1-norm $∥·∥$ and ∞-norm ${∥·∥}_{\infty }$, respectively. If $\nu \left(B\right):=\frac{∥B∥+{∥B∥}_{\infty }}{2}\le -\mu \left(A\right)$, then there exists a decreasing function ${\epsilon }_{\tau }=\epsilon \left(\tau \right)\in (0,1)$ for $\tau >0$ such that any solution of system (6) satisfies

$$∥z\left(t\right)∥\le C_0{\mathrm{e}}^{-{\epsilon }_{\tau }\sigma t},t>0,$$

(8)

where $C_0$ is a positive constant depending on initial data $z_0\left(s\right),s\in [-\tau ,0]$, and $\sigma =\left|\mu \right(A\left)\right|-\nu \left(B\right)$. In particular,

$$∥{\mathrm{e}}^{At}{\mathrm{e}}_{\tau }^{B_{1}t}∥\le C_0{\mathrm{e}}^{-{\epsilon }_{\tau }\sigma t},t>0.$$

(9)

where ${\mathrm{e}}_{\tau }^{B_{1}t}$ is defined in (7).

Next, we introduce the following scaling transformation for system (5),

$$\begin{array}{c}{\tilde{u}}_1=a_{11}{u}_1,{\tilde{u}}_2=a_{22}{u}_2,{\tilde{u}}_3=a_{33}{u}_3,\\ {\tilde{a}}_{12}=\frac{a_{12}}{a_{22}},{\tilde{a}}_{13}=\frac{a_{13}}{a_{33}},{\tilde{a}}_{21}=\frac{a_{21}}{a_{11}},{\tilde{a}}_{23}=\frac{a_{23}}{a_{33}},{\tilde{a}}_{31}=\frac{a_{31}}{a_{11}},{\tilde{a}}_{32}=\frac{a_{32}}{a_{22}}.\end{array}$$

(10)

Setting diffusion rates to $d_1=d_2=d_3=1$ and removing tildes for the sake of convenience, we note that system (5) can be changed into

$$\left\{\begin{array}{c}\frac{\partial u_1(t,x)}{\partial t}=\Delta u_1(t,x)+r_1{u}_1(t,x)[1-u_1(t,x)+a_{12}{u}_2(t-\tau ,x)-a_{13}{u}_3(t-\tau ,x)],\hfill \\ \frac{\partial u_2(t,x)}{\partial t}=\Delta u_2(t,x)+r_2{u}_2(t,x)[1+a_{21}{u}_1(t-\tau ,x)-u_2(t,x)-a_{23}{u}_3(t-\tau ,x)],\hfill \\ \frac{\partial u_3(t,x)}{\partial t}=\Delta u_3(t,x)+r_3{u}_3(t,x)[1-a_{31}{u}_1(t-\tau ,x)-a_{32}{u}_2(t-\tau ,x)-u_3(t,x)],\hfill \end{array}\right.$$

(11)

for $(t,x)\in {\mathbb{R}}_{+}\times {\mathbb{R}}^n$. Assume that $0\le a_{ij}<1(i,j=1,2,3,i\ne j)$. Then, the system (11) has the following equilibria:

$$\begin{array}{c}{E}_0=(0,0,0),E_1=(1,0,0),E_2=(0,1,0),E_3=(0,0,1),\hfill \\ E_4=\left(\frac{1+a_{12}}{1-a_{12}{a}_{21}},\frac{1+a_{21}}{1-a_{12}{a}_{21}},0\right):=(k_1,k_2,0),\hfill \\ E_5=\left(\frac{1-a_{13}}{1-a_{13}{a}_{31}},0,\frac{1-a_{31}}{1-a_{13}{a}_{31}}\right),\hfill \\ E_6=\left(0,\frac{1-a_{23}}{1-a_{23}{a}_{32}},\frac{1-a_{32}}{1-a_{23}{a}_{32}}\right).\hfill \end{array}$$

A planar traveling wave is denoted by $u_i(x,t)={\varphi }_i\left(\xi \right),$ ($\xi =\nu ·x+ct$, where $\nu \in {\mathbb{R}}^n$ is a fixed unit vector, here we set $\nu =e_1=(1,0,\cdots ,0)$ for simplicity) of Equation (13), then ${\varphi }_i$ has to satisfy the following system

$$\left\{\begin{array}{c}c{\varphi }_1^{\prime }\left(\xi \right)={\varphi }_1^{″}\left(\xi \right)+r_1{\varphi }_1\left(\xi \right)[1-{\varphi }_1\left(\xi \right)+a_{12}{\varphi }_2(\xi -c\tau )-a_{13}{\varphi }_3(\xi -c\tau )],\hfill \\ c{\varphi }_2^{\prime }\left(\xi \right)={\varphi }_2^{″}\left(\xi \right)+r_2{\varphi }_2\left(\xi \right)[1+a_{21}{\varphi }_1(\xi -c\tau )-{\varphi }_2\left(\xi \right)-a_{23}{\varphi }_3(\xi -c\tau )],\hfill \\ c{\varphi }_3^{\prime }\left(\xi \right)={\varphi }_3^{″}\left(\xi \right)+r_3{\varphi }_3\left(\xi \right)[1-a_{31}{\varphi }_1(\xi -c\tau )-a_{32}{\varphi }_2(\xi -c\tau )-{\varphi }_3\left(\xi \right)],\hfill \end{array}\right.\xi \in \mathbb{R}.$$

(12)

Huang and Weng ([52] Theorem 4.1) proved the existence of the solution of system (12) by using the upper–lower solutions and Schauder’s fixed point theorem.

Proposition 1. 

Let $m_i\in (0,k_i),i=1,2$, be small enough and satisfy

Hypothesis 1 (H1). 

$c_0:=a_{31}{k}_1+a_{32}{k}_2-1>0$;

Hypothesis 2 (H2). 

$1-a_{13}>m_1,1-a_{23}>m_2,a_{31}{m}_1+a_{32}{m}_2>max\left\{\frac{m_1}{k_1},\frac{m_2}{k_2}\right\}$.

Then, for $c>c_{*}:=max\left\{2\sqrt{r_1(1-a_{13})},2\sqrt{r_2(1-a_{23})},2\sqrt{r_3{c}_0}\right\}$, system (12) has a solution $\mathrm{\Phi }\left(\xi \right)=({\varphi }_1,{\varphi }_2,{\varphi }_3)\left(\xi \right)$ connecting $E_3$ and $E_4$, which satisfies

$${\varphi }_1^{\prime }\left(\xi \right)>0,{\varphi }_2^{\prime }\left(\xi \right)>0,{\varphi }_3^{\prime }\left(\xi \right)<0,\xi \in \mathbb{R}.$$

Remark 1. 

Here, we assume that the diffusion rates $d_1=d_2=d_3=1$ in system (5). It is easy to verify that the conditions in ([52] Theorem 4.1) for the existence of solution of system (12) can be reduced to the assumptions (H1) and (H2).

To utilize the comparison principle, we convert the competitive–cooperative system (11) into a cooperative system by introducing the variable transformation ${\tilde{u}}_1=u_1$, ${\tilde{u}}_2=u_2$, and ${\tilde{u}}_3=1-u_3$. For simplicity, we omit the tildes in the resulting system, which becomes as follows

$$\left\{\begin{array}{c}\frac{\partial u_1(t,x)}{\partial t}=\Delta u_1(t,x)+r_1{u}_1(t,x)[1-a_{13}-u_1(t,x)+a_{12}{u}_2(t-\tau ,x)+a_{13}{u}_3(t-\tau ,x)],\hfill \\ \frac{\partial u_2(t,x)}{\partial t}=\Delta u_2(t,x)+r_2{u}_2(t,x)[1-a_{23}+a_{21}{u}_1(t-\tau ,x)-u_2(t,x)+a_{23}{u}_3(t-\tau ,x)],\hfill \\ \frac{\partial u_3(t,x)}{\partial t}=\Delta u_3(t,x)+r_3(1-u_3(t,x))[a_{31}{u}_1(t-\tau ,x)+a_{32}{u}_2(t-\tau ,x)-u_3(t,x)]\hfill \end{array}\right.$$

(13)

with the initial conditions

$$u_i(s,x)=u_{i0}(s,x),(s,x)\in [-\tau ,0]\times {\mathbb{R}}^n,i=1,2,3.$$

(14)

According to the properties of the monotone semiflows [57], we have the following comparison principle.

Lemma 3 

(Comparison Principle). Let ${\mathit{u}}^{\pm }(t,x)=(u_1^{\pm },u_2^{\pm },u_3^{\pm })(t,x)$ be the solution of (13) with the initial data ${\mathit{u}}_0^{\pm }(s,x)=(u_1^{\pm },u_2^{\pm },u_3^{\pm })(s,x)$, respectively. If

$${\mathit{E}}_{-}:=(0,0,0)\le {\mathit{u}}_0^{-}(s,x)\le {\mathit{u}}_0^{+}(s,x)\le (k_1,k_2,1):={\mathit{E}}_{+},$$

for $(s,x)\in [-\tau ,0]\times {\mathbb{R}}^n$, then

$$(0,0,0)\le {\mathit{u}}^{-}(t,x)\le {\mathit{u}}^{+}(t,x)\le (k_1,k_2,1),$$

for $(t,x)\in {\mathbb{R}}_{+}\times {\mathbb{R}}^n$.

We obtain the following existence of planar traveling waves of system (13) by Proposition 1.

Theorem 1. 

Provided that conditions in Proposition 1 hold. For any $c>c_{*}$, there is a planar traveling wave $\mathrm{\Phi }\left(\xi \right)=({\varphi }_1,{\varphi }_2,{\varphi }_3)\left(\xi \right)$ of the system (13) which connects the equilibria $(0,0,0)$ and $(k_1,k_2,1)$, with the wave profile component ${\varphi }_i$ increasing.

To obtain the stability of planar traveling waves, we introduce the following assumption.

Hypothesis 3 (H3). 

$1-(a_{12}+a_{13})>0$, $1-(a_{21}+a_{23})>0$.

Remark 2. 

In the case of one-dimensional space, Ma et al. ([49] Theorem 2.2) obtained the stability of the traveling wave of system (5) without time delays by setting the following assumptions,

(P1) 

$1-(a_{12}+a_{13})>\frac{r_2{k}_2}{r_1{k}_1}{a}_{21}$, $1-(a_{21}+a_{23})>\frac{r_1{k}_1}{r_2{k}_2}{a}_{12}$,

(P2) 

$a_{31}{k}_1+a_{32}{k}_2-1>\frac{r_1}{r_3}{k}_1{a}_{13}+\frac{r_2}{r_3}{k}_2{a}_{23}$.

Here, we do not need the conditions of the ratios $\frac{r_i}{r_j}(i,j=1,2,3,i\ne j)$ in (P1)–(P2) and reduce the assumptions (P1) and (P2) to the assumptions (H3) and (H1), respectively. Thus, if we denote the range of parameters satisfying (H1)–(H3) by $\mathfrak{R}$, the set $\mathfrak{R}$ can contain a very large range of parameters corresponding to assumptions in ([49] Theorem 2.2). Indeed, the Hypotheses (H1)–(H3) can be ensured if we assume

$$\begin{array}{c}{a}_{12}=a_{21}=\frac{3}{5},a_{13}=a_{23}=\frac{1}{8},\hfill \\ a_{31}=a_{32}=\frac{4}{5},m_1=m_2\in (0,\frac{7}{8}).\hfill \end{array}$$

Set

$$c^{*}=2\sqrt{\nu \left(A\right)+\nu \left(B\right)},$$

(15)

where $\nu (·)=\frac{∥·∥+{∥·∥}_{\infty }}{2}$ and

$$A=\left(\begin{array}{ccc}{r}_1{k}_1& 0& 0\\ 0& r_2{k}_2& 0\\ 0& 0& r_3\end{array}\right),B=\left(\begin{array}{ccc}0& r_1{k}_1{a}_{12}& r_1{k}_1{a}_{13}\\ r_2{k}_2{a}_{21}& 0& r_2{k}_2{a}_{23}\\ r_3{a}_{31}& r_3{a}_{32}& 0\end{array}\right).$$

Let $c\ge c^{*}$ and $({\varphi }_1(x·e_1+ct),{\varphi }_2(x·e_1+ct),{\varphi }_3(x·e_1+ct))$ be the planar traveling wave of (13) with the speed c connecting ${\mathit{E}}_{-}$ and ${\mathit{E}}_{+}$. For such a speed c and large constant ${\zeta }_0$, we define a weighted function as

$$w\left(x\right)=\left\{\begin{array}{cc}{\mathrm{e}}^{-\lambda (x_1-{\zeta }_0)},\hfill & x_1\le {\zeta }_0,\hfill \\ 1,\hfill & x_1\ge {\zeta }_0,\hfill \end{array}\right.$$

(16)

where $\lambda =\frac{1}{2}c$.

Here, we present the main results in this paper.

Theorem 2 

(Stability). Supposed that (H1)–(H3) hold. For any given planar traveling wave $({\varphi }_1(x·e_1+ct),{\varphi }_2(x·e_1+ct),{\varphi }_3(x·e_1+ct))$ of the system (13) with speed $c\ge c^{*}$ connecting ${\mathit{E}}_{-}$ and ${\mathit{E}}_{+}$, if the initial data satisfy

$$0\le u_{i0}(s,x)\le k_i,(s,x)\in [-\tau ,0]\times {\mathbb{R}}^n,i=1,2,3,k_3=1,$$

(17)

and the initial perturbation

$$u_{i0}-{\varphi }_i\in C^1([-\tau ,0],W_w^{2,1}\left({\mathbb{R}}^n\right)\cap L^2\left({\mathbb{R}}^n\right)),i=1,2,3,$$

(18)

then a nonnegative solution of the Cauchy problem (13) and (14) uniquely exists and satisfies

$$(0,0,0)\le (u_1(t,x),u_2(t,x),u_3(t,x))\le (k_1,k_2,1),x\in {\mathbb{R}}^n,t>0,$$

and

$$u_i-{\varphi }_i\in C^1([0,+\infty ),W_w^{2,1}\left({\mathbb{R}}^n\right)\cap L^2\left({\mathbb{R}}^n\right)).$$

Furthermore, we obtain

(i) 

when $c>c^{*}$, the solution $(u_1(t,x),u_2(t,x),u_3(t,x))$ converges to the planar traveling wave $({\varphi }_1(x·e_1+ct),{\varphi }_2(x·e_1+ct),{\varphi }_3(x·e_1+ct))$ exponentially in time, i.e.,

$$\underset{x\in {\mathbb{R}}^n}{sup}|u_i(t,x)-{\varphi }_i(x·e_1+ct)|\le C{t}^{-\frac{n}{2}}{\mathrm{e}}^{-{\epsilon }_{\tau }\sigma t},t>0,i=1,2,3,$$

for some constant $\sigma >0$, where ${\epsilon }_{\tau }=\epsilon \left(\tau \right)\in (0,1)$ is a decreasing function for $\tau >0$.

(ii) 

When $c=c^{*}$, the solution $(u_1(t,x),u_2(t,x),u_3(t,x))$ converges to the planar traveling wave $({\varphi }_1(x·e_1+c^{*}t),{\varphi }_2(x·e_1+c^{*}t),{\varphi }_3(x·e_1+c^{*}t))$ algebraically in time, i.e.,

$$\underset{x\in {\mathbb{R}}^n}{sup}|u_i(t,x)-{\varphi }_i(x·e_1+ct)|\le C{t}^{-\frac{n}{2}},t>0,i=1,2,3.$$

Remark 3. 

We refer to the method of [25,49] to prove the Theorem 2, but there is a significant difference. Comparing the weighted function $w\left(x\right)$ in Theorem 2 with Theorem 1.2 in [25], it can be observed that the weight function $w\left(x\right)$ in Theorem 2 is unrelated to the eigenvalue, which serves as the root of the characteristic equation corresponding to the traveling wavefronts. In one-dimensional space, Ma et al. [49] only proved the stability of traveling waves with speed $c>c^{*}$ under the decay rate ${\mathrm{e}}^{-\mu t}$. In this work, we not only provide a more precise decay rate $(t^{-\frac{n}{2}}{\mathrm{e}}^{-{\epsilon }_{\tau }\sigma t})$ for the stability of planar traveling waves with speed $c>c^{*}$ in high dimensional space, but also establish the algebraic stability of planar traveling waves with speed $c=c^{*}$.

3. Stability

In this section, we focus on the proof of stability of planar traveling waves. To begin with, we present the global existence and uniqueness of the solution to the Cauchy problem (13) and (14), which can be demonstrated by the standard energy method and continuity extension method [58] or the theory of abstract functional differential equation [59].

Proposition 2 

(Global Existence and Uniqueness). Assume that the initial data satisfy

$$0\le u_{i0}(s,x)\le k_i,x\in {\mathbb{R}}^n,i=1,2,3.$$

for any given planar traveling wave $({\varphi }_1(x·e_1+ct),{\varphi }_2(x·e_1+ct),{\varphi }_3(x·e_1+ct))$ of (13) with speed $c\ge c^{*}$ connecting the equilibrium ${\mathit{E}}_{-}$ and ${\mathit{E}}_{+}$. If the initial perturbation satisfies

$$u_{i0}-{\varphi }_i\in C^1([-\tau ,0],W_w^{2,1}\left({\mathbb{R}}^n\right)\cap L^2\left({\mathbb{R}}^n\right)),i=1,2,3,$$

then there exists a unique solution $(u_1(t,x),u_2(t,x),u_3(t,x))$ of Cauchy problem (13) globally and (14) such that

$$(0,0,0)\le (u_1(t,x),u_2(t,x),u_3(t,x))\le (k_1,k_2,1),x\in {\mathbb{R}}^n,t>0,$$

and

$$u_i-{\varphi }_i\in C^1([0,+\infty ),W_w^{2,1}\left({\mathbb{R}}^n\right)\cap L^2\left({\mathbb{R}}^n\right)).$$

where the function $w\left(x\right)$ is defined by (16).

Let $c\ge c^{*}$ and initial data $u_{i0}(s,x)$ satisfy $0\le u_{i0}\left(x\right)\le k_i(k_3=1)$, and define

$$\begin{array}{c}{\tilde{u}}_{i0}^{+}(s,x)=max\{u_{i0}(s,x),{\varphi }_i(x·e_1+cs)\},\\ {\tilde{u}}_{i0}^{-}(s,x)=min\{u_{i0}(s,x),{\varphi }_i(x·e_1+cs)\},\end{array}$$

where $x\in {\mathbb{R}}^n,i=1,2,3$. It is apparent that

$$\begin{array}{c}0\le {\tilde{u}}_{i0}^{-}(s,x)\le u_{i0}(s,x)\le {\tilde{u}}_{i0}^{+}(s,x)\le k_i,\\ 0\le {\tilde{u}}_{i0}^{-}(s,x)\le {\varphi }_i(x·e_1+cs)\le {\tilde{u}}_{i0}^{+}(s,x)\le k_i.\end{array}$$

(19)

The piecewise continuity of the initial value ${\tilde{u}}_{i0}^{\pm }(s,x)$ can be guaranteed. However, the initial value exhibits low regularity, potentially resulting in a lack of regularity in the corresponding solutions. To address this limitation, we replace these initial data with smooth functions $u_{i0}^{\pm }(s,x)$ as the new initial data, ensuring that

$$0\le u_{i0}^{-}(s,x)\le {\tilde{u}}_{i0}^{-}(s,x)\le u_{i0}(s,x)\le {\tilde{u}}_{i0}^{+}(s,x)\le u_{i0}^{+}(s,x)\le k_i,$$

(20)

for $s\in [-\tau ,0],x\in {\mathbb{R}}^n,i=1,2,3$.

Let $(u_1^{\pm }(t,x),u_2^{\pm }(t,x),u_3^{\pm }(t,x))$ be the corresponding solutions of (13) with the initial data $(u_{10}^{\pm }(s,x),u_{20}^{\pm }(s,x),u_{30}^{\pm }(s,x))$. Thus, according to the comparison principle in Lemma 3, it follows that

$$\begin{array}{c}0\le u_i^{-}(t,x)\le u_i(t,x)\le u_i^{+}(t,x)\le k_i,\hfill \\ 0\le u_i^{-}(t,x)\le {\varphi }_i(x·e_1+ct)\le u_i^{+}(t,x)\le k_i,\hfill \end{array}$$

(21)

for $x\in {\mathbb{R}}^n,i=1,2,3$. Note that

$$\begin{array}{ccc}\hfill U_i(t,\xi )& =& u_i^{+}(t,x)-{\varphi }_i(x·e_1+ct),t>0,\hfill \\ \hfill U_i(s,\xi )& =& u_i^{+}(s,x)-{\varphi }_i(x·e_1+cs),s\in [-\tau ,0],\hfill \end{array}$$

for $x\in {\mathbb{R}}^n,t>0,i=1,2,3$, and $\xi =x+ct·e_1=(x_1+ct,x_2,\cdots ,x_n)(c\ge c^{*})$, it follows from (19) and (21) that

$$(0,0,0)\le (U_1(s,\xi ),U_2(s,\xi ),U_3(s,\xi ))\le (k_1,k_2,k_3),s\in [-\tau ,0],\xi \in {\mathbb{R}}^n,$$

(22)

and

$$(0,0,0)\le (U_1(t,\xi ),U_2(t,\xi ),U_3(t,\xi ))\le (k_1,k_2,k_3),t>0,\xi \in {\mathbb{R}}^n.$$

(23)

Then we obtain that $(U_1(t,\xi ),U_2(t,\xi ),U_3(t,\xi ))$ satisfies

$$\left\{\begin{array}{cc}{U}_{1t}+c{U}_{1{\xi }_1}=\hfill & \Delta U_1+r_1\left[1-a_{13}-(U_1+{\varphi }_1)+a_{12}(U_2^{\tau }+{\varphi }_2^{\tau })+a_{13}(U_3^{\tau }+{\varphi }_3^{\tau })\right]U_1\hfill \\ \hfill & +r_1{\varphi }_1\left[a_{12}{U}_2^{\tau }+a_{13}{U}_3^{\tau }-U_1\right],\hfill \\ U_{2t}+c{U}_{2{\xi }_1}=\hfill & \Delta U_2+r_2\left[1-a_{23}+a_{21}(U_1^{\tau }+{\varphi }_1^{\tau })-(U_2+{\varphi }_2)+a_{23}(U_3^{\tau }+{\varphi }_3^{\tau })\right]U_2\hfill \\ \hfill & +r_2{\varphi }_2\left[a_{21}{U}_1^{\tau }+a_{23}{U}_3^{\tau }-U_2\right],\hfill \\ U_{3t}+c{U}_{3{\xi }_1}=\hfill & \Delta U_3+r_3\left[({\varphi }_3+U_3)-a_{31}(U_1^{\tau }+{\varphi }_1^{\tau })-a_{32}(U_2^{\tau }+{\varphi }_2^{\tau })\right]U_3\hfill \\ \hfill & +r_3[1-{\varphi }_3][a_{31}{U}_1^{\tau }+a_{32}{U}_2^{\tau }-U_3],\hfill \end{array}\right.$$

(24)

where $U_i=U_i(t,\xi ),U_i^{\tau }=U_i(t-\tau ,\xi -c\tau ·e_1),{\varphi }_i={\varphi }_i(\xi ·e_1)$ and ${\varphi }_i^{\tau }={\varphi }_i(\xi ·e_1-c\tau )$ for $t>0,\xi \in {\mathbb{R}}^n,i=1,2,3$.

On the basis of conditions in Theorem 2, the following lemma indicates the decay rate of $U_i$ in ${\Omega }_{-}=(-\infty ,{\zeta }_0]\times {\mathbb{R}}^{n-1}$.

Lemma 4. 

There exists a decreasing function ${\epsilon }_{\tau }=\epsilon \left(\tau \right)\in (0,1)$ such that

$$\sum _{i=1}^3{\parallel U_i(t,·)\parallel }_{L^{\infty }\left({\Omega }_{-}\right)}\le C{t}^{-\frac{n}{2}}{\mathrm{e}}^{-{\epsilon }_{\tau }{\sigma }_{1}t},t>0,c>c^{*},$$

(25)

and

$$\sum _{i=1}^3{\parallel U_i(t,·)\parallel }_{L^{\infty }\left({\Omega }_{-}\right)}\le C{t}^{-\frac{n}{2}},t>0,c=c^{*},$$

(26)

where ${\sigma }_1={\sigma }_1\left(c\right)$ is defined in (35).

Proof. 

Owing to the nonnegative of $U_i$ and $U_i+{\varphi }_i=u_i^{+}\in (0,k_i),i=1,2,3$, we have

$$\begin{array}{cc}\hfill & r_1\left[1-a_{13}-(U_1+{\varphi }_1)+a_{12}(U_2^{\tau }+{\varphi }_2^{\tau })+a_{13}(U_3^{\tau }+{\varphi }_3^{\tau })\right]U_1+r_1{\varphi }_1\left[a_{12}{U}_2^{\tau }+a_{13}{U}_3^{\tau }-U_1\right]\hfill \\ \hfill & \le r_1(1+a_{12}{k}_2)U_1+r_1{k}_1\left[a_{12}{U}_2^{\tau }+a_{13}{U}_3^{\tau }\right]\hfill \\ \hfill & =r_1{k}_1\left[U_1+a_{12}{U}_2^{\tau }+a_{13}{U}_3^{\tau }\right],\hfill \end{array}$$

where we use the fact $1+a_{12}{k}_2=k_1$. Similarly, we have

$$\begin{array}{cc}\hfill & r_2\left[1-a_{23}+a_{21}(U_1^{\tau }+{\varphi }_1^{\tau })-(U_2+{\varphi }_2)+a_{23}(U_3^{\tau }+{\varphi }_3^{\tau })\right]U_2+r_2{\varphi }_2\left[a_{21}{U}_1^{\tau }+a_{23}{U}_3^{\tau }-U_2\right]\hfill \\ \hfill & \le r_2{k}_2[a_{21}{U}_1^{\tau }+U_2+a_{23}{U}_3^{\tau }],\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill & r_3\left[({\varphi }_3+U_3)-a_{31}(U_1^{\tau }+{\varphi }_1^{\tau })-a_{32}(U_2^{\tau }+{\varphi }_3^{\tau })\right]U_3+r_3[1-{\varphi }_3][a_{31}{U}_1^{\tau }+a_{32}{U}_2^{\tau }-U_3]\hfill \\ \hfill & \le r_3[a_{31}{U}_1^{\tau }+a_{32}{U}_2^{\tau }+U_3].\hfill \end{array}$$

Thus,

$$\left\{\begin{array}{c}{U}_{1t}+c{U}_{1{\xi }_1}\le \Delta U_1+r_1{k}_1\left[U_1+a_{12}{U}_2^{\tau }+a_{13}{U}_3^{\tau }\right],\hfill \\ U_{2t}+c{U}_{2{\xi }_1}\le \Delta U_2+r_2{k}_2[a_{21}{U}_1^{\tau }+U_2+a_{23}{U}_3^{\tau }],\hfill \\ U_{3t}+c{U}_{3{\xi }_1}\le \Delta U_3+r_3[a_{31}{U}_1^{\tau }+a_{32}{U}_2^{\tau }+U_3],\hfill \end{array}\right.$$

for $(t,\xi )\in {\mathbb{R}}_{+}\times {\mathbb{R}}^n$.

Let $v(t,\xi )={(v_1(t,\xi ),v_2(t,\xi ),v_3(t,\xi ))}^T$ be the solution of the following system with the same initial data $U_{i0}(s,\xi )$,

$$\left\{\begin{array}{c}{v}_{1t}(t,\xi )+c{v}_{1{\xi }_1}=\Delta v_1(t,\xi )+r_1{k}_1[v_1(t,\xi )+a_{12}{v}_2(t-\tau ,\xi -c\tau e_1)+a_{13}{v}_2(t-\tau ,\xi -c\tau e_1)],\hfill \\ v_{2t}(t,\xi )+c{v}_{2{\xi }_1}=\Delta v_2(t,\xi )+r_2{k}_2[a_{21}{v}_1(t-\tau ,\xi -c\tau e_1)+v_2(t,\xi )+a_{23}{v}_3(t-\tau ,\xi -c\tau e_1)],\hfill \\ v_{3t}(t,\xi )+c{v}_{3{\xi }_1}=\Delta v_3(t,\xi )+r_3[a_{31}{v}_1(t-\tau ,\xi -c\tau e_1)+a_{32}{v}_2(t-\tau ,\xi -c\tau e_1)+v_3(t,\xi )],\hfill \end{array}\right.$$

(27)

for $(t,\xi )\in {\mathbb{R}}_{+}\times {\mathbb{R}}^n$. Then, one obtains

$$(U_1(t,\xi ),U_2(t,\xi ),U_3(t,\xi ))\le (v_1(t,\xi ),v_2(t,\xi ),v_3(t,\xi )),(t,\xi )\in {\mathbb{R}}_{+}\times {\mathbb{R}}^n.$$

(28)

Letting $V_i(t,\xi )={\mathrm{e}}^{-\lambda ({\xi }_1-{\zeta }_0)}{v}_i(t,\xi )$. Then, we have

$$\left\{\begin{array}{c}{V}_{1t}+(c-2\lambda )V_{1{\xi }_1}=\Delta V_1+({\lambda }^2-c\lambda +r_1{k}_1)V_1+r_1{k}_1{\mathrm{e}}^{-c\lambda \tau }[a_{12}{V}_2^{\tau }+a_{13}{V}_3^{\tau }],\hfill \\ V_{2t}+(c-2\lambda )V_{2{\xi }_1}=\Delta V_2+({\lambda }^2-c\lambda +r_2{k}_2)V_2+r_2{k}_2{\mathrm{e}}^{-c\lambda \tau }[a_{21}{V}_1^{\tau }+a_{23}{V}_3^{\tau }],\hfill \\ V_{3t}+(c-2\lambda )V_{3{\xi }_1}=\Delta V_3+({\lambda }^2-c\lambda +r_3)V_3+r_3{\mathrm{e}}^{-c\lambda \tau }[a_{31}{V}_1^{\tau }+a_{32}{V}_2^{\tau }],\hfill \end{array}\right.$$

where $V_i=V_i(t,\xi ),V_i^{\tau }=V_i(t-\tau ,\xi -c\tau e_1)$ for $(t,\xi )\in {\mathbb{R}}_{+}\times {\mathbb{R}}^n,i=1,2,3$. Taking the Fourier transform to the above system, one has

$$\left\{\begin{array}{c}{\widehat{V}}_{1t}(t,\eta )={(-|\eta |}^2-\frac{1}{4}{c}^2+r_1{k}_1){\widehat{V}}_1(t,\eta )+r_1{k}_1{\mathrm{e}}^{-c\tau (\lambda +\mathrm{i}{\eta }_1)}[a_{12}{\widehat{V}}_2(t-\tau ,\eta )+a_{13}{\widehat{V}}_3(t-\tau ,\eta ),\hfill \\ {\widehat{V}}_{2t}(t,\eta )={(-|\eta |}^2-\frac{1}{4}{c}^2+r_2{k}_2){\widehat{V}}_2(t,\eta )+r_2{k}_2{\mathrm{e}}^{-c\tau (\lambda +\mathrm{i}{\eta }_1)}[a_{21}{\widehat{V}}_1(t-\tau ,\eta )+a_{23}{\widehat{V}}_3(t-\tau ,\eta )],\hfill \\ {\widehat{V}}_{3t}(t,\eta )={(-|\eta |}^2-\frac{1}{4}{c}^2+r_3){\widehat{V}}_3(t,\eta )+r_3{\mathrm{e}}^{-c\tau (\lambda +\mathrm{i}{\eta }_1)}[a_{31}{\widehat{V}}_1(t-\tau ,\eta )+a_{32}{\widehat{V}}_3(t-\tau ,\eta )],\hfill \end{array}\right.$$

(29)

where

$${\widehat{V}}_i(t,\eta )={\int }_{{\mathbb{R}}^n}{\mathrm{e}}^{-\mathrm{i}\xi ·\eta }{V}_i(t,\xi )\mathrm{d}\xi ,i=1,2,3.$$

Denote

$$A=\left(\begin{array}{ccc}{r}_1{k}_1& 0& 0\\ 0& r_2{k}_2& 0\\ 0& 0& r_3\end{array}\right),B_{\tau }\left(\eta \right)={\mathrm{e}}^{-c\tau (\lambda +\mathrm{i}{\eta }_1)}\left(\begin{array}{ccc}0& r_1{k}_1{a}_{12}& r_1{k}_1{a}_{13}\\ r_2{k}_2{a}_{21}& 0& r_2{k}_2{a}_{23}\\ r_3{a}_{31}& r_3{a}_{32}& 0\end{array}\right).$$

and $A\left(\eta \right)=-(|\eta {|}^2+\frac{1}{4}{c}^2)I+A$. Then, system (29) can be rewritten as

$${\widehat{V}}_t(t,\eta )=A\left(\eta \right)\widehat{V}(t,\eta )+B_{\tau }\left(\eta \right)\widehat{V}(t-\tau ,\eta ),$$

(30)

where $\widehat{V}(t,\eta )={({\widehat{V}}_1(t,\eta ),{\widehat{V}}_2(t,\eta ),{\widehat{V}}_3(t,\eta ))}^T$.

It follows form Lemma 1 that the solution of (30) has the following form

$$\begin{array}{cc}\hfill \widehat{V}(t,\eta )=& {\mathrm{e}}^{A\left(\eta \right)(t+\tau )}{\mathrm{e}}_{\tau }^{B_1\left(\eta \right)t}{\widehat{V}}_0(-\tau ,\eta )+{\int }_{-\tau }^0{\mathrm{e}}^{A\left(\eta \right)(t-s)}{\mathrm{e}}_{\tau }^{B_1\left(\eta \right)(t-s-\tau )}\left[\frac{\mathrm{d}{\widehat{V}}_0(s,\eta )}{\mathrm{d}s}-A\left(\eta \right){\widehat{V}}_0(s,\eta )\right]\mathrm{d}s,\hfill \\ \hfill :=& I_1(t,\eta )+{\int }_{-\tau }^0{I}_2(t-s,\eta )\mathrm{d}s,\hfill \end{array}$$

(31)

where $B_1\left(\eta \right)=B_{\tau }\left(\eta \right){\mathrm{e}}^{A\left(\eta \right)\tau }$. Therefore, by taking the inverse Fourier transform to (31), we obtain

$$\begin{array}{cc}\hfill V(t,\xi )=& {\mathcal{F}}^{-1}\left[I_1\right](t,\xi )+{\int }_{-\tau }^0{\mathcal{F}}^{-1}\left[I_2\right](t-s,\xi )\mathrm{d}s\hfill \\ \hfill =& \frac{1}{{\left(2\pi \right)}^n}{\int }_{{\mathbb{R}}^n}{\mathrm{e}}^{\mathrm{i}\xi ·\eta }{\mathrm{e}}^{A\left(\eta \right)(t+\tau )}{\mathrm{e}}_{\tau }^{B_1\left(\eta \right)t}{\widehat{V}}_0(-\tau ,\eta )\mathrm{d}\eta \hfill \\ \hfill & +{\int }_{-\tau }^0\frac{1}{{\left(2\pi \right)}^n}{\int }_{{\mathbb{R}}^n}{\mathrm{e}}^{\mathrm{i}\xi ·\eta }{\mathrm{e}}^{A\left(\eta \right)(t-s)}{\mathrm{e}}_{\tau }^{B_1\left(\eta \right)(t-s-\tau )}\left[\frac{\mathrm{d}{\widehat{V}}_0(s,\eta )}{\mathrm{d}s}-A\left(\eta \right){\widehat{V}}_0(s,\eta )\right]\mathrm{d}\eta \mathrm{d}s.\hfill \end{array}$$

(32)

For $c\ge c^{*}=2\sqrt{\nu \left(A\right)+\nu \left(B\right)}\ge 2\sqrt{\nu \left(A\right)+\nu \left(B_{\tau }\left(\eta \right)\right)}$, it is easy to verify that $\mu \left(A\right(\eta \left)\right)<0$ and

$$\nu \left(B_{\tau }\left(\eta \right)\right)\le -\mu \left(A\left(\eta \right)\right)={\left|\eta \right|}^2+\frac{1}{4}{c}^2-\nu \left(A\right),c\ge c^{*}.$$

(33)

By Lemma 2 and $\nu \left(B_{\tau }\left(\eta \right)\right)\le \nu \left(B\right)$, there exists a decreasing function ${\epsilon }_{\tau }=\epsilon \left(\tau \right)\in (0,1)$ such that

$$∥{\mathrm{e}}^{A\left(\eta \right)t}{\mathrm{e}}_{\tau }^{B_1\left(\eta \right)t}∥\le C_1{\mathrm{e}}^{-{\epsilon }_{\tau }\left(\right|\mu \left(A\left(\eta \right)\right)|-\nu \left(B_{\tau }\left(\eta \right)\right))t}\le C_1{\mathrm{e}}^{-{\epsilon }_{\tau }\left(\right|\mu \left(A\left(\eta \right)\right)|-\nu \left(B\right))t},$$

where $C_1$ is a positive constant depending on ${\widehat{V}}_0\left(s\right),s\in [-\tau ,0]$.

By the definition of the Fourier transform, we obtain

$$\underset{\eta \in {\mathbb{R}}^n}{sup}∥{\widehat{V}}_0(-\tau ,\eta )∥\le {\int }_{{\mathbb{R}}^n}∥V_0(-\tau ,x)∥\mathrm{d}x=\sum _{i=1}^3{∥V_{i0}(-\tau ,·)∥}_{L^1\left({\mathbb{R}}^n\right)}.$$

Thus,

$$\begin{array}{cc}\hfill \underset{\xi \in {\mathbb{R}}^n}{sup}∥{\mathcal{F}}^{-1}\left[I_1\right](t,\xi )∥=& \underset{\xi \in {\mathbb{R}}^n}{sup}∥\frac{1}{{\left(2\pi \right)}^n}{\int }_{{\mathbb{R}}^n}{\mathrm{e}}^{\mathrm{i}\xi ·\eta }{\mathrm{e}}^{A\left(\eta \right)(t+\tau )}{\mathrm{e}}_{\tau }^{B_1\left(\eta \right)t}{\widehat{V}}_0(-\tau ,\eta )\mathrm{d}\eta ∥\hfill \\ \hfill \le & C{\int }_{{\mathbb{R}}^n}{\mathrm{e}}^{-{\epsilon }_{\tau }\left(\right|\mu \left(A\left(\eta \right)\right)|-\nu \left(B\right))t}∥{\widehat{V}}_0(-\tau ,\eta )∥\mathrm{d}\eta \hfill \\ \hfill \le & C\underset{\eta \in {\mathbb{R}}^n}{sup}∥{\widehat{V}}_0(-\tau ,\eta )∥{\int }_{{\mathbb{R}}^n}|{\mathrm{e}}^{-{\epsilon }_{\tau }{\left(\right|\eta |}^2+\frac{1}{4}{c}^2-\left|\mu \left(A\right)\right|-\nu \left(B\right))t}|\mathrm{d}\eta \hfill \\ \hfill \le & C\sum _{i=1}^3{∥V_{i0}(-\tau ,·)∥}_{L^1\left({\mathbb{R}}^n\right)}{\mathrm{e}}^{-{\epsilon }_{\tau }{\sigma }_1\left(c\right)t}{\int }_{{\mathbb{R}}^n}{\mathrm{e}}^{-{\epsilon }_{\tau }{\left|\eta \right|}^{2}t}\mathrm{d}\eta \hfill \\ \hfill \le & C\sum _{i=1}^3{∥V_{i0}(-\tau ,·)∥}_{L^1\left({\mathbb{R}}^n\right)}{t}^{-\frac{n}{2}}{\mathrm{e}}^{-{\epsilon }_{\tau }{\sigma }_1\left(c\right)t},\hfill \end{array}$$

(34)

where we make use of the following derivation

$$\begin{array}{cc}\hfill {\int }_{{\mathbb{R}}^n}{\mathrm{e}}^{-{\epsilon }_{\tau }{\left|\eta \right|}^{2}t}\mathrm{d}\eta =& \prod _{i=1}^n{\int }_{\mathbb{R}}{\mathrm{e}}^{-{\epsilon }_{\tau }{\eta }_i^{2}t}\mathrm{d}{\eta }_i=\prod _{i=1}^n{\int }_{\mathbb{R}}{t}^{-\frac{1}{2}}{\mathrm{e}}^{-{\epsilon }_{\tau }{\zeta }_i^2}\mathrm{d}{\zeta }_i\le C{t}^{-\frac{n}{2}},t>0,\hfill \end{array}$$

and

$${\sigma }_1\left(c\right):=\frac{1}{4}{c}^2-\left(\right|\nu \left(A\right)|+\nu \left(B\right))\left\{\begin{array}{cc}=0,\hfill & c=c^{*},\hfill \\ >0,\hfill & c>c^{*}.\hfill \end{array}\right.$$

(35)

By using the property of the Fourier transform, we have

$$\left|{\left(\mathrm{i}{\eta }_i\right)}^{{\alpha }_i}{\widehat{V}}_j(t,\eta )\right|=\left|\mathcal{F}\left[{\partial }_{x_i}^{{\alpha }_i}{V}_j\right](t,\eta )\right|\le {\int }_{{\mathbb{R}}^n}\left|{\partial }_{x_i}^{{\alpha }_i}{V}_j(t,x)\right|\mathrm{d}x={∥{\partial }_{x_i}^{{\alpha }_i}{V}_j(t,·)∥}_{L^1\left({\mathbb{R}}^n\right)},$$

for $j=1,2,3$ and integers ${\alpha }_i\ge 0,i=1,2,\cdots n$. Then,

$$\underset{\eta \in {\mathbb{R}}^n}{sup}∥A\left(\eta \right){\widehat{V}}_0(s,\eta )∥=\underset{\eta \in {\mathbb{R}}^n}{sup}∥\left[-\left({\left|\eta \right|}^2+\frac{1}{4}{c}^2\right)I+A\right]{\widehat{V}}_0(s,\eta )∥\le C\sum _{i=1}^3{∥V_{i0}(s,·)∥}_{W^{2,1}\left({\mathbb{R}}^n\right)}.$$

Similarly, one has

$$\begin{array}{cc}\hfill & \underset{\xi \in {\mathbb{R}}^n}{sup}∥{\mathcal{F}}^{-1}\left[I_2\right](t-s,\xi )∥\hfill \\ \hfill =& \underset{\xi \in {\mathbb{R}}^n}{sup}∥\frac{1}{{\left(2\pi \right)}^n}{\int }_{{\mathbb{R}}^n}{\mathrm{e}}^{\mathrm{i}\xi ·\eta }{\mathrm{e}}^{A\left(\eta \right)(t-s)}{\mathrm{e}}_{\tau }^{B_1\left(\eta \right)(t-s-\tau )}\left[\frac{\mathrm{d}{\widehat{V}}_0(s,\eta )}{\mathrm{d}s}-A\left(\eta \right){\widehat{V}}_0(s,\eta )\right]\mathrm{d}\eta ∥\hfill \\ \hfill \le & C{\int }_{{\mathbb{R}}^n}{\mathrm{e}}^{-{\epsilon }_{\tau }\left(\right|\mu \left(A\left(\eta \right)\right)|-\nu \left(B\right))(t-s)}∥\frac{\mathrm{d}{\widehat{V}}_0(s,\eta )}{\mathrm{d}s}-A\left(\eta \right){\widehat{V}}_0(s,\eta )∥\mathrm{d}\eta \hfill \\ \hfill \le & C{\mathrm{e}}^{-{\epsilon }_{\tau }{\sigma }_1\left(c\right)(t-s)}{\int }_{{\mathbb{R}}^n}{\mathrm{e}}^{-{\epsilon }_{\tau }{\left|\eta \right|}^2(t-s)}∥\frac{\mathrm{d}{\widehat{V}}_0(s,\eta )}{\mathrm{d}s}-A\left(\eta \right){\widehat{V}}_0(s,\eta )∥\mathrm{d}\eta \hfill \\ \hfill \le & C{\mathrm{e}}^{-{\epsilon }_{\tau }{\sigma }_1\left(c\right)(t-s)}\underset{\eta \in {\mathbb{R}}^n}{sup}∥\frac{\mathrm{d}{\widehat{V}}_0(s,\eta )}{\mathrm{d}s}-A\left(\eta \right){\widehat{V}}_0(s,\eta )∥{\int }_{{\mathbb{R}}^n}{\mathrm{e}}^{-{\epsilon }_{\tau }{\left|\eta \right|}^2(t-s)}\mathrm{d}\eta \hfill \\ \hfill \le & C{(t-s)}^{-\frac{n}{2}}{\mathrm{e}}^{-{\epsilon }_{\tau }{\sigma }_1\left(c\right)(t-s)}\sum _{i=1}^3\left({∥\frac{\mathrm{d}{V}_{i0}(s,·)}{\mathrm{d}s}∥}_{L^1\left({\mathbb{R}}^n\right)}+{∥V_{i0}(s,·)∥}_{W^{2,1}\left({\mathbb{R}}^n\right)}\right).\hfill \end{array}$$

(36)

Substituting (34) and (36) to (32), we can obtain

$$\sum _{i=1}^3{∥V_i(t,·)∥}_{L^{\infty }\left({\mathbb{R}}^n\right)}\le C{t}^{-\frac{n}{2}}{\mathrm{e}}^{-\epsilon {\sigma }_1\left(c\right)t}.$$

(37)

Because of ${\mathrm{e}}^{\lambda ({\xi }_1-{\zeta }_0)}\le 1$ for ${\xi }_1\le {\zeta }_0$, we have

$$0\le U_i(t,\xi )\le v_i(t,\xi )={\mathrm{e}}^{\lambda ({\xi }_1-{\zeta }_0)}{V}_i(t,\xi )\le V_i(t,\xi ),i=1,2,3,$$

(38)

for $t>0,\xi \in {\Omega }_{-}\times {\mathbb{R}}^{n-1}$. Thus, (25) and (26) can be immediately deduced by (37) and (38). □

In what follows, we will derive the decay rate of $U_i$ in ${\Omega }_{+}=[{\zeta }_0,+\infty )\times {\mathbb{R}}^{n-1}$.

Lemma 5. 

It holds that

$$\sum _{i=1}^3{∥U_i(t,·)∥}_{L^{\infty }\left({\Omega }_{+}\right)}\le C{t}^{-\frac{n}{2}}{\mathrm{e}}^{-\gamma t},t>0,c>c^{*},$$

(39)

and

$$\sum _{i=1}^3{∥U_i(t,·)∥}_{L^{\infty }\left({\Omega }_{+}\right)}\le C{t}^{-\frac{n}{2}},t>0,c=c^{*},$$

(40)

where γ is a constant defined in (43).

Proof. 

It follows from (23) and Lemma 4 that $(U_1,U_2,U_3)$ satisfies, for $c>c^{*}$

$$\left\{\begin{array}{c}{U}_{1t}+c{U}_{1{\xi }_1}\le \Delta U_1+r_1\left(k_1-{\varphi }_1\right)U_1+r_1{\varphi }_1\left[a_{12}{U}_2^{\tau }+a_{13}{U}_3^{\tau }-U_1\right],t>0,\xi \in {\Omega }_{+},\hfill \\ U_{2t}+c{U}_{2{\xi }_1}\le \Delta U_2+r_2\left(k_2-{\varphi }_2\right)U_2+r_2{\varphi }_2\left[a_{21}{U}_1^{\tau }+a_{23}{U}_3^{\tau }-U_2\right],t>0,\xi \in {\Omega }_{+},\hfill \\ U_{3t}+c{U}_{3{\xi }_1}\le \Delta U_3+r_3\left(1-a_{31}{\varphi }_1^{\tau }-a_{32}{\varphi }_2^{\tau }\right)U_3+r_3[1-{\varphi }_3][a_{31}{U}_1^{\tau }+a_{32}{U}_2^{\tau }-U_3],t>0,\xi \in {\Omega }_{+},\hfill \\ U_i{|}_{{\xi }_1={\zeta }_0}\le C_2{(1+t)}^{-\frac{n}{2}}{\mathrm{e}}^{-{\epsilon }_{\tau }{\sigma }_1\left(c\right)t},t>0,({\xi }_2,\cdots ,{\xi }_n)\in {\mathbb{R}}^{n-1},\hfill \\ U_i{|}_{t=s}=U_{i0}(s,\xi ),s\in [-\tau ,0],\xi \in {\Omega }_{+},\hfill \end{array}\right.$$

(41)

and for $c=c^{*}$,

$$\left\{\begin{array}{c}{U}_{1t}+c{U}_{1{\xi }_1}\le \Delta U_1+r_1\left(k_1-{\varphi }_1\right)U_1+r_1{\varphi }_1\left[a_{12}{U}_2^{\tau }+a_{13}{U}_3^{\tau }-U_1\right],t>0,\xi \in {\Omega }_{+},\hfill \\ U_{2t}+c{U}_{2{\xi }_1}\le \Delta U_2+r_2\left(k_2-{\varphi }_2\right)U_2+r_2{\varphi }_2\left[a_{21}{U}_1^{\tau }+a_{23}{U}_3^{\tau }-U_2\right],t>0,\xi \in {\Omega }_{+},\hfill \\ U_{3t}+c{U}_{3{\xi }_1}\le \Delta U_3+r_3\left(1-a_{31}{\varphi }_1^{\tau }-a_{32}{\varphi }_2^{\tau }\right)U_3+r_3[1-{\varphi }_3][a_{31}{U}_1^{\tau }+a_{32}{U}_2^{\tau }-U_3],t>0,\xi \in {\Omega }_{+},\hfill \\ U_i{|}_{{\xi }_1={\zeta }_0}\le C_3{(1+t)}^{-\frac{n}{2}},t>0,({\xi }_2,\cdots ,{\xi }_n)\in {\mathbb{R}}^{n-1},\hfill \\ U_i{|}_{t=s}=U_{i0}(s,\xi ),s\in [-\tau ,0],\xi \in {\Omega }_{+}.\hfill \end{array}\right.$$

(42)

Let

$$0<\gamma <min\left\{{\epsilon }_{\tau }{\sigma }_1\left(c\right),\delta \right\},$$

(43)

where

$$\delta =min\left\{r_1{k}_1(1-a_{12}-a_{13}),r_2{k}_2(1-a_{21}-a_{23}),r_3(a_{31}{k}_1+a_{32}{k}_2-1)\right\}.$$

Because of ${lim}_{{\xi }_1\to +\infty }{\varphi }_i\left({\xi }_1\right)=k_i$ and assumption (P2), we can choose ${\zeta }_0$ and $t_{*}$ large enough to ensure that

$$\left\{\begin{array}{c}{r}_1{\varphi }_1\left({\xi }_1\right)(1-a_{12}-a_{13})+r_1[{\varphi }_1\left({\xi }_1\right)-k_1]-\frac{n}{2(1+t+\tau )}-\gamma >0,\hfill \\ r_2{\varphi }_2\left({\xi }_1\right)(1-a_{21}-a_{23})+r_2[{\varphi }_2\left({\xi }_1\right)-k_2]-\frac{n}{2(1+t+\tau )}-\gamma >0,\hfill \\ r_3(1-a_{31}-a_{32})[1-{\varphi }_3\left({\xi }_1\right)]-\frac{n}{2(1+t+\tau )}-\gamma +r_3\left[a_{31}{\varphi }_1({\xi }_1-c\tau )+a_{31}{\varphi }_2({\xi }_1-c\tau )-1\right]>0,\hfill \end{array}\right.$$

(44)

for ${\xi }_1>{\zeta }_0$ and $t>t_{*}$.

When $c>c^{*}$, let

$${\overline{U}}_i(t,\xi )=C_4{(1+t+\tau )}^{-\frac{n}{2}}{\mathrm{e}}^{-\gamma t},t>0,i=1,2,3,$$

where $C_4>C_2$ large enough such that ${\overline{U}}_i(t,\xi )\ge U_i(t,\xi ),t\in [0,t_{*}]\times {\mathbb{R}}^n$ and $\gamma \in (0,min\{{\epsilon }_{\tau }{\sigma }_1\left(c\right),\delta \})$. By direct computations, we can verify that $({\overline{U}}_1,{\overline{U}}_2,{\overline{U}}_3)$ is an upper solution to (41) in the form

$$\left\{\begin{array}{c}{\overline{U}}_{1t}+c{\overline{U}}_{1{\xi }_1}\ge \Delta {\overline{U}}_1+r_1\left(k_1-{\varphi }_1\right){\overline{U}}_1+r_1{\varphi }_1\left[a_{12}{\overline{U}}_2^{\tau }+a_{13}{\overline{U}}_3^{\tau }-{\overline{U}}_1\right],t>t_{*},\xi \in {\Omega }_{+},\hfill \\ {\overline{U}}_{2t}+c{\overline{U}}_{2{\xi }_1}\ge \Delta {\overline{U}}_2+r_2\left(k_2-{\varphi }_2\right){\overline{U}}_2+r_2{\varphi }_2\left[a_{21}{\overline{U}}_1^{\tau }+a_{23}{\overline{U}}_3^{\tau }-{\overline{U}}_2\right],t>t_{*},\xi \in {\Omega }_{+},\hfill \\ {\overline{U}}_{3t}+c{\overline{U}}_{3{\xi }_1}\ge \Delta {\overline{U}}_3+r_3\left(1-a_{31}{\varphi }_1^{\tau }-a_{32}{\varphi }_2^{\tau }\right){\overline{U}}_3+r_3[1-{\varphi }_3][a_{31}{\overline{U}}_1^{\tau }+a_{32}{\overline{U}}_2^{\tau }-{\overline{U}}_3],t>t_{*},\xi \in {\Omega }_{+},\hfill \\ {\overline{U}}_i{|}_{{\xi }_1={\zeta }_0}\ge C_2{(1+t)}^{-\frac{n}{2}}{\mathrm{e}}^{-{\epsilon }_{\tau }{\sigma }_1\left(c\right)t},t>t_{*},({\xi }_2,\cdots ,{\xi }_n)\in {\mathbb{R}}^{n-1},\hfill \\ {\overline{U}}_i{|}_{t=s}\ge {\overline{U}}_i(s,\xi ),s\in [-\tau ,0],\xi \in {\Omega }_{+}.\hfill \end{array}\right.$$

(45)

Hence, for $c>c^{*}$, we have

$$0\le U_i(t,\xi )\le {\overline{U}}_i(t,\xi )=C_4{(1+t+\tau )}^{-\frac{n}{2}}{\mathrm{e}}^{-\gamma t},t>0,\xi \in {\Omega }_{+},i=1,2,3.$$

(46)

When $c=c^{*}$, note that

$${\overline{U}}_i(t,\xi )=C_5{(1+t+\tau )}^{-\frac{n}{2}},t>0,i=1,2,3,$$

where $C_5>C_3$ large enough such that ${\overline{U}}_i(t,\xi )\ge U_i(t,\xi ),t\in [0,t_{*}]\times {\mathbb{R}}^n$. Similarly, $({\overline{U}}_1,{\overline{U}}_2,{\overline{U}}_3)$ satisfies

$$\left\{\begin{array}{c}{\overline{U}}_{1t}+c{\overline{U}}_{1{\xi }_1}\ge \Delta {\overline{U}}_1+r_1\left(k_1-{\varphi }_1\right){\overline{U}}_1+r_1{\varphi }_1\left[a_{12}{\overline{U}}_2^{\tau }+a_{13}{\overline{U}}_3^{\tau }-{\overline{U}}_1\right],t>t_{*},\xi \in {\Omega }_{+},\hfill \\ {\overline{U}}_{2t}+c{\overline{U}}_{2{\xi }_1}\ge \Delta {\overline{U}}_2+r_2\left(k_2-{\varphi }_2\right){\overline{U}}_2+r_2{\varphi }_2\left[a_{21}{\overline{U}}_1^{\tau }+a_{23}{\overline{U}}_3^{\tau }-{\overline{U}}_2\right],t>t_{*},\xi \in {\Omega }_{+},\hfill \\ {\overline{U}}_{3t}+c{\overline{U}}_{3{\xi }_1}\ge \Delta {\overline{U}}_3+r_3\left(1-a_{31}{\varphi }_1^{\tau }-a_{32}{\varphi }_2^{\tau }\right){\overline{U}}_3+r_3[1-{\varphi }_3][a_{31}{\overline{U}}_1^{\tau }+a_{32}{\overline{U}}_2^{\tau }-{\overline{U}}_3],t>t_{*},\xi \in {\Omega }_{+},\hfill \\ {\overline{U}}_i{|}_{{\xi }_1={\zeta }_0}\ge C_3{(1+t)}^{-\frac{n}{2}},t>t_{*},({\xi }_2,\cdots ,{\xi }_n)\in {\mathbb{R}}^{n-1},\hfill \\ {\overline{U}}_i{|}_{t=s}\ge {\overline{U}}_i(s,\xi ),s\in [-\tau ,0],\xi \in {\Omega }_{+}.\hfill \end{array}\right.$$

Thus, for $c=c^{*}$, we have

$$0\le U_i(t,\xi )\le {\overline{U}}_i(t,\xi )=C_5{(1+t+\tau )}^{-\frac{n}{2}},t>0,\xi \in {\Omega }_{+},i=1,2,3.$$

(47)

Then, (39) and (40) can be immediately derived by (46) and (47). □

It follows from Lemmas 4 and 5 that we can obtain decay rates of $U(t,\xi )$ in $L^{\infty }\left({\mathbb{R}}^n\right)$.

Lemma 6. 

It holds that

$$\sum _{i=1}^3{∥U_i(t,·)∥}_{L^{\infty }\left({\mathbb{R}}^n\right)}\le C{t}^{-\frac{n}{2}}{\mathrm{e}}^{-{\epsilon }_{\tau }\sigma t},t>0,c>c^{*},$$

and

$$\sum _{i=1}^3{∥U_i(t,·)∥}_{L^{\infty }\left({\mathbb{R}}^n\right)}\le C{t}^{-\frac{n}{2}},t>0,c=c^{*},$$

where $0<\sigma <min\{{\sigma }_1\left(c\right),\delta /{\epsilon }_{\tau }\}$.

Since $U_i(t,\xi )=u_i^{+}(t,x)-{\varphi }_i(x·e_1+ct),i=1,2,3$, we have the convergence result of the solution as follows.

Lemma 7. 

It holds that

$$\underset{x\in {\mathbb{R}}^n}{sup}\left|u_i^{+}(t,x)-{\varphi }_i(x·e_1+ct)\right|\le C{t}^{-\frac{n}{2}}{\mathrm{e}}^{-{\epsilon }_{\tau }\sigma t},c>c^{*},$$

and

$$\underset{x\in {\mathbb{R}}^n}{sup}\left|u_i^{+}(t,x)-{\varphi }_i(x·e_1+ct)\right|\le C{t}^{-\frac{n}{2}},c=c^{*},$$

for $t>0,i=1,2,3$, where $0<\sigma <min\{{\sigma }_1\left(c\right),\delta /{\epsilon }_{\tau }\}$.

In what follows, we prove our main result of Theorem 2.

Proof of Theorem 2. 

For $c\ge c^{*}$, let $\xi =x+ct·e_1$ and

$$V_i(t,\xi )={\varphi }_i(x·e_1+ct)-u_i^{-}(t,x),V_{i0}(s,\xi )={\varphi }_i(x·e_1+cs)-u_{i0}^{-}(s,x),$$

for $t>0,s\in [-\tau ,0],x\in {\mathbb{R}}^n,i=1,2,3$. We can obtain that $u_i^{-}(t,x)$ converges to ${\varphi }_i(x·e_1+ct)$, namely,

$$\underset{x\in {\mathbb{R}}^n}{sup}\left|u_i^{-}(t,x)-{\varphi }_i(x·e_1+ct)\right|\le C{t}^{-\frac{n}{2}}{\mathrm{e}}^{-{\epsilon }_{\tau }\sigma t},c>c^{*},$$

and

$$\underset{x\in {\mathbb{R}}^n}{sup}\left|u_i^{-}(t,x)-{\varphi }_i(x·e_1+ct)\right|\le C{t}^{-\frac{n}{2}},c=c^{*},$$

for $t>0,i=1,2,3$, where $0<\sigma <min\{{\sigma }_1\left(c\right),\delta /{\epsilon }_{\tau }\}$. Since

$$\begin{array}{c}0\le u_i^{-}(t,x)\le u_i(t,x)\le u_i^{+}(t,x)\le k_i,\hfill \\ 0\le u_i^{-}(t,x)\le {\varphi }_i(x·e_1+ct)\le u_i^{+}(t,x)\le k_i,\hfill \end{array}$$

for $x\in {\mathbb{R}}^n,t>0,i=1,2,3$, by the squeeze argument, we have

$$\underset{x\in {\mathbb{R}}^n}{sup}\left|u_i(t,x)-{\varphi }_i(x·e_1+ct)\right|\le C{t}^{-\frac{n}{2}}{\mathrm{e}}^{-{\epsilon }_{\tau }\sigma t},c>c^{*},$$

and

$$\underset{x\in {\mathbb{R}}^n}{sup}\left|u_i(t,x)-{\varphi }_i(x·e_1+ct)\right|\le C{t}^{-\frac{n}{2}},c=c^{*},$$

for $t>0,i=1,2,3$. □

4. Numerical Simulations

We prove exact planar traveling waves without time delay in Appendix A, and it is helpful to construct the initial data in numerical simulations. In this section, we conduct some numerical simulations to support the main result.

Next, we consider the system (11) with the Neumann boundary conditions, which are described as

$$\frac{\partial u_1(t,x)}{\partial n}=\frac{\partial u_2(t,x)}{\partial n}=\frac{\partial u_3(t,x)}{\partial n}=0,t\ge 0,x\in \Omega ,$$

(48)

and the initial values satisfy

$$u_1(s,x)=u_{10}(s,x),u_2(s,x)=u_{20}(s,x),u_3(s,x)=u_{30}(s,x),s\in [-\tau ,0],x\in \Omega ,$$

(49)

where $\Omega$ is a bounded domain in ${\mathbb{R}}^n$ (here, we assume $n=1$ for simplicity) with smooth boundary $\partial \Omega$, $\frac{\partial }{\partial n}$ represents the outward normal derivative on $\partial \Omega$; the homogenous Neumann boundary conditions indicate that the species cannot move across the boundary $\partial \Omega$.

In system (11), we choose $r_1=r_2=r_3=1,a_{12}=\frac{1}{2},a_{21}=\frac{1}{4},a_{13}=a_{23}=\frac{1}{8},a_{31}=a_{32}=\frac{7}{8},\tau =1$. Then, system (11) with the above coefficients has two steady states $E_3=(0,0,1)$ and $E_4=(\frac{12}{7},\frac{10}{7},0)$. From Proposition 1, system (11) admits a planar traveling wave with speed $c>c_{*}=max\left\{2\sqrt{r_1(1-a_{13})},2\sqrt{r_2(1-a_{23})},2\sqrt{r_3{c}_0}\right\}=2.6458$. Additionally, when the initial data satisfy (17) and the initial perturbation satisfies (18), the planar traveling wave with speed $c>c^{*}=2\sqrt{\nu \left(A\right)+\nu \left(B\right)}=3.7177$ is exponentially stable in $L^{\infty }\left({\mathbb{R}}^n\right)$, and the planar traveling wave solution with speed $c=c^{*}$ is algebraically stable in $L^{\infty }\left({\mathbb{R}}^n\right)$.

Now, we choose $\Omega =[-60,60],\tau =1$, and the initial data

$$\left\{\begin{array}{c}{u}_{10}(s,x)=\frac{12}{7}q\left(x\right)\left[1-\frac{1}{4}{(1-tanh\left(x\right))}^2\right],\hfill \\ u_{20}(s,x)=\frac{10}{7}q\left(x\right)\left[1-\frac{1}{4}{(1-tanh\left(x\right))}^2\right],\hfill \\ u_{30}(s,x)=q\left(x\right)\left[1-\frac{1}{4}{(1+tanh\left(x\right))}^2\right],\hfill \end{array}x\in [-60,60],s\in \left[-1,0\right],\right.$$

(50)

where $q\left(x\right)$ is the mollification function of $p\left(x\right)$ and

$$q\left(x\right)=\left\{\begin{array}{cc}\frac{15}{16}\left[1-\frac{1}{31+\frac{1}{2}x}cos\left(x\right)\right],\hfill & x\in [-30,30],\hfill \\ 1,\hfill & x\in \Omega \setminus [-30,30].\hfill \end{array}\right.$$

It is obvious that the initial conditions in Theorem 2 can be ensured by the initial data (50). Using MATLAB 2020b, we compute the numerical solution of (13) (see Figure 1). Figure 1 shows that the solution of system (11) eventually converges to the equilibrium $E_4=\left(\frac{12}{7},\frac{10}{7},0\right)$, indicating that species $u_1$ and $u_2$ survive while $u_3$ becomes extinct. For $t\ge 4$, the solution of system (11) exhibits the behavior of a stable planar traveling wave, maintaining its shape over time. Additionally, Figure 1 reveals that the solution of system (11) propagates from the positive to the negative direction along the x-axis. These facts can be seen more clearly in Figure 2, Figure 3 and Figure 4.

An ecological phenomenon corresponding to these findings can be described as follows. Initially, system (11) contains only the native species $u_3$, with species $u_1$ and $u_2$ absent. Later, two exotic cooperative species, $u_1$ and $u_2$, invade the system and compete with $u_3$, as governed by system (11). This raises a natural question: can all three species coexist? Based on our results, the native species $u_3$ will ultimately go extinct, while the invasive species $u_1$ and $u_2$ will persist. Moreover, Theorem 2 indicates that this process remains stable under suitable initial perturbations.

5. Conclusions

In this paper, we are devoted to studying the three-species Lotka–Volterra model with competition and cooperation. The reaction term takes the form of nonlocal nonlinearity. We prove the multidimensional stability of traveling waves in n-dimensional space. Specifically, in the case of the speed $c>c^{*}$, we show the exponential stability of traveling waves, which is in the form of $t^{-\frac{n}{2}}{\mathrm{e}}^{-{\epsilon }_{\tau }\sigma t}$ in $L^{\infty }\left({\mathbb{R}}^n\right)$, where $\sigma >0$ represents a constant and ${\epsilon }_{\tau }\in (0,1)$ depends on the time delay $\tau >0$. Here, the time delay is a decreasing function ${\epsilon }_{\tau }=\epsilon \left(\tau \right)$, which makes the convergence rate of the plane-traveling wave slow down. Additionally, when the speed $c=c^{*}$, we obtain the algebraic stability in the form of $t^{-\frac{n}{2}}$. The adapted methods are the Fourier transform and the weighted energy method with a finer weight function. Finally, we carry out some numerical simulations to reveal the theoretical results. The exact planar traveling waves without delay play a key role in the construction of the initial data. In the following, it is interesting to consider the influence of fractional Laplacian diffusion on the stability of traveling wave solutions.