# Front propagation into unstable states for periodic monotone reaction-diffusion systems

Liangliang Deng, Arnaud Ducrot, Quentin Griette

<figure><img src="image_1.png" /></figure>

## To cite this version:

Liangliang Deng, Arnaud Ducrot, Quentin Griette. Front propagation into unstable states for periodic monotone reaction-diffusion systems. 2025. hal-04991522v2

HAL Id: hal-04991522

https://hal.science/hal-04991522v2

Preprint submitted on 13 Nov 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Front propagation into unstable states for periodic monotone reaction-diffusion systems

Liangliang Deng<sup>a,*</sup>, Arnaud Ducrot<sup>b</sup>, and Quentin Griette<sup>b</sup>

<sup>a</sup> School of Science and Engineering, The Chinese University of Hong Kong, Shenzhen, Guangdong 518172, P.R. China

& School of Mathematics Sciences, University of Science and Technology of China, Hefei, Anhui 230026, P.R. China,

<sup>b</sup> Université Le Havre Normandie, Normandie Univ., LMAH UR 3821, 76600 Le Havre, France.

November 13, 2025

## Abstract

This paper is concerned with the invasion fronts of spatially periodic monotone reaction-diffusion systems in a multi-dimensional setting. We study the pulsating travelling waves that connect the trivial equilibrium, for which all components of the state variable are identically equal to zero, to a uniformly persistent stationary state, for which all components are uniformly positive. When the trivial equilibrium is linearly unstable, we prove that any pulsating wave travels at a speed no less than that of the linearized system at the equilibrium, namely the minimal speed is linearly determined. It is also shown that pulsating waves are strictly monotonic in time whenever the nonlinearity is subhomogeneous. Beyond these general qualitative properties, the main focus of the paper is to derive sufficient conditions for the existence and nonexistence of pulsating waves propagating in any given direction. Our proof of the existence part relies upon a new level of understanding of the multi-dimensional pulsating waves observed from a direction-dependent coordinate system.

**Keywords:** monotone systems, semilinear parabolic equations, pulsating waves, rational direction of propagation.

**AMS Subject Classification (2020):** 35C07; 35K40, 35K57.

## 1 Introduction and main results

We consider the following quasimonotone system of semilinear parabolic equations

$$ \frac{\partial u}{\partial t} + Lu = f(x, u), \quad (t, x) \in \mathbb{R} \times \mathbb{R}^N, \qquad (1.1) $$

where $u = (u_1, \dots, u_d)^T$ for some $d > 1$ is the unknown of $t \in \mathbb{R}$ and $x \in \mathbb{R}^N$ for some $N \ge 1$, and $f = (f_1, \dots, f_d)^T: \mathbb{R}^N \times \mathbb{R}^d \to \mathbb{R}^d$ is a nonlinear vector-valued function satisfying

$$ \forall u \in \Sigma \subseteq \mathbb{R}^d, \quad \frac{\partial f_i}{\partial u_j}(\cdot, u) \ge 0 \quad \text{for all } i, j \in \{1, \dots, d\} \text{ whenever } i \ne j. \qquad (1.2) $$

<footnote>*Corresponding author.</footnote>

<footnote>E-mail addresses: dengll@cuhk.edu.cn (L. Deng), arnaud.ducrot@univ-lehavre.fr (A. Ducrot), quentin.griette@univ-lehavre.fr (Q. Griette).</footnote>
The symbol $L := \operatorname{diag}(L^1, \dots, L^d)$ in (1.1) denotes a diagonal matrix of second-order elliptic operators being either in divergence form

$$L^i(x)u_i := -\nabla \cdot (A^i(x)\nabla u_i) + q^i(x) \cdot \nabla u_i, \quad i = 1, \dots, d, \qquad (1.3)$$

or in nondivergence form

$$L^i(x)u_i := -\operatorname{tr}(A^i(x)D^2u_i) + q^i(x) \cdot \nabla u_i, \quad i = 1, \dots, d, \qquad (1.4)$$

where $\nabla$, $\nabla \cdot$, and $D^2$ denote the gradient, divergence, and Hessian with respect to $x \in \mathbb{R}^N$, respectively. Such types of parabolic systems are also referred to as (quasi-)monotone or cooperative reaction-diffusion systems. In this regard, we refer to [45, 58] for many aspects of modelling biomedical and ecological applications. An extensive review of earlier mathematical theories is provided in the monograph by Smith [55, Chap.7].

As is well known, the large time behaviour of the homogeneous Fisher-KPP equation with compactly supported initial data is governed by travelling wave solutions [24, 38]. Such solutions are a special class of entire solutions to nonlinear diffusion equations. These two pioneering papers have led to extensive research on the analysis of front propagation phenomena in reaction-diffusion equations and their systems. As far as travelling waves for systems are concerned, we refer to Volpert et al. [57, 58] for a comprehensive introduction and numerous self-contained treatments, and to Ogiwara and Matano [49, 50] for some qualitative analyses. Although a great deal of literature is devoted to homogeneous equations, the notion of travelling waves has been extended to non-autonomous heterogeneous equations over the past two decades, see Berestycki and Hamel [5] for a systematic survey. In this work, we shall study travelling waves for the monotone system (1.1) with spatially periodic heterogeneities in arbitrary dimensions.

## 1.1 Setting of the problem

Before formulating the problem, let us fix some notations to prevent any future confusion. For any two vectors $u, v \in \mathbb{R}^d$, we write $u \le v$ (resp. $u \ll v$) to mean $u_i \le v_i$ (resp. $u_i < v_i$) for all $i \in \{1, \dots, d\}$, while $u < v$ if $u \le v$ and $u_i < v_i$ for some $i$. The maximum between the two is understood componentwise, i.e. $\max(u, v) = (\max(u_1, v_1), \dots, \max(u_d, v_d))^T$. The minimum is defined similarly. With this definition, the maximum between two vectors coincides with the least upper bound, and the minimum is the greatest lower bound.

We denote the maximum (or uniform) norm of a vector $u \in \mathbb{R}^d$ by

$$|u|_{\infty} := \max_{1 \le i \le d} (|u_i|) = \lim_{p \to \infty} |u|_p$$

where $|u|_p := (\sum_{i=1}^d |u_i|^p)^{1/p}$ with $p \ge 1$ is the $p$-norm and by convention $| \cdot | := | \cdot |_2$ denotes the usual Euclidean norm. The function spaces $C^{1+\alpha/2,2+\alpha}(\mathbb{R} \times \mathbb{R}^N)^d$ with $\alpha \ge 0$ and $L^p(\mathbb{R} \times \mathbb{R}^N)^d$ with $p \ge 1$ are equipped with the following norms:

$$\|w\|_{C^{1+\alpha/2,2+\alpha}(\mathbb{R} \times \mathbb{R}^N)^d} = \max_{1 \le i \le d} \|w_i\|_{C^{1+\alpha/2,2+\alpha}(\mathbb{R} \times \mathbb{R}^N)}, \quad \|w\|_{L^p(\mathbb{R} \times \mathbb{R}^N)^d} = |(\|w_i\|_{L^p(\mathbb{R} \times \mathbb{R}^N)})_{1 \le i \le d}|_p.$$

Let $BUC(\mathbb{R} \times \mathbb{R}^N)^d$ be the Banach space of all bounded and uniformly continuous functions from $\mathbb{R} \times \mathbb{R}^N$ to $\mathbb{R}^d$ endowed with the usual supremum norm

$$\|u\|_{\infty} = \max_{1 \le i \le d} \sup_{(t,x) \in \mathbb{R} \times \mathbb{R}^N} |u_i(t,x)|.$$

With a slight abuse of notation, let $\mathbf{1} := (1, \dots, 1)^T$, $\mathbf{0} := (0, \dots, 0)^T \in \mathbb{R}^d$ for future convenience. A vector $u \in \mathbb{R}^d$ is said to be nonnegative, positive, and strictly positive, if it satisfies $u \ge \mathbf{0}$, $u > \mathbf{0}$, and $u \gg \mathbf{0}$, respectively.

Let us now elaborate on the general hypotheses. Denote by $\mathbb{T}^N = (\mathbb{R}/\mathbb{Z})^N$ the $N$-dimensional unit torus. We assume that the coefficients involved in $L$ given by (1.3) or (1.4) satisfy:
(C) The matrix fields $A^i : \mathbb{T}^N \to \mathcal{S}_N(\mathbb{R})$ (the space of real symmetric $N \times N$ matrices) with $A^i(x) = (a_{lm}^i(x))_{1\le l,m\le N}$ are of class $C^{1,\alpha}$ for some $\alpha > 0$, and uniformly elliptic in the sense that there exist positive constants $\underline{\gamma} \le \bar{\gamma}$ such that

$$ \forall (x, \xi) \in \mathbb{T}^N \times \mathbb{R}^N, \quad \underline{\gamma}|\xi|^2 \le \xi^T A^i(x) \xi \le \bar{\gamma}|\xi|^2 \quad \text{for all } i \in \{1, \dots, d\}, \qquad (1.5) $$

and the vector fields $q^i : \mathbb{T}^N \to \mathbb{R}^N$ with $q^i(x) = (q_l^i(x))_{1\le l\le N}$ are of class $C^\alpha$.

To formulate the conditions for $f(x, u)$, we first recall some concepts from [10, 56].

**Definition 1.1.** A matrix-valued function $H(z) = (h_{ij}(z))_{1\le i,j\le d}$ from an arbitrary domain $\Omega$ to $\mathcal{M}_d(\mathbb{R})$ (the space of real $d \times d$ matrices) is said to be:

* cooperative if $h_{ij}(z) \ge 0$ whenever $i \ne j$ for all $i, j \in \{1, \dots, d\}$ and $z \in \Omega$;

* fully coupled if for any nontrivial partitions $I, J \subset \{1, \dots, d\}$, that is, $I \cap J = \emptyset$ and $I \cup J = \{1, \dots, d\}$, there exist $i_0 \in I$ and $j_0 \in J$ such that $h_{i_0j_0}(z) \ne 0$.

Equivalently, the matrix field $H$ with bounded continuous entries is cooperative if and only if

$$ \underline{H} := \left( \inf_{z \in \Omega} h_{ij}(z) \right)_{1 \le i, j \le d} \quad \text{is essentially nonnegative (or called quasi-positive),} $$

while it is fully coupled if and only if

$$ \bar{H} := (\|h_{ij}\|_{\infty})_{1 \le i, j \le d} \quad \text{is irreducible.} $$

We are now in a position to state the set of hypotheses about the nonlinearity.

(N1) The nonlinearity $f = (f_1, \dots, f_d)^T : \mathbb{T}^N \times \mathbb{R}^d \to \mathbb{R}^d$ is continuous and of class $C^\alpha$ in $x$, locally uniformly for $u \in \mathbb{R}^d$. The partial derivatives $\partial f_i / \partial u_j$ exist and are continuous in $\mathbb{T}^N \times \mathbb{R}^d$ for all $i, j \in \{1, \dots, d\}$. Moreover, $f(x, 0) = 0$ and $f$ is of class $C^{1,\beta}(\mathbb{T}^N \times B_\infty(0, \sigma))$ for some $\beta \in (0, 1]$, where $B_\infty(0, \sigma) := \{u \in \mathbb{R}^d \mid |u|_\infty \le \sigma\}$ with a small $\sigma > 0$.

(N2) The Jacobian matrix field of $f(x, u)$ with respect to $u$, denoted as $D_u f(x, u)$, is cooperative for all $(x, u) \in \mathbb{T}^N \times \Sigma$ with some $\Sigma \subseteq \mathbb{R}^d$ in the sense of Definition 1.1.

(N3) There exists a constant $\hat{\eta} > 0$ such that $f(x, \hat{\eta}1) \le 0$ for all $x \in \mathbb{T}^N$.

(N4) The matrix field $D_u f(\cdot, u) : \mathbb{T}^N \to \mathcal{M}_d(\mathbb{R})$ is fully coupled for every $u \in \mathbb{R}^d$ in the sense of Definition 1.1.

Assumption (N2) is equivalent to condition (1.2). It is also possible to weaken (N2) such that $D_u f(\cdot, u)$ is cooperative only for $u \in B_\infty^+(0, \hat{\eta}) := \{u \in [0, \infty)^d \mid |u|_\infty \le \hat{\eta}\}$. We mention that a weak condition of local monotonicity was formulated by Volpert [57, 58]. Under (N2), a sufficient condition for (N3) to hold is that for each $i \in \{1, \dots, d\}$,

$$ \forall u > 0, \quad \lim_{l \to +\infty} \frac{f_i(x, lu_1, \dots, lu_d)}{l|u|_\infty} < 0 \quad \text{uniformly for } x \in \mathbb{T}^N. $$

Assumption (N4) means that the system (1.1) links all the unknowns $u_i$ in a strong sense, i.e. it cannot be separated into several subsystems independent of each other. Unless otherwise specified, the assumptions of (C) and (N1)-(N4) are made throughout the paper.

Let us now introduce the object we want to study. In the periodic framework, the notion of travelling periodic/pulsating waves was first proposed by Shigesada, Kawasaki and Teramoto [54] for a one-dimensional equation of Fisher type, and further generalized by Berestycki and Hamel [4] in a quite general setting, see also Xin [62], Hudson and Zinner [36], Ogiwara and Matano [49, 50], Heinze [34] for earlier works. As far as the monotone system (1.1) is concerned, such solutions are defined as follows.
**Definition 1.2.** Let $e \in \mathbb{S}^{N-1} := \{\zeta \in \mathbb{R}^N | |\zeta| = 1\}$ be given. A pulsating travelling wave propagating in the direction $e$ with the average (or effective) speed $c \neq 0$ is a classical bounded entire solution $u \equiv u(t, x)$ to (1.1) which is strictly positive and satisfies

$$ \forall k \in \mathbb{Z}^N, \forall (t, x) \in \mathbb{R} \times \mathbb{R}^N, \quad u\left(t + \frac{k \cdot e}{c}, x\right) = u(t, x - k), \qquad (1.6) $$

$$ \liminf_{x \cdot e \to -\infty} u(t, x) \gg 0, \quad \lim_{x \cdot e \to +\infty} u(t, x) = 0, \qquad (1.7) $$

where the limits are understood to hold locally uniformly for $t \in \mathbb{R}$ and uniformly with respect to the vectors of $e^{\perp} := \{y \in \mathbb{R}^N | y \cdot e = 0\}$.

Note that, since we are interested only in invasion fronts propagating into the unstable state $0$, a persistence-like property is introduced to describe one of the limiting states of pulsating waves, rather than imposing an *a priori* given steady state for (1.1). The goal of this paper is threefold. When the trivial steady state of (1.1) is unstable, we first establish some qualitative estimates for pulsating waves. In particular, we prove that any pulsating wave given by Definition 1.2 is strictly monotonic in time if the nonlinearity is subhomogeneous. The second part is devoted to the existence of pulsating waves. In the third part, we prove the nonexistence of pulsating waves when the trivial steady state of (1.1) is stable. Throughout the last two parts, we also discuss the asymptotic persistence (known as the *hair-trigger effect*) and extinction for the Cauchy problem associated with (1.1).

Let us add some other related works on the pulsating wave-like solutions for spatially periodic reaction-diffusion equations and their systems. As far as scalar equations are concerned, we refer to [2, 8, 16, 17, 20, 21, 26, 32, 33]. Recent advances in non-monotone systems can also be found in [14, 29, 30] and the references therein. From a different perspective, Weinberger [61] developed a general theory of spreading properties for a discrete-time recursion governed by an order preserving periodic operator, see [60] for his earliest work in this direction. This theory introduces a profound understanding of the relationship between travelling waves and the asymptotic spreading of solutions for a wide class of reaction-diffusion models. The results in [61] have further been extended in the past few years, see Liang and Zhao [42, 43] for monostable case. We emphasize that applying the abstract results of [19, 42, 43] to monotone systems does require stronger assumptions on (1.1) than those made here. This is also one of our primary motivations for the present work. A detailed discussion will follow the statement of our main results. In addition, Fang and Zhao [23] conducted a deep study of travelling waves for bistable monotone semiflows, which includes results on one-dimensional bistable pulsating waves for monotone systems. Giletti and Rossi [26] extended the results of [21, 23] to multidimensional scalar equations of both bistable and multistable types. The aforementioned works contributed with different methodologies for studying travelling waves of monotone systems and generalized several results of Volpert et al. [57, 58]. It should also be noted that the presence of multiple steady states in such systems gives rise to numerous other types of travelling waves or pulses [57, 58], including bistable ones, which are not addressed here within the periodic framework adopted in this work.

Let us also mention that several works have focused on pulsating waves in space-time periodic media. We refer to Nadin [47] for scalar equations, and to Fang et al. [22] for monostable monotone semiflows with one spatial dimension, as well as Du et al. [19] with arbitrary space dimensions. Furthermore, the notion of generalized travelling waves in random media was introduced by Matano [44] (see also Shen [53]), and transition fronts under a fairly general setting, by Berestycki and Hamel [5]. The literature on this topic is now extensive and a complete review is beyond the scope of this work.

To conclude this section, let us point out that our motivation for this work is to understand the kind of additional conditions to be imposed on the system (1.1) which allow front propagation into the unstable trivial equilibrium.
## 1.2 Statement of the results

To state our main results, several notions related to linear problems are required. For this purpose, for each $e \in \mathbb{S}^{N-1}$ and each $\lambda \in \mathbb{R}$, we first introduce the operators $L_{\lambda e}^i$ as follows:

$$ 
\begin{aligned}
L_{\lambda e}^i \psi &:= e^{\lambda e \cdot x} L^i \left( e^{-\lambda e \cdot x} \psi \right) \\
&= -\nabla \cdot (A^i \nabla \psi) + 2\lambda e A^i \nabla \psi + q^i \cdot \nabla \psi - [\lambda^2 e A^i e - \lambda \nabla \cdot (A^i e) + \lambda q^i \cdot e] \psi
\end{aligned} 
\quad (1.8) $$

or

$$ L_{\lambda e}^i \psi := -\operatorname{tr}(A^i D^2 \psi) + 2\lambda e A^i \nabla \psi + q^i \cdot \nabla \psi - (\lambda^2 e A^i e + \lambda q^i \cdot e) \psi, $$

where we recall that the operators $L^i$ are given by (1.3) or (1.4).

Set $L_{\lambda e} := \operatorname{diag}(L_{\lambda e}^1, \dots, L_{\lambda e}^d)$ and let $k(\lambda, e)$ denote the periodic principal eigenvalue of the system of linear elliptic operators

$$ \phi \mapsto L_{\lambda e} \phi - D_u f(x, 0) \phi $$

posed for $\phi \in C^2(\mathbb{T}^N)^d$ with periodicity conditions. In other words, $k(\lambda, e)$ is the unique real number such that there exists a function $\phi_{\lambda e} \in C^2(\mathbb{T}^N)^d$ satisfying

