Introduction

In nature, all organisms migrate or disperse. This can take diverse forms such as walking, swimming, flying, or being transported by wind or flowing water, see Shigesada and Kawasaki [44]. Such a migration, or dispersion, can be to some extent related to human activities that bring drastic changes in the global environment. According to [44], dispersive movements become noticeable when an offspring or a seed leaves its natal site, or when an organism’s habitat deteriorates from overcrowding. The spatially explicit ecological theories of such events have been made possible thanks to successful development of mathematical models that have played a central role in the description of migrations, see Friedman [18], Shigesada and Kawasaki [44], Okubo and Levin [37], Cantrell and Cosner [12], Volpert [52], Logan [33, 34], Perthame [39], and the references therein.

Mathematical literature dealing with the species’ spread mostly relies on reaction–diffusion equations that assume that the dispersal is governed by random diffusion and that it, along with the growth processes, take place continuously in time and space, (Cantrell and Cosner [12], Lewis and Li [27]). This approach has had a remarkable success in explaining the rates at which species have invaded large open environments, see Shigesada and Kawasaki [44], Okubo and Levin [37], Cantrell and Cosner [12], Lewis and Li [27], Volpert [52], Logan [33, 34] and Perthame [39]. However, it is well-known that ecological species may experience several phenomena that, depending on circumstances, can be either time-continuous (growth, death, birth, release, etc.), or time-discrete (harvest, birth, death, release, etc.), see also Ma and Li [35], Dumont and Tchuenche [13], Yatat Djeumen et al. [60], Yatat Djeumen [58] and the references therein. In the case of time-discrete perturbations, whose duration is negligible in comparison with the duration of the process, it is natural to assume that these perturbations act instantaneously; that is, in the form of impulses (Lakshmikantham et al. [26], Bainov and Simeonov [6]). Hence, there is a need to create a meaningful mathematical framework to analyze models leading to, say, impulsive reaction–diffusion equations.

There are several works that considered the impact of impulsive perturbations on the dynamics of a system, both in the space-implicit and space-explicit case (see also Sects. 1.1 and 1.2). Very often, impulsive perturbations in space-implicit mathematical models result in the occurrence of periodic solutions in the model (e.g. Ma and Li [35], Yatat Djeumen et al. [60] and references therein). For space-explicit impulsive mathematical models in bounded domains, in addition to periodic solutions that may occur, the issue of minimal domain has been also addressed (e.g. Lewis and Li [27], Yatat Djeumen and Dumont [61]). On the other hand, in unbounded domains the problems of the existence of traveling wave solutions and the computation of the spreading speeds are rarely addressed (see below).

The study of traveling waves, as well as the computation of spreading speeds for monotone systems of reaction–diffusion equations, have been done by several authors (e.g. Weinberger et al. [56], Weinberger [54, 55], Lewis et al. [28], Li et al. [29], Volpert [52], Yatat Djeumen et al. [59] and references therein). However, little is known about the case of impulsive reaction–diffusion equations (see e.g. Yatat Djeumen and Dumont [61], Fazly et al. [15]). For a scalar impulsive reaction–diffusion equation, Yatat Djeumen and Dumont [61] considered the Fisher-Kolmogorov-Petrowsky-Piscounov (FKPP) equation and obtained conditions under which an invasive traveling wave, connecting the extinction steady state and the positive steady state may exist (see also Lewis and Li [27]). To the best of our knowledge, the existence of traveling waves for system of impulsive reaction–diffusion equation was studied only in Huang et al. [23], and only in a particular case. Precisely, the authors considered a stage-structured population model, were only one stage (or state variable) experiences a spatial diffusion, while the others are stationary. This assumption leads to a partially degenerate system of impulsive reaction–diffusion equations. Moreover, they also assumed the reaction term for the diffusing state variable to be linear.

In the following, we provide a brief literature review of both space-implicit and space-explicit impulsive mathematical models.

Space-Implicit Impulsive Models

Several types of perturbations, or instantaneous phenomena, have been considered as pulse events in mathematical models. These include events such as birth, vaccination, release, harvest, or fire events : see Mailleret and Lemesle [36] where various examples are given and references therein; see also White et al. [57], Dumont and Tchuenche [13], Strugarek et al. [45], Bliman et al. [9], Anguelov et al. [4], Liu et al. [31], Zhao et al. [65], Yatat Djeumen and Dumont [61], Yatat Djeumen et al. [63] and references therein. The resulting impulsive models were rigorously analyzed by the well-known theory due to Lakshmikantham et al. [26], Bainov and Simeonov [6] and [5], or Lakmeche and Arino [25]. We note that, by using a suitable comparison argument, the standard theory of ordinary differential equations (Hale [19, 20]) can also be used.

Modelling Releases as Pulse Events

In the framework of biological control of pests or vectors of infectious diseases, the sterile insect technique (SIT) is one of the promising ones. SIT control generally consists in massive releases of sterile insects in the targeted area in order to eliminate, or at least to lower the pest population under a certain threshold (Anguelov et al. [4]). Generally, SIT releases are done periodically. That is why several authors modelled the release process as a periodic impulsive event, while keeping the continuous-time differential equation framework for the birth, growth, death, or mating process (e.g. White et al. [57], Dumont and Tchuenche [13], Strugarek et al. [45], Bliman et al. [9] and the references therein). Based on the qualitative analysis of their impulsive models, the authors were able to derive meaningful relations between the period and the size of the releases in order to achieve the elimination of the vector or pest population in the long term.

In the context of interacting species such as prey-predator interactions, there exist mathematical models that tackled periodic releases of one of the interacting species (e.g. prey only, predator only). They considered periodic impulses as to model periodic release events. The authors found relations involving the pulse time period and the amount of released species that precluded extinction, in the long term dynamics, of interacting populations (see for instance Zhang et al. [64], Zhao et al. [65]).

Modelling Fires as Pulse Events in Tree–Grass Interactions in Fire-prone Savannas

Maintaining the balance between the grass and the trees in savanna is of utmost importance for both human and animal populations living in such areas. The problem is that in typical circumstances the trees encroach on the grassland making the environment inhabitable for many species. It turns out that periodic fires, either natural or man-made, are one of the way to maintain an acceptable equilibrium. Thus, several mathematical models have been developed to study tree–grass interactions in fire-prone savannas (see the review of Yatat Djeumen et al. [63]). It is well-known that in such savannas, grass biomass that dries up during dry seasons is used as fuel for fires (see op. cit. and references therein). Hence, low grass biomass will result in fires of lower intensity. We emphasize that here we are not dealing with ‘forest fires’, like that usually occurring in Europe, North America, or, sometimes, in Australia, but only with grass-induced fires. Thus, in particular, a dense tree population results in a lower fire intensity and consequently, in a lower tree destruction. Some of the models take into account fire as a time-continuous forcing in tree–grass interactions. However, as pointed out in Yatat Djeumen et al. [63], it is questionable whether it makes sense to model fire as a permanent forcing that continuously removes a fraction of the fire sensitive biomass. Indeed, since several months and even years can pass between two successive fires, they can be rather considered as instantaneous perturbations of the savanna ecosystem (see also Yatat Djeumen [58], Yatat Djeumen et al. [60], Tchuinté et al. [47, 48]). Several recent papers have proposed to model fires either as stochastic events, while keeping the continuous-time differential equation framework (Baudena et al. [7], Beckkage et al. [8], Synodinos et al. [46]), or by using a time-discrete model (Higgins et al. [22], Accatino et al. [1, 2], Klimasara and Tyran-Kamińska [24]). However, a drawback of many of the aforementioned recent time-discrete stochastic models (Higgins et al. [22], Baudena et al. [7], Beckage et al. [8]) is that they hardly lend themselves to analytical treatment. Thus, on the basis of recent publications (see for instance Yatat Djeumen et al. [63] and the references therein), we consider fires as impulsive time-periodic events. While certainly an approximation, such an approach results in impulsive differential equation models which are a good compromise combining the impact of time-discrete fires with a time-continuous process of the vegetation growth. As such, they are analytically tractable, while at the same time remain reasonably realistic.

Now we recall the minimalistic tree–grass interactions model with pulse fires that will be used later in the paper (see Sect. 3). We assume that the trees and grass form an amensalistic system in which grass is harmed by trees (which, for instance, block the sunlight) but itself does not affect them. The fires occur periodically every \(\tilde{\tau }\) units of time. We denote by \(T_n\) (resp. \(G_n\)) the tree (resp. grass) biomass during the inter-fire season number \(n\in \mathbb {N}^*\), where \(\mathbb {N}^*=\{1,2,\ldots \}\). Following the formalism of the recursion equations (see Weinberger et al. [56], Yatat Djeumen and Dumont [61]), the resulting minimalistic system of equations governing the trees–grass interactions with periodic fires is given by (see also Yatat Djeumen et al. [63] and references therein)