$$ 
\left\{ 
\begin{array}{l}
L_{\lambda e} \phi_{\lambda e} - D_u f(x, 0) \phi_{\lambda e} = k(\lambda, e) \phi_{\lambda e} \quad \text{in } \mathbb{T}^N, \\
\phi_{\lambda e} \gg 0 \quad \text{in } \mathbb{T}^N, \quad \|\phi_{\lambda e}\|_{\infty} = 1.
\end{array}
\right. 
\quad (1.9) $$

When $\lambda = 0$, we set

$$ \lambda_1^p := k(0, e). \quad (1.10) $$

Observe from (1.8) that $\lambda_1^p$ does not depend on $e$. In what follows, we say that $0$ is a linearly unstable solution of (1.1) if $\lambda_1^p < 0$, while it is stable if $\lambda_1^p \ge 0$.

We now introduce the Dirichlet principal eigenvalue for the operator $L_{\lambda e} - D_u f(x, 0)$ with $\lambda = 0$. Let $B_R$ be the open ball of $\mathbb{R}^N$ with centre $0$ and radius $R$. For any $R > 0$, we consider the unique $\lambda_1^R \in \mathbb{R}$ such that there exists a function $\varphi^R \in C^2(B_R)^d \cap C^1(\overline{B_R})^d$ which satisfies

$$ 
\left\{ 
\begin{array}{ll}
L\varphi^R - D_u f(x, 0)\varphi^R = \lambda_1^R \varphi^R & \text{in } B_R, \\
\varphi^R = 0 & \text{on } \partial B_R, \\
\varphi^R \gg 0 \text{ in } B_R, & \|\varphi^R\|_{\infty} = 1.
\end{array}
\right. 
\quad (1.11) $$

The existence and uniqueness of the principal eigenpairs follow from the Krein-Rutman theory and the fact that under assumptions (N2) and (N4), $D_u f(x, 0)$ is cooperative and fully coupled in the sense of Definition 1.1, see [10, 40, 56]. With the notion of the generalized principal eigenvalue for scalar elliptic equations (see [7, 9] for instance), we also introduce such a quantity for the operator $L - D_u f(x, 0)$ in $\mathbb{R}^N$. This definition here reads

$$ \lambda_1 := \sup \left\{ \lambda \in \mathbb{R} \mid \exists \varphi \in C^2(\mathbb{R}^N)^d, \varphi \gg 0 \text{ and } (L - D_u f(x, 0))\varphi \ge \lambda \varphi \text{ in } \mathbb{R}^N \right\}. \quad (1.12) $$

Assume now $\lambda_1^p < 0$. Then the following quantity turns out to be well defined:

$$ c^*(e) := \inf_{\lambda > 0} \frac{-k(\lambda, e)}{\lambda}. \quad (1.13) $$

Moreover, the infimum in (1.13) can also be reached on $(0, \infty)$ and there holds

$$ \lambda_1^\infty := \lim_{R \to +\infty} \lambda_1^R \in \mathbb{R}, \quad \lambda_1^p \le \lambda_1^\infty = \lambda_1 = \max_{\lambda e \in \mathbb{R}^N} k(\lambda, e). \quad (1.14) $$

Further properties for the quantities introduced above will be stated in Section 2.

With these quantities defined, our main results are stated as the following three theorems.
**Theorem 1.1.** Assume $\lambda_1^p < 0$. Let $(c, u)$ be a pulsating travelling wave solution to (1.1). Then the wave speed $c$ satisfies

$$c \geq c^*(e).$$

If $f$ is further assumed to be subhomogeneous, i.e.

$$\forall \theta \in (0, 1), \quad \forall (x, u) \in \mathbb{T}^N \times [0, \infty)^d, \quad \theta f(x, u) \leq f(x, \theta u), \qquad (1.15)$$

then the function $u = u(t, x)$ is componentwise increasing in $t$ if $c > 0$ and decreasing if $c < 0$. Moreover, one has $c\partial_t u(t, x) \gg 0$ for all $(t, x) \in \mathbb{R} \times \mathbb{R}^N$.

**Theorem 1.2.** Assume $\lambda_1 < 0$. Suppose the nonlinearity $f$ satisfies

$$\forall (x, u) \in \mathbb{T}^N \times [0, \infty)^d, \quad f(x, u) \leq D_u f(x, 0)u. \qquad (1.16)$$

Then, for every $e \in \mathbb{S}^{N-1}$ and for all $c \geq c^*(e)$, the system (1.1) admits a pulsating travelling wave in direction $e$ with speed $c$.

**Theorem 1.3.** Assume $\lambda_1^p \geq 0$. Let condition (1.16) be fulfilled, and assume further that the inequality (1.16) is strict whenever $u \gg 0$, i.e.

$$\forall (x, u) \in \mathbb{T}^N \times (0, \infty)^d, \quad f(x, u) < D_u f(x, 0)u. \qquad (1.17)$$

Then there exists no pulsating travelling wave of (1.1) in any given direction.

**Remark 1.4.** The set of our hypotheses (N1)-(N4) together with (1.16) is actually an analogue of (F1)-(F5) formulated by Barles et al. [3], in which the authors investigated front propagation in homogeneous media by developing a viscosity solution approach for Hamilton-Jacobi equations. In fact, under assumptions (N1) and (1.15), the sublinear condition (1.16) holds. Furthermore, if $f$ is strictly subhomogeneous, i.e.

$$\forall \theta \in (0, 1), \quad \forall (x, u) \in \mathbb{T}^N \times (0, \infty)^d, \quad \theta f(x, u) < f(x, \theta u), \qquad (1.18)$$

then the strictly sublinear condition (1.17) is also fulfilled.

**Remark 1.5.** If $\lambda_1 < 0$, it follows from (1.13) and (1.14) that $c^*(e) > 0$ for every $e \in \mathbb{S}^{N-1}$. Furthermore, if $\lambda_1^p = \lambda_1$, then the first two theorems show that the system (1.1) admits a pulsating wave if and only if $c \geq c^*(e)$. Nevertheless, these two quantities are not necessarily equal. Depending on the assumptions on the linearized operator, a strict inequality $\lambda_1^p < \lambda_1$ is possible, see Remark 2.4 below.

**Remark 1.6.** It is worth pointing out that the proof of the linear determinacy for the minimal wave speed does not require the limiting condition of pulsating waves at $x \cdot e = -\infty$ in (1.7). Moreover, the existence of pulsating waves is guaranteed under assumption (1.16) even if the monotonicity does not hold in general.

Theorem 1.3 asserts that the linearly stable equilibrium precludes invasion by any uniformly positive state. Indeed, under assumptions of Theorem 1.3, Proposition 5.1 in Sect. 5 implies that there is no nonnegative bounded entire solution to (1.1) besides the trivial one. It should also be stressed that condition (1.17) is essential to the validity of Theorem 1.3. Recall that in the case of scalar reaction-diffusion equations, pulsating waves are known to exist for both ignition-type nonlinearities (satisfying $\lambda_1^p = 0$, see [4]) and bistable ones (satisfying $\lambda_1^p > 0$, see [17, 20]).

Before proceeding, let us further comment on the above results and recall some earlier works in the literature. Theorem 1.1 provides a sufficient condition for pulsating waves to be monotonic in time. As a consequence, if (1.15) holds, then any pulsating wave $(u, c)$ converges, as $t \to +\infty$ if
$c > 0$ or as $t \to -\infty$ if $c < 0$, to a bounded strictly positive steady state for (1.1), which can also be verified to be $\mathbb{Z}^N$-periodic because of (1.6). This observation naturally motivates us to seek pulsating waves that connect the unstable state $0$ and a stable steady state $p$ for (1.1). To this end, if there exists a minimal periodic steady state $p$ in the class of all strictly positive bounded steady states for (1.1), and that it is asymptotically stable for the semiflow generated by (1.1) in the set $\{\varphi \in C(\mathbb{T}^N)^d | 0 < \varphi \le p\}$, then the limiting condition (1.7) could be refined to

$$ \lim_{x \cdot e \to -\infty} |u(t, x) - p(x)|_{\infty} = 0, \quad \lim_{x \cdot e \to +\infty} u(t, x) = 0, \qquad (1.19) $$

locally uniformly for $t \in \mathbb{R}$ and uniformly in $e^\perp$. Under such a monostable setting, the existence of pulsating waves satisfying (1.19) is a direct consequence of Weinberger's result in [61] for one component case and Liang-Zhao's result in [43] for one-dimensional case. Nevertheless, besides contributing additional qualitative properties, the present paper provides a feasible and alternative approach to establishing the existence of pulsating waves for monostable monotone systems. Furthermore, our work extends some of the results on planar travelling waves in homogeneous media from Volpert et al. [57] and Li et al. [41], as well as those on curved travelling waves in straight cylindrical domains from Volpert [58].

However, it is technically challenging to find some sharp sufficient conditions under which the system (1.1) goes under the above-mentioned framework. A special case worth noting is that by strengthening (1.15) to (1.18), one can apply the theory of strongly monotone (guaranteed here by assumption (N4)) and strictly subhomogeneous maps, see [64, Sect 2.3] for instance, to yield a threshold-type result for the semiflow generated by (1.1) in terms of the stability properties of the trivial equilibrium, see [15, 35, 59] for some examples of such arguments. This also suggests us to establish the Liouville-type results for monotone system (1.1) as preparation for studying invasion dynamics. Observe that under assumption (1.16), the system (1.1) shares structural similarities with scalar Fisher-KPP equations. According to the results of Berestycki et al. [7, 9], the strong-KPP condition is almost sharp for the uniqueness and stability of bounded positive steady states for scalar semilinear parabolic equations in the whole space. In fact, under assumptions $\lambda_1^p < 0$ and (1.16), one can readily prove that there exists at least a strictly positive periodic steady states for the monotone system (1.1). But the uniqueness (or even the minimality) and the stability within the class of all nonnegative bounded steady states, namely without assuming a priori that any steady state for (1.1) is uniformly positive and periodic, is far from straightforward.

On the other hand, in order to establish monotonicity results, we assume in Theorem 1.1 that $f$ is subhomogeneous in the sense of (1.15). This condition plays a role similar to that of the strong-KPP assumption made by Berestycki et al. [8] for scalar equations. But a weaker assumption was formulated by Hamel [31], which also guarantees the monotonicity of monostable pulsating waves. A generalization to monotone systems becomes then natural yet more complicated. To sum up, it remains an interesting and open problems to derive sharp conditions that yield Liouville-type results for the stationary system of (1.1) in the whole space. For these reasons, we will not here formulate any assertion leading to (1.19) but adopt the weaker condition (1.7).

Let us finally mention that Du et al. [18, 19] extended the results of [22, 63] concerning a class of monostable monotone systems to a multi-dimensional periodic medium in which the system admits equilibria lying on the boundary of the positive cone in a Banach lattice. Typically, the nonlinearity considered in [18] takes the form

$$ \begin{aligned} f_1(x, u_1, \dots, u_d) &= u_1 h_1(x, u_1, \dots, u_d), \\ f_i(x, u_1, \dots, u_d) &= \sum_{j=1}^{i-1} c_{ij}(x) u_j + u_i h_i(x, u_1, \dots, u_d), \quad 2 \le i \le d, \end{aligned} \qquad (1.20) $$

under certain conditions imposed on the functions $h_i$ and $c_{ij}$. A classical model that fits their setting is a cooperative system derived from the Lotka-Volterra competitive system after performing
a well known change of variables. However, there are many examples of $f$ which could not be
unified as (1.20), see Sect. 1.4 below for details. In this regard, our setting thus covers a large
class of other types of nonlinearities not previously treated, thereby contributing new results.

1.3 Elements of the proof

The proof of Theorem 1.1 adapts the arguments of Berestycki et al. [8] and Hamel [31] who investigated scalar equations of strong-KPP and monostable types. The proof of Theorem 1.2 uses a general methodology developed by the first two authors of the present paper in [13, 14]. Specifically, the first step constructs pulsating waves in each direction of unit vectors with rational coordinates (denoted as $Q^N \cap S^{N-1}$), and the second step performs an approximation procedure for a general direction in $S^{N-1}$. For the latter purpose and to prove Theorem 1.3, we adapt some arguments developed by Berestycki et al. for scalar equations [7, 9]. Observe that the profile of one-dimensional pulsating waves at time $T = t+1/c$ is an exact copy of the profile at $T = t$, shifted in space by one unit of length. This observation provides a key motivation for formulating the pulsating wave as a fixed-point problem. However, in the multi-dimensional setting, the problem (that is, to verify (1.6)) becomes much more involved.

We now highlight the main new ingredients in the present paper for constructing a pulsating
wave propagating along any rational direction, which differ from those in [13, 14]. Firstly, we
construct sub- and supersolutions for all admissible speeds, especially for the critical case, inspired
again by [31]. Secondly, by slightly modifying the coordinate transformation introduced in [13],
pulsating waves in a rational direction, observed from a new system of coordinates, become
(roughly speaking) pulsating waves propagating inside a cylinder with a periodic cross-section.
Based on that and inspired by Ogiwara and Matano [50, Sect.7], we reformulate a pulsating wave
with speed $c > 0$ as a fixed point of the solution map at $t = \tau_1/c$ composed with a shift by $\tau_1$ along
the cylinder axis, where $\tau_1 > 0$ is a parameter determined by the period of the medium and the
rational direction of propagation. Lastly, we address this issue by considering an initial/boundary-
value problem posed in a semi-infinite periodic cylinder with a moving (at velocity $c$) boundary
subject to the homogeneous Neumann data, which is then solved by the Schauder fixed-point
theorem (see Proposition 4.5). In particular, to verify the limiting conditions of the solution, a
novel lower barrier is also constructed (see (4.42)-(4.45)). This approach was developed by Griette
and Matano [29] who investigated pulsating one-dimensional waves for other models.

Let us provide further remarks on our method. For scalar equations, a common and powerful
approach to constructing pulsating waves is to solve a semilinear degenerate elliptic equation via
a regularization procedure [4, 8, 20]. This approach requires refined uniform estimates of the
solution with respect to a vanishing viscosity parameter. However, it would certainly require
more PDE techniques and possibly make additional assumptions on the coefficients in (1.4) to
extend such estimates, especially gradient estimates, to general systems as (1.1). In contrast, a
key contribution of our method, as already demonstrated for the scalar Fisher-KPP equation in
[13], lies in overcoming the difficulty caused by elliptic degeneracy. The present work generalizes
this method to the monotone system (1.1) and provides a new approach to constructing pulsating
waves in rational directions. In [13], the construction of pulsating waves in a rational direction
relies on the approximation of the equation in bounded domains. To do so, we first solved a
space-time periodic boundary value problem and then performed an iteration procedure. Here
we devise a subtle fixed-point problem as outlined above to construct pulsating waves in rational
directions, which actually leads to considerable simplifications compared to [2, 13]. For a non-
rational direction of propagation, we recover existence by passing the limit as in [13, 14].
## 1.4 Examples

Monotone reaction-diffusion systems arise in diverse applications [58]. As a topic in theoretical ecology, such systems describe the dynamics of mutualistic symbiosis between two or more species [45]. Here we briefly present several examples to illustrate the applications and motivations of the present work. We provide only the explicit form of the nonlinearities that fall within our settings, while the corresponding modelling mechanism refers to the references cited.

The first example is a very classical model (known as a “positive feedback loop” of d species) for controlling protein synthesis in the cell [55, Sect. 4.2 and Sect. 7.7] (see also [45, Sect. 7.2]) in which the nonlinearity is defined as

$$ 
\begin{aligned}
f_1(u_1, \dots, u_d) &= g(u_d) - \alpha_1 u_1, \\
f_i(u_1, \dots, u_d) &= u_{i-1} - \alpha_i u_i, \quad 2 \le i \le d,
\end{aligned} 
\qquad (1.21) 
$$

with positive coefficients $\alpha_i$ and $d > 2$. The function $g: \mathbb{R}_+ \to \mathbb{R}_+$ is assumed to be:

$$ \forall s > 0, \quad 0 < g(s) < M, \quad g'(s) > 0 \text{ and } g(0) = 0. \qquad (1.22) $$

A typical example of $g$ from [55] is given by

$$ g(s) = \frac{s^p}{1 + s^p} \quad \text{for some integer } p \ge 1. $$

With such a choice and for any $p > 1$, the function $f = (f_1, \dots, f_d)$ defined in (1.21) satisfies assumptions (N1)-(N4) but verifies neither subhomogeneity (1.15) nor sublinearity (1.16). When $p = 1$, then (1.15) is fulfilled. Furthermore, a slight modification of the assumptions on $g$ satisfying (1.22) can yield further insights. We consider two types for the function $g$ as follows:

* For all $s > 0$, $0 < g(s) \le g'(0)s$ (weak-KPP hypothesis),

* The map $s \mapsto \frac{g(s)}{s}$ is strictly decreasing in $(0, \infty)$ (strong-KPP hypothesis).

It is straightforward to verify that if $g$ is weak-KPP type, then the function $f$ defined in (1.21) only satisfies the sublinear condition (1.16). However, the strictly subhomogeneous condition (1.18) is verified provided $g$ is strong-KPP type. In the latter case, condition (1.17) is thus guaranteed from Remark 1.4.

A second example is an SIR-type model for rabies circulation between two populations [15], in which the model assumes that susceptible individuals remain motionless to capture the disease dynamics. After performing a change of variables, the authors turned such an epidemic model into a monotone system from which the nonlinearity is of the form

$$ 
\begin{aligned}
f_1(x, u_1, u_2) &= S_1^0(x) \left(1 - e^{-(\beta_{11}(x)u_1 + \beta_{12}(x)u_2)}\right) - \delta_1(x)u_1, \\
f_2(x, u_1, u_2) &= S_2^0(x) \left(1 - e^{-(\beta_{21}(x)u_1 + \beta_{22}(x)u_2)}\right) - \delta_2(x)u_2.
\end{aligned} 
$$

With all parameters being periodic in $x$ and positive everywhere, the function $f = (f_1, f_2)$ satisfies not only assumptions (N1)-(N4) but also the strictly (in fact, strongly) subhomogeneous condition (1.18). Hence all the aforementioned conditions are fulfilled in this example.

The last example is drawn from an ecological model of population migration among distinct patches (such as cities) [55, Sect. 4.4] in which the nonlinearity is defined as

$$ 
\begin{aligned}
f_1(u_1, \dots, u_d) &= \epsilon(u_2 - u_1) + r_1 u_1 \left(1 - \frac{u_1}{K_1}\right), \\
f_2(u_1, \dots, u_d) &= \epsilon(u_1 + u_3 - 2u_2) + r_2 u_2 \left(1 - \frac{u_2}{K_2}\right), \\
&\vdots \\
f_d(u_1, \dots, u_d) &= \epsilon(u_{d-1} - u_d) + r_d u_d \left(1 - \frac{u_d}{K_d}\right).
\end{aligned} 
$$
One can easily check that the function $f = (f_1, \dots, f_d)$ satisfies assumptions (N1)-(N4) and it is strongly subhomogeneous since the population growth obeys the Logistic law.

*Outline of the paper.* The remainder of the paper is as follows: Sect. 2 contains some preliminaries and the proof of which is provided in Appendix A. In Sect. 3, we derive the lower bound for the wave speeds, its linear determinacy, and the monotonicity results. Sect. 4 and Sect. 5 prove the existence and nonexistence of pulsating waves, respectively. Appendix B contains two Harnack-type estimates that are used several times in our proofs.

## 2 Preliminaries

In this section, we collect some important properties of the principal eigenpairs for a cooperative fully coupled system of elliptic operators with periodic coefficients. The relationship among the principal eigenvalues in several senses is also analyzed.

**Proposition 2.1.** Let $H \in C^\alpha(\mathbb{T}^N, \mathcal{M}_d(\mathbb{R}))$ be a cooperative and fully coupled matrix field according to Definition 1.1. Let $(\lambda_1^R, \varphi^R)$ be the principal eigenpair of the problem

$$ \begin{cases} (L-H)\varphi^R = \lambda_1^R \varphi^R & \text{in } B_R, \\ \varphi^R = \mathbf{0} & \text{on } \partial B_R, \\ \varphi^R \gg \mathbf{0} & \text{in } B_R, \end{cases} $$

where $L$ denotes a diagonal matrix of linear elliptic operators given by (1.4). Let

$$ \lambda_1 = \sup \left\{ \lambda \in \mathbb{R} \mid \exists \psi \in C^2(\mathbb{R}^N)^d, \psi \gg \mathbf{0} \text{ and } (L-H)\psi \ge \lambda \psi \text{ in } \mathbb{R}^N \right\} \quad (2.1) $$

denote the generalized principal eigenvalue of the operator $L-H$ in $\mathbb{R}^N$.

Then $\lambda_1^R \searrow \lambda_1$ as $R \to +\infty$. Moreover, there exists a generalized principal eigenfunction $\varphi \in C^2(\mathbb{R}^N)^d$ associated with $\lambda_1$, i.e. it satisfies $\varphi \gg \mathbf{0}$ and

$$ (L-H)\varphi = \lambda_1\varphi \quad \text{in } \mathbb{R}^N. $$

**Proposition 2.2.** Let $L_{\lambda e}$ be a diagonal matrix of linear elliptic operators given by (1.8) and $H \in C^\alpha(\mathbb{T}^N, \mathcal{M}_d(\mathbb{R}))$ be cooperative and fully coupled in the sense of Definition 1.1. Then the following statements hold.

(i) For all $\lambda \in \mathbb{R}$ and for each $e \in \mathbb{S}^{N-1}$, there exists a unique $k = k(\lambda, e) \in \mathbb{R}$ (principal eigenvalue) and a function $\phi = \phi_{\lambda e} \in C^2(\mathbb{T}^N)^d$ with $\phi_{\lambda e} \gg \mathbf{0}$, unique up to positive scalar multiple, which solve the system of periodic elliptic equations

$$ (L_{\lambda e} - H)\phi = k\phi \quad \text{in } \mathbb{T}^N. \quad (2.2) $$

(ii) The principal eigenvalue $k(\lambda, e)$ can be characterized by

$$ k(\lambda, e) = \max_{\phi \gg \mathbf{0}, \phi \in C^2(\mathbb{T}^N)^d} \min_{1 \le i \le d} \inf_{x \in \mathbb{T}^N} \frac{(e^i(L_{\lambda e} - H)\phi)(x)}{\phi_i(x)}, \quad (2.3) $$

where $e^i \in \mathbb{S}^{d-1}$ denotes the row vector with only the $i$-th coordinate equal to 1.

(iii) For any given $e \in \mathbb{S}^{N-1}$, the function $\lambda \mapsto k(\lambda, e)$ is analytic and strictly concave in $\mathbb{R}$. Furthermore, there exists $\beta \in \mathbb{R}$ such that

$$ \forall \lambda \in \mathbb{R}, \quad k(\lambda, e) \le -\underline{\gamma}\lambda^2 + \beta\lambda + \max_{1 \le i \le d} \|h_{ii}\|_{\infty}, \quad (2.4) $$

where $\underline{\gamma}$ has been defined in (1.5).
(iv) Let $\phi_{\lambda e}$ be the unique principal eigenfunction such that

$$ \phi_{\lambda e} \gg 0 \quad \text{and} \quad \|\phi_{\lambda e}\|_{\infty} = 1. \qquad (2.5) $$

Then for all $\lambda \in \mathbb{R}$, the map $\mathbb{S}^{N-1} \ni e \mapsto k(\lambda, e)$ is continuous and the normalized principal eigenfunction $\phi_{\lambda e}$ also depends continuously on $\lambda e \in \mathbb{R}^N$ with respect to the uniform topology.

(v) One has

$$ \lambda_1 = \max_{\lambda e \in \mathbb{R}^N} k(\lambda, e), \qquad (2.6) $$

where $\lambda_1$ is defined by (2.1). Furthermore, for any given $e \in \mathbb{S}^{N-1}$, the maximum in (2.6) is attained at a unique $\bar{\lambda} \in \mathbb{R}$ such that the function $\mathbb{R}^N \ni x \mapsto \phi_{\bar{\lambda}e}(x)e^{\bar{\lambda}e \cdot x}$ is a generalized principal eigenfunction associated with $\lambda_1$.

Before proceeding further, we mention a recent work by Girardin and Mazari [28] who made a deep study of principal eigenvalues for a linear cooperative system of space-time periodic parabolic operators. In particular, their results can be applied to explain in various ways how spatial heterogeneity affects $k(\lambda, e)$. The aforementioned results in the elliptic setting are analogues to their results. For instance, Propositions 3.2-3.3 in [28] are more general than Proposition 2.1. But we include in Appendix A the proof of partial results of Proposition 2.2 for completeness of the present paper. Most of the proofs are direct adaptations of the arguments developed by Berestycki et al. [4, 6, 8] and Nadin [46] for scalar equations.

The next result is a direct consequence of Proposition 2.2.

**Proposition 2.3.** Assume $k(0, e) < 0$. Then the following conclusions hold true.

(i) The quantity $c^*(e)$ defined by (1.13) is a real number. Moreover, for every fixed $c \ge c^*(e)$, the following quantity

$$ \lambda_c(e) := \min\{\lambda > 0 \mid k(\lambda, e) + c\lambda = 0\} \qquad (2.7) $$

is well-defined and there holds

$$ c^*(e) = \min_{\lambda > 0} \frac{-k(\lambda, e)}{\lambda}. \qquad (2.8) $$

(ii) Both $\lambda_c(e)$ and $c^*(e)$ are continuous in $e \in \mathbb{S}^{N-1}$.

(iii) The equation

$$ k(\lambda, e) + c\lambda = 0 \quad \text{for } \lambda > 0 \qquad (2.9) $$

has solutions if and only if $c \ge c^*(e)$. For $c = c^*(e)$, there exists a unique solution $\lambda^* := \lambda_{c^*}(e) > 0$, while for all $c > c^*(e)$, there are exactly two solutions $\lambda_c^- := \lambda_c(e)$ and $\lambda_c^+ = \lambda_c^+(e) > 0$.

**Remark 2.4.** From [29, Proposition 2.2], two sufficient conditions for $\lambda_1 = \lambda_1^p$ to hold can be summarized as follows:

* For each $i \in \{1, \dots, d\}$, the operator $L^i$ is in divergence form (1.3) and has no advection, that is, $q^i(x) \equiv 0$. Moreover, the matrix field $H$ in (2.2) is symmetric.

* For each $i \in \{1, \dots, d\}$, the operator $L^i$ is in nondivergence form (1.4). Furthermore, the diffusion matrix field $x \mapsto A^i(x)$ in $L^i$ and the function $x \mapsto H(x)$ in (2.2) are even, while the advection vector field $x \mapsto q^i(x)$ in $L^i$ is odd.
Indeed, one can prove that the function $\lambda \mapsto k(\lambda, e)$ is even provided one of the above conditions is fulfilled. Hence its global maximum over $\mathbb{R}$ is necessarily attained at $\lambda = 0$ due to the strict concavity of $\lambda \mapsto k(\lambda, e)$. We then have $\lambda_1 = k(0, \cdot) = \lambda_1^p$ by (2.6). However, a very interesting counter-example was constructed by Griette and Matano [29] for which $\lambda_1^p < \lambda_1$ holds, see Remark 4.2 and Appendix of [29] for details. Nevertheless, due to the smoothness assumptions of (C), our proof only needs to handle the operators defined by (1.4).

# 3 Qualitative properties of pulsating waves

This section is devoted to the proof of Theorem 1.1.

## 3.1 Lower bound on the wave speed

We prove here that the speed $c$ of pulsating waves for (1.1) is always bounded from below by the quantity $c^*(e)$, which especially implies that the system (1.1) has no such solutions for all $c < c^*(e)$. The proof of the first part of Theorem 1.1 is given by the following proposition.

**Proposition 3.1.** Let $u = (u_1, \dots, u_d)^T(t, x)$ be a pulsating travelling wave of (1.1) propagating in the direction $e \in \mathbb{S}^{N-1}$ with speed $c \neq 0$. Using the change of variables

$$ U(s, x) = u \left( \frac{x \cdot e - s}{c}, x \right) \quad \text{for } s \in \mathbb{R} \text{ and } x \in \mathbb{R}^N, \qquad (3.1) $$

then each component $U_i$ of $U = (U_1, \dots, U_d)^T$ satisfies

$$ -\infty < \underline{\mu} \le \liminf_{s \to +\infty} \left( \min_{x \in \mathbb{R}^N} \frac{-\partial_s U_i(s, x)}{U_i(s, x)} \right) \le \limsup_{s \to +\infty} \left( \max_{x \in \mathbb{R}^N} \frac{-\partial_s U_i(s, x)}{U_i(s, x)} \right) \le \bar{\Lambda} < +\infty, $$

where $\underline{\mu}$ and $\bar{\Lambda}$ are defined as

$$ \underline{\mu} = \min_{1 \le i \le d} \mu_i \quad \text{and} \quad \bar{\Lambda} = \max_{1 \le i \le d} \Lambda_i \qquad (3.2) $$

with

$$ \mu_i := \liminf_{s \to +\infty} \left( \min_{x \in \mathbb{R}^N} \frac{-\partial_s U_i(s, x)}{U_i(s, x)} \right) \quad \text{and} \quad \Lambda_i := \limsup_{s \to +\infty} \left( \max_{x \in \mathbb{R}^N} \frac{-\partial_s U_i(s, x)}{U_i(s, x)} \right). \qquad (3.3) $$

Assume further that $\lambda_1^p < 0$. Then both $\underline{\mu}$ and $\bar{\Lambda}$ are positive and solve (2.9). In particular, one has

$$ c \ge c^*(e), $$

where $\lambda_1^p$ and $c^*(e)$ are given by (1.10) and (2.8), respectively.

*Proof.* The proof relies on arguments that are close to the ones used by Hamel [31, Proposition 2.2]. Here is actually a generalization towards reaction-diffusion systems.

We first prove that two quantities $\underline{\mu}$ and $\bar{\Lambda}$ in (3.2) are well-defined. To this end, let us rewrite the system (1.1) as

$$ \partial_t u(t, x) + Lu(t, x) = H(t, x)u(t, x), \quad t \in \mathbb{R}, x \in \mathbb{R}^N, \qquad (3.4) $$

where the matrix-valued function $H = (h_{ij})_{1 \le i,j \le d} : \mathbb{R} \times \mathbb{R}^N \to \mathcal{M}_d(\mathbb{R})$ is defined by

$$ h_{ij}(t, x) = \int_0^1 \frac{\partial f_i}{\partial u_j}(x, s u_1(t, x), \dots, s u_d(t, x)) \, ds. \qquad (3.5) $$
By assumption (N1) on $f$ coupled with the boundedness and regularity of $u$ as a pulsating wave for (1.1), it is clearly true that all entries of $H$ are globally bounded and of class $C^{\alpha/2,\alpha}$. Fix any $(t_0, x_0) \in \mathbb{R} \times \mathbb{R}^N$ and let $P_r(t_0, x_0) := [t_0 - r^2, t_0] \times \bar{B}_r(x_0)$ be the closed parabolic cylinder of radius $r$, height $r^2$, and top center $(t_0, x_0)$. From the Schauder interior estimates for parabolic systems (see [39, Chap.VII] or [52]), there exists $C_r > 0$ independent of $(t_0, x_0)$ such that

$$ \forall r \in (0, 1), \quad \|u\|_{C^{1+\alpha/2,2+\alpha}(P_r)^d} \le C_r \|u\|_{L^\infty(P_1)^d}. $$

In particular, one has, for all $(t_0, x_0) \in \mathbb{R} \times \mathbb{R}^N$ and for each $i \in \{1, \dots, d\}$,

$$ |\partial_t u_i(t_0, x_0)| + |\nabla u_i(t_0, x_0)| \le C \max_{1 \le j \le d} \max_{(t,x) \in P_1} u_j(t, x). \quad (3.6) $$

Let now $\theta_c > 1$ be given where $c \ne 0$ is the speed of the pulsating wave $u$ of (1.1). By our assumptions (N2) and (N4), the matrix $H$ given by (3.5) is cooperative and fully coupled. Applying the Harnack inequality from Theorem B.2 to (3.4), we then deduce that there exists $\hat{\kappa}_{\theta_c} > 0$ independent of $(t_0, x_0)$ such that

$$ \begin{aligned} \max_{1 \le j \le d} \max_{(t,x) \in P_1} u_j(t, x) &\le \max_{1 \le j \le d} \max_{(t,x) \in [t_0 - \theta_c, t_0] \times [x_0 - \theta_c, x_0 + \theta_c]^N} u_j(t, x) \\ &\le \hat{\kappa}_{\theta_c} \min_{1 \le j \le d} \min_{(t,x) \in [t_0 + \frac{3\theta_c}{2}, t_0 + 2\theta_c] \times [x_0 - \theta_c, x_0 + \theta_c]^N} u_j(t, x). \end{aligned} $$

Choose a vector $k_0 \in \mathbb{Z}^N$ such that $3\theta_c/2 \le (k_0 \cdot e)/c \le 2\theta_c$ and $|k_0| \le \theta_c$ (increasing the value of $\theta_c$ if necessary, such a vector always exists). Then we have that

$$ \forall (t_0, x_0) \in \mathbb{R} \times \mathbb{R}^N, \quad \max_{1 \le j \le d} \sup_{(t,x) \in P_1} u_j(t, x) \le \hat{\kappa}_{\theta_c} u_i \left( t_0 + \frac{k_0 \cdot e}{c}, x_0 + k_0 \right), \quad i = 1, \dots, d. $$

Using the pulsating properties (1.6) of the wave $u$ together with (3.6), one gets

$$ \sup_{(t_0, x_0) \in \mathbb{R} \times \mathbb{R}^N} \left( \frac{|\partial_t u_i(t_0, x_0)|}{u_i(t_0, x_0)} + \frac{|\nabla u_i(t_0, x_0)|}{u_i(t_0, x_0)} \right) \le \hat{\kappa}_{\theta_c} C < +\infty, \quad i = 1, \dots, d. \quad (3.7) $$

Observe from (3.1) that

$$ \frac{-\partial_s U_i(s, x)}{U_i(s, x)} = \frac{\partial_t u_i((x \cdot e - s)/c, x)}{c u_i((x \cdot e - s)/c, x)}, \quad i = 1, \dots, d. $$

Hence (3.7) ensures that the function $\partial_s U_i/U_i$ for $i \in \{1, \dots, d\}$ is globally bounded in $\mathbb{R} \times \mathbb{R}^N$. As a consequence, all the quantities $\mu_i$ and $\Lambda_i$, defined in (3.3), are real numbers.

Our next aim is to prove that both $\mu$ and $\Lambda$, defined in (3.2), are positive. By (1.6), the function $U: \mathbb{R} \times \mathbb{R}^N \to \mathbb{R}^d$ defined by (3.1) is $\mathbb{Z}^N$-periodic with respect to the second variable. Moreover, it satisfies

$$ U(+\infty, \cdot) = 0 \quad \text{and} \quad U(s, x) \gg 0 \quad \text{for all } (s, x) \in \mathbb{R} \times \mathbb{R}^N. \quad (3.8) $$

From the definition of $\mu$ in (3.2) and using the periodicity of $U$, there exist $i_0 \in \{1, \dots, d\}$ and a sequence $\{(s_n, x_n)\}_{n \in \mathbb{N}}$ such that

$$ x_n \in \mathbb{T}^N, s_n \to +\infty \quad \text{and} \quad \frac{-\partial_s U_{i_0}(s_n, x_n)}{U_{i_0}(s_n, x_n)} \to \underline{\mu} \quad \text{as } n \to \infty. \quad (3.9) $$

By the compactness of $\mathbb{T}^N$, one may assume, up to extracting a subsequence, that $x_n \to x_\infty \in \mathbb{T}^N$ as $n \to \infty$. Define

$$ v^n(t, x) = \frac{u(t + t_n, x)}{u_{i_0}(t_n, x_n)} \quad \text{with} \quad t_n = \frac{x_n \cdot e - s_n}{c}. \quad (3.10) $$
For each $n \in \mathbb{N}$, each component of $v^n$ satisfies

$$ \partial_t v_i^n(t, x) + L^i v_i^n(t, x) - \frac{f_i(x, u(t+t_n, x))}{u_i(t+t_n, x)} v_i^n(t, x) = 0, \quad i = 1, \dots, d, \qquad (3.11) $$

where the operator $L^i$ is given by (1.4). Moreover, since $ct_n = x_n \cdot e - s_n \to -\infty$ as $n \to \infty$ by (3.9), it follows from (1.6) and (1.7) that

$$ \lim_{n \to \infty} u(t+t_n, x) = \mathbf{0} \quad \text{locally uniformly for } (t, x) \in \mathbb{R} \times \mathbb{R}^N. \qquad (3.12) $$

Now, using Taylor's expansion for each $f_i$ and by assumption (N1), one has

$$ \frac{f_i(x, u(t+t_n, x))}{u_i(t+t_n, x)} = \frac{\partial f_i(x, \mathbf{0})}{\partial u_1} \frac{u_1(t+t_n, x)}{u_i(t+t_n, x)} + \dots + \frac{\partial f_i(x, \mathbf{0})}{\partial u_i} + \dots + \frac{\partial f_i(x, \mathbf{0})}{\partial u_d} \frac{u_d(t+t_n, x)}{u_i(t+t_n, x)} \\ + \frac{1}{2} \frac{u^T(t+t_n, x) D_u^2 f_i(x, \mathbf{0}) u(t+t_n, x)}{u_i(t+t_n, x)} + \frac{o(|u(t+t_n, x)|^2)}{u_i(t+t_n, x)}. $$

Recalling the definition (3.10) of $v_n$, the above formula applied in (3.11) then yields

$$ \partial_t v_i^n + L^i v_i^n - \sum_{j=1}^{d} \frac{\partial f_i(x, \mathbf{0})}{\partial u_j} v_j^n = \frac{1}{2} \frac{u^T(t+t_n, x) D_u^2 f_i(x, \mathbf{0}) u(t+t_n, x)}{u_{i0}(t_n, x_n)} + \frac{o(|u(t+t_n, x)|^2)}{u_{i0}(t_n, x_n)} $$

for all $i \in \{1, \dots, d\}$. Note also that from (3.7), the functions $v^n(t, x) = u(t+t_n, x)/u_{i0}(t_n, x_n)$ are locally bounded in $\mathbb{R} \times \mathbb{R}^N$. Combining this fact with (3.12), and from standard parabolic estimates, $v^n$ converges in $C_{\text{loc}}^{1,2}(\mathbb{R} \times \mathbb{R}^N)^d$, up to a diagonal extraction of a subsequence, to a nonnegative entire solution $v^\infty$ to the system

$$ \partial_t v^\infty(t, x) + L v^\infty(t, x) - D_u f(x, \mathbf{0}) v^\infty(t, x) = \mathbf{0} \quad \text{for } (t, x) \in \mathbb{R} \times \mathbb{R}^N. \qquad (3.13) $$

Since $u$ is a pulsating wave according to Definition 1.2, the function $v^n$ defined in (3.10) satisfies the same pulsating relations as in (1.6) and $v^n \gg \mathbf{0}$. Therefore, the above convergence also yields $v^\infty \ge \mathbf{0}$ and

$$ \forall (t, x, k) \in \mathbb{R} \times \mathbb{R}^N \times \mathbb{Z}^N, \quad v^\infty \left( t + \frac{k \cdot e}{c}, x \right) = v^\infty(t, x - k). \qquad (3.14) $$

In particular, since $D_u f(x, \mathbf{0})$ is cooperative, we have

$$ \left\{ \begin{array}{l} \partial_t v_{i0}^\infty + L^{i0} v_{i0}^\infty - \frac{\partial f_{i0}(x, \mathbf{0})}{\partial u_{i0}} v_{i0}^\infty = \sum_{j=1, j \neq i0}^{d} \frac{\partial f_{i0}(x, \mathbf{0})}{\partial u_j} v_j^\infty \ge 0 \quad \text{in } \mathbb{R} \times \mathbb{R}^N, \\ \forall (t, x) \in \mathbb{R} \times \mathbb{R}^N, \quad v_{i0}^\infty(t, x) \ge 0, \quad v_{i0}^\infty(0, x_\infty) = 1, \end{array} \right. $$

where the operator $L^{i0}$ is given by (1.4) with $i = i_0$. It then follows from the strong maximum principle that $v_{i0}^\infty$ is positive everywhere in $\mathbb{R} \times \mathbb{R}^N$. Since $D_u f(x, \mathbf{0})$ is also fully coupled, applying the Harnack inequality from Theorem B.2 together with the pulsating property of $v^\infty$, one can easily deduce that all other components of $v^\infty$ are also positive everywhere, that is, $v^\infty(t, x) \gg \mathbf{0}$ for all $(t, x) \in \mathbb{R} \times \mathbb{R}^N$.

On the other hand, from (3.10), one has, for each $i \in \{1, \dots, d\}$,

$$ \frac{\partial_t v_i^n(t, x)}{v_i^n(t, x)} = \frac{\partial_t u_i(t+t_n, x)}{u_i(t+t_n, x)} = c \times \frac{-\partial_s U_i(x \cdot e - c(t+t_n), x)}{U_i(x \cdot e - c(t+t_n), x)}, \quad (t, x) \in \mathbb{R} \times \mathbb{R}^N. $$
Since $ct_n = x_n \cdot e - s_n \to -\infty$ as $n \to \infty$, we infer from the definition of the quantity $\underline{\mu}$ in (3.2) and (3.3) that for all $(t, x) \in \mathbb{R} \times \mathbb{R}^N$,

$$ \left( \frac{\partial_t v_1^\infty(t, x)}{v_1^\infty(t, x)}, \dots, \frac{\partial_t v_d^\infty(t, x)}{v_d^\infty(t, x)} \right) \geq c\underline{\mu} \mathbf{1} \quad \text{if } c > 0, \qquad (3.15) \\ \left( \frac{\partial_t v_1^\infty(t, x)}{v_1^\infty(t, x)}, \dots, \frac{\partial_t v_d^\infty(t, x)}{v_d^\infty(t, x)} \right) \leq c\underline{\mu} \mathbf{1} \quad \text{if } c < 0, $$

while $\partial_t v_{i_0}^\infty(0, x_\infty)/v_{i_0}^\infty(0, x_\infty) = c\underline{\mu}$ due to (3.9). Let us now prove that $\partial_t v^\infty = c\underline{\mu}v^\infty$. For that purpose, we define a function $w = (w_1, \dots, w_d)^T$ for $(t, x) \in \mathbb{R} \times \mathbb{R}^N$ as

$$ w(t, x) = \partial_t v^\infty(t, x) - c\underline{\mu}v^\infty(t, x). $$

It is straightforward to check that $w$ satisfies the same equation (3.13) as $v^\infty$. Due to $v^\infty \gg 0$ and (3.15), one has that

$$ \forall (t, x) \in \mathbb{R} \times \mathbb{R}^N, \quad w(t, x) \geq \mathbf{0} \quad (\text{resp. } \leq \mathbf{0}) \quad \text{if } c > 0 \quad (\text{resp. if } c < 0). \qquad (3.16) $$

Note also that

$$ w_{i_0}(0, x_\infty) = \partial_t v_{i_0}^\infty(0, x_\infty) - c\underline{\mu}v_{i_0}^\infty(0, x_\infty) = 0. $$

Therefore, $w_{i_0}$ reaches its minimum or maximum (depending on the sign of $c$) 0 at $(0, x_\infty)$. From the strong maximum principle and using the pulsating relations of $v_{i_0}^\infty$, one can deduce that $w_{i_0}(t, x) = 0$ for all $(t, x) \in \mathbb{R} \times \mathbb{R}^N$. As argued above, the fact that $D_u f(x, \mathbf{0})$ is fully coupled further yields $w = \mathbf{0}$.

In other words, one gets $\partial_t v^\infty = c\underline{\mu}v^\infty$. Since $v^\infty$ is a solution of (3.13), it can be written as

$$ v^\infty(t, x) = e^{-\underline{\mu}(x \cdot e - ct)} \psi(x). $$

where $\psi := (\psi_1, \dots, \psi_d)$ satisfies

$$ -L_{\underline{\mu}e}\psi + D_u f(x, \mathbf{0})\psi = c\underline{\mu}\psi \quad \text{in } \mathbb{T}^N. $$

Here the operator $L_{\underline{\mu}e} := \text{diag}(L_{\underline{\mu}e}^1, \dots, L_{\underline{\mu}e}^d)$ is given by (1.8) with $\lambda = \underline{\mu}$. By (3.14) and $v^\infty \gg \mathbf{0}$, it is easy to verify that $\psi$ is $\mathbb{Z}^N$-periodic and $\psi \gg \mathbf{0}$. We then infer from Proposition 2.2 that

$$ -k(\underline{\mu}, e) = c\underline{\mu}. $$

As for the quantity $\bar{\Lambda}$ defined in (3.2), replacing $\underline{\mu}$ by $\bar{\Lambda}$ and reversing inequalities in (3.15), one can similarly repeat the above proof to get

$$ -k(\bar{\Lambda}, e) = c\bar{\Lambda}. $$

Since $\lambda \mapsto k(\lambda, e)$ is strictly concave, $k(0, e) = \lambda_1^p < 0$ and $c \neq 0$, it follows that $\underline{\mu}$ and $\bar{\Lambda}$ cannot be zero and have the same sign. From the definitions of $\Lambda_i$ in (3.3) and $\bar{\Lambda}$ in (3.2) together with (3.8), we conclude that both $\underline{\mu}$ and $\bar{\Lambda}$ are positive.

Finally, from Proposition 2.3, one has

$$ c = \frac{-k(\bar{\Lambda}, e)}{\bar{\Lambda}} = \frac{-k(\underline{\mu}, e)}{\underline{\mu}} \geq \min_{\lambda > 0} \frac{-k(\lambda, e)}{\lambda} = c^*(e). $$

This ends the proof of Proposition 3.1. $\square$
## 3.2 Monotonicity of the wave with respect to time

We prove here that if the nonlinearity is subhomogeneous, then any pulsating travelling wave of (1.1) is strictly monotonic in time.

The following result corresponds to the second part of Theorem 1.1.

**Proposition 3.2.** Assume that (1.15) holds. Then any pulsating wave $(c, u)$ of (1.1) is componentwise increasing in $t$ if $c > 0$ and decreasing if $c < 0$. Moreover, one has $c \partial_t u(t, x) \gg 0$ for all $(t, x) \in \mathbb{R} \times \mathbb{R}^N$.

*Proof.* The proof follows from the arguments of Berestycki et al. [8, Theorem 1.2] in which the monotonicity of pulsating fronts for strong KPP equations was proved.

Let $u(t, x)$ be a pulsating travelling wave of (1.1) propagating along the direction $e \in \mathbb{S}^{N-1}$ with speed $c \neq 0$. With the change of variables of Proposition 3.1, we consider the function $U = (U_1, \dots, U_d)^T(s, x)$ for $(s, x) \in \mathbb{R} \times \mathbb{T}^N$ defined by

$$U(s, x) = u(t, x) \quad \text{with} \quad t = \frac{x \cdot e - s}{c} \in \mathbb{R}.$$

Then it is continuous and positive. By Proposition 3.1, each component of $U$ satisfies

$$\liminf_{s \to +\infty} \inf_{x \in \mathbb{T}^N} \frac{-\partial_s U_i(s, x)}{U_i(s, x)} > 0, \quad i = 1, \dots, d,$$

and thus there exists $\bar{s} \in \mathbb{R}$ such that

$$\forall s \ge \bar{s}, \forall x \in \mathbb{T}^N, \partial_s U_i(s, x) < 0, \quad i = 1, \dots, d.$$

Moreover, since the limiting condition (1.7) for $u(t, x)$ is equivalent to

$$\lim_{s \to +\infty} U(s, x) = \mathbf{0} \quad \text{and} \quad \liminf_{s \to -\infty} U(s, x) \gg \mathbf{0} \quad \text{uniformly in } \mathbb{T}^N, \qquad (3.17)$$

then by continuity, we have

$$\inf_{s \le \bar{s}, x \in \mathbb{T}^N} U(s, x) \gg \mathbf{0}. \qquad (3.18)$$

As a consequence, there exists some $R \in \mathbb{R}$ such that $R \ge \bar{s}$ and

$$\forall \sigma \ge 0, \forall s \ge R, \forall x \in \mathbb{T}^N, U^\sigma(s, x) := U(s + \sigma, x) \le U(s, x). \qquad (3.19)$$

One can assume $R > 0$. In the sequel, fix any $\sigma \ge 0$ and set

$$\gamma^* := \sup \{\gamma \ge 0 \mid \forall (s, x) \in (-\infty, R] \times \mathbb{T}^N, \gamma U^\sigma(s, x) < U(s, x)\}.$$

Since $U$ is bounded and $\inf_{(-\infty, R) \times \mathbb{T}^N} U^\sigma \gg \mathbf{0}$ by (3.18), the quantity $\gamma^*$ is finite and positive.

One has to prove that $\gamma^* \ge 1$. To this end, we argue by contradiction and assume $\gamma^* < 1$. Since $U(-\infty, \cdot) \gg \mathbf{0}$ by (3.17) and from the definition of $\gamma^*$, there exist $(s_0, x_0) \in (-\infty, R] \times \mathbb{T}^N$ and $i_0 \in \{1, \dots, d\}$ such that

$$\gamma^* U^\sigma(s_0, x_0) \le U(s_0, x_0) \quad \text{and} \quad \gamma^* U_{i_0}^\sigma(s_0, x_0) = U_{i_0}(s_0, x_0). \qquad (3.20)$$

By continuity, we have

$$\gamma^* U^\sigma(s, x) \le U(s, x) \quad \text{for all } (s, x) \in (-\infty, R] \times \mathbb{T}^N,$$

and

$$\gamma^* U^\sigma(s, x) \le U(s, x) \quad \text{for all } (s, x) \in [R, +\infty) \times \mathbb{T}^N$$
due to $\gamma^* < 1$ and (3.19). Coming back to the original variables $(t, x) \in \mathbb{R} \times \mathbb{R}^N$, we define a function $w = (w_1, \dots, w_d)^T(t, x)$ as

$$w(t, x) = U(x \cdot e - ct, x) - \gamma^* U^\sigma(x \cdot e - ct, x). \quad (3.21)$$

Then, $w(t, x) \ge 0$ for all $(t, x) \in \mathbb{R} \times \mathbb{R}^N$. Moreover, $w$ satisfies

$$\partial_t w(t, x) + L w(t, x) = f(x, U(x \cdot e - ct, x)) - \gamma^* f(x, U^\sigma(x \cdot e - ct, x)).$$

From assumption (1.15), one has $\gamma^* f(x, U^\sigma) \le f(x, \gamma^* U^\sigma)$ due to $0 < \gamma^* < 1$ and thus

$$\partial_t w + L w = f(x, U) - \gamma^* f(x, U^\sigma) \ge f(x, U) - f(x, \gamma^* U^\sigma).$$

Let us now examine the differential inequality satisfied by the $i_0$-th component of $w$:

$$\partial_t w_{i_0} + L^{i_0} w_{i_0} \ge f_{i_0}(x, U) - f_{i_0}(x, \gamma^* U^\sigma),$$

where $L^{i_0}$ is given by (1.4) with $i = i_0$. Since $f$ is quasimonotone and $w \ge 0$, then

$$f_{i_0}(x, U) - f_{i_0}(x, \gamma^* U^\sigma) = \sum_{j=1}^{d} h_{i_0 j}(t, x) w_j \ge h_{i_0 i_0}(t, x) w_{i_0},$$

where one has set

$$h_{i_0 j}(t, x) := \int_0^1 \frac{\partial f_{i_0}}{\partial u_j}(x, \gamma^* U_1^\sigma(x \cdot e - ct, x) + s w_1(t, x), \dots, \gamma^* U_d^\sigma(x \cdot e - ct, x) + s w_d(t, x)) ds.$$

Therefore, we obtain that

$$\forall (t, x) \in \mathbb{R} \times \mathbb{R}^N, \quad \partial_t w_{i_0}(t, x) + L^{i_0} w_{i_0}(t, x) - h_{i_0 i_0}(t, x) w_{i_0}(t, x) \ge 0.$$

Setting $t_0 := (x_0 \cdot e - s_0)/c$, we then infer from (3.20) and (3.21) that $w_{i_0}(t_0, x_0) = 0$. As $w_{i_0}$ is nonnegative in $\mathbb{R} \times \mathbb{R}^N$, it follows from the strong maximum principle that

$$\forall t \le t_0, \forall x \in \mathbb{R}^N, \quad w_{i_0}(t, x) = 0.$$

Furthermore, since $U^\sigma(x \cdot e - ct, x)$ and $U(x \cdot e - ct, x)$ are $\mathbb{Z}^N$-periodic in the second variable, one can easily check that the function $w$ defined by (3.21) satisfies

$$w\left(t + \frac{k \cdot e}{c}, x\right) = w(t, x - k) \quad \text{for all } (t, x, k) \in \mathbb{R} \times \mathbb{R}^N \times \mathbb{Z}^N$$

and thus $w_{i_0}(t, x) = 0$ for all $(t, x) \in \mathbb{R} \times \mathbb{R}^N$. However, since $U^\sigma$ has the same asymptotics as in (3.17) and due to $\gamma^* < 1$, then there exists $\epsilon > 0$ independent of $\sigma$ such that

$$\forall x \in \mathbb{R}^N, \quad \liminf_{t \to +\infty} w_{i_0}(t, x) = \liminf_{t \to +\infty} (U_{i_0} - \gamma^* U_{i_0}^\sigma)(x \cdot e - ct, x) = (1 - \gamma^*) \epsilon > 0 \quad \text{if } c > 0,$$

$$\forall x \in \mathbb{R}^N, \quad \liminf_{t \to -\infty} w_{i_0}(t, x) = \liminf_{t \to -\infty} (U_{i_0} - \gamma^* U_{i_0}^\sigma)(x \cdot e - ct, x) = (1 - \gamma^*) \epsilon > 0 \quad \text{if } c < 0.$$

This leads to a contradiction. Consequently, one gets $\gamma^* \ge 1$. From the definition of $\gamma^*$ and by (3.19), one has that

$$\forall \sigma \ge 0, \quad U(s, x) \ge U(s + \sigma, x) \quad \text{for all } (s, x) \in \mathbb{R} \times \mathbb{T}^N.$$

Coming back to the original variables $(t, x)$, this means that

$$\forall \sigma \ge 0, \quad u(t, x) \ge u\left(t - \frac{\sigma}{c}, x\right) \quad \text{for all } (t, x) \in \mathbb{R} \times \mathbb{R}^N.$$
We further claim that the inequality is componentwise strict as soon as $\sigma > 0$. To see this, for any given $\sigma > 0$ and $c \neq 0$, we set

$$v(t, x) := u(t, x) - u\left(t - \frac{\sigma}{c}, x\right).$$

From the uniqueness of the corresponding Cauchy problem, applying the strong maximum principle together with assumption (N4) to the linear parabolic system satisfied by $v$, we can deduce that $v \equiv 0$ or $v \gg 0$ in $\mathbb{R} \times \mathbb{R}^N$. If $v \equiv 0$, we obtain in particular that

$$\forall x \in \mathbb{R}^N, \quad u(0, x) = U(x \cdot e, x) = U(x \cdot e + \sigma, x) = u(-\sigma/c, x).$$

which is impossible because of (3.17). We thereby conclude that $u(t, x)$ is componentwise increasing in the $t$ variable if $c > 0$ and decreasing if $c < 0$.

Lastly, by (3.7), the function $w(t, x) := \partial_t u(t, x)$ is globally bounded. If $c > 0$, then $w$ is nonnegative and $w \neq 0$. Since $f$ is quasimonotone, each component of $w$ satisfies

$$\partial_t w_i + L^i w_i - \frac{\partial f_i(x, u)}{\partial u_i} w_i = \sum_{j=1, j \neq i}^{d} \frac{\partial f_i(x, u)}{\partial u_j} w_j \geq 0, \quad i = 1, \dots, d,$$

Note that the coefficients of the zero-order term in the above equation are bounded. With the same arguments as before, one can prove that all components of $w$ are positive everywhere in $\mathbb{R} \times \mathbb{R}^N$, that is, $\partial_t u \gg 0$ if $c > 0$. The case $c < 0$ can be treated similarly. That completes the proof of Proposition 3.2. $\square$

# 4 Existence of pulsating waves

This section is devoted to the proof of Theorem 1.2.

## 4.1 Construction of two sub- and supersolutions

In order to prove the existence results, we shall use appropriate sub- and supersolutions. Let $h = (h_1, \dots, h_d)^T$ be the vector-valued function defined by

$$h(t, x) = e^{-\lambda_c(e)(x \cdot e - ct)} \phi_{\lambda_c(e)}(x) \quad (4.1)$$

where $\phi_{\lambda_c(e)} \in C^2(\mathbb{T}^N)^d$ is the unique principal eigenfunction of (2.2) in the sense of (2.5) corresponding to $\lambda = \lambda_c(e)$, and $\lambda_c(e)$ is defined by (2.7). By Propositions 2.2-2.3 and from (2.2) applied with $H \equiv D_u f(x, 0)$, it is straightforward to check that $h$ is a supersolution of (1.1) for all $c \geq c^*(e)$ due to assumption (1.16). In what follows, we shall construct another supersolution of (1.1) for $c = c^*(e)$ and two subsolutions of (1.1) for $c > c^*(e)$ and $c = c^*(e)$, respectively.

For notational concision, let us temporarily forget the dependence of $k(\lambda, e)$, $\phi_{\lambda_c(e)}$, $\lambda_c(e)$ and $c^*(e)$ on $e \in \mathbb{S}^{N-1}$ in the next three lemmas. Since $f(x, 0) \equiv 0$ and $f \in C^{1, \beta}(\mathbb{T}^N \times B_\infty(0, \sigma))^d$ with $\sigma > 0$ by assumption (N1), then there exists $M > 0$ such that

$$|f(x, u) - D_u f(x, 0)u|_{\infty} \leq M|u|_{\infty}^{1+\beta} \quad \text{for all } (x, u) \in \mathbb{T}^N \times \mathbb{R}^d \text{ with } |u|_{\infty} \leq \sigma. \quad (4.2)$$

Even if it means decreasing $\sigma > 0$, one can assume without loss of generality that

$$0 < \sigma < \hat{\eta} \quad (4.3)$$

where $\hat{\eta}$ is given in assumption (N3). Let us mention that conditions similar to (4.2) have early been formulated for scalar nonlinear functions in [31], but there is a slight difference: We allow in (4.2) that all coordinates of the vector $u$ may not have the same sign. What it comes into play will be clearly understood in the next two lemmas.
**Lemma 4.1.** Let $\omega = (\omega_1, \dots, \omega_d)^T$ be the vector-valued function defined by

$$ \omega(t, x) = e^{-\lambda_c(x \cdot e - ct)} \phi_{\lambda_c}(x) - Ke^{-(\lambda_c+\delta)(x \cdot e - ct)} \phi_{\lambda_c+\delta}(x) \quad (4.4) $$

where $\phi_{\lambda_c}$ and $\phi_{\lambda_c+\delta}$ are the unique principal eigenfunctions of (2.2) in the sense of (2.5) corresponding to $\lambda = \lambda_c$ and $\lambda_c+\delta$ with $\delta > 0$, respectively. Then for any $c > c^*$, there exist $K > 0$ and $\delta \in (0, \min(\beta\lambda_c, \lambda_c^+ - \lambda_c))$, where $\lambda_c^+$ and $\beta$ are given in Proposition 2.3 and (4.2), such that:

* $\sup_{(t,x) \in \mathbb{R} \times \mathbb{R}^N} \omega(t,x) \ll \hat{\eta} \mathbf{1},$

* $\omega(t,x) \ll \mathbf{0}$ for all $(t,x) \in \mathbb{R} \times \mathbb{R}^N$ satisfying $x \cdot e - ct \le 0,$

* there holds

$$ \partial_t \omega + L\omega \le f(x, \omega) $$

on the set

$$ \Omega_+ := \{ (t, x) \in \mathbb{R} \times \mathbb{R}^N \mid \max_{1 \le i \le d} \omega_i(t, x) > 0 \}, \quad (4.5) $$

where $\hat{\eta}$ is given in (N3) and $L = \text{diag}(L^1, \dots, L^d)$ denotes a diagonal matrix of the operators $L^i$ given by (1.4).

*Proof.* Let

$$ a_0 := \min_{1 \le i \le d} \min_{x \in \mathbb{T}^N} (\phi_{\lambda_c})_i(x) \quad \text{and} \quad b_0 := \min_{1 \le i \le d} \min_{x \in \mathbb{T}^N} (\phi_{\lambda_c+\delta})_i(x). \quad (4.6) $$

Since $\mathbf{0} \ll \phi_{\lambda} \le \mathbf{1}$ by (2.5), one has

$$ \begin{aligned} \max_{1 \le i \le d} \omega_i(t, x) &\le e^{-\lambda_c(x \cdot e - ct)} - b_0 K e^{-(\lambda_c + \delta)(x \cdot e - ct)} \\ &\le \left(\frac{1}{b_0 K}\right)^{\frac{\lambda_c}{\delta}} \times \left[ \left(\frac{\lambda_c}{\lambda_c + \delta}\right)^{\frac{\lambda_c}{\delta}} - \left(\frac{\lambda_c}{\lambda_c + \delta}\right)^{\frac{\lambda_c + \delta}{\delta}} \right] \end{aligned} $$

where the second inequality holds since the resulting upper bound is the global maximum of the middle term and it is reached on the line $x \cdot e - ct = s_0$ with $s_0 := \delta^{-1} \ln[b_0 K(\lambda_c + \delta)/\lambda_c]$. By (4.3), we then obtain that

$$ \sup_{(t,x) \in \mathbb{R} \times \mathbb{R}^N} \omega(t,x) \le \sigma \mathbf{1} \ll \hat{\eta} \mathbf{1} \quad \text{if } K > K_1 $$

where

$$ K_1 := \frac{1}{b_0} \left(\frac{1}{\sigma}\right)^{\frac{\delta}{\lambda_c}} \times \left[ \left(\frac{\lambda_c}{\lambda_c + \delta}\right)^{\frac{\lambda_c}{\delta}} - \left(\frac{\lambda_c}{\lambda_c + \delta}\right)^{\frac{\lambda_c + \delta}{\delta}} \right]^{\frac{\delta}{\lambda_c}} > 0. \quad (4.7) $$

Furthermore, in order to match the condition (4.2), we need to search for $K$ such that

$$ \forall (t, x) \in \Omega_+, \quad |\omega(t, x)|_{\infty} = \max_{1 \le i \le d} |\omega_i(t, x)| \le \sigma, \quad (4.8) $$

where $\Omega_+$ is defined by (4.5) and $\sigma$ is given in (4.2). It suffices to choose some $K > K_1$ such that

$$ \forall (t, x) \in \Omega_+, \quad a_0 e^{-\lambda_c(x \cdot e - ct)} - K e^{-(\lambda_c + \delta)(x \cdot e - ct)} \ge -\sigma, $$

where $a_0$ has been defined in (4.6). Observe that $\Omega_+ \subset \{(t, x) | x \cdot e - ct \ge \delta^{-1} \ln(b_0 K)\}$. Then

$$ a_0 e^{-\lambda_c(x \cdot e - ct)} - K e^{-(\lambda_c + \delta)(x \cdot e - ct)} \ge a_0 (b_0 K)^{-\frac{\lambda_c}{\delta}} - K (b_0 K)^{-\frac{\lambda_c + \delta}{\delta}} \quad \text{in } \Omega_+, $$
where $b_0$ has been defined in (4.6). Hence (4.8) holds as long as $K \ge \max(K_1, K_2)$ where $K_1$ is given by (4.7) and one has set

$$K_2 := \left( \frac{1 - a_0 b_0}{\sigma b_0} \right)^{\frac{\delta}{\lambda_c}} > 0. \qquad (4.9)$$

Next, one can easily check that $\omega(t, x) \ll 0$ for any $\delta > 0$ and for all $(t, x) \in \mathbb{R} \times \mathbb{R}^N$ satisfying $x \cdot e - ct \le 0$ as long as $K > 1/b_0$ where $b_0$ is defined in (4.6). On the other hand, direct computations show

$$\partial_t \omega + L\omega = D_u f(x, \mathbf{0})\omega - (c\lambda_c + c\delta + k(\lambda_c + \delta)) K \phi_{\lambda_c + \delta}(x) e^{-(\lambda_c + \delta)(x \cdot e - ct)}.$$

From Propositions 2.2 and 2.3, if $c > c^*$, one has

$$r_{\delta} := c\lambda_{c} + c\delta + k(\lambda_{c} + \delta) > 0 \quad \text{for all } \delta \in (0, \lambda_{c}^{+} - \lambda_{c}).$$

To conclude one has to check that

$$\forall (t, x) \in \Omega_{+}, \quad D_{u}f(x, \mathbf{0})\omega(t, x) - f(x, \omega(t, x)) \le r_{\delta}K\phi_{\lambda_{c}+\delta}(x)e^{-(\lambda_{c}+\delta)(x \cdot e - ct)}. \qquad (4.10)$$

Fix $\delta \in (0, \min\{\beta\lambda_c, \lambda_c^+ - \lambda_c\})$ and $K \ge \max(K_1, K_2, b_0^{-1}, Mb_0^{-1}/r_\delta)$. Due to (4.8), one can then infer from (4.2) and the expression (4.4) of $\omega$ that for all $(t, x) \in \Omega_+$,

$$|f(x, \omega(t, x)) - D_u f(x, \mathbf{0})\omega(t, x)|_{\infty} \le M|\omega(t, x)|_{\infty}^{1+\beta} \le M \left( \max_{1 \le i \le d} (\phi_{\lambda_c})_i(x) \right)^{1+\beta} e^{-\lambda_c(1+\beta)(x \cdot e - ct)}.$$

As a consequence, (4.10) holds if

$$M \left( \max_{1 \le i \le d} (\phi_{\lambda_c})_i(x) \right)^{1+\beta} e^{-\lambda_c(1+\beta)(x \cdot e - ct)} \le r_{\delta} K \min_{1 \le i \le d} (\phi_{\lambda_c+\delta})_i(x) e^{-(\lambda_c+\delta)(x \cdot e - ct)}.$$

With the above choice for the parameters $K$ and $\delta$, we further obtain that all points $(t, x) \in \mathbb{R} \times \mathbb{R}^N$ contained in $\Omega_+$ satisfy $x \cdot e - ct > 0$ and that for all $(t, x) \in \Omega_+$,

$$\begin{align*}
& M \left( \max_{1 \le i \le d} (\phi_{\lambda_c})_i(x) \right)^{1+\beta} e^{-(1+\beta)\lambda_c(x \cdot e - ct)} - r_{\delta} K \min_{1 \le i \le d} (\phi_{\lambda_c+\delta})_i(x) e^{-(\lambda_c+\delta)(x \cdot e - ct)} \\
\le & \left[ M \left( \max_{1 \le i \le d} (\phi_{\lambda_c})_i(x) \right)^{1+\beta} - r_{\delta} K \min_{1 \le i \le d} (\phi_{\lambda_c+\delta})_i(x) \right] e^{-(\lambda_c+\delta)(x \cdot e - ct)} \\
\le & (M - r_{\delta} K b_0) e^{-(\lambda_c+\delta)(x \cdot e - ct)} \\
\le & 0.
\end{align*}$$

This ends the proof of Lemma 4.1. $\square$

Inspired by [31, Proposition 4.5], we now construct a delicate subsolution of (1.1) for $c = c^*$. Recall first that from Proposition 2.3 and the formula (2.8) of $c^*$, we know that $\lambda^*$ is a unique minimizer of the function $(0, \infty) \ni \lambda \mapsto -k(\lambda)/\lambda$. Therefore, by analyticity and strict concavity of $\lambda \mapsto k(\lambda)$, we have

$$k(\lambda^*) + c^*\lambda^* = 0, \quad k'(\lambda^*) + c^* = 0 \quad \text{and} \quad k''(\lambda^*) \le 0. \qquad (4.11)$$

Furthermore, it also holds that

$$r_{\delta}^{*} := (\lambda^{*} + \delta)c^{*} + k(\lambda^{*} + \delta) < 0 \quad \text{for all } \delta > 0. \qquad (4.12)$$
To proceed we define

$$v^*(t, x; \lambda) := e^{-\lambda(x \cdot e - c^*t)} \phi_\lambda(x) \quad \text{for } (t, x) \in \mathbb{R} \times \mathbb{R}^N \text{ and } \lambda \in \mathbb{R}. \qquad (4.13)$$

Since the normalized principal eigenfunction $\phi_\lambda$ of (2.2) in the sense of (2.5) is also analytic with respect to $\lambda \in \mathbb{R}$, we have that

$$\frac{\partial v^*(t, x; \lambda)}{\partial \lambda} = - \left[ (x \cdot e - c^*t) \phi_\lambda(x) - \frac{\partial \phi_\lambda(x)}{\partial \lambda} \right] e^{-\lambda(x \cdot e - c^*t)}, \qquad (4.14)$$

and that $\frac{\partial \phi_\lambda}{\partial \lambda} \in C^2(\mathbb{T}^N)^d$ satisfies

$$L_\lambda \frac{\partial \phi_\lambda}{\partial \lambda} - D_u f(x, \mathbf{0}) \frac{\partial \phi_\lambda}{\partial \lambda} - k(\lambda) \frac{\partial \phi_\lambda}{\partial \lambda} = L'_\lambda \phi_\lambda - k'(\lambda) \phi_\lambda \quad \text{in } \mathbb{T}^N,$$

where $L'_\lambda$ denotes the derivative of the operator $L_\lambda$ given by (1.8) with respect to $\lambda \in \mathbb{R}$.

**Lemma 4.2.** Let $\omega^* = (\omega_1^*, \dots, \omega_d^*)^T$ be the vector-valued function defined by

$$\omega^*(t, x) = \begin{cases} \kappa_1 v^*(t, x; \lambda^* + \delta) - \kappa_2 \left( \left. \frac{\partial v^*(t, x; \lambda)}{\partial \lambda} \right|_{\lambda=\lambda^*} + K v^*(t, x; \lambda^*) \right) & \text{if } x \cdot e - c^*t \ge 0, \\ \mathbf{0} & \text{if } x \cdot e - c^*t < 0, \end{cases} \qquad (4.15)$$

where $v^*$ is defined by (4.13). Then for any fixed $\delta \in (0, \beta\lambda^*)$ where $\beta$ is given in (4.2), there exist $K > 0$, $s_0 > 0$ and $0 < \kappa_1 \le \kappa_2$ (depending on $s_0$) such that:

* $\omega^*(t, x) < \mathbf{0}$ for all $(t, x) \in \mathbb{R} \times \mathbb{R}^N$ with $x \cdot e - c^*t = 0$,

* $\omega^*(t, x) \gg \mathbf{0}$ for all $(t, x) \in \mathbb{R} \times \mathbb{R}^N$ with $x \cdot e - c^*t \ge s_0$,

* $\sup_{(t,x) \in \mathbb{R} \times \mathbb{R}^N} \omega(t, x) \ll \hat{\eta} \mathbf{1}$,

* there holds

$$\partial_t \omega^* + L \omega^* \le f(x, \omega^*)$$

on the set

$$\Omega_+^* := \left\{ (t, x) \in \mathbb{R} \times \mathbb{R}^N \mid \max_{1 \le i \le d} \omega_i^*(t, x) > 0 \right\}, \qquad (4.16)$$

where $\hat{\eta}$ is given in (N3) and $L = \text{diag}(L^1, \dots, L^d)$ denotes a diagonal matrix of the operators $L^i$ given by (1.4).

**Proof.** Since $\phi_{\lambda^*}$ and $\phi_{\lambda^*+\delta}$ are positive, one may fix some $K > 0$ such that

$$\forall x \in \mathbb{T}^N, \quad \max_{1 \le i \le d} \max_{x \in \mathbb{T}^N} (\phi_{\lambda^*+\delta})_i(x) - \left( \left. \frac{\partial \phi_{\lambda^*}(x)}{\partial \lambda} \right|_{\lambda=\lambda^*} \right)_i \le K \min_{1 \le i \le d} \min_{x \in \mathbb{T}^N} (\phi_{\lambda^*})_i(x). \qquad (4.17)$$

From (4.14) applied at $\lambda = \lambda^*$, such a choice for the parameter $K$ in (4.17) then ensures that the function $\omega^*$ given by (4.2) is nonpositive and nonzero whenever $x \cdot e - c^*t = 0$ as soon as $\kappa_2 \ge \kappa_1 > 0$. To proceed let us set, for notational concision,

$$W(t, x) := - \left. \frac{\partial v^*(t, x; \lambda)}{\partial \lambda} \right|_{\lambda=\lambda^*} - K v^*(t, x; \lambda^*) \qquad (4.18)$$

so that

$$\omega^*(t, x) = \kappa_1 v^*(t, x; \lambda^* + \delta) + \kappa_2 W(t, x) \quad \text{if } x \cdot e - c^*t \ge 0. \qquad (4.19)$$
Notice also that for each $\lambda > 0$,

$$ \frac{\partial v^*(t, x; \lambda)}{\partial \lambda} = -(x \cdot e - c^*t)v^*(t, x; \lambda) + o((x \cdot e - c^*t)v^*(t, x; \lambda)) \quad \text{as } x \cdot e - c^*t \to +\infty. \quad (4.20) $$

Applying this estimate at $\lambda = \lambda^*$, we find that the leading term of $W$ is still given by

$$ (x \cdot e - c^*t)e^{-\lambda^*(x \cdot e - c^*t)}\phi_{\lambda^*}(x) \quad \text{as } x \cdot e - c^*t \to +\infty. $$

Thus, as $\delta > 0$ and $\mathbf{0} \ll \phi_{\lambda} \le \mathbf{1}$ (by (2.5)), then there exists $s_0 > 0$ large enough such that

$$ \begin{gathered} v^*(t, x; \lambda^* + \delta) \le e^{-(\lambda^*+\delta)(x \cdot e - c^*t)} \mathbf{1} \le W(t, x) \le 2(x \cdot e - c^*t)e^{-\lambda^*(x \cdot e - c^*t)} \mathbf{1}, \quad (4.21) \\ \text{for all } (t, x) \in \mathbb{R} \times \mathbb{R}^N \text{ satisfying } x \cdot e - c^*t \ge s_0, \end{gathered} $$

whence it holds that

$$ \forall \kappa_1 > 0, \kappa_2 > 0, \quad \omega^*(t, x) = \kappa_1 v^*(t, x; \lambda^* + \delta) + \kappa_2 W(t, x) \gg \mathbf{0} $$

for all $(t, x) \in \mathbb{R} \times \mathbb{R}^N$ satisfying $x \cdot e - c^*t \ge s_0$.

On the other hand, direct computations show that the function $v^*$ given by (4.13) satisfies

$$ \partial_t v^*(t, x; \lambda) + Lv^*(t, x; \lambda) - D_u f(x, \mathbf{0}) v^*(t, x; \lambda) = (c^*\lambda + k(\lambda)) v^*(t, x; \lambda). \quad (4.22) $$

Differentiating both sides of (4.22) with respect to $\lambda$ yields

$$ \begin{aligned} & \partial_t \left( \frac{\partial v^*(t, x; \lambda)}{\partial \lambda} \right) + L \left( \frac{\partial v^*(t, x; \lambda)}{\partial \lambda} \right) - D_u f(x, \mathbf{0}) \frac{\partial v^*(t, x; \lambda)}{\partial \lambda} \\ & = (c^* + k'(\lambda)) v^*(t, x; \lambda) + (k(\lambda) + c^*\lambda) \frac{\partial v^*(t, x; \lambda)}{\partial \lambda}. \end{aligned} $$

It then follows from (4.11) that

$$ \partial_t \left( \left. \frac{\partial v^*(t, x; \lambda)}{\partial \lambda} \right|_{\lambda=\lambda^*} \right) + L \left( \left. \frac{\partial v^*(t, x; \lambda)}{\partial \lambda} \right|_{\lambda=\lambda^*} \right) = D_u f(x, \mathbf{0}) \left. \frac{\partial v^*(t, x; \lambda)}{\partial \lambda} \right|_{\lambda=\lambda^*}. \quad (4.23) $$

Note also that the function $v^*$ given by (4.13) with $\lambda = \lambda^* + \delta$ satisfies

$$ \partial_t v^*(t, x; \lambda^* + \delta) + Lv^*(t, x; \lambda^* + \delta) - D_u f(x, \mathbf{0}) v^*(t, x; \lambda^* + \delta) = r_{\delta}^* v^*(t, x; \lambda^* + \delta) $$

where $r_{\delta}^*$ has been defined in (4.12). As a result, the function $\omega^*$ given by (4.15) when $x \cdot e - c^*t \ge 0$ then satisfies

$$ \partial_t \omega^*(t, x) + L\omega^*(t, x) = D_u f(x, \mathbf{0})\omega^*(t, x) + \kappa_1 r_{\delta}^* v^*(t, x; \lambda^* + \delta). $$

Our next aim is to prove that

$$ \forall (t, x) \in \Omega_+^*, \quad D_u f(x, \mathbf{0})\omega^*(t, x) + \kappa_1 r_{\delta}^* v^*(t, x; \lambda^* + \delta) \le f(x, \omega^*(t, x)), \quad (4.24) $$

where $\Omega_+^*$ is defined by (4.16). Let

$$ \Omega' := \{(t, x) \in \mathbb{R} \times \mathbb{R}^N \mid x \cdot e - c^*t \ge 0\} \quad (4.25) $$

and $\overline{W} := \|W\|_{L^{\infty}(\Omega')}$ where $W$ is given by (4.18). From (4.19), we then find that

$$ \forall (t, x) \in \Omega', \quad \omega^*(t, x) \le (\kappa_1 + \kappa_2 \overline{W}) \mathbf{1}. \quad (4.26) $$
To match the condition (4.2), we take now any

$$0 < \kappa_1 \le \kappa_2 \le \frac{\sigma}{1 + \bar{W}}$$

so that $|\omega^*(t, x)|_\infty \le \sigma$ for all $(t, x) \in \Omega'$ where $\sigma$ is given in (4.2). In particular, due to (4.3), it is straightforward to see from the expression (4.15) of $\omega^*$ that

$$\sup_{(t,x) \in \mathbb{R} \times \mathbb{R}^N} \omega^*(t, x) \ll \eta \mathbf{1}.$$

By (4.2), one also has

$$|f(x, \omega^*(t, x)) - D_u f(x, \mathbf{0}) \omega^*(t, x)|_\infty \le M |\omega^*(t, x)|_\infty^{1+\beta} \quad \text{for all } (t, x) \in \Omega'$$

On the other hand, due to (4.17), all points $(t, x) \in \mathbb{R} \times \mathbb{R}^N$ contained in $\Omega_+^*$ necessarily satisfy $x \cdot e - c^*t > 0$ and thus $\Omega_+^* \subset \Omega'$, where $\Omega_+^*$ and $\Omega'$ are defined by (4.16) and (4.25), respectively. Note furthermore that the constant $r_\delta^*$, given in (4.12), is negative. Consequently, to prove (4.24) it suffices to show that

$$\forall (t, x) \in \Omega', \quad \kappa_1 r_\delta^* v^*(t, x; \lambda^* + \delta) + M \mathbf{1} |\omega^*(t, x)|_\infty^{1+\beta} \le \mathbf{0}.$$

To that aim, we distinguish between two cases. If $x \cdot e - c^*t \ge s_0$, it follows from $\mathbf{0} \ll \phi_{\lambda^*+\delta} \le \mathbf{1}$, (4.12), (4.19)-(4.21) and $\kappa_2 \ge \kappa_1 > 0$ that

$$\begin{aligned} & \kappa_1 r_\delta^* v^*(t, x; \lambda^* + \delta) + M \mathbf{1} |\omega^*(t, x)|_\infty^{1+\beta} \\ & \le \kappa_1 r_\delta^* e^{-(\lambda^*+\delta)(x \cdot e - c^*t)} \mathbf{1} + M \mathbf{1} \left[ 2\kappa_2 \times 2(x \cdot e - c^*t) e^{-\lambda^*(x \cdot e - c^*t)} \right]^{1+\beta} \\ & = \left[ \kappa_1 r_\delta^* e^{-\delta(x \cdot e - c^*t)} + M \kappa_2^{1+\beta} 4^{1+\beta} (x \cdot e - c^*t)^{(1+\beta)} e^{-\beta \lambda^*(x \cdot e - c^*t)} \right] \mathbf{1} e^{-\lambda^*(x \cdot e - c^*t)}, \end{aligned}$$

where $s_0$ is given in (4.21). Since $0 < \delta < \beta\lambda^*$ and $r_\delta^* < 0$ by (4.12), then there exists a constant $C_1 = C_1(s_0) > 0$ such that

$$0 < (-r_\delta^{*-1}) 4^{1+\beta} M (x \cdot e - c^*t)^{(1+\beta)} e^{\delta(x \cdot e - c^*t)} \le C_1 e^{\beta \lambda^* (x \cdot e - c^*t)} \\ \text{for all } (t, x) \in \mathbb{R} \times \mathbb{R}^N \text{ satisfying } x \cdot e - c^*t \ge s_0,$$

whence there holds

$$\kappa_1 r_\delta^* v^*(t, x; \lambda^* + \delta) + M \mathbf{1} |\omega^*(t, x)|_\infty^{1+\beta} \le r_\delta^* (\kappa_1 - C_1 \kappa_2^{1+\beta}) \mathbf{1} e^{-(\lambda^*+\delta)(x \cdot e - c^*t)}.$$

If $0 \le x \cdot e - c^*t < s_0$, it follows from $\mathbf{0} \ll \phi_{\lambda^*+\delta} \le \mathbf{1}$, $\kappa_2 \ge \kappa_1 > 0$, (4.12) and (4.26) that

$$\begin{aligned} & \kappa_1 r_\delta^* v^*(t, x; \lambda^* + \delta) + M \mathbf{1} |\omega^*(t, x)|_\infty^{1+\beta} \\ & \le \kappa_1 r_\delta^* e^{-(\lambda^*+\delta)s_0} \mathbf{1} + M \mathbf{1} [\kappa_2(1 + \bar{W})]^{1+\beta} \\ & \le r_\delta^* (\kappa_1 - C_2 \kappa_2^{1+\beta}) e^{-(\lambda^*+\delta)s_0} \mathbf{1}, \end{aligned}$$

where one has set

$$C_2 := (-r_\delta^{*-1}) M (1 + \bar{W})^{1+\beta} e^{(\lambda^*+\delta)s_0} > 0.$$

In both cases, (4.24) then follows.

Finally, to conclude we set $C(s_0) := \max(C_1, C_2)$ and fix $\kappa_1, \kappa_2$ such that

$$0 < C \kappa_2^{1+\beta} \le \kappa_1 \le \kappa_2 \le \frac{\sigma}{1 + \bar{W}}.$$

Under such a choice and with the parameter $K$ satisfying (4.17), the function $\omega^*$ given by (4.15) then satisfies all the properties stated in Lemma 4.2 and thereby the proof is done. $\square$
**Lemma 4.3.** With the notations of Lemma 4.2, let $h^* = (h_1^*, \dots, h_d^*)^T$ be the vector-valued function given by

$$h^*(t, x) = \begin{cases} Mv^*(t, x; \lambda^*) - \kappa_2 \left. \frac{\partial v^*(t, x; \lambda)}{\partial \lambda} \right|_{\lambda=\lambda^*} & \text{if } x \cdot e - c^*t > 0, \\ \hat{\eta} \mathbf{1} & \text{if } x \cdot e - c^*t \le 0, \end{cases} \quad (4.27)$$

Then there exists a constant $M > 0$ large enough such that

$$\forall (t, x) \in \mathbb{R} \times \mathbb{R}^N, \quad h^*(t, x) \gg \omega^*(t, x), \quad h^*(t, x) \gg \mathbf{0},$$

and $h^* \to \mathbf{0}$ uniformly as $x \cdot e - c^*t \to +\infty$. Moreover, $h^*$ is a classical supersolution of (1.1) whenever $x \cdot e - c^*t > 0$.

*Proof.* Due to $\phi_{\lambda^*} \gg \mathbf{0}$ and $\phi_{\lambda^*+\delta} \gg \mathbf{0}$, under the fixed parameters in (4.15), one can choose $M_1 > 0$ such that

$$(M_1 + \kappa_2 K) \min_{1 \le i \le d} \min_{x \in \mathbb{T}^N} (\phi_{\lambda^*})_i(x) > \kappa_1 \max_{1 \le i \le d} \max_{x \in \mathbb{T}^N} (\phi_{\lambda^*+\delta})_i(x).$$

By Lemma 4.2, one has $|\omega^*(t, x)|_\infty < \hat{\eta}$ for all $(t, x) \in \mathbb{R} \times \mathbb{R}^N$ where $\omega^*$ is given by (4.15). Thus, for any $M > M_1$, one can easily check that $h^*(t, x) \gg \omega^*(t, x)$ for all $(t, x) \in \mathbb{R} \times \mathbb{R}^N$. Thanks to (4.20), one also has $h^*(t, x) \to \mathbf{0}$ uniformly as $x \cdot e - c^*t \to +\infty$.

Now, using (4.22) with $\lambda = \lambda^*$ and (4.23), we obtain that

$$\partial_t h^* + L h^* = D_u f(x, \mathbf{0}) h^*.$$

It then follows from (1.16) that $h^*(t, x)$ is a supersolution of (1.1) for all $(t, x) \in \mathbb{R} \times \mathbb{R}^N$ satisfying $x \cdot e - c^*t > 0$ as soon as $h^* \ge \mathbf{0}$.

It remains to fix a suitable $M > 0$ such that $h^* \gg \mathbf{0}$. As in the proof of Lemma 4.2, we distinguish between two cases. If $x \cdot e - c^*t \ge s_0$ ($s_0 > 0$ is as in Lemma 4.2), then the estimate (4.20) applied at $\lambda = \lambda^*$ yields $h^* \gg \mathbf{0}$ for any $M > 0$. Now, recall from (4.14) that

$$h^*(t, x) = e^{-\lambda^*(x \cdot e - c^*t)} \left[ M \phi_{\lambda^*}(x) + \kappa_2(x \cdot e - c^*t) - \kappa_2 \left. \left( \frac{\partial \phi_{\lambda^*}(x)}{\partial \lambda} \right) \right|_{\lambda=\lambda^*} \right].$$

where $\partial_\lambda \phi_\lambda|_{\lambda=\lambda^*}$ is bounded and $\phi_{\lambda^*} \gg \mathbf{0}$. Thus, if $0 < x \cdot e - c^*t < s_0$, then there exists $M_2 > 0$ large enough such that

$$M_2 \min_{1 \le i \le d} \min_{x \in \mathbb{T}^N} (\phi_{\lambda^*})_i(x) + \kappa_2(x \cdot e - c^*t) \ge \kappa_2 \max_{1 \le i \le d} \max_{x \in \mathbb{T}^N} \left( \left. \frac{\partial \phi_\lambda}{\partial \lambda} \right|_{\lambda=\lambda^*} \right)_i(x)$$

for all $(t, x) \in \mathbb{R} \times \mathbb{R}^N$ satisfying $0 < x \cdot e - c^*t < s_0$. To conclude we fix some $M > \max(M_1, M_2)$ and then Lemma 4.3 follows. $\square$

## 4.2 Description of a direction-dependent coordinate system

Here we introduce a coordinate system that depends on the direction of propagation for pulsating waves.

Let $e \in \mathbb{S}^{N-1}$ be given and denote by $\{e^1, \dots, e^N\}$ the canonical basis of $\mathbb{R}^N$. Consider a linear orthogonal transformation $R \in \mathcal{O}(\mathbb{R}^N)$ such that

$$Re^1 = e \quad \text{and} \quad e^\perp = \operatorname{span}\{Re^i\}_{2}^N = \{x \in \mathbb{R}^N \mid \forall y \in \mathbb{R}^{N-1}, x = R(0, y)^T\}. \quad (4.28)$$
For any $k \in \mathbb{Z}^N$, we write $k = (k \cdot e)e + k_{\perp}$ with $k_{\perp} := k - (k \cdot e)e \in e^{\perp}$. By (4.28), we then have

$$R^{-1}k = \begin{pmatrix} k \cdot e \\ \mathbf{0} \end{pmatrix} + R^{-1}k_{\perp} \quad \text{and} \\ R^{-1}k_{\perp} = \begin{pmatrix} 0 \\ \hat{R}k_{\perp} \end{pmatrix} \in \bigoplus_{i=2}^{N} \mathbb{R}e^i = \{0\} \times \mathbb{R}^{N-1} \tag{4.29}$$

for a linear map $\hat{R} : e^{\perp} \to \mathbb{R}^{N-1}$. Let $u(t, x)$ be a pulsating travelling wave of (1.1) propagating in the direction $e$ with speed $c \neq 0$. We define a function $\tilde{u} = (\tilde{u}_1, \dots, \tilde{u}_d)^T(t, r, y)$ for $(t, r, y) \in \mathbb{R} \times \mathbb{R} \times \mathbb{R}^{N-1}$ as

$$\tilde{u}(t, r, y) = u(t, R(r, y)^T) \quad \text{with } r = x \cdot e. \tag{4.30}$$

From (1.6) and (4.29), we obtain

$$\begin{aligned} \tilde{u} \left( t + \frac{k \cdot e}{c}, r, y \right) &= u \left( t, R(r, y)^T - k \right) \\ &= u \left( t, R(r - k \cdot e, r)^T - R^{-1}k_{\perp} \right) \\ &= \tilde{u} \left( t, r - k \cdot e, y - \hat{R}k_{\perp} \right). \end{aligned}$$

Furthermore, $\tilde{u}$ satisfies the system of parabolic equations

$$\partial_t \tilde{u} + \mathcal{L} \tilde{u} = \tilde{f}(r, y, \tilde{u}), \quad (t, r, y) \in \mathbb{R} \times \mathbb{R} \times \mathbb{R}^{N-1}, \tag{4.31}$$

where the nonlinearity $\tilde{f} = (\tilde{f}_1, \dots, \tilde{f}_d)^T : \mathbb{R} \times \mathbb{R}^{N-1} \times \mathbb{R}^d \to \mathbb{R}^d$ is defined by

$$\tilde{f}(r, y, \tilde{u}) = f(R(r, y)^T, \tilde{u}),$$

and $\mathcal{L} := \text{diag}(\mathcal{L}^1, \dots, \mathcal{L}^d)$ denotes a diagonal matrix of operators with $\mathcal{L}^i$ being a second-order elliptic operator of the form

$$\mathcal{L}^i(r, y) \tilde{u}_i := - \text{tr} \left( \tilde{A}^i(r, y) D_{(r,y)}^2 \tilde{u}_i \right) + \tilde{q}^i(r, y) \cdot \nabla_{(r,y)} \tilde{u}_i, \quad i = 1, \dots, d. \tag{4.32}$$

Here the matrix-valued function $\tilde{A}^i : \mathbb{R} \times \mathbb{R}^{N-1} \to S_N(\mathbb{R})$ is given by

$$\tilde{A}^i(r, y) = R A^i(R(r, y)^T) R^T, \quad i = 1, \dots, d,$$

and the vector-valued function $\tilde{q}^i : \mathbb{R} \times \mathbb{R}^{N-1} \to \mathbb{R}^N$ is given by

$$\tilde{q}^i(r, y) = q^i(R(r, y)^T)R, \quad i = 1, \dots, d.$$

As a conclusion, for each given $e \in S^{N-1}$, the pulsating condition (1.6) for the wave $(c, u)$ of (1.1) can be rewritten as the following condition for the function $\tilde{u}$ given by (4.30):

$$\begin{gathered} \forall k = (k \cdot e)e + k_{\perp} \in \mathbb{Z}^N, \quad \forall (t, r, y) \in \mathbb{R} \times \mathbb{R} \times \mathbb{R}^{N-1}, \\ \tilde{u} \left( t + \frac{k \cdot e}{c}, r, y \right) = \tilde{u} \left( t, r - k \cdot e, y - \hat{R}k_{\perp} \right), \end{gathered} \tag{4.33}$$

and the limiting condition (1.7) is now equivalent to

$$\liminf_{r \to -\infty} \tilde{u}(t, r, y) \gg \mathbf{0} \quad \text{and} \quad \lim_{r \to +\infty} \tilde{u}(t, r, y) = \mathbf{0} \tag{4.34}$$

locally uniformly for $t \in \mathbb{R}$ and uniformly for $y \in \mathbb{R}^{N-1}$, in which $\tilde{u}$ becomes an entire solution to the parabolic system (4.31). As in [13], we call (4.33) the *R-pulsating condition* in the sequel.

Since $f(\cdot, u)$ and all coefficients of (1.3) are $\mathbb{Z}^N$-periodic, it is straightforward to check from (4.28)-(4.29) that $\tilde{f}(\cdot, \cdot, \tilde{u})$ and all coefficients of the operator $\mathcal{L}$ involved in (4.31) satisfy

$$\begin{gathered} \forall k = (k \cdot e)e + k_{\perp} \in \mathbb{Z}^N, \quad \forall (r, y) \in \mathbb{R} \times \mathbb{R}^{N-1}, \\ \tilde{a} \left( r + k \cdot e, y + \hat{R}k_{\perp} \right) = \tilde{a}(r, y), \end{gathered} \tag{4.35}$$

where $\tilde{a} \equiv \tilde{A}^i(r, y), \tilde{q}^i(r, y), \tilde{f}(r, y, \cdot)$ and one has set $\tilde{a}(r, y) = a(x)$ with $x = R(r, y)^T$.
## 4.3 Construction of the wave in each rational direction

Under the framework of section 4.2, we prove in this part the existence of pulsating waves for (1.1) propagating in the direction of any unit vector with rational coordinates.

In this and the next subsection we always assume $\lambda_1 < 0$, where $\lambda_1$ is defined by (1.12), from which the formulae (2.6) and (2.8) imply that

$$c^*(e) > 0 \text{ for all } e \in \mathbb{S}^{N-1}.$$

We also require a key lemma, which has been proved in [13].

**Lemma 4.4.** Let $\zeta \in \mathbb{Q}^N \cap \mathbb{S}^{N-1}$ for some $N \ge 2$ be given. Then one can find an orthogonal basis $\{\zeta^1, \dots, \zeta^N\}$ of $\mathbb{R}^N$ and $N$ constants $\tau_i > 0$ such that

$$\zeta^1 = \zeta, \quad \tau_1 = \min_{k \cdot \zeta > 0, k \in \mathbb{Z}^N} k \cdot \zeta, \qquad (4.36)$$

and

$$\mathbb{Z}^N = \bigoplus_{i=1}^{N} \tau_i \mathbb{Z}^{\zeta^i}, \quad \zeta^i \perp \zeta^j \ (i \neq j), \quad i, j = 1, \dots, N.$$

Now, let $\zeta \in \mathbb{Q}^N \cap \mathbb{S}^{N-1}$ be given. For any $k \in \mathbb{Z}^N$, we decompose $k$ along $\zeta$ and $\zeta^\perp$, that is, $k = (k \cdot \zeta)\zeta + k_\perp$ with $k_\perp \in \zeta^\perp$. By Lemma 4.4, we have

$$k \cdot \zeta = \tau_1 p_1, \quad k_\perp = \sum_{i=2}^{N} \tau_i p_i \zeta^i,$$

where $p_i \in \mathbb{Z}$ for all $1 \le i \le N$. The R-pulsating condition (4.33) then reduces to

$$\left\{ \begin{array}{l} \forall (t, r) \in \mathbb{R}^2, \quad \forall p_1 \in \mathbb{Z}, \quad \tilde{u}\left(t + \frac{p_1 \tau_1}{c}, r, y\right) = \tilde{u}(t, r - p_1 \tau_1, y), \\ \forall y \in \mathbb{R}^{N-1}, \quad \forall l \in \tau_2 \mathbb{Z} \times \dots \times \tau_N \mathbb{Z}, \quad \tilde{u}(t, r, y + l) = \tilde{u}(t, r, y), \end{array} \right. \qquad (4.37)$$

where $r = x \cdot \zeta$ by (4.30). Moreover, the property (4.35) satisfied by $\tilde{f}$ and all coefficients involved in (4.32) reduces to the following periodicity:

$$\begin{array}{l} \forall (r, y) \in \mathbb{R} \times \mathbb{R}^{N-1}, \quad \forall l \in \tau_2 \mathbb{Z} \times \dots \times \tau_N \mathbb{Z}, \\ \tilde{a}(r + \tau_1, y) = \tilde{a}(r, y + l) = \tilde{a}(r, y). \end{array} \qquad (4.38)$$

For notational simplicity, set $\tau' := (\tau_2, \dots, \tau_N)$ and define an $(N-1)$-dimensional torus as

$$\mathbb{T}_{\tau'}^{N-1} = \mathbb{R}/\tau_2\mathbb{Z} \times \dots \times \mathbb{R}/\tau_N\mathbb{Z}.$$

We then conclude from (4.31)-(4.35) that a pulsating travelling wave of (1.1) propagating along the direction $\zeta \in \mathbb{Q}^N \cap \mathbb{S}^{N-1}$ with speed $c \ne 0$ is a solution to the problem

$$\left\{ \begin{array}{ll} \partial_t \tilde{u} + \mathcal{L} \tilde{u} = \tilde{f}(r, y, \tilde{u}) & \text{for } (t, r, y) \in \mathbb{R} \times \mathbb{R} \times \mathbb{T}_{\tau'}^{N-1}, \\ \forall (t, r) \in \mathbb{R}^2, \quad y \in \mathbb{T}_{\tau'}^{N-1}, \quad \tilde{u}\left(t + \frac{\tau_1}{c}, r, y\right) = \tilde{u}(t, r - \tau_1, y), & \\ \mathbf{0} \ll \tilde{u} \le \hat{\eta} \mathbf{1}, \quad \lim_{r \to +\infty} \tilde{u}(t, r, y) = \mathbf{0}, \quad \liminf_{r \to -\infty} \tilde{u}(t, r, y) \gg \mathbf{0}, & \end{array} \right. \qquad (4.39)$$

where the convergence holds locally uniformly for $t \in \mathbb{R}$ and uniformly for $y \in \mathbb{T}_{\tau'}^{N-1}$.

To address (4.39), we follow and extend some ideas taken from [29]. Roughly speaking, we first consider a similar problem posed in a semi-infinite periodic cylinder with a boundary moving at a velocity $c > 0$ equipped with homogeneous Neumann data. This approximated problem will be solved via a fixed-point argument. Afterwards, we pass to the limit to find a one-dimensional front propagation inside the whole cylinder.

The main result of this subsection reads as follows.
**Proposition 4.5.** Let $\zeta \in \mathbb{Q}^N \cap \mathbb{S}^{N-1}$ be given and assume (1.16) holds. Then for every $c \ge c^*(\zeta)$, the problem (4.39) has a classical solution $\tilde{u}_c$.

*Proof.* The proof is divided into three steps. We shall treat two cases $c > c^*(\zeta)$ and $c = c^*(\zeta)$ in parallel, rather than using a usual limiting argument to deal with the latter. In addition, the uniform positivity of solutions as $r \to -\infty$ will be directly derived from a novel lower barrier.

**Step 1.** *Some preliminaries.*

For any $c > c^*(\zeta)$, we define for all $(t, r, y) \in \mathbb{R} \times \mathbb{R} \times \mathbb{R}^{N-1}$:

$$ \tilde{h}(t, r, y) := h(t, R(r, y)^T) = e^{-\lambda_c(r-ct)} \phi_{\lambda_c}(R(r, y)^T), \quad (4.40) $$

$$ \tilde{\omega}(t, r, y) := \omega(t, R(r, y)^T) = h(t, R(r, y)^T) - K e^{-(\lambda_c+\delta)(r-ct)} \phi_{\lambda_c+\delta}(R(r, y)^T), \quad (4.41) $$

where $h$ and $\omega$ are defined by (4.1) and (4.4) (recall that $x = R(r, y)^T$ under the change of variables (4.30)). Similarly, for $c = c^*(\zeta)$, we define

$$ \tilde{h}^*(t, r, y) := h^*(t, R(r, y)^T), \quad \tilde{\omega}^*(t, r, y) := \omega^*(t, R(r, y)^T), $$

where $h^*$ and $\omega^*$ are defined by (4.27) and (4.15). It is clear to see that $\tilde{h}$, $\tilde{\omega}$, $\tilde{h}^*$ and $\tilde{\omega}^*$ satisfy the same pulsating properties of $(t, r)$ and periodicity of $y$ as in (4.37). Observe also that

$$ \forall n \in \mathbb{N}^*, \quad \forall (t, r, y) \in \mathbb{R} \times \mathbb{R} \times \mathbb{T}_{\tau'}^{N-1}, \\ \tilde{\omega}(t, r + n\tau_1, y) \ll \tilde{h}(t, r + n\tau_1, y) = e^{-n\tau_1\lambda_c} \tilde{h}(t, r, y) \ll \tilde{h}(t, r, y), $$

where equality holds since the principal eigenfunction $\phi_\lambda(R(r, y)^T)$ in the direction $\zeta \in \mathbb{Q}^N \cap \mathbb{S}^{N-1}$ satisfies the same periodicity of $(r, y)$ as in (4.38). As for the case $c = c^*(e)$, from (4.13), (4.14) and (4.27), we find that

$$ h^*(t, R(r + n\tau_1, y)^T) = e^{-n\tau_1\lambda^*} h^*(t, R(r, y)^T) + \kappa_2 n\tau_1 e^{-\lambda^*(r+n\tau_1-c^*(\zeta)t)} \phi_{\lambda^*}(R(r, y)^T). $$

Thus, if $r + n\tau_1 - c^*(\zeta)t > 0$, one can choose a constant $\hat{\gamma} = \hat{\gamma}_{\kappa_2, \tau_1} > 0$ such that

$$ \forall n \in \mathbb{N}^*, \quad \forall (t, r, y) \in \mathbb{R}^2 \times \mathbb{T}_{\tau'}^{N-1}, \quad \tilde{h}^*(t, r + n\tau_1, y) \le \hat{\gamma} \tilde{h}^*(t, r, y). $$

By Lemma 4.3, one also has

$$ \forall n \in \mathbb{N}^*, \quad \forall (t, r, y) \in \mathbb{R}^2 \times \mathbb{T}_{\tau'}^{N-1}, \quad \tilde{\omega}^*(t, r + n\tau_1, y) \ll \hat{\gamma} \tilde{h}^*(t, r, y). $$

Consider now the following functions defined for $(t, r, y) \in \mathbb{R} \times \mathbb{R} \times \mathbb{T}_{\tau'}^{N-1}$ as

$$ \bar{\tilde{u}}_c(t, r, y) = \min(\hat{\eta}\mathbf{1}, \tilde{h}(t, r, y)), \quad \tilde{u}_c(t, r, y) = \max\left(\mathbf{0}, \sup_{n \in \mathbb{N}} \tilde{\omega}(t, r + n\tau_1, y)\right), \quad (4.42) $$

$$ \bar{\tilde{u}}_{c^*}(t, r, y) = \min(\hat{\eta}\mathbf{1}, \hat{\gamma}\tilde{h}^*(t, r, y)), \quad \tilde{u}_{c^*}(t, r, y) = \max\left(\mathbf{0}, \sup_{n \in \mathbb{N}} \tilde{\omega}^*(t, r + n\tau_1, y)\right). \quad (4.43) $$

Then, $\bar{\tilde{u}}_c \gg \tilde{u}_c \ge \mathbf{0}$ and $\bar{\tilde{u}}_{c^*} \gg \tilde{u}_{c^*} \ge \mathbf{0}$. Moreover, $\tilde{u}_c$ and $\tilde{u}_{c^*}$ are uniformly positive (that is, bounded from below by a positive constant) as $r \to -\infty$.

**Step 2.** *An approximated problem posed in semi-infinite periodic cylinders.*

Denote by $\lfloor x \rfloor$ the largest integer below $x$ and fix $a < 0$ such that

$$ \forall t \ge 0, \ r \le \lfloor -a/\tau_1 \rfloor + ct, \ y \in \mathbb{T}_{\tau'}^{N-1}, \ \min_{1 \le i \le d} \tilde{h}_i(t, r, y) \ge \hat{\eta} \quad \text{if } c > c^*(\zeta). $$
Note that by Lemma 4.3, the above holds for any $a < 0$ when $c = c^*(\zeta)$. Define

$$ \underline{\tilde{u}}_c^a(t, r, y) = \max \left( \mathbf{0}, \max_{0 \le n \le \lfloor -a/\tau_1 \rfloor} \tilde{\omega}(t, r + n\tau_1, y) \right), \quad (4.44) $$

$$ \underline{\tilde{u}}_{c^*}^a(t, r, y) = \max \left( \mathbf{0}, \max_{0 \le n \le \lfloor -a/\tau_1 \rfloor} \tilde{\omega}^*(t, r + n\tau_1, y) \right). \quad (4.45) $$

For all $c \ge c^*(\zeta)$, it is clear that $\underline{\tilde{u}}_c^a \to \tilde{u}_c$ as $a \to -\infty$ where $\tilde{u}_c$ and $\underline{\tilde{u}}_{c^*}^a$ are defined in (4.42) and (4.43). For brevity we drop the subscripts $c$ and $c^*$ in (4.42)-(4.45) when $t = 0$.

Set $\Sigma_a^+ := (a, \infty) \times \mathbb{T}_{\tau'}^{N-1}$. From the above construction and Lemmas 4.1-4.2, we find that $\bar{u}(0, r, y) \gg \underline{\tilde{u}}_c^a(0, r, y)$, and that $\bar{u}(0, r, y)$ and $\underline{\tilde{u}}_c^a(0, r, y)$ are identically equal to $\hat{\eta}_1$ and $\mathbf{0}$ near the bottom of $\Sigma_a^+$, respectively. Let $E_a$ be the set defined by

$$ E_a = \{ \varphi \in \mathrm{BUC}(\overline{\Sigma_a^+})^d \mid \underline{\tilde{u}}_c^a(0, r, y) \le \varphi(r, y) \le \bar{u}(0, r, y) \}, $$

which is convex and compactly included in $\mathrm{BUC}(\overline{\Sigma_a^+})^d$. We now equip $E_a$ with the uniform convergence topology and set up a fixed-point problem for an approximation of the solution to (4.39) in $E_a$. To do so, fix $c \ge c^*(\zeta)$ and define a map as

$$ \begin{gathered} Q^a: E_a \to \mathrm{BUC}(\overline{\Sigma_a^+})^d, \\ \varphi \mapsto \tilde{u} \left( \frac{\tau_1}{c}, \cdot + \tau_1, \cdot \right), \end{gathered} \quad (4.46) $$

where $\tilde{u}$ is the solution to the initial-boundary value problem

$$ \left\{ \begin{array}{ll} \partial_t \tilde{u} + \mathcal{L} \tilde{u} = \tilde{f}(r, y, \tilde{u}) & \text{for } (t, r, y) \in (0, \infty) \times (a + ct, \infty) \times \mathbb{T}_{\tau'}^{N-1}, \\ \partial_r \tilde{u}(t, a + ct, y) = \mathbf{0} & \text{for } (t, y) \in (0, \infty) \times \mathbb{T}_{\tau'}^{N-1}, \\ \tilde{u}(0, r, y) = \varphi(r, y) & \text{for } (r, y) \in \overline{\Sigma_a^+}. \end{array} \right. \quad (4.47) $$

It is known that for any given $\varphi \in \mathrm{BUC}(\overline{\Sigma_a^+})^d$, problem (4.47) admits a unique solution $\tilde{u}(t; \varphi) \in \mathrm{BUC}([a + ct, \infty) \times \mathbb{T}_{\tau'}^{N-1})^d$ defined on $\mathbb{R}_+$, and it further becomes smooth as $t$ runs.

Our next aim is to check that the set $E_a$ is stable under $Q^a$, and that the map $Q^a$ is compact from $E_a$ into itself. The former follows from the following claim.

**Claim 4.6.** The solution $\tilde{u} = \tilde{u}(t, r, y)$ of (4.47) satisfies

$$ \underline{\tilde{u}}_c^a(t, r, y) \le \tilde{u}(t, r, y) \le \bar{u}_c(t, r, y) \quad \text{in } \mathbb{R}_+ \times [a + ct, \infty) \times \mathbb{T}_{\tau'}^{N-1}, $$

where $\bar{u}_c$ and $\underline{\tilde{u}}_c^a$ for all $c \ge c^*(\zeta)$ are defined in (4.42)-(4.45), respectively.

Before proving this claim let us proceed with our argument. For any given $\varphi \in E_a$, we infer from Claim 4.6 and the pulsating properties of $\underline{\tilde{u}}_c^a$ and $\bar{u}_c$ that the map $Q^a$, defined by (4.46), satisfies

$$ \underline{\tilde{u}}_c^a(0, r, y) = \underline{\tilde{u}}_c^a \left( \frac{\tau_1}{c}, r + \tau_1, y \right) \le Q^a(\varphi)(r, y) \le \bar{u}_c \left( \frac{\tau_1}{c}, r + \tau_1, y \right) = \bar{u}(0, r, y) $$

for all $(r, y) \in \overline{\Sigma_a^+}$, that is, $Q^a(\varphi) \in E_a$. The latter compactness property is standard, owing to the smoothing effect of the uniform parabolic operators. For completeness we include its proof. Take $\{\varphi^n\}_{n \in \mathbb{N}} \subset E_a$ and set $\tilde{u}^n := Q^a(\varphi^n)$. Due to $\tau_1/c > 0$, the classical parabolic estimates yield that for every fixed $M > 0$, $\{\tilde{u}^n(\tau_1/c, \cdot + \tau_1, \cdot)\}_{n \in \mathbb{N}}$ is uniformly bounded in $C^{2+\alpha}([a, M] \times \mathbb{T}_{\tau'}^{N-1})^d$. Hence one may extract (using a diagonal process) a subsequence of $\{\tilde{u}^n\}$, not relabelled, which converges in $C_{\mathrm{loc}}^2([a, \infty) \times \mathbb{T}_{\tau'}^{N-1})^d$ to a vector-valued function $\tilde{u}$ satisfying (4.46) and remaining
in $E_a$ due to Claim 4.6. Furthermore, the convergence is uniform for $r \in [a, \infty)$. Indeed, since $\tilde{u}(0, r, \cdot) \to \mathbf{0}$ as $r \to +\infty$, by continuity, for any given $\epsilon > 0$ and $M$ sufficiently large, there exists $n_0 \in \mathbb{N}$ independent of $M$ such that

$$ \forall n \ge n_0, \forall r \ge M, \forall y \in \mathbb{T}_{\tau'}^{N-1}, |\tilde{u}^n(\tau_1/c, r + \tau_1, y)|_{\infty} \le \epsilon. $$

Thus, the map $Q^a$ is compact from $E_a$ into itself.

As a consequence, the set $\overline{Q^a(E_a)}$ is compact in $E_a$. We then conclude from the Schauder fixed-point theorem that $Q^a$ has a fixed point in $E_a$, denoted as $\tilde{u}^a$ satisfying

$$ Q^a(\tilde{u}^a) = \tilde{u}^a \left( \frac{\tau_1}{c}, \cdot + \tau_1, \cdot \right) \in E_a. $$

**Step 3.** *Passage to the limit in the infinite cylinders.*

Take any sequence $\{a_n\}_{n \in \mathbb{N}} \subset (-\infty, 0)$ such that $a_n \to -\infty$ as $n \to \infty$. Let

$$ Q^{a_n}(\tilde{u}^{a_n})(r, y) = \tilde{u}^{a_n} \left( \frac{\tau_1}{c}, r + \tau_1, y \right) \quad \text{for } (r, y) \in \overline{\Sigma_{a_n}^+} := [a_n, \infty) \times \mathbb{T}_{\tau'}^{N-1} $$

be the solutions to the problem (4.46). Due to parabolic regularity, $Q^{a_n}(\tilde{u}^{a_n})$ for each $n \in \mathbb{N}$ is bounded in $C^{2+\alpha}(\overline{\Sigma_{a_n}^+})^d$. Up to extracting a subsequence, one may assume that

$$ \tilde{u}^{a_n} \left( \frac{\tau_1}{c}, r + \tau_1, \cdot \right) \to \tilde{u}^{\infty} \left( \frac{\tau_1}{c}, r + \tau_1, \cdot \right) \text{ locally uniformly for } r \in \mathbb{R}. $$

Moreover, the positive invariance of $E_{a_n}$ ensures that

$$ \tilde{u}^{\infty} \in E_{\infty} := \left\{ \varphi \in \mathrm{BUC}(\mathbb{R} \times \mathbb{T}_{\tau'}^{N-1})^d \mid \tilde{u}(0, r, y) \le \varphi(r, y) \le \tilde{u}(0, r, y) \right\}, $$

whence $\tilde{u}^{\infty}$ becomes a fixed point of the map $Q^{\infty}$ given by (4.46) with $a = -\infty$, that is, $Q^{\infty}(\tilde{u}^{\infty}) = \tilde{u}^{\infty}$. As a consequence, for each given $c \ge c^*(\zeta)$, the problem (4.39) admits a classical entire solution $\tilde{u} \equiv \tilde{u}^{\infty}(t, r, y)$ satisfying

$$ \left\{ \begin{array}{ll} \tilde{u}^{\infty} \left( t + \frac{\tau_1}{c}, r, y \right) = \tilde{u}^{\infty} \left( t, r - \tau_1, y \right) & \text{for all } (t, r, y) \in \mathbb{R} \times \mathbb{R} \times \mathbb{T}_{\tau'}^{N-1}. \\ \tilde{u}_c \left( t, r, y \right) \le \tilde{u}^{\infty} \left( t, r, y \right) \le \tilde{u}_c \left( t, r, y \right) & \end{array} \right. $$

Furthermore, from our construction of the lower and upper bounds in (4.42)-(4.43), it is straightforward to see that

$$ \lim_{r \to +\infty} \tilde{u}^{\infty}(t, r, y) = \mathbf{0} \quad \text{and} \quad \liminf_{r \to -\infty} \tilde{u}^{\infty}(t, r, y) \gg \mathbf{0} $$

locally uniformly in $r \in \mathbb{R}$ and uniformly in $y \in \mathbb{T}_{\tau'}^{N-1}$.

Eventually, we conclude from the above three steps that for each $c \ge c^*(\zeta)$, the function $\tilde{u}_c := \tilde{u}^{\infty}(t, r, y)$ is the desired solution to (4.39). This ends the proof of Proposition 4.5. $\square$

*Proof of Claim 4.6.* We only deal with the case $c > c^*(\zeta)$. The case $c = c^*(\zeta)$ can be handled similarly. Moreover, it is straightforward to derive the upper estimates for the solution $\tilde{u}$ of problem (4.47) because $\hat{\eta}\mathbf{1}$ and $\hat{h}$ are two global supersolutions to (4.31).

Let us now turn to the proof of the lower estimates. Due to $\tilde{u}(0, r, y) \ge \underline{u}^a(0, r, y) \not\equiv \mathbf{0}$, the strong maximum principle and Hopf lemma together with assumption (N4) yields

$$ \tilde{u}(t, r, y) \gg \mathbf{0} \quad \text{for all } t > 0 \text{ and } (r, y) \in [a + ct, \infty) \times \mathbb{T}_{\tau'}^{N-1}. \quad (4.48) $$
Since $\tilde{f}$ and all coefficients of $\mathcal{L}$ in (4.47) are periodic in $(r, y)$ as stated in (4.38), we infer from Lemma 4.1 that

$$ \partial_t \tilde{u} + \mathcal{L} \tilde{u} - \tilde{f}(r, y, \tilde{u}) \geq \partial_t \tilde{v}^l + \mathcal{L} \tilde{v}^l - \tilde{f}(r, y, \tilde{v}^l) \quad \text{in } \{t > 0\} \cap \Omega_+^l, $$

where one has set (remember that $r = x \cdot \zeta$ by (4.30))

$$ \tilde{v}^l(t, r, y) := \tilde{\omega}(t, r + l\tau_1, y) \quad \text{for each } l \in \{0, 1, \dots, \lfloor -a/\tau_1 \rfloor\}, \text{ and} \\ \Omega_+^l := \{ (t, r, y) \in \mathbb{R} \times \mathbb{R} \times \mathbb{T}_{\tau'}^{N-1} \mid \max_{1 \le i \le d} \tilde{v}_i^l(t, r, y) > 0 \}. \tag{4.49} $$

Notice that the map $l \mapsto \tilde{v}^l$ is decreasing. Due to $a < 0$ and $l \le \lfloor -a/\tau_1 \rfloor$, one also has

$$ \Omega_+^l \subset \Omega^a := \{ (t, r, y) \in \mathbb{R} \times \mathbb{R} \times \mathbb{T}_{\tau'}^{N-1} \mid r - ct \ge a \}, \tag{4.50} $$

which implies that

$$ \tilde{v}^l < 0 \quad \text{in } \Omega^a \setminus \Omega_+^l. \tag{4.51} $$

On the other hand, the above construction (4.49) implies that $\tilde{v}^l \le 0$ on $\partial \Omega_+^l$. Combining (4.48) with (4.50), it then follows that $\tilde{u} \gg \tilde{v}^l$ on $\{t > 0\} \cap \partial \Omega_+^l$. Recalling the choice of initial data in (4.47), we know that for all $(r, y) \in \mathbb{R} \times \mathbb{T}_{\tau'}^{N-1}$ and for each $l \in \mathbb{N}$, $(\tilde{u} - \tilde{v}^l)(0, r, y) \ge (\tilde{u}^a - \tilde{v}^l)(0, r, y) > 0$, where $\tilde{u}^a$ is defined by (4.44). To summarize, we have

$$ \left\{ \begin{array}{l} \partial_t(\tilde{u}_i - \tilde{v}_i^l) + \mathcal{L}^i(\tilde{u}_i - \tilde{v}_i^l) \ge \sum_{j=1}^d \tilde{h}_{ij}^l (\tilde{u}_j - \tilde{v}_j^l) \quad \text{in } \{t > 0\} \cap \Omega_+^l, \\ \tilde{u}_i - \tilde{v}_i^l > 0 \quad \text{on } \{t > 0\} \cap \partial \Omega_+^l, \quad i = 1, \dots, d, \\ \forall (r, y) \in \mathbb{R} \times \mathbb{T}_{\tau'}^{N-1}, \quad (\tilde{u}_i - \tilde{v}_i^l)(0, r, y) \ge (\ne) 0, \quad i = 1, \dots, d, \end{array} \right. \tag{4.52} $$

where the operator $\mathcal{L}^i$ is defined by (4.32) and one has set

$$ \tilde{h}_{ij}^l(t, r, y) := \int_0^1 \frac{\partial \tilde{f}_i}{\partial \tilde{u}_j}(r, y, \tilde{v}_1^l + s\tilde{u}_1, \dots, \tilde{v}_d^l + s\tilde{u}_d)ds, \quad l \in \mathbb{N}, i, j = 1, \dots, d. $$

Note that the functions $\tilde{h}_{ij}^l$ are bounded. Following [51, page 189 in Chap.3] and [6, Lemma 2.2], we next prove that

$$ \forall (t, r, y) \in \{t > 0\} \cap \Omega_+^l, \quad \tilde{u}(t, r, y) \gg \tilde{v}^l(t, r, y) \quad \text{for each } l \in \mathbb{N}. \tag{4.53} $$

To do so, fix $\epsilon > 0$ and define

$$ w^\epsilon(t, r, y) = \tilde{u}(t, r, y) - \tilde{v}^l(t, r, y) + \epsilon \mathbf{1} e^{Mt}, \quad (t, r, y) \in \mathbb{R}_+ \times \mathbb{R} \times \mathbb{T}_{\tau'}^{N-1}, \tag{4.54} $$

where the constant $M > 0$ is chosen such that

$$ M - \sum_{j=1}^{d} \tilde{h}_{ij}^{l}(t, r, y) > 0 \quad \text{for all } i = 1, \dots, d \text{ and } l \in \mathbb{N}. $$

From (4.52) and (4.54), each component of $w^\epsilon$ then satisfies

$$ \left\{ \begin{array}{l} \partial_t w_i^\epsilon + \mathcal{L}^i w_i^\epsilon - \sum_{j=1}^d \tilde{h}_{ij}^l w_i^\epsilon > 0 \quad \text{in } \{t > 0\} \cap \Omega_+^l, \\ w_i^\epsilon \ge \epsilon \quad \text{on } \partial(\{t > 0\} \cap \Omega_+^l), \quad i = 1, \dots, d. \end{array} \right. \tag{4.55} $$
We now claim that $w^\epsilon \gg 0$ in $\{t>0\} \cap \Omega_+^l$ for any $\epsilon > 0$. To prove this we argue by contradiction. Set $\Omega_T := \{0 < t \le T\} \cap \Omega_+^l$ for $T > 0$. Suppose there exist $\epsilon > 0$, $(t, r, y) \in \bar{\Omega}_T$ and, without loss of generality, $i_0 \in \{1, \dots, d\}$ such that $w_{i_0}^\epsilon(t, r, y) = 0$. Notice first that, by (4.4), we can easily check that the function $\tilde{v}^l$ defined in (4.49) satisfies $\max_{1 \le i \le d} \tilde{v}_i^l(t, r, \cdot) \to -\infty$ as $r \to -\infty$ locally uniformly for $(t, l) \in \mathbb{R} \times \mathbb{Z}$, and $\max_{1 \le i \le d} \tilde{v}_i^l(t, r, \cdot) \to 0$ as $r \to +\infty$ locally uniformly in $t$ and uniformly in $l$. From (4.48) and (4.54), we find that the minimum of $w_{i_0}^\epsilon$ over $\bar{\Omega}_T$ cannot be attained at $|r| = \infty$. Consequently, there exists $(t_0, r_0, y_0) \in \bar{\Omega}_T$ such that

$$w_{i_0}^{\epsilon}(t_0, r_0, y_0) = \min_{\bar{\Omega}_T} w_{i_0}^{\epsilon} = 0 \quad \text{and} \quad w_j^{\epsilon} \ge 0 \text{ for } j \ne i_0. \qquad (4.56)$$

Observe from (4.55) that one necessarily has $(t_0, r_0, y_0) \in \Omega_T$. Therefore, $\nabla_{(r,y)} w_{i_0}^\epsilon(t_0, r_0, y_0) = 0$ and $D^2_{(r,y)} w_{i_0}^\epsilon(t_0, r_0, y_0)$ is nonnegative definite. Using the uniform ellipticity and symmetry involved in the operator $\mathcal{L}^{i_0}$, one can further deduce that

$$(\mathcal{L}^{i_0} w_{i_0}^{\epsilon})(t_0, r_0, y_0) \le 0.$$

Combining (4.55), (4.56) and $\tilde{h}_{ij} \ge 0$ for $i \ne j$, we obtain

$$\partial_t w_{i_0}^{\epsilon}(t_0, r_0, y_0) > 0.$$

However, the definition of the minimum in (4.56) yields that for all $t \in [0, t_0]$, if $(t, r_0, y_0) \in \Omega_+^l$, then $w_{i_0}^\epsilon(t, r_0, y_0) \ge w_{i_0}^\epsilon(t_0, r_0, y_0)$ by continuity. Since $t_0 > 0$, for $\sigma > 0$ sufficiently small, one has $(t_0 - \sigma, r_0, y_0) \in \{t > 0\} \cap \Omega_+^l$. Therefore, it is possible to differentiate the inequality to obtain $\partial_t w_{i_0}^\epsilon(t_0, r_0, y_0) \le 0$. This leads to a contradiction.

As a result, we obtain

$$\forall T > 0, \forall \epsilon > 0, \quad w^{\epsilon} = (\tilde{u} - \tilde{v}^{l}) + \epsilon \mathbf{1} e^{Mt} \gg 0 \quad \text{in } \Omega_T,$$

whence it also holds in $\{t > 0\} \cap \Omega_+^l$. Letting $\epsilon \to 0$, then $\tilde{u} - \tilde{v}^l \ge 0$ in $\{t > 0\} \cap \Omega_+^l$. Furthermore, the inequality is strict, that is, $\tilde{u} \gg \tilde{v}^l$ in $\{t > 0\} \cap \Omega_+^l$. Indeed, suppose not, then there exist $(t', r', y') \in \{t \ge 0\} \cap \bar{\Omega}_+$ and $i' \in \{1, \dots, d\}$ such that $(\tilde{u} - \tilde{v})_{i'}(t', r', y') = 0$. Applying the strong maximum principle and Hopf lemma [51, Theorem 7 in Chap.3] to the $i$-th differential inequality of (4.52), one gets $(\tilde{u} - \tilde{v})_{i'} = 0$ in the connected component of $\{0 < t \le t'\} \cap \Omega_+$ containing $(t', r', y')$. Since the matrix $(\tilde{h}_{ij})_{1 \le i, j \le d}$ is fully coupled, one also has $(\tilde{u} - \tilde{v})_j = 0$ for all $j \in \{1, \dots, d\}$ over the same domain and, by continuity, up to the boundary as well. Recalling (4.48) and the definition of $\Omega_+$ in (4.49), a contradiction has been reached.

Finally, we conclude from (4.48)-(4.51) and (4.53) that

$$\forall t \ge 0, \forall (r, y) \in [a + ct, \infty) \times \mathbb{T}_{\tau'}^{N-1}, \quad \tilde{u}(t, r, y) \ge \max \left( \mathbf{0}, \max_{0 \le n \le \lfloor -a/\tau_1 \rfloor} \tilde{\omega}(t, r + n\tau_1, y) \right).$$

Claim 4.6 then follows. $\square$

## 4.4 Completion of the proof of Theorem 1.2

In this subsection, we return to the original coordinate system with the aim of proving the existence of pulsating waves propagating in any given direction.

Proposition 4.5 equivalently says that for each $c \ge c^*(\zeta)$, system (1.1) admits a pulsating travelling wave $u = u(t, x; \zeta)$ propagating in the direction $\zeta \in \mathbb{Q}^N \cap \mathbb{S}^{N-1}$ with speed $c$. Moreover, from the proof of Proposition 4.5, one also has, for all $(t, x) \in \mathbb{R} \times \mathbb{R}^N$,

$$\max \left( \mathbf{0}, \sup_{n \in \mathbb{N}} \omega(t, x + n\tau_1 \zeta) \right) \le u(t, x; \zeta) \le \min \left( \hat{\eta} \mathbf{1}, h(t, x) \right) \quad \text{for } c > c^*(\zeta), \qquad (4.57)$$

$$\max \left( \mathbf{0}, \sup_{n \in \mathbb{N}} \omega^*(t, x + n\tau_1 \zeta) \right) \le u^*(t, x; \zeta) \le \min \left( \hat{\eta} \mathbf{1}, \hat{\gamma} h^*(t, x) \right) \quad \text{for } c = c^*(\zeta), \qquad (4.58)$$
where $h, \omega, \omega^*$ and $h^*$ are defined by (4.1), (4.4), (4.15) and (4.27), and the parameter $\hat{\gamma}$ in (4.58) is the same as in (4.43). Recall that the asymptotics of the waves at infinity were directly obtained from (4.57)-(4.58). However, such lower bounds no longer hold for a general direction $e \in \mathbb{S}^{N-1}$ (more on this later). Furthermore, since $f$ is not assumed to be subhomogeneous in Theorem 1.2, the monotonicity results derived in Proposition 3.2 may not be valid.

On this account, we prove a persistence property for the solution of the Cauchy problem associated with the monotone system (1.1) before proving the existence of pulsating waves propagating in any given direction. This kind of property is known as the *hair-trigger effect*.

**Proposition 4.7 (Asymptotic persistence).** Let $\psi = (\psi_1, \dots, \psi_d)^T \in BUC(\mathbb{R}^N)^d$ be such that $\psi \ge \mathbf{0}$ and $\psi_i \neq 0$ for some $i$. Denote by $u = u(t, x; \psi)$ the solution to (1.1) with the initial datum $u(0, x) = \psi(x)$. Let $\lambda_1$ be the generalized principal eigenvalue defined by (1.12).

Assume $\lambda_1 < 0$. Then there exists $\epsilon > 0$ independent of $\psi$ such that

$$ \liminf_{t \to +\infty} u(t, x; \psi) \geq \epsilon \mathbf{1} \quad \text{locally uniformly for } x \in \mathbb{R}^N. $$

*Proof.* Applying Proposition 2.1 to problem (1.11) yields $\lambda_1^R \to \lambda_1$ as $R \to +\infty$. Therefore, under the assumption $\lambda_1 < 0$, there exists $R_0 > 0$ large enough such that

$$ \forall R \geq R_0, \quad \lambda_1^R < 0. $$

In the sequel, fix some $R > R_0$ and let $\varphi^R$ be the corresponding principal eigenfunction of (1.11), normalized such that $\|\varphi^R\|_{L^\infty(B_R)^d} = 1$. Recall from our assumption (N1) that $f(x, \mathbf{0}) = \mathbf{0}$ and $f \in C^1(\mathbb{T}^N \times B_\infty(\mathbf{0}, \sigma))^d$ with $\sigma > 0$. Thus, one can choose $\kappa_0 \in (0, \sigma]$ small enough such that

$$ \forall (x, \kappa) \in \mathbb{T}^N \times (0, \kappa_0], \quad |D_u f(x, \mathbf{0})\kappa\mathbf{1} - f(x, \kappa\mathbf{1})|_{\infty} \leq -\kappa\lambda_1^R. \quad (4.59) $$

Due to $\mathbf{0} \ll \varphi^R \le \mathbf{1}$, we then have that

$$ \forall \kappa \in (0, \kappa_0], \quad L(\kappa\varphi^R) - f(x, \kappa\varphi^R) = D_u f(x, \mathbf{0})\kappa\varphi^R + \lambda_1^R \kappa\varphi^R - f(x, \kappa\varphi^R) \le \mathbf{0} \quad \text{in } B_R, \quad (4.60) $$

where $L$ denotes a diagonal matrix of the operators given by (1.4).

On the other hand, since the initial datum $\psi$ is assumed to be nonnegative and not identically equal to $\mathbf{0}$, we infer from the strong maximum principle coupled with assumption (N4) that $u(1, x; \psi) \gg \mathbf{0}$ for all $x \in \mathbb{R}^N$. Consequently, possibly by further decreasing $\kappa_0$ in (4.60), $u(1, x; \psi) \gg \kappa_0 \varphi^R(x)$ for all $x \in \bar{B}_R$. For any $\kappa \in (0, \kappa_0]$, let $\psi^\kappa$ be the function defined by

$$ \psi^\kappa(x) = \begin{cases} \kappa\varphi^R(x) & \text{if } x \in B_R, \\ \mathbf{0} & \text{otherwise.} \end{cases} $$

From (4.60), the function $\psi^\kappa$ then becomes a generalized subsolution of the stationary equation associated with (1.1) in $\mathbb{R}^N$. Moreover, one has $u(1, x; \psi) \ge \psi^\kappa(x)$ for all $x \in \mathbb{R}^N$. Consider now the solution $u = \underline{u}(t, x; \psi^\kappa)$ to (1.1) with the initial datum $\psi^\kappa$. By the parabolic comparison principle, we obtain that

$$ \forall (t, x) \in (0, \infty) \times \mathbb{R}^N, \quad \mathbf{0} \ll \underline{u}(t, x; \psi^\kappa) \le u(t+1, x; \psi), \quad (4.61) $$

and that $\underline{u}$ is nondecreasing in $t$. The classical parabolic estimates then imply that $\underline{u}$ converges in $C_{\text{loc}}^2(\mathbb{R}^N)^d$ as $t \to +\infty$ to a bounded stationary solution to (1.1), denoted as $p^\kappa$ satisfying $\mathbf{0} \le p^\kappa(x) \le \hat{\eta}\mathbf{1}$ for all $x \in \mathbb{R}^N$, where $\hat{\eta}$ is given in assumption (N3). Furthermore, due to

$$ p^\kappa(0) \ge \underline{u}(0, 0; \psi^\kappa) = \kappa \varphi^R(0) \gg \mathbf{0}, $$
it follows from the strong maximum principle that $p^\kappa \gg 0$ in $\mathbb{R}^N$ for any $\kappa \in (0, \kappa_0]$, and thus $\min_{x \in \overline{B_R}} p^\kappa(x) \gg 0$. Our next aim is to prove that

$$p^{\kappa}(x) \geq \kappa_0 \varphi^R(x) \quad \text{in } B_R. \qquad (4.62)$$

Note that one can define

$$\gamma^* := \sup \{ \gamma > 0 \mid p^{\kappa} > \gamma \varphi^R \text{ in } \overline{B_R} \}.$$

We claim that $\gamma^* \geq \kappa_0$. Indeed, if $\gamma^* < \kappa_0$, then there exist $x_0 \in \overline{B_R}$ and $i_0 \in \{1, \dots, d\}$ such that

$$\gamma^* \varphi_{i_0}^R(x_0) = p_{i_0}^\kappa(x_0) \quad \text{and} \quad \gamma^* \varphi^R(x) \leq p^\kappa(x) \text{ in } \overline{B_R}.$$

One necessarily has $x_0 \in B_R$ since $\varphi^R \equiv 0$ on $\partial B_R$ while $\min_{\overline{B_R}} p^\kappa \gg 0$. Observe also that due to $\gamma^* < \kappa_0$ and (4.60), $\gamma^* \varphi^R$ is still a subsolution of the stationary equations associated with (1.1) in $B_R$. From the strong elliptic maximum principle coupled with assumption (N4), we deduce that $\gamma^* \varphi^R \equiv p^\kappa$ in $\overline{B_R}$, which contradicts the values of both on $\partial B_R$. We thereby conclude from (4.61) and (4.62) that

$$\liminf_{t \to +\infty} u(t, x; \psi) \geq \kappa_0 \varphi^R(x) \quad \text{in } \overline{B_{R/2}}.$$

Now, recall that $R$ was chosen sufficiently large. Then, $[0, 1]^N \subset B_{R/2}$. Furthermore, since all coefficients of the operator $L$ and $f(\cdot, u)$ are $\mathbb{Z}^N$-periodic, it follows that for all $l \in \mathbb{Z}^N$ and any $\kappa \in (0, \kappa_0]$, $u(t, x-l; \psi^\kappa(\cdot-l))$ satisfies the same system of equations as $u(t, x; \psi^\kappa)$. Note also that the Dirichlet principal eigenvalue $\lambda_1^p$ is invariant under periodic translations of the coefficients. With the same arguments as above, we can obtain that

$$\forall l \in \mathbb{Z}^N, \quad \liminf_{t \to +\infty} u(t, x; \psi) \geq \kappa_0 \varphi^R(x-l) \quad \text{in } \overline{B_{R/2}(l)},$$

where $B_R(l)$ denotes the ball of radius $R$ centered at $l \in \mathbb{Z}^N$. Finally, since $\inf_{B_R} \varphi^R \gg 0$ for all $R > 0$, then there exists a constant $\epsilon > 0$ independent of $\psi$ such that

$$\forall x \in \mathbb{R}^N, \quad \liminf_{t \to +\infty} u(t, x; \psi) \geq \epsilon \mathbf{1},$$

which concludes the proof of Proposition 4.7. $\square$

**Remark 4.8.** Recall from Remark 2.4 that in general only $\lambda_1^p < \lambda_1$ holds, where both quantities are given by (1.10) and (1.12), respectively. It is worthy pointing out that in the case of $\lambda_1^p < 0 \le \lambda_1$, the persistence property relies strongly on the kind of additional conditions to be imposed on the initial data. We refer to [27, 48] for more details in a space-time periodic setting. In particular, if the initial datum is periodic, a proof similar to the one above yields the uniform persistence of solutions associated with (1.1).

Following [13, 14] for the rational approximation to any direction of propagation, we are now in a position to complete the:

*Proof of Theorem 1.2.* Fix any $e \in \mathbb{S}^{N-1}$. Since $\mathbb{Q}^N \cap \mathbb{S}^{N-1}$ is a dense subset of $\mathbb{S}^{N-1}$ (see [16, Lemma 4.1]), there exists a sequence $\{\zeta^m\}_{m \in \mathbb{N}} \subset \mathbb{Q}^N \cap \mathbb{S}^{N-1}$ such that $\zeta^m \to e$ as $m \to +\infty$. Moreover, for each given $c_e > c^*(e)$, we can choose a sequence $\{c^m\}_{m \in \mathbb{N}} \subset (c^*(\zeta^m), \infty)$ such that $c^m \to c_e$ as $m \to +\infty$ since the map $\mathbb{S}^{N-1} \ni e \mapsto c^*(e)$ is continuous. Consider now the sequence of pulsating waves propagating in the direction $\zeta^m$ with speed $c^m$ for (1.1):

$$c^m > c^*(\zeta^m), \quad u^m(t, x) := u(t, x; \zeta^m), \quad m \in \mathbb{N}.$$

For each $m \in \mathbb{N}$, the vector-valued function $u^m$ satisfies

$$u^m \left( t + \frac{k \cdot \zeta^m}{c^m}, x \right) = u^m(t, x-k) \quad \text{for all } (t, x, k) \in \mathbb{R} \times \mathbb{R}^N \times \mathbb{Z}^N,$$
and by (4.57), there holds

$$ \max \left( \mathbf{0}, \sup_{n \in \mathbb{N}} \omega(t, x + n \tau_1^m \zeta^m; \zeta^m) \right) \le u^m(t, x) \le \min(\hat{\eta} \mathbf{1}, h(t, x; \zeta^m)), \quad (4.63) $$

where the parameter $\tau_1^m$ is given by (4.36) with $\zeta = \zeta^m$. From standard parabolic estimates, $\{u^m\}$ converges (up to extracting a subsequence) in $C_{\text{loc}}^{1,2}(\mathbb{R} \times \mathbb{R}^N)^d$ as $m \to +\infty$ to a bounded entire solution $u^\infty$ to (1.1), which is nonnegative and verifies

$$ \forall (t, x, k) \in \mathbb{R} \times \mathbb{R}^N \times \mathbb{Z}^N, \quad u^\infty \left( t + \frac{k \cdot e}{c_e}, x \right) = u^\infty(t, x). \quad (4.64) $$

From statement (iv) of Proposition 2.2 and Proposition 2.3, we find that the functions $h$ and $\omega$ given by (4.1) and (4.4) are continuous in both $e \in \mathbb{S}^{N-1}$ and $c \in (c^*(e), \infty)$ with respect to the uniform topology. Note also that $u^\infty$ is bounded and $c_e > c^*(e)$. We then infer from (4.63) and the periodicity of equations that

$$ \forall (t, x) \in \mathbb{R} \times \mathbb{R}^N, \quad \max(\mathbf{0}, \omega(t, x; e)) \le u^\infty(t, x) \le \min(\hat{\eta} \mathbf{1}, h(t, x; e)), \quad (4.65) $$

where

$$ \begin{aligned} h(t, x; e) &= e^{-\lambda_{c_e}(e)(x \cdot e - c_e t)} \phi_{\lambda_{c_e}(e)}(x), \\ \omega(t, x; e) &= h(t, x; e) - K e^{-(\lambda_{c_e}(e)+\delta)(x \cdot e - c_e t)} \phi_{\lambda_{c_e}(e)+\delta}(x), \end{aligned} $$

with $\lambda_{c_e}(e) > 0$ given by (2.7). Here only such a lower bound stated in (4.65) holds since $\tau_1^m \to 0$ as $m \to +\infty$ by the definition of $\tau_1$ in (4.36). From the strong maximum principle coupled with assumption (N4), one further gets $u^\infty \gg \mathbf{0}$.

It remains to verify the limiting condition (1.7). Using the pulsating property of $u^\infty$ in (4.64) and the positivity of speed $c_e$ (remember that $c^*(e) > 0$ under assumption $\lambda_1 < 0$), it is equivalent to prove that

$$ \lim_{t \to -\infty} u^\infty(t, x) = \mathbf{0} \quad \text{and} \quad \liminf_{t \to +\infty} u^\infty(t, x) \gg \mathbf{0} $$

locally uniformly for $x \in \mathbb{R}^N$. The former follows immediately from (4.65) and the fact that the function $h$ is componentwise decreasing in $t$. To prove the latter limit, we consider a solution $\underline{u} : \mathbb{R}_+ \times \mathbb{R}^N \to \mathbb{R}^d$ to (1.1) starting from $\underline{u}(0, x) = \max(\mathbf{0}, \omega(0, x; e))$. Since $u^\infty$ is a solution to (1.1), the parabolic comparison principle then ensures that

$$ \forall (t, x) \in \mathbb{R}_+ \times \mathbb{R}^N, \quad \underline{u}(t, x) \le u^\infty(t, x) \le \hat{\eta} \mathbf{1}. $$

Applying Proposition 4.7 to $\underline{u}$, we obtain that there exists $T > 0$ large enough such that

$$ \forall t \ge T, \forall R > 0, \quad \inf_{|x| \le R} u^\infty(t, x) \ge \inf_{|x| \le R} \underline{u}(t, x) \ge \epsilon \mathbf{1}, $$

whence the latter limit follows. Thus, the function $u(t, x; e) = u^\infty(t, x)$ is the desired pulsating travelling wave for (1.1) propagating along direction $e$ with speed $c_e$.

For the critical case $c_e = c^*(e)$, let us set

$$ c^m := c^*(\zeta^m), \quad u^m(t, x) := u^*(t, x; \zeta^m), \quad m \in \mathbb{N}. $$

Passing to the limit $m \to +\infty$ as before, we can obtain the existence of a bounded entire solution of (1.1) verifying pulsating conditions as in (4.64) with $c_e = c^*(e)$, and similarly, by (4.58) and Lemma 4.3, there also holds

$$ \forall (t, x) \in \mathbb{R} \times \mathbb{R}^N, \quad \max(\mathbf{0}, \omega^*(t, x; e)) \le u^*(t, x; e) \le \min(\hat{\eta} \mathbf{1}, h^*(t, x; e)), $$
where $\omega^*$ and $h^*$ are given by (4.15) and (4.27). The strong maximum principle coupled with assumption (N4) further yields $u^*(t, x; e) \gg 0$ in $\mathbb{R} \times \mathbb{R}^N$. Furthermore, by Lemma 4.3 and $c^*(e) > 0$, it is straightforward to check that $u^*(t, x) \to 0$ as $t \to -\infty$ locally uniformly for $x \in \mathbb{R}^N$. With the same argument as before, by Proposition 4.7, we can also prove the other limit of $u^*$, that is, $\liminf_{t \to +\infty} u^*(t, x; e) \gg 0$ locally uniformly for $x \in \mathbb{R}^N$.

This ends the proof of Theorem 1.2. □

# 5 Nonexistence of pulsating waves

This section is devoted to the proof of Theorem 1.3. The non-existence results will rest on the following proposition.

**Proposition 5.1 (Extinction).** Let $\psi = (\psi_1, \dots, \psi_d)^T \in BUC(\mathbb{R}^N)^d$ be such that $\psi \ge 0$. Denote by $u = u(t, x; \psi)$ the solution to (1.1) with the initial datum $u(0, x) = \psi(x)$. Let $\lambda_1^p$ be the periodic principal eigenvalue given by (1.10).

Assume that (1.17) holds. If $\lambda_1^p \ge 0$, then

$$ \lim_{t \to +\infty} u(t, x; \psi) = \mathbf{0} \quad \text{uniformly for } x \in \mathbb{R}^N. $$

*Proof.* Consider the solution $\bar{u} = \bar{u}(t, x)$ of (1.1) with the initial datum given by

$$ \bar{u}(0, x) = \mathbf{1} \max(\hat{\eta}, \|\psi\|_{\infty}), $$

where the constant $\hat{\eta}$ is given in (N3). It then follows from the comparison principle that

$$ \forall (t, x) \in \mathbb{R}_+ \times \mathbb{R}^N, \quad \mathbf{0} \le u(t, x; \psi) \le \bar{u}(t, x), \qquad (5.1) $$

and that $\bar{u}$ is nonincreasing in $t$. Moreover, for all $t > 0$, $\bar{u}(t, x)$ is periodic in $x$ because $\bar{u}(0, x)$ is constant. As a consequence, the classical parabolic estimates imply that $\bar{u}$ converges as $t \to +\infty$ in $C^2(\mathbb{R}^N)^d$ to a $\mathbb{Z}^N$-periodic stationary solution to (1.1), denoted as $U$, that is,

$$ \left\{ \begin{array}{ll} LU(x) = f(x, U(x)) & \text{for } x \in \mathbb{R}^N, \\ \forall (x, k) \in \mathbb{R}^N \times \mathbb{Z}^N, & U(x + k) = U(x), \\ \forall x \in \mathbb{R}^N, & \mathbf{0} \le U(x) \le \bar{u}(0, x), \end{array} \right. \qquad (5.2) $$

where $L$ denotes a diagonal matrix of the operators given by (1.4).

Next, we prove that $U = \mathbf{0}$. To do so, let $\phi_0 \in C^2(\mathbb{T}^N)^d$ denote the periodic principal eigenfunction of (1.9) associated with $\lambda_1^p$. Set

$$ \gamma^* := \inf \{ \gamma > 0 \mid \gamma\phi_0 \gg U \text{ in } \mathbb{T}^N \} \ge 0. $$

This quantity is well-defined since $U$ is bounded and $\phi_0 \gg 0$. We claim that $\gamma^* = 0$. Suppose $\gamma^* > 0$. Then, $\gamma^*\phi_0 \ge U$ and there exists $\{x_n\}_{n \in \mathbb{N}} \subset \mathbb{T}^N$ such that, for at least one index $i \in \{1, \dots, d\}$, $(\gamma^*\phi_0 - U)_i(x_n) \to 0$ as $n \to \infty$. By the compactness of $\mathbb{T}^N$, one may assume that $x_n \to x_\infty$ as $n \to \infty$, whence there holds $(\gamma^*\phi_0 - U)_i(x_\infty) = 0$ by continuity. As $\lambda_1^p \ge 0$, we infer from (1.16) and (5.2) that

$$ L(\gamma^*\phi_0 - U) - D_u f(x, \mathbf{0})(\gamma^*\phi_0 - U) \ge \lambda_1^p \gamma^*\phi_0 \ge \mathbf{0} \quad \text{in } \mathbb{T}^N. $$

Hence the strong maximum principle coupled with assumption (N4) yields $\gamma^*\phi_0 \equiv U$. However, since $\phi_0 \gg \mathbf{0}$ and $\gamma^* > 0$, and by assumption (1.17), one has

$$ L(\gamma^*\phi_0) - f(x, \gamma^*\phi_0) = (D_u f(x, \mathbf{0}) + \lambda_1^p) \gamma^*\phi_0 - f(x, \gamma^*\phi_0) \ge (\mathfrak{a}) \mathbf{0} \quad \text{in } \mathbb{T}^N. $$
Recalling (5.2), a contradiction has been achieved. One gets $\gamma^* = 0$ and thereby $U = 0$.

As a consequence, we infer from (5.1) that

$$0 \le \lim_{t \to +\infty} u(t, x; \psi) \le \lim_{t \to +\infty} \bar{u}(t, x) = 0 \quad \text{uniformly for } x \in \mathbb{R}^N.$$

The proof of Proposition 5.1 is complete.

<figure><img src="image_3.png" /></figure>

*Proof of Theorem 1.3.* Assume that $u = u(t, x)$ is a pulsating travelling wave of (1.1) according to Definition 1.2. Consider the solution $u = u(t, x; \hat{\eta}1)$ to (1.1) starting from $u(0, x; \hat{\eta}1) = \hat{\eta}1$. Since $u(t, x)$ is bounded from above by $\hat{\eta}1$, the parabolic comparison principle ensures that for each $t \in \mathbb{R}$ and $s \in \mathbb{R}$, there holds

$$\forall t \ge s, \forall x \in \mathbb{R}^N, \quad 0 \le u(t, x) \le u(t - s, x; \hat{\eta}1).$$

By Proposition 5.1, $u(t - s, x; \hat{\eta}1)$ converges to $0$ as $s \to -\infty$ locally uniformly for $t \in \mathbb{R}$ and uniformly for $x \in \mathbb{R}^N$, which implies that $u(t, x) = 0$. This leads to a contradiction. □

<figure><img src="image_4.png" /></figure>

## Data availability statement

No new data were created or analysed in this study.

## Acknowledgements

L. Deng's research is supported by the China Postdoctoral Science Foundation (No. 2025M773053) and in part by the NSF of China (No.12471203). He would like to thank the last two authors for generous hospitality during his visit to the Université Le Havre Normandie, where this work was initiated. He would also like to express his sincere gratitude to Professors Xuefeng Wang and Xing Liang for supporting his postdoctoral position.

The authors are very grateful to the referees for their valuable comments and suggestions which helped to improve our original manuscript, as well as for pointing out some inspiring problems for further study.

# Appendix

## A Proof of Proposition 2.2

Statement (i) follows from the Krein-Rutman theory. The max-min characterization in (ii) is quite similar to the one by Sweers [56, Remark 3.1] (see also [10]). The analyticity of $\lambda \mapsto k(\lambda, e)$ ($\phi_{\lambda e}$ as well) is standard owing to Kato's perturbation theory (see [37, Section 7.2]). Estimate (2.4) is easy to derive by virtue of formula (2.3). Statement (iv) is a direct consequence of the foregoing properties. Let us now prove the rest of the results.

*Proof of Proposition 2.2 (iii).* Only the concavity remains to be proved. Following [8, Lemma 3.1], let us first show that for any given $e \in \mathbb{S}^{N-1}$ the function $\lambda \mapsto k(\lambda, e)$ is concave in $\mathbb{R}$. One has to check that

$$\forall \hat{\lambda}, \tilde{\lambda} \in \mathbb{R}, \forall r \in (0, 1), \quad k(r\hat{\lambda} + (1-r)\tilde{\lambda}, e) \ge rk(\hat{\lambda}, e) + (1-r)k(\tilde{\lambda}, e). \quad (A.1)$$

Let $E_{\lambda}$ be the set defined by

$$E_{\lambda} = \{ \psi \in C^2(\mathbb{R}^N)^d \mid \psi \gg 0 \text{ and } \psi e^{\lambda \cdot x} \text{ is } \mathbb{Z}^N\text{-periodic} \}. \quad (A.2)$$
Then, formula (2.3) can be rewritten as

$$k(\lambda, e) = \max_{\psi \in E_{\lambda}} \min_{1 \le i \le d} \inf_{x \in \mathbb{T}^N} \frac{e^{i(L-H)\psi}(x)}{\psi_i(x)}, \quad (A.3)$$

where $L = \text{diag}(L^1, \dots, L^d)$ denotes the operator given by (1.4). Let now $\hat{\psi} \in E_{\lambda}$ and $\tilde{\psi} \in E_{\lambda}$ be arbitrarily chosen, respectively. Define

$$\hat{z} = (\hat{z}_1, \dots, \hat{z}_d)^T = (\ln \hat{\psi}_1, \dots, \ln \hat{\psi}_d)^T \quad \text{and} \quad \tilde{z} = (\tilde{z}_1, \dots, \tilde{z}_d)^T = (\ln \tilde{\psi}_1, \dots, \ln \tilde{\psi}_d)^T.$$

Set $z := r\hat{z} + (1-r)\tilde{z}$, $\lambda := r\hat{\lambda} + (1-r)\tilde{\lambda}$, and $\psi(x) := (e^{z_1(x)}, \dots, e^{z_d(x)})^T$. One can easily verify that $\psi \in E_{\lambda}$. Hence, it can be used as a test function in (A.3) so that

$$k(\lambda, e) \ge \min_{1 \le i \le d} \inf_{x \in \mathbb{T}^N} \frac{e^{i(L-H)\psi}(x)}{\psi_i(x)}.$$

Direct computations show that for each $i \in \{1, \dots, d\}$,

$$\frac{L^i \psi_i}{\psi_i} = \frac{1}{\psi_i} [-\text{tr}(A^i D^2 \psi_i) + q^i \cdot \nabla \psi_i] = -\text{tr}(A^i D^2 z_i) - \nabla z_i A^i \nabla z_i + q^i \cdot \nabla z_i,$$

$$\frac{e^i B \psi}{\psi_i} = \frac{1}{\psi_i} \sum_{j=1}^{d} b_{ij} \psi_j = \sum_{j=1}^{d} b_{ij} e^{r(\hat{z}_j - \hat{z}_i) + (1-r)(\tilde{z}_j - \tilde{z}_i)}.$$

Moreover, due to $0 < r < 1$ and by (1.5), one has

$$\begin{aligned} \nabla z_i A^i \nabla z_i &= r \nabla \hat{z}_i A^i \nabla \hat{z}_i - r(1-r)(\nabla \hat{z}_i - \nabla \tilde{z}_i)A^i(\nabla \tilde{z}_i - \nabla \hat{z}_i) + (1-r)\nabla \tilde{z}_i A^i \nabla \tilde{z}_i \\ &\le r \nabla \hat{z}_i A^i \nabla \hat{z}_i + (1-r)\nabla \tilde{z}_i A^i \nabla \tilde{z}_i, \end{aligned} \quad (A.4)$$

and also

$$e^{r(\hat{z}_j - \hat{z}_i) + (1-r)(\tilde{z}_j - \tilde{z}_i)} \le r e^{\hat{z}_j - \hat{z}_i} + (1-r)e^{\tilde{z}_j - \tilde{z}_i}$$

since the function $x \mapsto e^x$ is convex. As a consequence, we obtain

$$\begin{aligned} \frac{e^i(L-H)\psi}{\psi_i} &\ge r \left( -\text{tr}(A^i D^2 \hat{z}_i) - \nabla \hat{z}_i A^i \nabla \hat{z}_i + q^i \cdot \nabla \hat{z}_i - \sum_{j=1}^d h_{ij} e^{\hat{z}_j - \hat{z}_i} \right) \\ &\quad + (1-r) \left( -\text{tr}(A^i D^2 \tilde{z}_i) - \nabla \tilde{z}_i A^i \nabla \tilde{z}_i + q^i \cdot \nabla \tilde{z}_i - \sum_{j=1}^d h_{ij} e^{\tilde{z}_j - \tilde{z}_i} \right) \\ &\ge r \left( \frac{-\text{tr}(A^i D^2 \hat{\psi}_i) + q^i \cdot \nabla \hat{\psi}_i}{\hat{\psi}_i} - \sum_{j=1}^d h_{ij} \frac{\hat{\psi}_j}{\hat{\psi}_i} \right) \\ &\quad + (1-r) \left( \frac{-\text{tr}(A^i D^2 \tilde{\psi}_i) + q^i \cdot \nabla \tilde{\psi}_i}{\tilde{\psi}_i} - \sum_{j=1}^d h_{ij} \frac{\tilde{\psi}_j}{\tilde{\psi}_i} \right) \\ &= r \frac{e^i(L-H)\hat{\psi}}{\hat{\psi}_i} + (1-r) \frac{e^i(L-H)\tilde{\psi}}{\tilde{\psi}_i}, \quad i=1, \dots, d. \end{aligned}$$

This implies that

$$k(\lambda, e) \ge \min_{1 \le i \le d} \inf_{x \in \mathbb{T}^N} \frac{e^i(L-H)\psi}{\psi_i} \ge r \min_{1 \le i \le d} \inf_{x \in \mathbb{T}^N} \frac{e^i(L-H)\hat{\psi}}{\hat{\psi}_i} + (1-r) \min_{1 \le i \le d} \inf_{x \in \mathbb{T}^N} \frac{e^i(L-H)\tilde{\psi}}{\tilde{\psi}_i}.$$
Eventually, since $\hat{\psi}$ and $\tilde{\psi}$ were arbitrarily chosen, (A.1) then follows.

Next, inspired by [29, Proposition 2.2], we prove the strict concavity of $\lambda \mapsto k(\lambda, e)$. Assume by contradiction that there exist two real numbers $\hat{\lambda} < \tilde{\lambda}$ such that (A.1) for $r = 1/2$ (without loss of generality) is an equality. Let $\phi_{\hat{\lambda}e}$ and $\phi_{\tilde{\lambda}e}$ be the unique principal eigenfunctions of (2.2) in the sense of (2.5) corresponding to $\lambda = \hat{\lambda}$ and $\tilde{\lambda}$, respectively. Set

$$ \hat{\psi} := \phi_{\hat{\lambda}e} e^{-\hat{\lambda}e \cdot x} \in E_{\hat{\lambda}}, \quad \tilde{\psi} := \phi_{\tilde{\lambda}e} e^{-\tilde{\lambda}e \cdot x} \in E_{\tilde{\lambda}}, \qquad (A.5) $$

and

$$ \phi = (\phi_1, \dots, \phi_d)^T = \left( \sqrt{(\phi_{\hat{\lambda}e})_1 (\phi_{\tilde{\lambda}e})_1}, \dots, \sqrt{(\phi_{\hat{\lambda}e})_d (\phi_{\tilde{\lambda}e})_d} \right)^T. \qquad (A.6) $$

Such a choice then yields

$$ \psi := \phi e^{-\frac{\hat{\lambda}+\tilde{\lambda}}{2}e \cdot x} \in E_{(\hat{\lambda}+\tilde{\lambda})/2}, $$

where $E_{(\hat{\lambda}+\tilde{\lambda})/2}$ is given by (A.2). Thus,

$$ k((\hat{\lambda} + \tilde{\lambda})/2, e) \geq \min_{1 \leq i \leq d} \inf_{\mathbb{T}^N} \frac{e^i(L-H)\psi}{\psi_i} \geq \frac{1}{2}(k(\hat{\lambda}, e) + k(\tilde{\lambda}, e)). $$

Combining our contradiction hypothesis with formula (2.3), one gets

$$ \min_{1 \le i \le d} \inf_{\mathbb{T}^N} \frac{e^i (L_{(\hat{\lambda}+\tilde{\lambda})/2} - H) \phi_i}{\phi_i} = \frac{1}{2} \min_{1 \le i \le d} \inf_{\mathbb{T}^N} \frac{e^i (L_{\hat{\lambda}e} - H) \phi_{\hat{\lambda}e}}{(\phi_{\hat{\lambda}e})_i} + \frac{1}{2} \min_{1 \le i \le d} \inf_{\mathbb{T}^N} \frac{e^i (L_{\tilde{\lambda}e} - H) \phi_{\tilde{\lambda}e}}{(\phi_{\tilde{\lambda}e})_i}. $$

This means that the function $\phi$, defined by (A.6), is a principal eigenfunction corresponding to $k((\hat{\lambda} + \tilde{\lambda})/2, e)$. From our preceding computations, and in particular (A.4), we then have that $\nabla \psi_i / \hat{\psi}_i = \nabla \tilde{\psi}_i / \tilde{\psi}_i$ for all $i \in \{1, \dots, d\}$. By (A.5), we further obtain

$$ \frac{\nabla(\phi_{\hat{\lambda}e})_i}{(\phi_{\hat{\lambda}e})_i} - \frac{\nabla(\phi_{\tilde{\lambda}e})_i}{(\phi_{\tilde{\lambda}e})_i} = (\hat{\lambda} - \tilde{\lambda})e, \quad i = 1, \dots, d. $$

By integrating these equations over $\mathbb{T}^N$, one can deduce that $\hat{\lambda} = \tilde{\lambda}$. This leads to a contradiction, thereby completing the proof. $\square$

**Proof of Proposition 2.2 (v).** The proof is inspired by [6, Proposition 4.4] (see also [46, Theorem 2.12]). Firstly, by Proposition 2.1, there always exists a generalized principal eigenfunction $\varphi \in C^2(\mathbb{R}^N)^d$ associated with $\lambda_1$, that is,

$$ (L - H)\varphi = \lambda_1 \varphi, \quad \varphi \gg 0 \text{ in } \mathbb{R}^N. \qquad (A.7) $$

Here $\varphi = (\varphi_1, \dots, \varphi_d)^T$ is not necessarily bounded, nor uniformly positive (that is, it is possible that $\inf_{\mathbb{R}^N} \min_{1 \le i \le d} \varphi_i = 0$). Moreover, if $\phi_{\lambda e} \in C^2(\mathbb{T}^N)^d$ is a principal eigenfunction of (2.2), then the function $\mathbb{R}^N \ni x \mapsto e^{\lambda e \cdot x} \phi_{\lambda e}$ satisfies (A.7). By (2.1), one has $\lambda_1 \ge k(\lambda, e)$ for all $\lambda \in \mathbb{R}$ and $e \in \mathbb{S}^{N-1}$ and thus $\lambda_1 \ge \sup_{\lambda e \in \mathbb{R}^N} k(\lambda, e)$. Note also that for any given $e \in \mathbb{S}^{N-1}$, $\lambda \mapsto k(\lambda, e)$ is strictly concave. Therefore, to conclude it suffices to show that there exists $\alpha \in \mathbb{R}^N$ such that $\lambda_1 = k(\alpha)$.

Let $e^1 = (1, 0, \dots, 0)^T \in \mathbb{R}^d$ and define a function $\psi = (\psi_1, \dots, \psi_d)^T(x)$ for $x \in \mathbb{R}^N$ as

$$ \psi_i(x) = \frac{\varphi_i(x + e^1)}{\varphi_i(x)}, \quad i = 1, \dots, d. \qquad (A.8) $$

Applying the Harnack inequality from Theorem B.1 to (A.7), we deduce that for any fixed $R > 0$, there exists a constant $C_R > 0$ such that for each $x \in \mathbb{R}^N$,

$$ \forall y \in B_{R+1}(x), \quad \max_{1 \le i \le d} \varphi_i(y) \le C_R \min_{1 \le i \le d} \varphi_i(x). $$
We then infer from the periodicity of equations that the function $\psi$ defined by (A.8) is globally bounded. Set

$$M_1 := \sup_{x \in \mathbb{R}^N} \max_{1 \le i \le d} \psi_i(x). \qquad (A.9)$$

Then there exist $\{x_n\}_{n \in \mathbb{N}} \subset \mathbb{R}^N$ and $i_0 \in \{1, \dots, d\}$ such that $\psi_{i_0}(x_n) \to M_1$ as $n \to \infty$. For each $n \in \mathbb{N}$, let $\bar{x}_n \in [0, 1]^N$ be such that $x_n - \bar{x}_n \in \mathbb{Z}^N$. Up to extracting a subsequence, $\bar{x}_n \to \bar{x}_\infty \in [0, 1]^N$ as $n \to \infty$. Now, we define two sequences $\{\tilde{\varphi}^n\}_{n \in \mathbb{N}}$ and $\{\hat{\varphi}^n\}_{n \in \mathbb{N}}$ of vector-valued functions as follows:

$$\begin{aligned}
\tilde{\varphi}_i^n(x) &:= \frac{\varphi_i(x + x_n)}{\varphi_{i_0}(x_n)} \quad \text{and} \\
\hat{\varphi}_i^n(x) &:= \frac{\varphi_i(x + e^1 + x_n) - M_1\varphi_i(x + x_n)}{\varphi_{i_0}(x_n)}, \quad i = 1, \dots, d.
\end{aligned} \qquad (A.10)$$

For each $n \in \mathbb{N}$, $\tilde{\varphi}^n$ satisfies (A.7) with a spatial shift of $\bar{x}_n$. Applying the Harnack inequality again, we deduce that $\tilde{\varphi}^n(x)$ is locally bounded in $\mathbb{R}^N$. Indeed, due to $\tilde{\varphi}_{i_0}^n(0) = 1$, for any given $R > 0$, there exists a constant $C_R > 0$ such that

$$\sup_{x \in B_R} \max_{1 \le i \le d} \tilde{\varphi}_i^n(x) \le C_R \inf_{x \in B_R} \min_{1 \le i \le d} \tilde{\varphi}_i^n(x) \le C_R, \quad n \in \mathbb{N}.$$

From Schauder interior estimates and Sobolev's injections, there exists a subsequence of $\{\tilde{\varphi}^n\}$ that converges in $C^{2+\alpha'}(\overline{B_{R/2}})^d$ for all $\alpha' \in [0, \alpha)$. Since the convergence holds for any $R > 0$, a diagonal argument allows us to extract a particular subsequence (not relabeled) such that $\tilde{\varphi}^n$ converges in $C_{\text{loc}}^{2+\alpha'}(\mathbb{R}^N)^d$ to a function $\tilde{\varphi}^\infty$ satisfying

$$\begin{cases} L[x + \bar{x}_{\infty}, \tilde{\varphi}^{\infty}] - H(x + \bar{x}_{\infty})\tilde{\varphi}^{\infty} = \lambda_1 \tilde{\varphi}^{\infty}, \\ \tilde{\varphi}^{\infty} \ge \mathbf{0} \text{ in } \mathbb{R}^N, \quad \tilde{\varphi}_{i_0}^{\infty}(0) = 1. \end{cases}$$

It then follows from the strong maximum principle, together with the fact that $H$ is cooperative and fully coupled, that $\tilde{\varphi}^\infty(x) \gg \mathbf{0}$ for all $x \in \mathbb{R}^N$.

On the other hand, due to the periodicity of equations, the function $\hat{\varphi}^n$ defined in (A.10) also satisfies (A.7) with a spatial shift of $\bar{x}_n$. Moreover, it follows from (A.8) and (A.9) that $\hat{\varphi}^n(x) \le \mathbf{0}$ for all $x \in \mathbb{R}^N$ and $\hat{\varphi}_{i_0}^n(0) = \psi_{i_0}(x_n) - M_1$. From standard elliptic estimates, we have that $\hat{\varphi}^n(x) \to \hat{\varphi}^\infty(x)$ as $n \to \infty$, after passing to a subsequence, locally uniformly for $x \in \mathbb{R}^N$, where $\hat{\varphi}^\infty$ satisfies

$$\begin{cases} L[x + \bar{x}_{\infty}, \hat{\varphi}^{\infty}] - H(x + \bar{x}_{\infty})\hat{\varphi}^{\infty} = \lambda_1 \hat{\varphi}^{\infty}, \\ \hat{\varphi}^{\infty} \le \mathbf{0} \text{ in } \mathbb{R}^N, \quad \hat{\varphi}_{i_0}^{\infty}(0) = 0. \end{cases}$$

By the strong maximum principle, we deduce that $\hat{\varphi}^\infty = \mathbf{0}$. Recalling the definition of $\hat{\varphi}^n$ in (A.10), this implies in particular that

$$\forall x \in \mathbb{R}^N, \quad \tilde{\varphi}^\infty(x + e^1) = M_1 \tilde{\varphi}^\infty(x) \gg \mathbf{0}.$$

As a consequence, if we set $\alpha_1 := \ln M_1$, then $\phi_{\alpha_1}(x) := \tilde{\varphi}^\infty(x)e^{\alpha_1 x_1}$ is 1-periodic in $x_1$. Since $\tilde{\varphi}^\infty$ is a generalized principal eigenfunction of the operator $L[\cdot + \bar{x}_\infty] - H(\cdot + \bar{x}_\infty)$ in $\mathbb{R}^N$ and $\tilde{\varphi}_{i_0}^\infty(0) = 1$, we then obtain

$$\begin{cases} e^{\alpha_1 x_1} L[x + \bar{x}_\infty, \phi_{\alpha_1} e^{-\alpha_1 x_1}] = H(x + \bar{x}_\infty)\phi_{\alpha_1} + \lambda_1 \phi_{\alpha_1}, \\ \phi_{\alpha_1} \gg \mathbf{0} \text{ in } \mathbb{R}^N, \quad (\phi_{\alpha_1})_{i_0}(0) = 1. \end{cases}$$

Now let $e^2 = (0, 1, 0, \dots, 0)^T \in \mathbb{R}^d$ and modify the definition (A.8) of $\psi$ as follows:

$$\psi_i(x) = \frac{(\phi_{\alpha_1})_i(x + e^2)}{(\phi_{\alpha_1})_i(x)}, \quad i = 1, \dots, d,$$
which is also globally bounded by the Harnack inequalities once again. Set

$$M_2 := \sup_{x \in \mathbb{R}^N} \max_{1 \le i \le d} \psi_i(x).$$

Repeating the above proof, we obtain a function $\tilde{\phi}_{\alpha_1}^{\infty}$ such that the following function

$$\phi_{\alpha_1 e^{1+\alpha_2 e^2}(x)} := \tilde{\phi}_{\alpha_1}^{\infty}(x) e^{\alpha_1 x_1 + \alpha_2 x_2} \quad \text{with } \alpha_2 := \ln M_2$$

is $\mathbb{Z}^2$-periodic in $(x_1, x_2)$. Going on the above construction until the N-th time, one can find $\alpha_i = \ln M_i$ for all $i \in \{3, \dots, N\}$ and then get a function $\varphi$ verifying

$$\left\{ \begin{array}{l} L[x + \bar{z}_{\infty}, \varphi] - H(x + \bar{z}_{\infty})\varphi = \lambda_1 \varphi \quad \text{for } (x, \bar{z}_{\infty}) \in \mathbb{R}^N \times \mathbb{R}^N, \\ \varphi \gg \mathbf{0}, \quad \varphi_{i_0}(0) = 1 \quad \text{and} \quad \varphi(x)e^{\alpha \cdot x} \text{ is } \mathbb{Z}^N\text{-periodic in } x, \end{array} \right.$$

where one has set $\alpha := (\alpha_1, \dots, \alpha_N)^T \in \mathbb{R}^N$. Define now $\phi_\alpha(x) := \varphi(x)e^{\alpha \cdot x}$ and recall that

$$e^{\alpha \cdot x} L^i(\phi_\alpha e^{-\alpha \cdot x}) = L_\alpha^i \phi_\alpha \quad \text{for } (\alpha, x) \in \mathbb{R}^N \times \mathbb{R}^N,$$

where $L_\alpha^i$ denotes the operator given by (1.8) with $\lambda e = \alpha$. Set $L_\alpha := \text{diag}(L_\alpha^1, \dots, L_\alpha^d)$. Then,

$$\left\{ \begin{array}{l} L_\alpha[x + \bar{z}_\infty, \phi_\alpha] - H(x + \bar{z}_\infty)\phi_\alpha = \lambda_1 \phi_\alpha \quad \text{for } (x, \bar{z}_\infty) \in \mathbb{R}^N \times \mathbb{R}^N, \\ \forall (x, k) \in \mathbb{R}^N \times \mathbb{Z}^N, \quad \phi_\alpha(x + k) = \phi_\alpha(x), \quad \phi_\alpha(x) \gg \mathbf{0}. \end{array} \right.$$

Since the periodic principal eigenvalue of (2.2) is invariant under a translation of the coefficients, it follows from the uniqueness of the principal eigenpairs that $C\phi_\alpha$ for some constant $C > 0$ is the principal eigenfunction of (2.2), and $\lambda_1 = k(\alpha)$. This concludes the proof. $\square$

## B Harnack-type estimates for linear cooperative and fully coupled systems

For the sake of convenience, we present the Harnack inequalities for linear cooperative and fully coupled systems in the context of the present work. We refer to [12, Theorem 1.2], [1, Theorem 2.2] and [11, Corollary 8.1] for results on elliptic systems.

**Theorem B.1 (Elliptic case).** Let $L = \text{diag}(L^1, \dots, L^d)$ denote a diagonal matrix of the linear elliptic operators given by (1.4). Let $H \in C^\alpha(\mathbb{T}^N, \mathcal{M}_d(\mathbb{R}))$ be cooperative and fully coupled in the sense of Definition 1.1. Consider the system of linear elliptic equations

$$(L - H)u = \mathbf{0} \quad \text{in } \mathbb{R}^N. \qquad (B.1)$$

Denote by $\bar{h}_0$ the smallest positive entry of $\bar{H} = (\bar{h}_{ij})_{1 \le i,j \le d}$ where $\bar{h}_{ij} := \max_{x \in \mathbb{T}^N} h_{ij}(x)$.

There exists a constant $C_R > 0$ determined only by $N$, $d$, the coefficients of $L$ and $\bar{h}_0$ such that if $u = (u_1, \dots, u_d)^T \in C(\bar{B}_{R+1})^d$ is a nonnegative solution to (B.1), then

$$\max_{1 \le i \le d} \sup_{x \in B_R} u_i(x) \le C_R \min_{1 \le i \le d} \inf_{x \in B_R} u_i(x).$$

In view of parabolic systems, we refer to [25, Theorem 3.9], which assumes pointwise irreducibility that is slightly stronger than the fully coupled property. Girardin and Mazari have recently revisited and refined such estimates in the space-time periodic setting, see [28, Proposition 2.4] for details.
**Theorem B.2 (Parabolic case).** Let $H(t,x) = (h_{ij}(t,x))_{1\le i,j\le d}$ be a cooperative and fully coupled matrix field according to Definition 1.1. Assume that each entry $h_{ij}: \mathbb{R} \times \mathbb{R}^N \to \mathbb{R}$ is bounded and continuous. Consider the system of linear parabolic equations

$$ \partial_t u_i - L^i u_i = \sum_{j=1}^{d} h_{ij}(t,x) u_j, \quad (t,x) \in \mathbb{R} \times \mathbb{R}^N, \quad i = 1, \dots, d, \qquad (B.2) $$

where $L^i$ is defined by (1.4). Denote by $\bar{h}_0$ be the smallest positive entry of $\bar{H} = (\bar{h}_{ij})_{1\le i,j\le d}$ where $\bar{h}_{ij} := \max_{(t,x)\in\mathbb{R}\times\mathbb{R}^N} h_{ij}(t,x)$.

Let $\theta > 1$ be given. There exists a constant $\hat{\kappa}_\theta > 0$ determined only by $N$, $d$, the coefficients of $L^i$ and $\bar{h}_0$ such that if $u = (u_1, \dots, u_d)^T \in C([-2\theta, 2\theta] \times [-\theta, \theta]^N)^d$ is a nonnegative solution to (B.2), then

$$ \min_{1 \le i \le d} \min_{(t,x) \in [\frac{3\theta}{2}, 2\theta] \times [-\theta, \theta]^N} u_i(t,x) \ge \hat{\kappa}_\theta \max_{1 \le i \le d} \max_{(t,x) \in [-\theta, 0] \times [-\theta, \theta]^N} u_i(t,x). $$

Note that in [28, Proposition 2.4], only the diagonal entries of the zero-order term matrix are not required to be periodic. But going their proof, the Harnack inequality from Theorem B.2 still holds as long as $H$ is globally bounded.

## References

[1] A. Arapostathis, M. K. Ghosh, and Marcus S. I., *Harnack's inequality for cooperative weakly coupled elliptic systems*, Comm. Partial Differential Equations **24** (1999), 1555-1571.

[2] M. Bages, P. Martinez, and J.-M. Roquejoffre, *How travelling waves attract the solutions of KPP-type equations*, Trans. Amer. Math. Soc. **364** (2012), 5415-5468.

[3] G. Barles, L. C. Evans, and P. E. Souganidis, *Wavefront propagation for reaction-diffusion systems of PDE*, Duke Math. J. **61** (1990), 835-858.

[4] H. Berestycki and F. Hamel, *Front propagation in periodic excitable media*, Comm. Pure Appl. Math. **55** (2002), 949-1032.

[5] H. Berestycki and F. Hamel, *Generalized travelling waves for reaction-diffusion equations, In: Perspectives in Nonlinear Partial Differential Equations. In honor of H. Brezis*, Contemp. Math., vol. 446, AMS, Providence, RI, 2007.

[6] H. Berestycki, F. Hamel, and G. Nadin, *Asymptotic spreading in heterogeneous diffusive excitable media*, J. Funct. Anal. **255** (2008), 2146-2189.

[7] H. Berestycki, F. Hamel, and L. Roques, *Analysis of the periodically fragmented environment model: I-species persistence*, J. Math. Biol. **51** (2005), 75-113.

[8] H. Berestycki, F. Hamel, and L. Roques, *Analysis of the periodically fragmented environment model: II-biological invasions and pulsating travelling fronts*, J. Math. Pures Appl. **84** (2005), 1101-1146.

[9] H. Berestycki, F. Hamel, and L. Rossi, *Liouville-type results for semilinear elliptic equations in unbounded domains*, Ann. Mat. Pura Appl. **186** (2007), 469-507.

[10] I. Birindelli, E. Mitidieri, and G. Sweers, *Existence of the principal eigenvalue for cooperative elliptic systems in a general domain*, Differ. Uravn. **35** (1999), 326-334.

[11] J. Busca and B. Sirakov, *Harnack type estimates for nonlinear elliptic systems and applications*, Ann. Inst. H. Poincaré C Anal. Non Lin. **21** (2004), 543-590.

[12] Z.-Q. Chen and Z. Zhao, *Harnack principle for weakly coupled elliptic systems*, J. Differential Equations **139** (1997), 261-282.

[13] L. Deng and A. Ducrot, *Existence of multi-dimensional pulsating fronts for KPP equations: A new formulation approach*, Calc. Var. Partial Differential Equations **62** (2023), No. 134.
[14] L. Deng and A. Ducrot, *Pulsating waves in a multidimensional reaction-diffusion system of epidemic type*, J. Eur. Math. Soc. (JEMS), online first, DOI: 10.4171/JEMS/1511 (2024).

[15] L. Deng and Z.-C. Wang, *Propagation phenomena for a criss-cross infection model with non-diffusive susceptible population in periodic media*, Discrete Contin. Dyn. Syst. Ser. B **26** (2021), 4789-4814.

[16] W. Ding and T. Giletti, *Admissible speeds in spatially periodic bistable reaction-diffusion equations*, Adv. Math. **389** (2021), No.107889.

[17] W. Ding, F. Hamel, and X.-Q. Zhao, *Bistable pulsating fronts for reaction-diffusion equations in a periodic habitat*, Indiana Univ. Math. J. **66** (2017), 1189-1265.

[18] L.-J. Du and X. Bao, *On the existence and monotonicity of pulsating traveling waves for reaction-diffusion-advection systems in high dimensional and periodic media*, Nonlinear Anal. Real World Appl. **64** (2022), No. 103452.

[19] L.-J. Du, W.-T. Li, and W. Shen, *Propagation phenomena for time-space periodic monotone semiflows and applications to cooperative systems in multi-dimensional media*, J. Funct. Anal. **282** (2022), No. 109415.

[20] A. Ducrot, *A multi-dimensional bistable nonlinear diffusion equation in a periodic medium*, Math. Ann. **366** (2016), 783-818.

[21] A. Ducrot, T. Giletti, and H. Matano, *Existence and convergence to a propagating terrace in one-dimensional reaction-diffusion equations*, Trans. Amer. Math. Soc. **366** (2014), 5541-5566.

[22] J. Fang, X. Yu, and X.-Q. Zhao, *Traveling waves and spreading speeds for time-space periodic monotone systems*, J. Funct. Anal. **272** (2017), 4222-4262.

[23] J. Fang and X.-Q. Zhao, *Bistable traveling waves for monotone semiflows with applications*, J. Eur. Math. Soc. (JEMS) **17** (2015), 2243-2288.

[24] R. A. Fisher, *The wave of advance of advantageous genes*, Ann. Eugenics **7** (1937), 335-369.

[25] J. Földes and P. Poláčik, *On cooperative parabolic systems: Harnack inequalities and asymptotic symmetry*, Discrete Contin. Dyn. Syst. **25** (2009), 133-157.

[26] T. Giletti and L. Rossi, *Pulsating solutions for multidimensional bistable and multistable equations*, Math. Ann. **378** (2019), 1555-1611.

[27] L. Girardin, *Persistence, extinction, and spreading properties of noncooperative Fisher-KPP systems in space-time periodic media*, SIAM J. Math. Anal. **57** (2025), 233-261.

[28] L. Girardin and I. Mazari, *Generalized principal eigenvalues of space-time periodic, weakly coupled, cooperative, parabolic systems*, To appear in *Mémoires de la Société Mathématique de France*.

[29] Q. Griette and H. Matano, *Propagation dynamics of solutions to spatially periodic reaction-diffusion systems with hybrid nonlinearity*, preprint, arXiv:2108.10862, 2021.

[30] Q. Griette and H. Matano, *Front propagation in hybrid reaction-diffusion epidemic models with spatial heterogeneity. Part I: Spreading speed and asymptotic behavior*, Asymptot. Anal. **144** (2025), 1291-1326.

[31] F. Hamel, *Qualitative properties of monostable pulsating fronts: exponential decay and monotonicity*, J. Math. Pures Appl. **89** (2008), 355-399.

[32] F. Hamel, F. Lutscher, and M. Zhang, *Propagation phenomena in periodic patchy landscapes with interface conditions*, J. Dynam. Differential Equations **36** (2024), 435-486.

[33] F. Hamel and L. Roques, *Uniqueness and stability properties of monostable pulsating fronts*, J. Eur. Math. Soc. (JEMS) **13** (2011), 345-390.

[34] S. Heinze, *Wave solutions to reaction-diffusion systems in perforated domains*, Z. Anal. Anwendungen **20** (2001), 661-676.

[35] M. Huang, S.-L. Wu, and X.-Q. Zhao, *Propagation dynamics for time-periodic and partially degenerate reaction-diffusion systems*, SIAM J. Math. Anal. **54** (2022), 1860-1897.
[36] W. Hudson and B. Zinner, *Existence of travelling waves for reaction-diffusion equations of Fisher type in periodic media. In: J. Henderson (ed.) Boundary Value Problems for Functional-Differential Equations*, World Scientific, 1995.

[37] T. Kato, *Perturbation Theory for Linear Operators*, 2nd ed., Springer, 1984.

[38] A. N. Kolmogorov, I. G. Petrovskii, and N. S. Piskunov, *Étude de l'équation de la diffusion avec croissance de la quantité de matière et son application à un problème biologique*, Bull. Univ. d'État à Moscou, Sér. Intern. A 1 (1937), 1–26.

[39] O. A. Ladyženskaja, V. A. Solonnikov, and N. N. Ural'ceva, *Linear and Quasi-linear Equations of Parabolic Type*, vol. 23, AMS, 1988.

[40] K.-Y. Lam and Y. Lou, *Asymptotic behavior of the principal eigenvalue for cooperative elliptic systems and applications*, J. Dynm. Differential Equations 28 (2016), 29–48.

[41] B. Li, H. F. Weinberger, and M. A. Lewis, *Spreading speeds as slowest wave speeds for cooperative systems*, Math. Biosci. 196 (2005), 82–98.

[42] X. Liang and X. Zhao, *Asymptotic speeds of spread and traveling waves for monotone semiflows with applications*, Comm. Pure Appl. Math. 60 (2007), 1–40.

[43] X. Liang and X.-Q. Zhao, *Spreading speeds and traveling waves for abstract monostable evolution systems*, J. Funct. Anal. 259 (2010), 857–903.

[44] H. Matano, *Traveling waves in spatially random media*, RIMS Kokyuroku 1337 (2003), 1–9.

[45] J.D. Murray, *Mathematical Biology*, I: An Introduction, 3rd ed., Springer, Berlin, 2002.

[46] G. Nadin, *The principal eigenvalue of a space-time periodic parabolic operator*, Ann. Mat. Pura Appl. 188 (2009), 269–295.

[47] G. Nadin, *Traveling fronts in space-time periodic media*, J. Math. Pures Appl. 9 (2009), 232–262.

[48] G. Nadin, *Existence and uniqueness of the solution of a space-time periodic reaction-diffusion equation*, J. Differential Equations 249 (2010), 1288–1304.

[49] T. Ogiwara and H. Matano, *Monotonicity and convergence results in order-preserving systems in the presence of symmetry*, Discrete Contin. Dyn. Syst. 5 (1998), 1–34.

[50] T. Ogiwara and H. Matano, *Stability analysis in order-preserving systems in the presence of symmetry*, Proc. Roy. Soc. Edinburgh Sect. A 129 (1999), 395–438.

[51] M. H. Protter and H. F. Weinberger, *Maximum Principles in Differential Equations*, Springer Verlag, New York, 1984.

[52] W. Schlag, *Schauder and $L^p$ estimates for parabolic systems via Campanato spaces*, Comm. Partial Differential Equations 21 (1996), 1141–1175.

[53] W. Shen, *Traveling waves in diffusive random media*, J. Dynam. Differential Equation, 16 (2004), 1011–1060.

[54] N. Shigesada, K. Kawasaki, and E. Teramoto, *Traveling periodic waves in heterogeneous environments*, Theor. Popul. Biol. 30 (1986), 143–160.

[55] H. L. Smith, *Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems*, vol. 41, AMS, Providence, RI, 1995.

[56] G. Sweers, *Strong positivity in $C(\bar{\Omega})$ for elliptic systems*, Math. Z. 209 (1992), 251–271.

[57] A.I. Volpert, V.A. Volpert, and V.A. Volpert, *Traveling Wave Solutions of Parabolic Systems*, Translations of Mathematical Monographs, vol. 140, AMS, Providence, RI, 1994.

[58] V.A. Volpert, *Elliptic Partial Differential Equations*, Vol. 2, Reaction-Diffusion Equations, Birkhäuser, Basel, 2014.

[59] X.-S. Wang and X.-Q. Zhao, *Pulsating waves of a partially degenerate reaction-diffusion system in a periodic habitat*, J. Differential Equations 259 (2015), 7238–7259.
[60] H. F. Weinberger, *Long-time behavior of a class of biological models*, SIAM J. Math. Anal. **13** (1982), 353–396.

[61] H. F. Weinberger, *On spreading speeds and traveling waves for growth and migration models in a periodic habitat*, J. Math. Biol. **45** (2002), 511–548.

[62] J. Xin, *Front propagation in heterogeneous media*, SIAM Rev. **42** (2000), 161–230.

[63] X. Yu and X.-Q. Zhao, *Propagation phenomena for a reaction-advection-diffusion competition model in a periodic habitat*, J. Dyn. Differential Equations **29** (2017), 41–66.

[64] X.-Q. Zhao, *Dynamical Systems in Population Biology*, 2nd ed., Springer, Cham, 2017.