$$\begin{aligned} \left\{ \begin{array}{l} \left. \begin{array}{lcl} \displaystyle \frac{d T_n}{d t}&{}=&{}\gamma _{T}(\mathbf{W })\left( 1-\displaystyle \frac{T_n}{K_{T}(\mathbf{W })}\right) T_n-\delta _{T}T_n,\\ &{}&{}\\ \displaystyle \frac{d G_n}{d t}&{}=&{}\gamma _{G}(\mathbf{W })\left( 1-\displaystyle \frac{G_n}{K_{G}(\mathbf{W })}\right) G_n-\delta _{G}G_n-\eta _{TG}T_nG_n,\\ \end{array} \right. {0< t\le \tilde{\tau }},\\ \\ \left. \begin{array}{l} T_{n+1}(0)={[1}-\psi (T_n(\tilde{\tau }))w_G(G_n(\tilde{\tau })){]}T_n(\tilde{\tau }),\\ \\ G_{n+1}(0)=(1-\eta )G_n(\tilde{\tau }),\\ \end{array} \right. \end{array} \right. \end{aligned}$$
(1)

with non negative initial conditions \((T_1(0), G_1(0))\).

Here, between two successive fires; that is, in the inter-fire season n, the dynamics of both the tree and grass biomasses is modelled by the first two equations of system (1), where \(\gamma _G(\mathbf{W })\) and \(\gamma _T(\mathbf{W })\) denote the unrestricted rates of growth of the grass and the tree biomass, respectively, while \(K_{G}(\mathbf{W })\) and \(K_{T}(\mathbf{W })\) are the carrying capacities for grass and the trees, respectively. All these functions are assumed to be increasing and bounded functions of the water availability \(\mathbf{W }\) which is supposed to be known: see [62] for further details about the definition of the growth rates and the carrying capacities with respect to \(\mathbf{W }\). Further, \(\delta _{G}\) and \(\delta _{T}\) denote, respectively, the rates of the grass and the tree biomass loss due to natural causes, herbivores (grazing and/or browsing) or human actions, while \(\eta _{TG}\) denotes rate of the loss of the grass biomass due to the existence of trees per units of the biomasses.

At the end of the inter-fire season a fire occurs and impacts both the tree and grass biomasses. Thus there is an update of the biomasses at the beginning of the next inter-fire season. This event is modelled by the last two equations of system (1). We model the fire intensity, denoted by \(w_G,\) using the Holling III functional response (see (66)) ; that is, we assume that \(w_G\) is an increasing function tending to 1 as the argument increases and satisfying \(w_G(0) = w'_G(0)=0\). This reflects the fact that the fire becomes significant only when the grass biomass becomes sufficiently large. Impulsive fire-induced tree/shrub mortality, denoted by \(\psi\), is assumed to be a positive, decreasing, and nonlinear function of the tree biomass (see (67)). Indeed, fires affect differently large and small trees since fires with high intensity (flames more than ca. 2 metres high) cause greater mortality of shrubs and topkill of trees while fires of lower intensity (flames less than ca. 2 metres high) kill only shrubs and subshrubs (see also [62, 63] and references therein). Further, \(\eta\) is the specific loss of the grass biomass due to the fire. To avoid the extinction of either \(T_n,\) or \(G_n\), we assume that (see also Yatat Djeumen et al. [63])

$$\begin{aligned} \gamma _{G}(\mathbf{W })-\delta _{G}>0\quad \text{ and }\quad \gamma _{T}(\mathbf{W })-\delta _{T}>0. \end{aligned}$$
(2)

System (1) will be studied in Sect. 3.

Readers are referred to Yatat Djeumen et al. [63] for the derivation and analysis of system (1) following the formalism of Lakshmikantham et al. [26].

Space-Explicit Impulsive Models

The formulation of space-explicit impulsive models generally consists in the addition of local or non-local spatial operators to a temporal impulsive model. In Akhmet et al. [3], Li et al. [30] and Liu et al. [32], the authors considered impulsive reaction–diffusion equations to model spatio-temporal dynamics of ecological species with prey-predator interactions and experiencing pulse and periodic perturbations like harvest, release, etc. The spatial movement of species was modelled by the Laplace operator with a constant diffusion rate. Qualitative analysis of these models was done by using the theory of sectorial operators (Henry [21], Rogovchenko [41, 42], Li et al. [30]) and comparison arguments (Rogovchenko [40], Walter [53], Liu et al. [32], Akhmet et al. [3]). More precisely, the authors obtained some conditions involving the pulse time period that ensured the permanence of the predator-prey system and the existence of a unique globally stable periodic solution (Akhmet et al. [3], Li et al. [30], Liu et al. [32]). Vasilyeva et al. [50] dealt with the question of persistence versus extinction in a single population model featuring a non-local impulsive reaction–advection–diffusion model for an insect population. The non-local term was used to describe the dispersal of the adult insects by flight. The authors employed a dispersal kernel that gave the probability density function of the signed dispersal distances.

We note that the study of the species spread and their wave speeds, when they experience impulsive and periodic perturbations, in the case of scalar equations was done in Vasilyeva et al. [50], Lewis and Li [27], or Yatat Djeumen and Dumont [61]. However, until recently, systems of impulsive reaction–diffusion have not received much attention (Huang et al. [23] and Fazly et al. [15]). We aim to address this question here by extending the minimalistic trees–grass interactions model (1). To this end, we assume that both the woody and herbaceous plants can propagate in space through diffusion; see Yatat Djeumen et al. [59] for a discussion of the construction of the trees–grass interactions partial differential equations models. The resulting minimalistic system of impulsive reaction–diffusion equations is then given by

$$\begin{aligned} \left\{ \begin{array}{l} \left. \begin{array}{l} \displaystyle \frac{\partial T_n}{\partial t}=d_T(\mathbf{W })\displaystyle \frac{\partial ^2 T_n}{\partial x^2}+\gamma _{T}(\mathbf{W })\left( 1-\displaystyle \frac{T_n}{K_{T}(\mathbf{W })}\right) T_n-\delta _{T}T_n,\\ \\ \displaystyle \frac{\partial G_n}{\partial t}=d_G(\mathbf{W })\displaystyle \frac{\partial ^2 G_n}{\partial x^2}+\gamma _{G}(\mathbf{W })\left( 1-\displaystyle \frac{G_n}{K_{G}(\mathbf{W })}\right) G_n-\delta _{G}G_n-\eta _{TG}T_nG_n,\\ \end{array} \right. {0< t\le \tilde{\tau }},\quad x\in {\mathbb {R}},\\ \\ \left. \begin{array}{l} T_{n+1}(x,0)={[}1-\psi (T_n(x,\tilde{\tau }))w_G(G_n(x,\tilde{\tau })){]}T_n(x,\tilde{\tau }),\\ \\ G_{n+1}(x,0)=(1-\eta )G_n(x,\tilde{\tau }),\\ \end{array} \right. \end{array} \right. \end{aligned}$$
(3)

with given non negative initial conditions

$$\begin{aligned} (T_1(x,0), G_1(x,0)). \end{aligned}$$
(4)

In system (3), \(d_T(\mathbf{W })\) and \(d_G(\mathbf{W })\) denote the woody, respectively, herbaceous biomass spatial vegetative diffusion coefficient, while the remaining coefficients and assumptions on them are as in (1).

The aim of this paper is to give some insights into the existence of traveling wave solutions for monotone systems of impulsive reaction–diffusion equations as well as the computation of their spreading speeds. Precisely, we use the vector-valued recursion theory proposed by Weinberger et al. [56] (see also Weinberger [54, 55], Lewis et al. [28], Li et al. [29], Lewis and Li [27], and references therein) to develop a framework in which we are able to deal with the existence of traveling wave solutions for monotone systems of impulsive reaction–diffusion equations and the computation of spreading speeds. After this long introduction, Sect. 2 presents the framework for the computation of the spreading speeds and the existence of traveling wave solutions for monotone systems of impulsive reaction–diffusion equations. Section 3 deals with the application of this framework to a system of two impulsive reaction–diffusion equations that models tree–grass interactions in fire-prone savannas. We also provide some numerical illustrations of our theoretical results and, in particular, we show an approximation of the spreading speeds. The paper ends with a conclusion, where further extensions are discussed.

Preliminaries and Hypotheses

In this section we introduce the basic notation, definition and results that allow us to deal with the issues of the existence of traveling wave solutions and/or computation of the spreading speeds for impulsive reaction–diffusion (IRD) systems. Let us denote:

\(\mathbf{u }=(u_1, u_2,...,u_N)\), \(\mathbf{F }=(F_1, F_2,...,F_N)\), \(\mathbf{H }=(H_1, H_2,...,H_N)\), \(D=diag(d_1, d_2,...,d_N)\) with \(d_i>0\), for \(i=1,2,...,N\) and \(\displaystyle \frac{\partial ^2 \mathbf{u }}{\partial x^2}:=\left( \displaystyle \frac{\partial ^2 u_1}{\partial x^2},\displaystyle \frac{\partial ^2 u_2}{\partial x^2},...,\displaystyle \frac{\partial ^2 u_N}{\partial x^2}\right)\). In the case when there are no impulsive perturbations, the reaction–diffusion system is written as

$$\begin{aligned} \displaystyle \frac{\partial \mathbf{u }(x,t)}{\partial t} = D\displaystyle \frac{\partial ^2 \mathbf{u }(x,t)}{\partial x^2}+\mathbf{F }(\mathbf{u }(x,t)), \quad t> 0,\quad x\in \mathbb {R}, \end{aligned}$$
(5)

together with sufficiently smooth and nonnegative initial condition

$$\begin{aligned} \mathbf{u }(x,0) = \mathbf{u }_0(x), \quad x\in \mathbb {R}. \end{aligned}$$
(6)

For \(\tau\)-periodic impulsive perturbations, we consider the whole time interval, \([0,+\infty )\), as a succession of inter-perturbation seasons of length \(\tau\). Let us denote the state variables at time \(t\in [0, \tau ]\) and location x during the inter-perturbation season \(n\in \mathbb {N}^*\) as \(\mathbf{u }_n(t, x)=(u_{n,1}, u_{n,2},...,u_{n,N})\). Following the recursion formalism (Weinberger [54, 55], Weinberger et al. [56], Lewis and Li [27], Vasilyeva et al. [50], Fazly et al. [16], Huang et al. [23], Yatat Djeumen and Dumont [61]), the impulsive reaction–diffusion system is written as

$$\begin{aligned} \displaystyle \frac{\partial \mathbf{u }_n(x,t)}{\partial t}= & {} D\displaystyle \frac{\partial ^2 \mathbf{u }_n(x,t)}{\partial x^2}+\mathbf{F }(\mathbf{u }_n(x,t)), \quad {0< t\le \tau },\quad x\in \mathbb {R}, \end{aligned}$$
(7)

together with the updating condition

$$\begin{aligned} \mathbf{u }_{n+1}(x,0) = \mathbf{H }(\mathbf{u }_n(x,\tau )) \end{aligned}$$
(8)

and with sufficiently smooth and nonnegative initial data \({\mathbf{u }}_1(x,0)= \mathbf{u }_{0}(x)\). We note that (3) is a special case of (8).

Our work is based on the results of Li et al. [29] (see also Weinberger et al. [56]) concerning the existence of monostable traveling wave solutions and the computation of spreading speeds for systems of reaction–diffusion equations. We note that they focused on the case, when \(\mathbf{H }\) was just the time-\(\tau\)-map operator solution, \(Q_\tau\), of system (7). Since, however, the reduction of an IRD system to the recursion form does not depend on the updating condition, the results of op. cit. on the existence of traveling waves for recursions and the spreading speeds determined by them can be used verbatim to the recursion obtained from (7)–(8). Thus we recall the relevant results from [29].

We first assume that \(\mathbf{F }\), \(\mathbf{H }\) and the initial data are such that the IRD system (7)–(8) admits a unique nonnegative classical solution for each \(n \in \mathbb N^*\) (Zheng [66], Volpert [52], Logan [33, 34], Perthame [39]).

We begin with some notation (see Weinberger et al. [56], Li et al. [29]). For two vector-valued functions \(\mathbf{u }(x)\) and \(\mathbf{v} (x)\), \(\mathbf{u }(x)\le \mathbf{v} (x)\) means that \(u_i(x)\le v_i(x)\) for all \(i=1,2,...,N\) and \(x\in {\mathbb {R}}\), \(\max \{{\mathbf{u }}(x), \mathbf{v }(x)\}\) means the vector-valued function whose \(i\) th component at x is \(\max \{u_i(x), v_i(x)\}\), and \(\limsup \nolimits _{n\rightarrow \infty }\mathbf{u }^{(n)}(x)\) is the function whose \(i\)th component at x is \(\limsup \nolimits _{n\rightarrow \infty }u^{(n)}(x)\). We shall, moreover, use the usual symbol \(\mathbf{u }\gg \mathbf{v}\) if \(u_i(x) > v_i(x)\) for all i and x. We use the notation 0 for the constant vector whose all components are 0. If \(\beta \gg \mathbf{0 }\) is a constant vector, we define the set of functions

$$\begin{aligned} \mathcal {C}_\beta :=\{\mathbf{u }(x):\quad \mathbf{u }\quad \text{ is } \text{ continuous } \text{ and }\quad \mathbf{0 }\le \mathbf{u }(x)\le \beta \}. \end{aligned}$$

Let \(Q_{\tau }\) be the time-\(\tau\)-map operator solution of system (7). Then (8) can be written as

$$\begin{aligned} \mathbf{u }_{n+1}(x,0)=\mathbf{H }(Q_{\tau }[\mathbf{u }_{n}(x,0)])=:Q[\mathbf{u }_{n}(x,0)], \quad n \ge 1, \quad x\in {\mathbb {R}}, \end{aligned}$$
(9)

with the initial condition \({\mathbf{u }}_0(x)\).

In the sequel, we recall the key assumptions of Li et al. [29], related to the operator Q defined in (9). For a fixed \(y\in {\mathbb {R}},\) the translation operator by y is defined by \(T_y[{\mathbf{u }}](x)={\mathbf{u }}(x-y)\), for all \(x\in {\mathbb {R}}\).

Hypothesis 2.1

  1. i.

    The operator Q is order preserving in the sense that if \(\mathbf{u }\) and \(\mathbf{v}\) are any two functions in \(\mathcal {C}_\beta\) with \(\mathbf{v} \ge \mathbf{u }\), then \(Q[\mathbf{v} ]\ge Q[\mathbf{u }]\). In biological terms, the dynamics is cooperative.

  2. ii.

    \(Q[\mathbf{0 }]=\mathbf{0 }\), there is a constant vector \(\beta \gg \mathbf{0 }\) such that \(Q[\beta ]=\beta\), and if \(\mathbf{u }_0\) is any constant vector with \(\mathbf{u }_0\gg \mathbf{0 }\), then the constant vector \(\mathbf{u }_n\), obtained from the recursion (9), converges to \(\beta\) as n approaches infinity. This hypothesis, together with (i), imply that Q takes \(\mathcal {C}_\beta\) into itself, and that the equilibrium \(\beta\) attracts all initial functions in \(\mathcal {C}_\beta\) with uniformly positive components. There may also be other equilibria lying between \(\beta\) and the extinction equilibrium 0, in each of which at least one of the species is extinct.

  3. iii.

    Q is translation invariant. In biological terms this means that the habitat is homogeneous, so that the growth and migration properties are independent of location.

  4. iv.

    For any \(\mathbf{v} , \mathbf{u }\in \mathcal {C}_\beta\) and any fixed y, \(|Q[\mathbf{v} ](y)-Q[\mathbf{u }](y)|\) is arbitrarily small, provided \(|\mathbf{v} (x)-\mathbf{u }(x)|\) is sufficiently small on a sufficiently long interval centered at y.

  5. v.

    Every sequence \(\mathbf{v} _n\) in \(\mathcal {C}_{\beta }\) has a subsequence \(\mathbf{v} _{n_l}\) such that \(Q[\mathbf{v} _{n_l}]\) converges uniformly on every bounded set.

We are now in position to recall results of Li et al. [29] that deal with spreading speeds as well as traveling wave solutions for the IRD system (7)–(8), rewritten following the recursion formalism (9). In the sequel, we assume that Hypothesis 2.1. holds for the recursion operator Q of system (9).

Following Li et al. [29] (see also Weinberger et al. [56]), we consider a continuous \({\mathbb {R}}^N\)-valued function \(\phi (x)\) with the properties

$$\begin{aligned} \begin{array}{l} i.\quad \phi (x)\quad \text{ is } \text{ non-increasing } \text{ in }\quad x;\\ ii.\quad \phi (x)=\mathbf{0 }\quad \text{ for } \text{ all }\quad x\ge 0;\\ iii.\quad \mathbf{0 }\ll \phi (-\infty )\ll \beta . \end{array} \end{aligned}$$
(10)

We let, for all fixed \(c\in {\mathbb {R}}\), \(\mathbf{a} _0(c;s)=\phi (s)\), and define the sequence \(\mathbf{a} _n(c;s)\) by the recursion

$$\begin{aligned} \mathbf{a} _{n+1}(c;s)=\max \{\phi (s),Q[\mathbf{a} _n(c;x)](s+c)\}. \end{aligned}$$
(11)

Li et al. [29] showed that the sequence \(\mathbf{a} _n\) converges to a limit function \(\mathbf{a} (c;s)\) such that \(\mathbf{a} (c;\pm \infty )\) are equilibria of Q and \(\mathbf{a} (c;\infty )\) is independent of the initial function \(\phi\). Following this, they defined the slowest spreading speed \(c^*\le \infty\) by the equation

$$\begin{aligned} c^*=\sup \{c: \mathbf{a} (c;\infty )=\beta \}. \end{aligned}$$
(12)

The following result holds.

Theorem 1

[29, Theorem 2.1] There is an index j for which the following statement is true: Suppose that the initial function \(\mathbf{u }_0(x)\) is 0 for all sufficiently large x, and that there are positive constants \(0<\rho \le \sigma <1\) such that \({\mathbf{0 }}\le \mathbf{u }_0\le \sigma \beta\) for all x and \(\mathbf{u }_0\ge \rho \beta\) for all sufficiently negative x. Then for any positive \(\varepsilon\) the solution \(\mathbf{u }_n\) of recursion (9) has the properties

$$\begin{aligned} \lim \limits _{n\rightarrow +\infty }\left[ \sup \limits _{x\ge n(c^*+\varepsilon )}\{\mathbf{u }_n\}_j(x)\right] =0 \end{aligned}$$
(13)

and

$$\begin{aligned} \lim \limits _{n\rightarrow +\infty }\left[ \sup \limits _{x\le n(c^*-\varepsilon )}\{\beta -\mathbf{u }_n(x)\}\right] ={\mathbf{0 }}. \end{aligned}$$
(14)

That is, the \(j^{th}\) component spreads at a speed no higher than \(c^*\), and no component spreads at a lower speed.

In order to define the fastest speed \(c^*_f\), we choose \(\phi\) with the properties (10), and let \(\mathbf{b} _n(x)\) be the solution of the recursion (9) with \(\mathbf{b} _0(x)=\phi (x)\). Following Li et al. [29], we define the function

$$\begin{aligned} \mathbf{B} (c;x)=\limsup \limits _{n\rightarrow +\infty }{} \mathbf{b} _n(x+nc). \end{aligned}$$

Li et al. [29] showed that \(\mathbf{B} (c;\infty )\) is independent of the choice of the initial function \(\phi\) as long as \(\phi\) has the properties (10). We therefore can define the fastest spreading speed \(c^*_f\ge c^*\) by the formula

$$\begin{aligned} c^*_f=\sup \{c: \mathbf{B} (c;\infty )\ne \mathbf{0 }\}. \end{aligned}$$
(15)

The following result holds.

Theorem 2

[29, Theorem 2.2] There is an index i for which the following statement is true: Suppose that the initial function \(\mathbf{u }_0(x)\) is 0 for all sufficiently large x, and that there are positive constants \(0<\rho \le \sigma <1\) such that \({\mathbf{0 }}\le \mathbf{u }_0\le \sigma \beta\) for all x and \(\mathbf{u }_0\ge \rho \beta\) for all sufficiently negative x. Then for any positive \(\varepsilon\) the solution \(\mathbf{u }_n\) of the recursion (9) has the properties

$$\begin{aligned} \limsup \limits _{n\rightarrow +\infty }\left[ \inf \limits _{x\le n(c^*_f-\varepsilon )}\{\mathbf{u }_n\}_i(x)\right] >0 \end{aligned}$$
(16)

and

$$\begin{aligned} \lim \limits _{n\rightarrow +\infty }\left[ \sup \limits _{x\ge n(c^*_f+\varepsilon )}\mathbf{u }_n(x)\right] ={\mathbf{0 }}. \end{aligned}$$
(17)

That is, the ith component spreads at a speed no less than \(c^*_f\), and no component spreads at a higher speed.

Let \({\hat{\mathbf{u }}}_n, n\ge 0,\) be the solution to the recursion

$$\begin{aligned} {\hat{\mathbf{u }}}_{n+1}(x)= Q[{\hat{\mathbf{u }}}_{n}(x)],\quad n\ge 1,\quad x\in {\mathbb {R}}, \end{aligned}$$
(18)

with \({\hat{\mathbf{u }}}_{0}(x) = {\mathbf{u }}_0(x)\). Recall that a traveling wave of speed c is a solution of the recursion (9) which has the form \(\mathbf{u }_n(x,0)=\mathbf{Z} (x-nc)\) with \(\mathbf{Z} (s)\) a function in \(\mathcal {C}_{\beta }\); that is, the solution at time \(n + 1\) is simply the translate by c of its value at n. Then such a traveling wave defines a traveling wave solution for (7)–(8) in the following sense. By (9) we have

$$\begin{aligned} {\mathbf{u }}_n(x,0) = \hat{{\mathbf{u }}}_n(x) = \mathbf{Z} (x-nc) ={\mathbf{u }}_0(x-nc) \end{aligned}$$

and thus, by Hypothesis 2.1. (iii)., for \(n\tau \le t<(n+1)\tau\) we have \({\mathbf{u }}_{n+1}(x,t) = {\mathbf{u }}_{1}(x-nc,t)\). We observe that since the model is translation invariant, we obtain a traveling wave for the system (7) without the updating conditions.

Using the definitions of \(c^*\) and \(c^*_f\), we have the following result that deals with the existence of traveling wave solutions for the IRD systems (7)–(8).

Theorem 3

[29, Theorem 3.1] If \(c\ge c^*\), then there is a non-increasing traveling wave solution \(\mathbf{Z} (x-nc)\) of speed c with \(\mathbf{Z} (-\infty )=\beta\) and \(\mathbf{Z} (+\infty )\) an equilibrium other than \(\beta\).

If there is a traveling wave \(\mathbf{Z} (x-nc)\) with \(\mathbf{Z} (-\infty )=\beta\) such that for at least one component i

$$\begin{aligned} \liminf \limits _{x\rightarrow \infty }Z_i(x)=0, \end{aligned}$$

then \(c\ge c^*\). If this property is valid for all components of \(\mathbf{Z}\), then \(c\ge c^*_f\).

In practice, assumptions (iii), (iv). and (v) are typically satisfied for biologically reasonable (impulsive) models. The most challenging assumptions for IRD systems (7)–(8) are (i) and (ii).

Application to a Minimalistic Trees–Grass Interactions IRD System

In this section, we consider the minimalistic tree–grass interactions IRD system (3)–(4). Using a similar normalization procedure to Yatat Djeumen et al. [59] (see also Appendix A), system (3)–(4) becomes

$$\begin{aligned} \left\{ \begin{array}{lcl} \displaystyle \frac{\partial U_n}{\partial t} &{}=&{} U_n(1-U_n)+d_u\displaystyle \frac{\partial ^2 U_n}{\partial x^2}, \quad {0< t\le \tau },\quad x\in \mathbb {R},\\ \displaystyle \frac{\partial V_n}{\partial t} &{}=&{} \lambda V_n(1-V_n- \gamma U_n)+d_v\displaystyle \frac{\partial ^2 V_n}{\partial x^2}, \end{array} \right. \end{aligned}$$
(19)

together with the updating conditions

$$\begin{aligned} \left\{ \begin{array}{lcl} U_{n+1}(x,0) &{}=&{} {[}1- w_V(V_{n}(x,\tau ))\psi (U_n(x,\tau )){]}U_n(x,\tau ),\\ V_{n+1}(x,0) &{}=&{} (1- \eta )V_{n}(x,\tau ),\\ \end{array} \right. \end{aligned}$$
(20)

and sufficiently smooth and nonnegative initial data \(U_{1}(x,0), V_{1}(x,0)\). We are now looking for traveling wave solutions as well as the spreading speeds involving semi-trivial steady states; that is, steady states where either \(U_n=0,\) or \(V_n=0\) but not simultaneously \(U_n=0 =V_n=0.\)

Basic Properties of (19)–(20)

Let \(C_{ub}({\mathbb {R}})\) be the Banach space of bounded, uniformly continuous function on \({\mathbb {R}}\) and

$$\begin{aligned} C^{2}_b({\mathbb {R}}):=\{f\in C_{ub}({\mathbb {R}}): f''\in C_{ub}({\mathbb {R}})\}. \end{aligned}$$

\(C_{ub}({\mathbb {R}})\) and \(C^{2}_b({\mathbb {R}})\) are endowed with the following (sup) norms

$$\begin{aligned} \Vert f\Vert _{C_{ub}({\mathbb {R}})}=\Vert f\Vert _{\infty }=\sup \limits _{x\in {\mathbb {R}}}|f(x)| \end{aligned}$$
(21)

and

$$\begin{aligned} \Vert f\Vert _{C_{b}^2({\mathbb {R}})}=\Vert f\Vert _{C_{ub}({\mathbb {R}})}+\Vert f''\Vert _{C_{ub}({\mathbb {R}})}. \end{aligned}$$
(22)

\(C_{b}^2({\mathbb {R}})\) endowed with the norm \(\Vert \cdot \Vert _{C_{b}^2({\mathbb {R}})}\) is a Banach space.

We recall that we assumed that \(w_V\) was an increasing \(\mathcal {C}^1({\mathbb {R}})\) function such that for all \(V\in {\mathbb {R}}\),

$$\begin{aligned} w_V(0)=0,\quad w_V'(0)=0,\quad 0\le w_V(V)<1. \end{aligned}$$
(23)

Similarly, \(\psi\) is a decreasing \(\mathcal {C}^1({\mathbb {R}})\) function such that for all \(U\in C_{ub}({\mathbb {R}})\),

$$\begin{aligned} \psi (0)>0,\quad \psi '(0)<0,\quad 0< \psi (U)\le 1. \end{aligned}$$
(24)

For simplicity we note \((U_{n}(x,0), V_{n}(x,0))=(U_{n,0}(x),V_{n,0}(x)).\) In the sequel, we first address the question of the existence and uniqueness of solutions of the reaction–diffusion (RD) system (19) in unbounded domains.

For fixed \(n\in \mathbb {N}^*\), we set \({\mathbf {w}}=(w_1,w_2) :=(U_n,V_n)\). System (19) can be written as the abstract Cauchy problem

$$\begin{aligned} \left\{ \begin{array}{l} \displaystyle \frac{d {\mathbf {w}}}{d t} +A {\mathbf {w}}= {\mathbf {F}}({\mathbf {w}}),\\ {\mathbf {w}}(0)={\mathbf {w}}_0, \end{array} \right. \end{aligned}$$
(25)

where in the Banach space \(B=C_{ub}({\mathbb {R}})\times C_{ub}({\mathbb {R}})\) we have

$$\begin{aligned} \left\{ \begin{array}{l} D(A) = C^{2}_b({\mathbb {R}})\times C^{2}_b({\mathbb {R}}), \\ a=diag(d_u,d_v),\\ A{\mathbf {w}} = -a{\mathbf {w}}'',\\ {\mathbf {F}}: D(A) \rightarrow D(A), {\mathbf {F}}({\mathbf {w}})=(U_n(1-U_n),\lambda V_n(1-V_n-\gamma U_n)). \end{array} \right. \end{aligned}$$
(26)

For \(X\in \{C_{ub}({\mathbb {R}}),C^{2}_b({\mathbb {R}})\}\) and \((a,b)\in X\times X\) we define

$$\begin{aligned} \Vert (a,b)\Vert _{X\times X}=\Vert a\Vert _{X}+\Vert b\Vert _{X}. \end{aligned}$$

We shall consider (25) as a nonlinear perturbation of the linear part that, in this case, consists of two uncoupled diffusion equations. Thus, the corresponding semigroup is the diagonal semigroup consisting of the Gauss semigroups

$$\begin{aligned} \mathbf {S}(t) {\mathbf {w}} = diag (G_{d_u}(t)\star w_1, G_{d_v}(t)\star w_2), \quad t>0, \qquad {\mathbf {S}}(0){\mathbf {w}} ={\mathbf {w}}, \end{aligned}$$
(27)

where, for \(d=d_u,d_v\) and \(f\in X\),

$$\begin{aligned} (S_d(t)f)(x) = [G_{d}\star f](x,t) = \displaystyle \frac{1}{\sqrt{4\pi d t}}\displaystyle \int _{{\mathbb {R}}}\exp \left( -\displaystyle \frac{(x-y)^2}{{4d t}}\right) f(y)dy, \end{aligned}$$
(28)

and

$$\begin{aligned} G_d(x,t)=\displaystyle \frac{1}{\sqrt{4\pi d t}}\exp \left( -\displaystyle \frac{x^2}{{4d t}}\right) . \end{aligned}$$

Then, by e.g. [10, Section 7.3.10], the family \(\{{\mathbf {S}}(t)\}_{t\ge 0}\) is a \(C_0-\)semigroup of contractions (even analytic) on B,  with the generator (AD(A)). Furthermore, since \({\mathbf {F}}\) is a quadratic function, it is continuously Fréchet differentiable in B and therefore (25) has a unique local in time (defined on \([0,t_{max})\)) classical solution, provided \({\mathbf {w}}(0)\in D(A)\) (due to the analyticity, there is a local classical solution with \({\mathbf {w}}(0)\in B\) on \((0,t_{max})\), see e.g. [66, Theorem 2.3.5]).

Our problem is posed on the whole line and thus comparison theorems for the solutions are a little more delicate. Though in various forms they appear in many papers, see e.g. [17, 38, 49, 51] and references therein, a comprehensive proof of them, starting from the first principles, is difficult to find. Therefore, for pedagogical reasons as well as because we shall need some intermediate estimates later on, we decided to provide a simple proof for the problem at hand that uses the positivity of the semigroup \(\{S_d(t)\}_{t\ge 0}\) and the triangular structure of the nonlinearity in (19). In fact, the semigroup for the scalar problem,

$$\begin{aligned} \phi _t= & {} d\phi _{xx} + c(x,t)\phi ,\quad x \in {\mathbb {R}},\quad t\ge 0,\nonumber \\ \phi (x,0)= & {} \mathring{\phi }(x), \end{aligned}$$
(29)

where \(|c(x,t)|\le L\) on \({\mathbb {R}}\times {\mathbb {R}}_+\), is positive. Indeed, the equation can be re-written as

$$\begin{aligned} \varPhi _t= & {} d\varPhi _{xx} +C(x,t) \varPhi ,\quad x \in {\mathbb {R}}, \quad t\ge 0,\nonumber \\ \varPhi (x,0)= & {} \mathring{\phi }(x), \end{aligned}$$
(30)

where \(C(x,t)=c(x,t)+L\ge 0\) and \(\varPhi (x,t)=e^{Lt}\phi (x,t)\) and the positivity of the semigroup solving (30) follows from the Dyson-Phillips expansion [14, Theorem III.1.10]. Then, considering two solutions \(u_1\) and \(u_2\) with \(u_1(x,0)\le u_2(x,0)\) to the scalar nonlinear problem

$$\begin{aligned} u_t = d u_{xx} + F(u,t),\quad x \in {\mathbb {R}}, \quad t>0, \end{aligned}$$
(31)

on a common interval of existence \([0,t']\), where F is a differentiable function on \({\mathbb {R}}\times {\mathbb {R}}_+\), we find that \(z= u_2-u_1\) satisfies

$$\begin{aligned} z_t = d z_{xx} + c(x,t) z, \quad x \in {\mathbb {R}},\quad 0< t\le t' \end{aligned}$$
(32)

where \(c = F'((1-\theta )u_1 + \theta u_2), 0< \theta < 1,\) is bounded on \({\mathbb {R}}\times [0,t']\). By the above linear result, \(u_2 -u_1 = z\ge 0\). Returning now to (19), we see that the first equation is the Fisher equation and functions identically equal to 0 and to 1 are its solutions defined globally in time. Thus for any \(0\le U(x,0)\le 1\) we obtain \(0\le U(x,t)\le 1\) on \([0,t_{max})\). Hence U is defined globally in t and satisfies \(0\le U(x,t)\le 1\) for all \((x,t)\in {\mathbb {R}}\times {\mathbb {R}}_+\). Now, let V be the solution of the second equation in (19) on the maximum interval of existence \([0,t_{max})\),

$$\begin{aligned} V_t = d_v V_{xx} + V(1-U -V), \end{aligned}$$

with \(0\le V(x,0)\le 1\). Since the function identically equal to zero solves the above equation, as before we get \(V(x,t)\ge 0\) on \([0,t_{max})\) as long as \(V(x,0)\ge 0\). But then, using \(V(1-U-V) \le V(1-V)\) on account of \(U\ge 0\) we see, by e.g. Picard iterates, that V is dominated by the solution of the Fisher equation with the same initial condition and so, in particular, by 1. This gives the global in time existence of V and the bound \(0\le V\le 1\), and hence global in time solvability of the system (19) with initial conditions between 0 and 1.

Since the updating conditions (20) are non-increasing, if the initial data \(U_{1}(\cdot ,0), V_{1}(\cdot ,0)\) satisfy

$$\begin{aligned} \Vert U_1(\cdot ,0)\Vert _\infty \le 1 \quad \text{ and } \quad \Vert V_1(\cdot ,0)\Vert _\infty \le 1, \end{aligned}$$
(33)

then, for each \(n\in \mathbb {N}^*\), the solutions \((U_n,V_n)\) of system (25) satisfy

$$\begin{aligned} \Vert U_n\Vert _\infty \le 1\quad \text{ and }\quad \Vert V_n\Vert _\infty \le 1. \end{aligned}$$
(34)

The Existence of Equlibria of (19)–(20)

The First Coordinate Change

System (19) is monotone competitive and system (20) is not monotone. Therefore, the full system (19)–(20) is not monotone. Hence, in order to be able to apply results of Li et al. [29], we first apply a coordinate change in order to obtain a monotone cooperative system. We set

$$\begin{aligned} \left\{ \begin{array}{ccl} u_n &{} = &{} U_n, \\ v_n &{} = &{} 1-V_n, \end{array} \right. \end{aligned}$$
(35)

so that system (19)–(20) is transformed to

$$\begin{aligned} \left\{ \begin{array}{lcl} \displaystyle \frac{\partial u_n}{\partial t} &{}=&{} u_n(1-u_n)+d_u\displaystyle \frac{\partial ^2 u_n}{\partial x^2}, \quad {0< t\le \tau },\quad x\in \mathbb {R},\\ \displaystyle \frac{\partial v_n}{\partial t} &{}=&{} - \lambda v_n(1-v_n) + \lambda \gamma u_n(1-v_n)+d_v\displaystyle \frac{\partial ^2 v_n}{\partial x^2}, \end{array} \right. \end{aligned}$$
(36)

together with the updating conditions

$$\begin{aligned} \left\{ \begin{array}{lcl} u_{n+1}(x,0) &{}=&{} (1- w_v(v_{n}(x,\tau ))\psi (u_n(x,\tau )))u_n(x,\tau ),\\ v_{n+1}(x,0) &{}=&{} (1- \eta )v_{n}(x,\tau )+\eta . \end{array} \right. \end{aligned}$$
(37)

Properties (34), (33) and (23) imply that \(w_v\) is a decreasing \(\mathcal {C}^1({\mathbb {R}})\) function such that

$$\begin{aligned} w_v(1)=0,\quad w_v'(1)=0,\quad 0\le w_v<1. \end{aligned}$$
(38)

We also deduce that system (36) is monotone cooperative and the sequence defined in (37) is monotone increasing. Hence system (36)–(37) is monotone cooperative as long as the initial conditions belong to [0, 1].

Space Implicit Model

Space homogeneous solutions of system (36)–(37), during the \(n+1\) inter-fire season, satisfy

$$\begin{aligned} \left\{ \begin{array}{lcl} \displaystyle \frac{d u_{n+1}}{d t} &{}=&{} u_{n+1}(1-u_{n+1}), \quad 0< t\le \tau ,\quad n\in \mathbb {N}^*,\\ \displaystyle \frac{d v_{n+1}}{d t} &{}=&{} - \lambda v_{n+1}(1-v_{n+1})+ \lambda \gamma (1-v_{n+1})u_{n+1}, \end{array} \right. \end{aligned}$$
(39)

together with the updating conditions

$$\begin{aligned} \left\{ \begin{array}{lcl} u_{n+1}(0) &{}=&{} (1- w_v(v_{n}(\tau ))\psi (u_n(\tau )))u_n(\tau ),\\ v_{n+1}(0) &{}=&{} (1- \eta )v_{n}(\tau )+\eta . \end{array} \right. \end{aligned}$$
(40)

Solving the logistic equation (39)\(_1\) leads to

$$\begin{aligned} u_{n+1}(t) = \displaystyle \frac{u_{n+1}(0)}{u_{n+1}(0)+(1-u_{n+1}(0))e^{-t}},\quad 0\le t\le \tau . \end{aligned}$$
(41)

In addition, direct computations give

$$\begin{aligned} \displaystyle \int _0^tu_{n+1}(s)ds= & {} \displaystyle \int _0^t\frac{u_{n+1}(0)}{u_{n+1}(0)+(1-u_{n+1}(0))e^{-s}}ds \nonumber \\= & {} \displaystyle \int _0^t\frac{u_{n+1}(0)e^{s}}{u_{n+1}(0)e^{s}+(1-u_{n+1}(0))}ds\nonumber \\= & {} \ln \left( 1+u_{n+1}(0)(e^{t}-1)\right) \nonumber \\=: & {} \ln I_u(t). \end{aligned}$$
(42)

Now, returning to (39)\(_2\) and setting \(z=1/(1-v_{n+1})\), we get

$$\begin{aligned} \dot{z} = \lambda (1-z+\gamma u_{n+1}z) = \lambda -\lambda (1-\gamma u_{n+1})z. \end{aligned}$$

Using the integrating factor \(e^{\lambda \int _0^t(1-\gamma u_{n+1}(s))ds } = e^{\lambda t} [I_u(t)]^{-\lambda \gamma }\), we get

$$\begin{aligned} z(t) = e^{-\lambda t} [I_u(t)]^{\lambda \gamma } \frac{1}{1-v_{n+1}(0)} + \lambda e^{-\lambda t} [I_u(t)]^{\lambda \gamma } \int _0^t e^{\lambda s} [I_u(s)]^{-\lambda \gamma }ds \end{aligned}$$

so that

$$\begin{aligned} v_{n+1}(t)=1-\displaystyle \frac{(1-v_{n+1}(0))e^{\lambda t}[I_u(t)]^{-\lambda \gamma }}{1+\lambda (1-v_{n+1}(0))\displaystyle \int _0^te^{\lambda s}[I_u(s)]^{-\lambda \gamma }ds}. \end{aligned}$$

Using the updating condition (40) leads to

$$\begin{aligned} \left\{ \begin{array}{lcl} u_{n+1}(\tau ) &{} = &{} \displaystyle \frac{(1-w_v(v_n(\tau ))\psi (u_n(\tau )))u_{n}(\tau )}{e^{-\tau }+(1-e^{-\tau })(1-w_v(v_n(\tau ))\psi (u_n(\tau )))u_{n}(\tau )}=:\bar{F}_1(u_n(\tau ),v_n(\tau )), \\ v_{n+1}(\tau ) &{}=&{} 1-\displaystyle \frac{(1-\eta )(1-v_{n}(\tau )))e^{\lambda \tau }[I_u(\tau )]^{-\lambda \gamma }}{1+(1-\eta )(1-v_{n}(\tau ))\lambda \displaystyle \int _0^\tau e^{\lambda s}[I_u(s)]^{-\lambda \gamma }ds}=:\bar{F}_2(u_n(\tau ),v_n(\tau )). \end{array} \right. \end{aligned}$$
(43)

Thus, the solution of system (39)–(40) given by system (43) generates a discrete dynamical system. Space homogeneous steady states of system (36)–(37) are steady states of model (43).

Space Homogeneous Steady States of System (36)–(37)

In this section we compute the space homogeneous steady states of system (36)–(37) by solving the fixed point problem associated to system (43).

\(\bar{F}_1(u,v)=u\) implies \(u=0\) or

$$\begin{aligned} 1-w_v(v)\psi (u)=e^{-\tau }+(1-e^{-\tau })(1-w_v(v)\psi (u))u. \end{aligned}$$
(44)

Similarly, \(\bar{F}_2(u,v)=v\) implies \(v=1\) or

$$\begin{aligned} 1+(1-\eta )(1-v)\lambda \displaystyle \int _0^\tau e^{\lambda s}[I_u(s)]^{-\lambda \gamma }ds=(1-\eta )e^{-\lambda \tau }[I_u(\tau )]^{\lambda \gamma }. \end{aligned}$$
(45)

We therefore deduce the first space homogeneous steady state \({\mathbf {E}}_0=(0,1)\). Substituting \(u=0\) (i.e. \(I_0(t)=1\)) in (45) implies

$$\begin{aligned} \bar{v}=\displaystyle \frac{\eta }{(1-\eta )(e^{\lambda \tau }-1)}>0. \end{aligned}$$

Note that

$$\begin{aligned} \bar{v}<1\quad \mathrm {if\;and\;only\;if} \quad \mathcal {R}_0>1, \end{aligned}$$

where \(\mathcal {R}_0:=(1-\eta )e^{\lambda \tau }\). Substituting \(v=1\) in (44) implies \(u=1\). Hence, we obtain Lemma 1.

Lemma 1

System (36)–(37) admits as trivial and semi-trivial space homogeneous steady states in the feasible region:

  • \({\mathbf {E}}_0=(0,1)\) and \({\mathbf {E}}_u=(1,1)\) that always exist;

  • \({\mathbf {E}}_v=(0,\bar{v})=\left( 0,\displaystyle \frac{\eta }{(1-\eta )(e^{\lambda \tau }-1)}\right)\) if and only if

    $$\begin{aligned} \mathcal {R}_0=(1-\eta )\exp (\lambda \tau )>1. \end{aligned}$$

In this study, we are mainly concerned with the existence of traveling wave solutions of system (36)–(37) involving steady states \({\mathbf {E}}_u\) and \({\mathbf {E}}_v\) computed in Lemma 1. In order to use the results of Li et al. [29], we translate the steady state \({\mathbf {E}}_v\) to 0 through another coordinates change.

The Existence of Traveling Waves

The Second Coordinate Change

Recall that

$$\begin{aligned} \bar{v}=\displaystyle \frac{\eta }{(1-\eta )(e^{\lambda \tau }-1)} \end{aligned}$$

and

$$\begin{aligned} \mathcal {R}_0=(1-\eta )e^{\lambda \tau }. \end{aligned}$$

Recall also that \(0\le \bar{v}<1\) is equivalent to \(\mathcal {R}_0>1.\) Therefore, in this section we assume that

$$\begin{aligned} \mathcal {R}_0>1. \end{aligned}$$

We set

$$\begin{aligned} \left\{ \begin{array}{ccl} u_n &{} = &{} u_n, \\ q_n &{} = &{} v_n-\bar{v}. \end{array} \right. \end{aligned}$$
(46)

Hence, system (36)–(37) becomes

$$\begin{aligned} \left\{ \begin{array}{lcl} \displaystyle \frac{\partial u_n}{\partial t} &{}=&{} u_n(1-u_n)+d_u\displaystyle \frac{\partial ^2 u_n}{\partial x^2}, \quad 0 < t\le \tau ,\quad x\in \mathbb {R},\\ \displaystyle \frac{\partial q_n}{\partial t} &{}=&{} - \lambda (q_n+\bar{v})(1-q_n-\bar{v}) + \lambda \gamma u_n(1-q_n-\bar{v})+{d_v\displaystyle \frac{\partial ^2 q_n}{\partial x^2}}, \end{array} \right. \end{aligned}$$
(47)

together with the updating conditions

$$\begin{aligned} \left\{ \begin{array}{lcl} u_{n+1}(x,0) &{}=&{} (1- w(q_{n}(x,\tau ))\psi (u_n(x,\tau )))u_n(x,\tau )=:H_1(u_n(x,\tau ), q_n(x,\tau )),\\ q_{n+1}(x,0) &{}=&{} (1- \eta )q_{n}(x,\tau )+\eta (1-\bar{v})=:H_2(u_n(x,\tau ), q_n(x,\tau )). \end{array} \right. \end{aligned}$$
(48)

As previously, we deduce from properties (38) that w is a decreasing \(\mathcal {C}^1({\mathbb {R}})\) function such that

$$\begin{aligned} w(0)>0,\quad w(1-\bar{v})=0,\quad w'(1-\bar{v})=0,\quad 0\le w(q)<1. \end{aligned}$$
(49)

For simplicity, we set \((u_{n,0},q_{n,0})=(u_{n}(x,0),q_{n}(x,0))\). We also set \({\mathbf {P}}_{n,0}=(u_{n,0},q_{n,0})\) and let \(Q_{\tau }\) denote the time-\(\tau\)-map operator solution of system (47). Then

$$\begin{aligned} {\mathbf {P}}_{n+1,0}={\mathbf {H}}(Q_{\tau }[{\mathbf {P}}_{n,0}])=:Q[{\mathbf {P}}_{n,0}], \end{aligned}$$
(50)

where \({\mathbf {H}}(u,q)=(H_1(u,q), H_2(u,q))\).

Using the coordinate change (46), system (43) becomes

$$\begin{aligned} \left\{ \begin{array}{lcl} u_{n+1}(\tau ) &{} = &{} \displaystyle \frac{(1-w(q_n(\tau ))\psi (u_n(\tau )))u_{n}(\tau )}{e^{-\tau }+(1-e^{-\tau })(1-w(q_n(\tau ))\psi (u_n(\tau )))u_{n}(\tau )}=:F_1(u_{n}(\tau ),q_{n}(\tau )), \\ q_{n+1}(\tau ) &{}=&{} 1-\bar{v}-\displaystyle \frac{(1-\eta )(1-\bar{v}-q_{n}(\tau ))e^{\lambda \tau }[I_u(\tau )]^{-\lambda \gamma }}{1+\lambda (1-\eta )(1-\bar{v}-q_{n}(\tau ))\displaystyle \int _0^\tau e^{\lambda s}[I_u(s)]^{-\lambda \gamma }ds}=:{F_2(u_{n}(\tau ),q_{n}(\tau ))} \end{array} \right. \end{aligned}$$
(51)

and, using Lemma 1, it is straightforward to deduce that system (47)–(48) admits as space homogeneous steady states

$$\begin{aligned} {\mathbf {e}}_0=(0,1-\bar{v}),\quad {\mathbf {e}}_u=(1,1-\bar{v})\quad \text{ and }\quad {\mathbf {e}}_v=(0,0). \end{aligned}$$

Stability Analysis of Space Homogeneous Steady States of System (47)–(48)

We first focus on the integral term that appears in \(F_2\) (see Eq. (51)\(_2\)). Recalling (42), we consider

$$\begin{aligned} Y(u,q):=\int _0^\tau e^{\lambda s}\left( 1+u(e^s-1)(1-w(q)\psi (u))\right) ^{-\lambda \gamma }ds \end{aligned}$$

and an auxiliary function \(Z:\mathbb {R}^+\times \mathbb {R}^+\times [0, \tau ]\rightarrow \mathbb {R}\) defined by

$$\begin{aligned} Z(u,q,s)=e^{\lambda s}\left( 1+u(e^s-1)(1-w(q)\psi (u))\right) ^{-\lambda \gamma }. \end{aligned}$$

For every \(u,q\in \mathbb {R}^+\), the function \(s\mapsto Z(u,q,s)\) is continuous on the interval \([0,\tau ]\). In addition,

$$\begin{aligned} \displaystyle \frac{\partial Z}{\partial u}=-\lambda \gamma e^{\lambda s}(e^s-1)\left( 1-w(q)(\psi (u)-u\psi '(u))\right) \left( 1+u(e^s-1)(1-w(q)\psi (u))\right) ^{-\lambda \gamma -1} \end{aligned}$$

exists and is continuous for all \((u,q,s)\in \mathbb {R}^+\times \mathbb {R}^+\times [0, \tau ]\). Consequently, \(\displaystyle \frac{\partial Y}{\partial u}=\displaystyle \int _0^\tau \frac{\partial Z}{\partial u}ds\). Similarly, \(\displaystyle \frac{\partial Y}{\partial q}=\displaystyle \int _0^\tau \frac{\partial Z}{\partial q}ds\). For convenience, we set (see (51))

$$\begin{aligned} F_1(u,q)=\displaystyle \frac{A_1(u,q)}{A_2(u,q)} \quad \text{ and } \quad F_2(u,q)=1-\bar{v}+\displaystyle \frac{B_1(u,q)}{B_2(u,q)}. \end{aligned}$$

Computing the partial derivatives of \(A_1, A_2,B_1\) and \(B_2\) defined in (51) gives

$$\begin{aligned} \displaystyle \frac{\partial A_1}{\partial u}= & {} 1-w(q)(\psi (u)+u\psi '(u)), \nonumber \\ \displaystyle \frac{\partial A_1}{\partial q}= & {} -w'(q)\psi (u)u,\nonumber \\ \displaystyle \frac{\partial A_2}{\partial u}= & {} (1-e^{-\tau })\displaystyle \frac{A_1}{\partial u},\nonumber \\ \displaystyle \frac{\partial A_2}{\partial q}= & {} (1-e^{-\tau })\displaystyle \frac{A_1}{\partial q},\nonumber \\ \displaystyle \frac{\partial B_1}{\partial u}= & {} \gamma \lambda (1-\eta )(1-q-\bar{v})e^{\lambda \tau }(1+(e^\tau -1)u(1-w(q)\psi (u)))^{-\gamma \lambda -1} (e^\tau -1)(1-w(q)(\psi (u)+u\psi '(u))),\nonumber \\ \displaystyle \frac{\partial B_1}{\partial q}= & {} (1-\eta )e^{\lambda \tau }(1+(e^\tau -1)u(1-w(q)\psi (u)))^{-\gamma \lambda } \nonumber \\&-\gamma \lambda (1-\eta )(1-q-\bar{v})e^{\lambda \tau }(e^\tau -1)uw'(q)\psi (u)(1+(e^\tau -1)u(1-w(q)\psi (u)))^{-\gamma \lambda -1},\nonumber \\ \displaystyle \frac{\partial B_2}{\partial u}= & {} -\gamma \lambda ^2 (1-\eta )(1-q-\bar{v})\nonumber \\&\times \displaystyle \int _0^\tau e^{\lambda s}(e^s-1)(1-w(q)(\psi (u)+u\psi '(u))(1+(e^s-1)u(1-w(q)\psi (u)))^{-\gamma \lambda -1}ds,\nonumber \\ \displaystyle \frac{\partial B_2}{\partial q}= & {} -(1-\eta )\lambda \displaystyle \int _0^\tau e^{\lambda s}(1+(e^s-1)u(1-w(q)\psi (u)))^{-\gamma \lambda }ds\nonumber \\&+(1-\eta )(1-q-\bar{v})\gamma \lambda ^2\displaystyle \int _0^\tau e^{\lambda s}(e^s-1)uw'(q)\psi (u) (1+(e^s-1)u(1-w(q)\psi (u)))^{-\gamma \lambda -1}ds. \end{aligned}$$
(52)

Let \({\mathcal {J}}= \{J_{ij}\}_{1\le i,j\le 2}\) denote the Jacobian matrix of (51). Using (52), the quotient rule and the properties of w (see (49)) and \(\psi\) (see (24)), we obtain the following results:

  • Local stability of \({\mathbf {e}}_0\). At the steady state \({\mathbf {e}}_0=(0,1-\bar{v})\) the matrix \({\mathcal {J}}\) has the following entries:

    $$\begin{aligned} J_{11}= & {} e^\tau , \\ J_{12}= & {} 0, \\ J_{21}= & {} 0, \\ J_{22}= & {} (1-\eta )e^{\lambda \tau }. \\ \end{aligned}$$

    Indeed, for example,

    $$\begin{aligned} J_{12}=\displaystyle \frac{\partial F_1(0,1-\bar{v})}{\partial q}=\left. \displaystyle \frac{\displaystyle \frac{\partial A_1}{\partial q}A_2-A_1\displaystyle \frac{\partial A_2}{\partial q}}{(A_2)^2}\right| _{(u,q)=(0,1-\bar{v})}=0, \end{aligned}$$

    because

    $$\begin{aligned} \displaystyle \frac{\partial A_1(0,1-\bar{v})}{\partial q}=0 \end{aligned}$$

    thanks to \(w'(1-\bar{v})=0\) (see (49)). Since \(e^\tau\) is an eigenvalue of \({\mathcal {J}}\) at \({\mathbf {e}}_0\) and \(e^\tau >1,\) the steady state \({\mathbf {e}}_0\) is unstable.

  • Local stability of \({\mathbf {e}}_u\). At the steady state \({\mathbf {e}}_u=(1,1-\bar{v})\) the matrix \({\mathcal {J}}\) has the following entries:

    $$\begin{aligned} J_{11}= & {} e^{-\tau }, \\ J_{12}= & {} 0,\\ J_{21}= & {} 0, \\ J_{22}= & {} (1-\eta )e^{\lambda \tau (1-\gamma )}=:\mathcal {R}_1. \\ \end{aligned}$$

    Eigenvalues of the Jacobian matrix at \({\mathbf {e}}_u\) are \(e^{-\tau }\) and \(\mathcal {R}_1\) with \(e^{-\tau }<1\). Therefore, \({\mathbf {e}}_u\) is locally asymptotically stable (LAS) whenever \(\mathcal {R}_1<1\).

  • Local stability of \({\mathbf {e}}_v\). \({\mathcal {J}}\) at the steady state \({\mathbf {e}}_v=(0,0)\) has the following entries:

    $$\begin{aligned} J_{11}= & {} \left( 1-w(0)\psi (0)\right) e^\tau , \\ J_{12}= & {} 0, \\ J_{22}= & {} \displaystyle \frac{1}{(\mathcal {R}_0)^2}.\\ \end{aligned}$$

    Knowing the explicit value of \(J_{21}\) is not necessary since \(J_{12}=0\). Eigenvalues of the Jacobian matrix at \({\mathbf {e}}_v\) are \(J_{11}\) and \(J_{22}\). Recall that we assumed \(\mathcal {R}_0>1\). Hence, \(J_{22}<1\). Therefore, the steady state \({\mathbf {e}}_v\) is LAS whenever

    $$\begin{aligned} \mathcal {R}_2=\left( 1-w(0)\psi (0)\right) e^\tau <1. \end{aligned}$$

Hence, the following lemma holds true.

Lemma 2

The space homogeneous steady states of system (47)–(48) have the following stability properties.

  1. 1.

    The steady state \({\mathbf {e}}_0=(0,1-\bar{v})\) is unstable.

  2. 2.

    The steady state \({\mathbf {e}}_u=(1,1-\bar{v})\) is LAS whenever \(\mathcal {R}_1<1\).

  3. 3.

    The steady state \({\mathbf {e}}_v=(0,0)\) is LAS whenever \(\mathcal {R}_2<1\).

Application of the Results of [29]

In the sequel, we study the recursion operator Q defined in Eq. (50) and we check if it satisfies Hypotheses 2.1 of Li et al. [29]. Recall that \({\mathbf {P}}_{n,0}=(u_{n,0},v_{n,0})\) and \(Q_{\tau }\) is the time-\(\tau\)-map solution operator of reaction–diffusion system (47). We consider the order interval \(\mathcal {C}_{{\mathbf {e}}_{u}}=[{\mathbf {e}}_v, {\mathbf {e}}_{u}]\), where \({\mathbf {e}}_v=\mathbf{0 }\) and \({\mathbf {e}}_{u}=(1,1-\bar{v})\) is the positive coexistence steady state defined in Lemma 2.

Lemma 3

(Some properties of \(Q_\tau\))

  1. 1.

    The operator \(Q_\tau\) is order preserving in the sense that if \(\mathbf{u }\) and \(\mathbf{v}\) are any two functions in \(\mathcal {C}_{{\mathbf {e}}_{u}}\) with \(\mathbf{v} \ge \mathbf{u }\), then \(Q_\tau [\mathbf{v} ]\ge Q_\tau [\mathbf{u }]\).

  2. 2.

    \(Q_\tau\) is translation invariant.

  3. 3.

    For any \(\mathbf{v} , \mathbf{u }\in \mathcal {C}_{{\mathbf {e}}_u}\) and fixed x, \(|Q_\tau [\mathbf{v} ](x)-Q_\tau [\mathbf{u }](x)|\) is arbitrarily small, provided \(|\mathbf{v} (y)-\mathbf{u }(y)|\) is sufficiently small on a sufficiently long interval centered at x.

  4. 4.

    Every sequence \(\mathbf{v} _n(x)\) in \(\mathcal {C}_{{\mathbf {e}}_{u}}\) has a subsequence \(\mathbf{v} _{n_l}\) such that \(Q_\tau [\mathbf{v} _{n_l}]\) converges uniformly on every bounded set.

Proof

  1. 1.

    The reaction–diffusion system (47) is a cooperative system. Hence, following the analysis at the beginning of this section, we deduce that the time-\(\tau\)-map solution operator of system (47) is order preserving.

  2. 2.

    Let \({\mathbf {u}}\) be the solution of system (47) initiated at \({\mathbf {u}}_0\). For \(y\in {\mathbb {R}}\), we set \({\mathbf {v}}=T_y[{\mathbf {u}}]\), where, as defined earlier, \(T_y\) is the translation operator. In particular \({\mathbf {v}}_0=T_y[{\mathbf {u}}_0]\). We have \({\mathbf {v}}_t=(T_y[{\mathbf {u}}])_t=T_y[{\mathbf {u}}_t]\), \(A{\mathbf {v}}=AT_y[{\mathbf {u}}]=T_y[A{\mathbf {u}}]\) and \({\mathbf {F}}({\mathbf {v}})={\mathbf {F}}(T_y[{\mathbf {u}}])=T_y[F({\mathbf {u}})]\) since \({\mathbf {F}}\) does not explicitly depend on \(x\in {\mathbb {R}}\). Therefore, \({\mathbf {v}}_t+A{\mathbf {v}}-{\mathbf {F}}({\mathbf {v}})=T_y[{\mathbf {u}}_t+A{\mathbf {u}}-{\mathbf {F}}({\mathbf {u}})]=0\) and, by the uniqueness of solutions, we have

    $$\begin{aligned} T_y[Q_\tau [{\mathbf {u}}_0]](x) = T_y[{\mathbf {u}}](x) = {\mathbf {v}}(x) = Q_\tau [T_y[{\mathbf {u}}_0]](x). \end{aligned}$$

    Hence the time-\(\tau\)-map solution operator of system (47), \(Q_\tau\), is translation invariant.

    To prove 3. and 4. we write the solution, see e.g. [11, page 95], as

    $$\begin{aligned} \left\{ \begin{array}{l} u=G_{d_u}\star u_0+G_{d_u}\star \star f_{d_u}(u,q),\\ q=G_{d_v}\star q_0+G_{d_v}\star \star f_{d_v}(u,q),\\ \end{array} \right. \end{aligned}$$
    (53)

    where

    $$\begin{aligned} f_{d_u}(u,q)=u(1-u), \qquad f_{d_v}(u,q)=-\lambda (q+\bar{v})(1-q-\bar{v})+\lambda \gamma u(1-q-\bar{v}), \end{aligned}$$

    \(G_d\), \(d=d_u,d_v\), as well as the convolution \(\star\), were defined in (28) and the spatio-temporal convolution is given by

    $$\begin{aligned} G_d\star \star f_d =\displaystyle \int _{0}^t\frac{1}{\sqrt{4\pi d (t-s)}}\int _{{\mathbb {R}}}\exp \left( -\displaystyle \frac{(x-y)^2}{4d (t-s)}\right) f_d(u(y,s),q(y,s))dyds. \end{aligned}$$
  3. 3.

    Let (uq) and (vp) be two solutions of system (47) initiated at \((u_0,q_0)\) and \((v_0,p_0)\) respectively. We assume that \((u_0,q_0)\), \((v_0,p_0)\in \mathcal {C}_{{\mathbf {e}}_{u}}\), hence (uq) and (vp) are also uniformly bounded. Let \(\{S_d^c(t)\}_{t\ge 0}\) denotes the (positive) semigroup solving (29) for some function c satisfying \(|c(x,t)|\le L\). Using \(|S_d(t)u_0|\le S_d(t)|u_0|, u_0\in C_{ub}({\mathbb {R}}),\) where \(\{S_d(t)\}_{t\ge 0}\) is the diffusion semigroup (28), and the Phillips-Dyson expansion to (29) we ascertain that for any \(u_0\)

    $$\begin{aligned} |S^c_d(t)u_0| \le S^L_d(t)|u_0| = e^{Lt}S_d(t)|u_0|. \end{aligned}$$

    As before, we begin with solutions u and v to (47)\(_1\). Repeating the argument leading to (32), we see that \(z(x,t) = u(x,t)-v(x,t)\) can be estimated as

    $$\begin{aligned} |z(x,t)| \le e^{Lt}[S_{d_u}(t)|u_0-v_0|](x) + \displaystyle \frac{e^{Lt}}{\sqrt{4\pi d_u t}}\displaystyle \int _{{\mathbb {R}}}\exp \left( -\displaystyle \frac{(x-y)^2}{{4d_u t}}\right) |u_0(y)-v_0(y)| dy. \end{aligned}$$

    Let, for \(\epsilon >0\), \(r>0\) be such that

    $$\begin{aligned} \left( \int _{-\infty }^\frac{-r}{2\sqrt{d_u\tau }} + \int ^{\infty }_\frac{r}{2\sqrt{d_u\tau }}\right) e^{-z^2}dz \le \frac{\epsilon \sqrt{\pi }}{4e^{L\tau }}. \end{aligned}$$

    Then let us fix x and let \(|u_0(x)-v_0(x)| \le \delta \le \epsilon /2e^{L\tau }\) on \((x-r,x+r)\) so that we obtain for \(0<t\le \tau\)

    $$\begin{aligned} |z(x,t)|\le & {} \displaystyle \frac{e^{Lt}}{\sqrt{4\pi d_u t}}\displaystyle \int _{x-r}^{x+r}\exp \left( -\displaystyle \frac{(x-y)^2}{{4d_u t}}\right) |u_0(y)-v_0(y)| dy \nonumber \\&+ \displaystyle \frac{e^{Lt}}{\sqrt{4\pi d_u t}}\displaystyle \left( \int _{-\infty }^{x-r}+\int _{x+r}^\infty \right) \exp \left( -\displaystyle \frac{(x-y)^2}{{4d_u t}}\right) |u_0(y)-v_0(y)| dy \nonumber \\\le & {} e^{Lt} \delta + \displaystyle \frac{2 e^{Lt}}{\sqrt{\pi }}\displaystyle \left( \int _{-\infty }^\frac{-r}{2\sqrt{d_ut}}+\int _\frac{r}{2\sqrt{d_ut}}^\infty \right) e^{-z^2} dz\le e^{L\tau } \delta + \displaystyle \frac{2 e^{L\tau }}{\sqrt{\pi }}\displaystyle \left( \int _{-\infty }^\frac{-r}{2\sqrt{d_u\tau }}+\int _\frac{r}{2\sqrt{d_u\tau }}^\infty \right) e^{-z^2} dz\le \epsilon . \end{aligned}$$
    (54)

    By choosing appropriate r, we see that the estimate is valid for x in any given bounded subset of \({\mathbb {R}}\).

    Considering now (47)\(_2\), we see that \(Z(x,t) = q(x,t)-p(x,t)\) is a solution to

    $$\begin{aligned} Z_t= & {} d_v Z_{xx} + (q-qu-q^2-p+pv + p^2) = d_v Z_{xx} + Z(1-(p+q) - u) - p(u-v),\nonumber \\ Z(x,0)= & {} q(x,0)-p(x,0) =: Z_0(x,0) \end{aligned}$$
    (55)

    and considerations as above show that \(|Z(x,t)| \le e^{L_1t}\varPsi (x,t)\), where

    $$\begin{aligned} \varPsi _t= & {} d_v \varPsi _{xx}+ e^{-L_1 t}z,\nonumber \\ \varPsi (x,0)= & {} |Z_0(x,0)|. \end{aligned}$$
    (56)

    In the above, \(L_1\) is a constant bounding \(|1-(p+q)-u|, 0\le p,q,u \le 1\), and we used \(0\le v\le 1\). Hence

    $$\begin{aligned} |Z(x,\tau )|\le & {} \displaystyle \frac{e^{L_1\tau }}{\sqrt{4\pi d_v \tau }}\displaystyle \int _{{\mathbb {R}}}\exp \left( -\displaystyle \frac{(x-y)^2}{{4d_v \tau }}\right) |q_0(y)-p_0(y)| dy\\&+ \displaystyle \int _{0}^\tau \frac{e^{L_1(\tau -s)}}{\sqrt{4\pi d_v (\tau -s)}}\int _{{\mathbb {R}}}\exp \left( -\displaystyle \frac{(x-y)^2}{4d_v (\tau -s)}\right) z(y,s)dyds \end{aligned}$$

    and the estimates follow as above where, in the second term, we use the fact that (54) is uniform on \([0,\tau ]\) and any bounded subset of \({\mathbb {R}}\).

  4. 4.

    For each \(t\in (0, \tau ]\), the functions \(Q_t[{\mathbf {w}}_0]\) with \({\mathbf {w}}_0=(u_0,q_0)\in \mathcal {C}_{{\mathbf {e}}_{u}}\) form an equicontinuous family. Indeed, for \(0<t\le \tau\), \({\mathbf {w}}_0\in \mathcal {C}_{{\mathbf {e}}_{u}}\) and \(x\in {\mathbb {R}}\), \(Q_{t}[{\mathbf {w}}_0(x)]=:{\mathbf {w}}(t,x)=(u(t,x),q(t,x))\), following (53) and by using the property of the spatial convolution, we obtain

    $$\begin{aligned} \left\{ \begin{array}{l} \displaystyle \frac{\partial u}{\partial x}=\displaystyle \frac{\partial G_{d_u}}{\partial x}\star u_0+\displaystyle \frac{\partial G_{d_u}}{\partial x}\star \star f_{d_u}(u,q),\\ \displaystyle \frac{\partial q}{\partial x}=\displaystyle \frac{\partial G_{d_v}}{\partial x}\star q_0+\displaystyle \frac{\partial G_{d_v}}{\partial x}\star \star f_{d_v}(u,q).\\ \end{array} \right. \end{aligned}$$
    (57)

    Since \((u_0, q_0)\in \mathcal {C}_{{\mathbf {e}}_{u}}\) i.e. \(0\le u_0\le 1\) and \(0\le q_0\le 1-\bar{v}\), we have \(\Vert f_{d_u}(u,q)\Vert _{\infty }\le 1\) and \(\Vert f_{d_v}(u,q)\Vert _{\infty }\le M\) for some M. In addition, by direct calculation or, more generally, by [18, Theorem 11], for \(d=d_u,d_v\), there exist positive constants \(\alpha _d\) and \(\beta _d\) such that for \(t>0\)

    $$\begin{aligned} \left| \displaystyle \frac{\partial G_d(x,t)}{\partial x}\right| \le \frac{\alpha _de^{-\beta _d\frac{x^2}{t}}}{t} \end{aligned}$$

    Hence,

    $$\begin{aligned} \left\{ \begin{array}{l} \displaystyle \left| \frac{\partial u}{\partial x}(x,t)\right| \le \displaystyle \int _{{\mathbb {R}}}\frac{\alpha _{d_u}e^{-\beta _{d_u}\frac{(x-y)^2}{t}}}{t}dy+\int _0^t\int _{{\mathbb {R}}}\frac{\alpha _{d_u}e^{-\beta _{d_u}\frac{(x-y)^2}{t-s}}}{t-s}dyds,\\ \displaystyle \left| \frac{\partial q}{\partial x}(x,t)\right| \le \displaystyle (1-\bar{v})\int _{{\mathbb {R}}}\frac{\alpha _{d_v}e^{-\beta _{d_v}\frac{(x-y)^2}{t}}}{t}dy+M\int _0^t\int _{{\mathbb {R}}}\frac{\alpha _{d_v}e^{-\beta _{d_v}\frac{(x-y)^2}{t-s}}}{t-s}dyds. \end{array} \right. \end{aligned}$$
    (58)

    Evaluating the integrals in (58) we obtain

    $$\begin{aligned} \left\{ \begin{array}{l} \left| \displaystyle \frac{\partial u}{\partial x}(x,t)\right| \le \alpha _{d_u}\sqrt{\displaystyle \frac{\pi }{\beta _{d_u}t}}+\alpha _{d_u}\sqrt{\displaystyle \frac{\pi t}{\beta _{d_u}}}=:\delta _{d_u}(t),\\ \left| \displaystyle \frac{\partial q}{\partial x}(x,t)\right| \le (1-\bar{v})\alpha _{d_v}\sqrt{\displaystyle \frac{\pi }{\beta _{d_v}t}}+M\alpha _{d_v}\sqrt{\displaystyle \frac{\pi t}{\beta _{d_v}}}=:\delta _{d_v}(t), \end{array} \right. \end{aligned}$$
    (59)

    where \(\delta _{d_u}(t)\) and \(\delta _{d_v}(t)\) do not depend on \(u_0\), \(q_0\) and \(x\in {\mathbb {R}}\). Thus the first spatial derivative of the solution \({\mathbf {w}}=(u,q)\) is uniformly bounded. Hence, using the Mean Value Theorem, we deduce that the family of solutions of system (47) is equicontinuous. Then, part 4 of Lemma 3 follows from the Arzela-Ascoli’s theorem, see e.g. [43, Corollary 41] (i.e. any bounded and equicontinuous sequence of continuous functions on a separable metric space contains a uniformly convergent subsequence on every bounded subset).

\(\square\)

In order to obtain properties of the operator Q defined in (50), we first formulate results for the nonlinear operator \({\mathbf {H}}: B\rightarrow B\). Recall that the Banach space considered here is \(B=C_{ub}({\mathbb {R}})\times C_{ub}({\mathbb {R}})\) endowed with the sup-norm. For \(\mathbf{u }=(u,q)\in B\),

$$\begin{aligned} {\mathbf {H}}(\mathbf{u })=((1-w(q)\psi (u))u; \;(1-\eta )q+\eta (1-\bar{v})). \end{aligned}$$

Lemma 4

(Some properties of \({\mathbf {H}}\))

  1. 1.

    The nonlinear operator \({\mathbf {H}}\) is order preserving in the sense that if \(\mathbf{u }\) and \(\mathbf{v}\) are any two functions with \(\mathbf{v} \ge \mathbf{u }\), then \({\mathbf {H}}(\mathbf{v} )\ge {\mathbf {H}}(\mathbf{u })\).

  2. 2.

    \({\mathbf {H}}\) is translation invariant.

  3. 3.

    For any two functions \(\mathbf{u }\) and \(\mathbf{v}\), \(\Vert {\mathbf {H}}(\mathbf{v} )-{\mathbf {H}}(\mathbf{u })\Vert _B \le C \Vert \mathbf{v} -\mathbf{u }\Vert _B\), where \(C\in {\mathbb {R}}\) and depends on \(\Vert \mathbf{u }\Vert _B\), \(\Vert \mathbf{v} \Vert _B\).

  4. 4.

    If a sequence \(\mathbf{v} _n(x)\) converges uniformly on every bounded set, then \({\mathbf {H}}(\mathbf{v} _n(x))\) also has the same property.

Proof

Let \(\mathbf{u }=(u,q)\) and \(\mathbf{v} =(v,p)\) be such that \(\mathbf{u }\le \mathbf{v}\). Hence \((1-w(q)\psi (u))u\le (1-w(p)\psi (v))v\) since w (resp. \(\psi\)) is increasing (resp. decreasing) and \((1-\eta )q+\eta (1-\bar{v})\le (1-\eta )p+\eta (1-\bar{v})\). Thus \({\mathbf {H}}(\mathbf{u })\le H(\mathbf{v} )\) and part 1 of Lemma 4 holds. Since \({\mathbf {H}}\) does not explicitly depend on \(x\in {\mathbb {R}}\), then \({\mathbf {H}}\) is translation invariant and part 2 of Lemma 4 is valid. Part 3 follows from the local Lipschitz property of \({\mathbf {H}},\) while part 4 follows from the continuity of \({\mathbf {H}}\). This ends the proof. \(\square\)

Combining Lemmas 3 and 4, we deduce the following result for the recursion operator \(Q:={\mathbf {H}}\circ Q_\tau\), defined in (50).

Lemma 5

(Some properties of Q)

  1. 1.

    The operator Q is order preserving in the sense that if \(\mathbf{u }\) and \(\mathbf{v}\) are any two functions in \(\mathcal {C}_{{\mathbf {e}}_{u}}\) with \(\mathbf{v} \ge \mathbf{u }\), then \(Q[\mathbf{v} ]\ge Q[\mathbf{u }]\).

  2. 2.

    Q is translation invariant.

  3. 3.

    For any fixed x, \(|Q[\mathbf{v} ](x)-Q[\mathbf{u }](x)|\) is arbitrarily small, provided \(|\mathbf{v} (y)-\mathbf{u }(y)|\) is sufficiently small on a sufficiently long interval centered at x.

  4. 4.

    Every sequence \(\mathbf{v} _n(x)\) in \(\mathcal {C}_{{\mathbf {e}}_{u}}\) has a subsequence \(\mathbf{v} _{n_l}\) such that \(Q[\mathbf{v} _{n_l}]\) converges uniformly on every bounded set.

In the sequel, we assume that \(\mathcal {R}_0=(1-\eta )e^{\lambda \tau }>1\) and \(\mathcal {R}_1=(1-\eta )e^{\lambda \tau (1-\gamma )}<1\); that is, the coexistence steady state \({\mathbf {e}}_{u}=(1,1-\bar{v})\) exists and is stable. We also assume that \(\mathcal {R}_2=\left( 1-w(0)\psi (0)\right) e^\tau >1\); that is, \({\mathbf {e}}_v=(0,0)\) is unstable.

Taking into account Lemma 5, we deduce that the recursion operator Q defined in (50) verifies all conditions of Hypothesis 2.1. Consequently, we can apply the results of Li et al. [29] that deal with the spreading speeds and existence of traveling wave solutions for systems (47)–(48). Recall that a traveling wave of speed c is a solution of the recursion (50) which has the form \(\mathbf{u }_n(x,0)=\mathbf{Z} (x-nc)\) with \(\mathbf{Z} (s)\) being a function in \(\mathcal {C}_{{\mathbf {e}}_{u}}\). That is, the solution at time \(n + 1\) is simply the translate by c of its value at n. Using the definition of \(c^*\) (see (12)) and \(c^*_f\) (see (15)), the following result holds true.

Theorem 4

(Spreading speeds and traveling waves)

  1. (a).

    Slowest spreading speed: There is an index \(j\in \{1,2\}\) for which the following statement is true: Suppose that the initial function \(\mathbf{u }_0(x)\) is 0 for all sufficiently large x, and that there are positive constants \(0<\rho \le \sigma <1\) such that \({\mathbf{0 }}\le \mathbf{u }_0\le \sigma {\mathbf {e}}_{u}\) for all x and \(\mathbf{u }_0\ge \rho {\mathbf {e}}_{u}\) for all sufficiently negative x. Then for any positive \(\varepsilon\) the solution \(\mathbf{u }_n\) of the recursion (50) has the properties

    $$\begin{aligned} \lim \limits _{n\rightarrow +\infty }\left[ \sup \limits _{x\ge n(c^*+\varepsilon )}\{\mathbf{u }_n\}_j(x)\right] =0 \end{aligned}$$
    (60)

    and

    $$\begin{aligned} \lim \limits _{n\rightarrow +\infty }\left[ \sup \limits _{x\le n(c^*-\varepsilon )}\{{\mathbf {e}}_{u}-\mathbf{u }_n(x)\}\right] ={\mathbf{0 }}; \end{aligned}$$
    (61)

    that is, the j th component spreads at a speed no higher than \(c^*\), and no component spreads at a lower speed.

  2. (b).

    Fastest spreading speed: There is an index \(i\in \{1,2\}\) for which the following statement is true: Suppose that the initial function \(\mathbf{u }_0(x)\) is 0 for all sufficiently large x, and that there are positive constants \(0<\rho \le \sigma <1\) such that \({\mathbf{0 }}\le \mathbf{u }_0\le \sigma {\mathbf {e}}_{u}\) for all x and \(\mathbf{u }_0\ge \rho {\mathbf {e}}_{u}\) for all sufficiently negative x. Then for any positive \(\varepsilon\) the solution \(\mathbf{u }_n\) of the recursion (50) has the properties

    $$\begin{aligned} \limsup \limits _{n\rightarrow +\infty }\left[ \inf \limits _{x\le n(c^*_f-\varepsilon )}\{\mathbf{u }_n\}_i(x)\right] >0 \end{aligned}$$
    (62)

    and

    $$\begin{aligned} \lim \limits _{n\rightarrow +\infty }\left[ \sup \limits _{x\ge n(c^*_f+\varepsilon )}\mathbf{u }_n(x)\right] ={\mathbf{0 }}; \end{aligned}$$
    (63)

    that is, the i th component spreads at a speed no less than \(c^*_f\), and no component spreads at a higher speed.

  3. (c).

    Monostable traveling wave: If \(c\ge c^*\), there is a non-increasing monostable traveling wave solution \(\mathbf{Z} (x-nc)\) of speed c with \(\mathbf{Z} (-\infty )={\mathbf {e}}_{u}\) and \(\mathbf{Z} (+\infty )\) a steady state other than \({\mathbf {e}}_{u}\).

    If there is a traveling wave \(\mathbf{Z} (x-nc)\) with \(\mathbf{Z} (-\infty )=e_{u}\) such that for at least one component \(i\in \{1,2\}\)

    $$\begin{aligned} \liminf \limits _{x\rightarrow \infty }Z_i(x)=0, \end{aligned}$$

    then \(c\ge c^*\). If this property is valid for all components of \(\mathbf{Z}\), then \(c\ge c^*_f\).

Proof

Since the recursion operator Q defined in (50) verifies Hypothesis 2.1., the proof of Theorem 4 follows directly from Theorems 2.1, 2.2 and 3.1 of Li et al. [29]. \(\square\)

Let us point out that if, instead of the first coordinates change (35), we considered

$$\begin{aligned} \left\{ \begin{array}{ccl} u_n &{} = &{} 1-U_n, \\ v_n &{} = &{} V_n, \end{array} \right. \end{aligned}$$
(64)

then we would obtain a monotone increasing system (see Appendix B). Hence, by reasoning as before, one can study the case where the steady state \(\mathbf{0 }\) is stable and the steady state \({\mathbf {e}}_{u}\) is unstable. However, the bistable case, i.e. when both 0 and \({\mathbf {e}}_{u}\) are simultaneously stable, remains an open problem.

Numerical Simulations

In this section we provide numerical simulations of the impulsive tree–grass reaction–diffusion model (3)–(4). We note that the parameters with \(\tilde{}\) below refer to this system and can be derived from the corresponding parameters related to the normalized system (47)–(48), see Appendix A. Thus, we consider fire events as periodic and pulse perturbations with the time period \(\tilde{\tau }\). The form of the functions \(\gamma _{G}(\mathbf{W })\), \(\gamma _{T}(\mathbf{W })\), \(K_{G}(\mathbf{W })\), \(K_{T}(\mathbf{W })\), \(\psi\) and \(w_G\) is considered following Yatat Djeumen et al. [63]. The readers are referred to Appendix A for their definition and parametrization. The parameter values used in the following simulations are also given in Appendix A: see Tables 1 and 2.

Fig. 1
figure 1

Illustration of the spreading of both tree (see panel a) and grass (see panel b) biomasses toward the stable forest steady state \({\mathbf {E}}_T=(273.3955,0)\) with system (3)–(4). \(\mathbf{W }=1200\) mm.year\(^{-1}\), \(\tilde{\tau }=2\) year, \(d_T\)=0.001 and \(d_G\)=0.002. Remaining parameters are in Table 1

Using the parameters values given in Table 1, Fig. 1 depicts the spreading of tree and grass biomasses toward the stable forest homogeneous steady state \({\mathbf {E}}_T=(273.3955,0)\). In this case, \(\tilde{\mathcal {R}}_0=1.3667\), \(\tilde{\mathcal {R}}_1=0.0058\) and \(\tilde{\mathcal {R}}_2=3.3801\). Recall that \({\mathbf {E}}_T\) is LAS whenever \(\tilde{\mathcal {R}}_1<1\), while the grassland homogeneous steady state \({\mathbf {E}}_G=(0,3.3888)\) exists when \(\tilde{\mathcal {R}}_0>1\) and is LAS whenever \(\tilde{\mathcal {R}}_2<1\). In terms of tree–grass interactions, Fig. 1 illustrates the spreading of forest or the so-called ’forest encroachment’ phenomenon (Yatat Djeumen et al. [59]).

In the setting of the forest encroachment phenomenon, we carry out numerical simulations to compute the spreading speed of forest biomass. We investigate the relationship between the tree biomass diffusion coefficient and its spreading speed. To estimate the spreading speed of the tree biomass that undergoes a forest encroachment, grass biomass diffusion coefficient \(d_G\) is kept constant and equal to 0.002 while tree biomass diffusion coefficient \(d_T\) varies in the range [0.001, 0.9]. In the diffusive logistic equation, a linear relationship is obtained between the wave speed and the square root of the diffusion coefficient (e.g. Volpert [52], Yatat Djeumen et al. [59], Yatat Djeumen and Dumont [61]). Hence, for the tree biomass, we consider an equation of the form \(c_T(d_T)=a_1d_T^{a_2}\) to be fitted for the data shown in Fig. 2(b), where \(c_T\in [0.0865, 1.6899]\). We found that \(a_1\in (1.7624, 1.7861)\) and \(a_2\in (0.4816, 0.4911)\) with 95% confidence. In fact, \(a_1=1.7743\) and \(a_2=0.48634\), with \(r^2=1\), indicating that 100% of the variance of the data is explained by the equation.

With the parameter values given in Table 2, Fig. 3 illustrates the spreading of both tree and grass biomasses toward the grassland homogeneous steady state \({\mathbf {E}}_G=(0, 2.1096)\). In this case, \(\tilde{\mathcal {R}}_0=2.0425\), \(\tilde{\mathcal {R}}_1=1.2541\) and \(\tilde{\mathcal {R}}_2=0.9932\). Recall that \({\mathbf {E}}_G\) exists when \(\tilde{\mathcal {R}}_0>1\) and is LAS whenever \(\tilde{\mathcal {R}}_2<1\). We further investigate the relationship between the diffusion coefficient and the spreading speed of the grass biomass. We assume that \(d_T=0.001\) and \(d_G\in [0.001, 1]\). Motivated by the linear relationship obtained between the wave speed and the square root of the diffusion coefficient in the diffusive logistic equation, an equation like \(c_G(d_G)=a_1d_G^{a_2}\) was fitted to the data shown in Fig. 4b, where \(c_G\in [0.0465, 1.0877]\). We found that \(a_1\in (1.0829, 1.0935)\) and \(a_2\in (0.4839, 0.4915)\) with 95% confidence. In fact \(a_1=1.0882\) and \(a_2=0.48769\), with \(r^2=1\), indicating that 100% of the variance of the data is explained by the equation.

Fig. 2
figure 2

The forest homogeneous steady state \({\mathbf {E}}_T=(273.3955,0)\) is stable, while the grassland homogeneous steady state \({\mathbf {E}}_G=(0,3.3888)\) is unstable. We illustrate the forest encroachment phenomenon (see panel a) with \(\mathbf{W }=1200\) mm.year\(^{-1}\), \(\tilde{\tau }=2\) year, \(d_{T}\)=0.001, \(d_G\)=0.002, and we fit the spreading speeds for the tree biomass (see panel b). Remaining parameters are in Table 1

Fig. 3
figure 3

Illustration of the spreading of both tree (panel a) and grass (panel b) biomasses toward the grassland steady state \({\mathbf {E}}_G=(0, 2.1096)\). \(\mathbf{W }=450\) mm.year\(^{-1}\), \(\tilde{\tau }=2\) year, \(d_{T}\)=0.001, \(d_G\)=0.002. Remaining parameters are in Table 2

Fig. 4
figure 4

The forest homogeneous steady state \({\mathbf {E}}_T=(24.3905,0)\) is unstable while the grassland homogeneous steady state \({\mathbf {E}}_G=(0, 2.1096)\) is stable. We illustrate the grassland encroachment phenomenon (see panel a) and fit the spreading speeds for the grass biomass (see panel b). In panel a, \(\mathbf{W }=450\) mm.year\(^{-1}\), \(\tilde{\tau }=2\) year, \(d_{T}\)=0.001, \(d_G\)=0.002. Remaining parameters are in Table 2

Conclusion

In this paper, we used the vector-valued recursion equation theory (e.g. Weinberger et al. [56], Weinberger [54, 55], Lewis et al. [28], Li et al. [29]) to propose a framework in which we can systematically address the questions of the existence of traveling waves for monotone systems of impulsive reaction–diffusion equations and the computation of spreading speeds. This study extends the previous one that dealt with the impulsive Fisher-Kolmogorov-Petrowsky-Piscounov (FKPP) equation (Yatat Djeumen and Dumont [61]) but is restricted to monostable situations; that is, when only one of the steady states is stable. However, the bistable case (i.e. when two steady states are simultaneously stable) is also meaningful and needs to be studied for monotone systems of impulsive reaction–diffusion equations. Traveling waves in bistable reaction–diffusion systems without impulsive perturbations are treated in Volpert [52] (see also Yatat Djeumen et al. [60] for application in the context of bistable tree–grass reaction–diffusion model).

The computation of spreading speeds and the existence of traveling waves for bistable monotone systems of impulsive reaction–diffusion equations will be the aim of future studies. It first requires to elaborate a recursion equations theory that includes bistable cases and thus extending the results of Li et al. [29] that also deal only with monostable cases. Last but not least, in a recent work (Yatat Djeumen et al. [62]) we studied a space-implicit model that gives quite satisfactory results, in term of the model’s predictions, across the rainfall gradient and improving already published results obtained with tree–grass interactions models of similar complexity or less. Hence, it is desirable to also study its space-explicit counterpart acknowledging impulsive fire events. The aim will be to analyse the impact of impulsive fires on the dynamics of the vegetation mosaics like forest-grassland and forest-savanna in humid tropical savannas (see also Yatat Djeumen et al. [59]).