Acoustic Wind in a Hyperbolic Predator—Prey System

Abstract

We address a hyperbolic model for prey-sensitive predators interacting with purely diffusive prey. We adopt the Cattaneo formulation for describing the predators’ transport. Given the hyperbolicity, the long-lived short-wave patterns occur for sufficiently weak prey diffusivities. The main result is that the non-linear interplay of the short waves generically excites the slowly growing amplitude modulation for wide ranges of the model parameters. We have observed such a feature in the numerical experiments and support our conclusions with a short-wave asymptotic solution in the limit of vanishing prey diffusivity. Our reasoning relies on the so-called homogenized system that governs slow evolutions of the amplitudes of the short-wave parcels. It includes a term (called wind) which is absent in the original model and only comes from averaging over the short waves. It is the wind that (unlike any of the other terms!) is capable of exciting the instability and pumping the growth of solutions. There is quite a definite relationship between the predators’ transport coefficients to be held for getting rid of the wind. Interestingly, this relationship had been introduced in prior studies of small-scale mosaics in the spatial distributions of some real-life populations.

1. Introduction

We will be considering two species, one of which (say, predator) is sensitive to the concentration gradient of the other one (say, prey). Prey spread in a purely diffusive way. For the predators, we employ the Cattaneo model [1]. The resulting system of equations belongs to the family of hyperbolic models [2]. The last reference alludes to a review that delivers much helpful information on the subject.

The present article considers the slow evolution of short-wave packages for the outlined PDE system. In many respects, it continues article [3] that regards the effect of small-scale fluctuations of the environment on the stability of the homogeneous equilibria for one parabolic model of the Patlak–Keller–Segel family. Currently, we address a similar problem for a hyperbolic model. There are at least two reasons for doing so. The first one is the completeness of the study. Indeed, the hyperbolic and parabolic models, such as the Cattaneo and the Patlak–Keller–Segel formulations, arise from the same kinetic model under different hypotheses. Both models are natural but very different. For instance, the parabolic model must perceive a short wavelength external signal constantly to produce a short wave pattern in response [3]. In great contrast, the short waves can propagate due to an instantaneous perturbation in the hyperbolic case. The second reason regards the so-called aggregation that review [2] distinguishes as an important application of the hyperbolic models. Since the aggregation means the spontaneous forming of large-scale patterns from small-scale ones, there should be an intrinsic mechanism of destabilization and destruction of the latter. At the same time, there should be some exclusions as real-life populations sometimes display the small-scale mosaic patterns of spatial distributions [4].

Of course, such topics as the relevance of the short-wavelength limit to a real life population or the wave solutions of the hyperbolic models (which have got much less attentions than their parabolic counterparts as yet) are worth continuing discussion. However, we skip it for the sake of brevity and allude to freely accessible papers [3,5] for these matters and for more references to the subject.

In contrast with the Patlak–Keller–Segel parabolic model, which we addressed in article [3], the short-wave limit for the Cattaneo model, due to its hyperbolicity, leads to the elementary wave equation, the solutions to which form a kind of short-wave ether. We apply the homogenization over this ether to derive equations for the slow modulation of the short-wave packages. To some extent, this procedure follows that proposed in article [6]. The main result in the present article is the pumping of the slow amplitudes by the non-linear interactions of the short waves. This feature manifests itself in the excitement of wind, by which we mean some terms that enter the flux equation due to the homogenization. It turns out that this wind is capable of destabilizing the equilibria of the averaged amplitude equation and supporting the growth of their finite perturbations.

The article is organized as follows. In Section 2, we formulate the governing equation. In Section 3, we present several motivating examples. In Section 4, we derive the system that governs the slow evolution of the amplitudes of the short-wave parcels by the homogenization over the short waves. In Section 5, we consider the simplest model that invokes identically constant transport coefficients of the predator and the Lotka–Volterra kinetics. In Section 6, we unfold the destabilizing effect of wind through linear stability analysis. In Section 7, we present numerical solutions of the homogenized system that support our outcomes regarding the crucial part of the wind in supporting the instability and perturbation growth. In Section 8, we list several conclusions from the research presented.

2. The Governing Equations and Scaling

As we have mentioned above, we invoke the Cattaneo model [1] for the predators’ transport and the common reaction–diffusion equation for the prey. The former involves the equations for the balance of the predators and the equation for their flux that, in contrast with the Patlak–Keller–Segel law, takes into account the inertia of their responding. The dimensionless form of the governing equations reads as

$$\begin{array}{ccc}\hfill p_t+q_x& =& F(p,s),\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill q_t+\nu q& =& \kappa \chi (p,s)p{s}_x-\mu (p,s)p_x,\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill s_t& =& G(p,s)+\delta s_{xx}.\hfill \end{array}$$

(3)

In article [5], there are more details on scaling the dimensional governing equations for bringing them into the above dimensionless form.

In system (1)–(3), the independent variables, x and t, stand for the spatial and temporal coordinates, correspondingly, $x\in \mathbb{R}$, $t>0$. The dependent variables are $p=p(x,t)$, $q=q(x,t)$ and $s=s(x,t)$. The first and the last ones play the parts of the densities of the species. The first two equations constitute the Cattaneo model for the prey-sensitive movement of the predators, where the remaining dependent variable, q, stands for predator’s flux. Coefficient $\kappa >0$ stands as the sensitivity factor. Moreover, Equation (2) involves term $\nu q$ with constant factor $\nu >0$, which resembles the resistance force in mechanics. We will be calling the corresponding coefficient, $\nu$, resistivity. The notation of $\delta$ stands for the prey diffusivity factor that is constant by assumption.

By default, the predator diffusivity, $\mu$, the sensitivity, $\chi$, and the reactions terms, F and G, are smooth functions in variables $p>0, s>0$; at that

$$\mu (p,s)>0\forall p>0,s>0,p\chi (p,s)\to 0,p\to +0,\forall s\ge 0.$$

(4)

3. Motivating Examples

Throughout this section, we will be discussing three numerical solutions to the Cauchy problem for Equations (1)–(3) with small diffusivity of prey and identically constant coefficients $\mu$ and $\chi$. The initial data and other settings are as follows

$${p|}_{t=0}=p_e\left(1+ϵ{\mathrm{e}}^{-x^2}+a{sin}^N\left(2\pi nx\right)\right),q{{|}_{t=0}=\sqrt{\mu }{p}_{e}a{sin}^N\left(2\pi nx\right),s|}_{t=0}=s_e,$$

(5)

$$F=p(\gamma s-\beta ),G=s(\alpha -p-s),s_e=\beta /\gamma ,p_e=\alpha -s_e,\beta /\left(\gamma \alpha \right)=0.4,\alpha =5/4,\gamma =2/3,$$

(6)

$$\mu \equiv const=4,\nu =1,\delta =0.01,\chi \equiv const=\mu \nu /\left(p_e{s}_e\right),\kappa =1.$$

(7)

Quantities $s_e$ and $p_e$ are the equilibrium species densities in the sense that

$$F(p_e,s_e)=0,G(p_e,s_e)=0,$$

i.e., the system (1)–(3) allows the homogeneous equilibrium in which $p\equiv p_e$, $s\equiv s_e$ and $q\equiv 0$. The value assigned to the parameter $\chi$ is a threshold for the linear instability of such an equilibrium in the case of $\delta =0$ (see Section 6).

Figure 1 illustrates the numerical solutions, which cover the spatial range of $(-35,35)$. For more details of the implementation, see Section 7. Each row of snapshots shows the instantaneous graphs of the deviation of the species densities from their equilibrium values for one concrete set of initial data (5). The top row shows those patterns that arise from a short-wave modulation of the homogeneous equilibrium. The middle row demonstrates an interplay of the same modulation with a small smooth localized perturbation. The bottom row corresponds to the same perturbation but with no underlying short wave.

Table 1. This table states the correspondence between patterns shown in Figure 1 and concrete versions of the initial data (5).

Row	Settings of the Initial Data
Top	$N=5,n=3,ϵ=0,a=1/2$
Middle	$N=5,n=3,ϵ=0.1,a=1/2$
Bottom	$N=5,n=3,ϵ=0.1,a=0$

states the correspondence between the rows and the versions of initial data (5).

Overall, the top snapshots row displays a short-wave profile of the predator density amid almost homogeneous prey distribution. The former resembles a linear travelling wave subject to slow amplitude damping. It is not entirely unexpected since the short-wave limit of the governing equations leads to the linear acoustic system with the sound speed equal to $\sqrt{\mu }$, namely,

$$p_t+q_x=0,q_t+\mu p_x=0.$$

The amplitude damping is due to the resistance term, $\nu q$, which is not entirely negligible for the wavelengths adopted for computing.

The middle row shows that the small smooth localized initial perturbation persists and even gets amplified despite the existing resistance. Notice the asymmetry too. It cannot be explained using the above acoustic approximation only. These features sharply contrast with the bottom row, which displays a symmetrical picture of decaying perturbation. This symmetry is due to the mirror invariance of Equations (1)–(3) and initial data (5) (the latter takes place for $a=0$ only).

The listed observations suggest that short waves are capable of pumping the perturbations. In the rest of the article, we will be addressing this conjecture with the use of homogenization.

4. Homogenization and Wind

Let $\delta \to +0$ in Equation (3). We seek the short wavelength solutions that read as

$$p=p(x,t,\xi ,\tau ,\delta ),q=q(x,t,\xi ,\tau ,\delta ),s=s(x,t,\xi ,\tau ,\delta ),\mathrm{where}\tau =t/\delta ,\xi =x/\delta .$$

For calculating such a solution we employ the formal power expansions as follows

$$(p,q,s)(x,t,\xi ,\tau ,\delta )=\sum _{k=0}{(p,q,s)}_k(x,t,\xi ,\tau ){\delta }^k.$$

We assume that for every coefficient in these series there exists an averaged value defined as follows

$$\langle u\rangle (x,t)\stackrel{def}{=}\underset{L\to \infty }{lim}\frac{1}{4L^2}\underset{{[-L,L]}^2}{\int }u(x,t,\xi ,\tau )d\xi d\tau .$$

(8)

Hence, for every $k=0,1,\dots$ there is a decomposition

$${(p,q,s)}_k(x,t,\xi ,\tau )={(\overline{p},\overline{q},\overline{s})}_k(x,t)+{(\tilde{p},\tilde{q},\tilde{s})}_k(x,t,\xi ,\tau ),\mathrm{where}{(\overline{p},\overline{q},\overline{s})}_k={(〈p〉,〈q〉,〈s〉)}_k,$$

so that the tilde terms vanish on average. To simplify the notation, we will be omitting the zero-valued indices here and there below—that is, the following agreement is valid from here on:

$${\overline{p}}_0=\overline{p},{\tilde{p}}_0=\tilde{p},{\overline{q}}_0=\overline{q},{\tilde{q}}_0=\tilde{q},{\overline{s}}_0=\overline{s},{\tilde{s}}_0=\tilde{s}.$$

Finally, we assume that

$$\tilde{s}={\tilde{s}}_0=0.$$

We write the governing equations with the use of the slow variables, $x,t$, and the fast variables, $\xi ,\tau$. For this purpose, we substitute ${\partial }_t$ and ${\partial }_x$ with ${\partial }_t+{\delta }^{-1}{\partial }_{\tau }$ and ${\partial }_x+{\delta }^{-1}{\partial }_{\xi }$. Then replacing the dependent variables, $p,q,s$, with the above expansions gives a chain of equations for the approximations of successive orders. Specifically, collecting the terms of order ${\delta }^{-1}$ yields equations:

$${\tilde{p}}_{\tau }+{\tilde{q}}_{\xi }=0;{\tilde{q}}_{\tau }+\mu (\overline{p}+\tilde{p},\overline{s}){\tilde{p}}_{\xi }=0.$$

(9)

This system is quasi-linear and hyperbolic with respect to unknowns $\tilde{p},\tilde{q}$ since function $\mu$ is positive by assumption. The value of $\overline{s}$ stands as a constant parameter. Let

$${\partial }_{p}M(p,\overline{s})=\mu (p,\overline{s}).$$

Given this, the second equation in system (9) takes the form of a conservation law that reads as

$${\tilde{q}}_{\tau }+{\partial }_{\xi }M(\overline{p}+\tilde{p},\overline{s})=0,$$

Eliminating of variable q reduces system (9) to a single quasi-linear wave equation that reads as

$${\tilde{p}}_{\tau \tau }={\left(\mu (\overline{p}+\tilde{p},\overline{s}){\tilde{p}}_{\xi }\right)}_{\xi }.$$

(10)

The next step after resolving Equation (10) is recovering the other unknown, $\tilde{q}$, from system (9). Some compatibility conditions can arise at that. One of them is Equation (10) itself and it can seem that one more results from system (9) upon averaging. However, the latter holds identically due to the aforementioned conservatism.

Collecting the terms of order 1 yields equations:

$${\tilde{p}}_{1\tau }+{\tilde{q}}_{1\xi }=F(p_0,\overline{s})-p_{0t}-q_{0x},$$

(11)

$$\begin{array}{c}{\tilde{q}}_{1\tau }+{\left(\mu (p_0,\overline{s})p_1\right)}_{\xi }+s_1{\left(M_s(p_0,\overline{s})\right)}_{\xi }=\\ \kappa p_0\chi (p_0,\overline{s})({\overline{s}}_x+{\tilde{s}}_{1\xi })-\mu (p_0,\overline{s})p_{0x}-q_{0t}-\nu q_0,\end{array}$$

(12)

$${\tilde{s}}_{1\tau }-{\tilde{s}}_{1\xi \xi }=G(p_0,\overline{s})-{\overline{s}}_t.$$

(13)

Regarding Equation (12), it is worth noting that the collection of the terms of order 1 resulting from the term ${\delta }^{-1}\mu (p,s)p_{\xi }$ reads as $({\mu }_s(p_0,\overline{s})s_1+{\mu }_p(p_0,\overline{s})p_1){\tilde{p}}_{0\xi }+\mu (p_0,\overline{s}){\tilde{p}}_{1\xi }$ and this, by definition of function M, is equal to ${\left(\mu (p_0,\overline{s})p_1\right)}_{\xi }+s_1{\left(M_s(p_0,\overline{s})\right)}_{\xi }.$

Again, there are several compatibility conditions to hold for isolating unknown ${\tilde{q}}_1$ from Equations (11) and (12). These are

$$\langle F(p_0,\overline{s})-p_{0t}-q_{0x}\rangle =0,$$

(14)

$$\langle \kappa p_0\chi (p_0,\overline{s}){\overline{s}}_x+{\tilde{s}}_{1\xi }\left(\kappa p_0\chi (p_0,\overline{s})+M_s(p_0,\overline{s})\right)-\nu q_0-\mu (p_0,\overline{s})p_{0x}-q_{0t}\rangle =0,$$

(15)

$$\begin{array}{c}{\tilde{p}}_{1\tau \tau }-{\left((\mu (p_0,\overline{s}){({\tilde{p}}_1+{\overline{p}}_1)}_{\xi }+({\tilde{s}}_1+{\overline{s}}_1){\left(M_s(p_0,\overline{s})\right)}_{\xi }\right)}_{\xi }=\\ {({\mu }_0(p_0,\overline{s})p_{0x}+q_{0t}+\nu q_0-\kappa p_0\chi (p_0,\overline{s})({\overline{s}}_x+{\tilde{s}}_{1\xi }))}_{\xi }+{(F(p_0,\overline{s})-p_{0t}-q_{0x})}_{\tau }.\end{array}$$

(16)

Equation (16) is a linear wave equation in fast variables with respect to function ${\tilde{p}}_1$ provided that we have determined $\overline{p}$, $\tilde{p}$, ${\tilde{s}}_1$ and $\overline{s}$. More accurately, this is a two-parametric family of the linear wave equations where the values of ${\overline{p}}_1$ and ${\overline{s}}_1$ stand as parameters.

Assume there exists a two-parametric family of solutions $\tilde{p}=P(\overline{p},\overline{s})$, $\tilde{q}=Q(\overline{p},\overline{s})$ to system (9). Equation (13) entails the following compatibility condition

$${\overline{s}}_t=\langle G(\overline{p}+\tilde{p},\overline{s})\rangle .$$

(17)

This equation is nothing else than Equation (13) on average. Assume that it holds. Then Equation (13) defines a two-parametric family of operators

$$S_1=S_1(\overline{p},\overline{s}):\tilde{p}\mapsto {\tilde{s}}_1,S_1=\mathcal{L}(G(\overline{p}+·,\overline{s})-\langle G(\overline{p}+·,\overline{s})\rangle ),$$

where $\mathcal{L}$ is the inverse to the heat operator on a suitable space of functions vanishing on average.

Further, let us write conditions (14) and (15) as equations for $\overline{p}$ and $\overline{q}$:

$$\begin{array}{ccc}& {\overline{p}}_t+{\overline{q}}_x=\langle F(\overline{p}+\tilde{p},\overline{s})\rangle ,& \hfill \end{array}$$

(18)

$$\begin{array}{ccc}& {\overline{q}}_t+\nu \overline{q}+\overline{\mu }(\overline{p},\overline{s}){\overline{p}}_x=\kappa \overline{p}\overline{\chi }(\overline{p},\overline{s}){\overline{s}}_x+W(\overline{p},\overline{s}),& \hfill \end{array}$$

(19)

where $\tilde{p}=P(\overline{p},\overline{s})$ and ${\tilde{s}}_1=S_1(\overline{p},\overline{s})\tilde{p}$, and

$$\overline{\mu }(\overline{p},\overline{s})=\langle \overline{\mu }(\overline{p}+\tilde{p},\overline{s})\rangle ,\overline{\chi }(\overline{p},\overline{s})=\langle \overline{\chi }(\overline{p}+\tilde{p},\overline{s})\rangle ,$$

(20)

$$\begin{array}{c}W(\overline{p},\overline{s})=\langle \kappa ((\overline{p}+\tilde{p})\chi (\overline{p}+\tilde{p},\overline{s})-\overline{\chi }(\overline{p},\overline{s})\overline{p}){\overline{s}}_x+\\ {\tilde{s}}_{1\xi }\left(\kappa (\overline{p}+\tilde{p})\chi (\overline{p}+\tilde{p},\overline{s})+M_s(\overline{p}+\tilde{p},\overline{s})\right)-\mu (\overline{p}+\tilde{p},\overline{s}){(\overline{p}+\tilde{p})}_x+\overline{\mu }(\overline{p},\overline{s}){\overline{p}}_x\rangle .\end{array}$$

(21)

Equations (17)–(19) constitute a closed system regarding the averaged unknowns, $\overline{p},\overline{q},\overline{s}$, which we will call a homogenized or averaged system.

The homogenized system is similar to the original governing equations but not the same. In the generic case, the effect of the short waves is threefold. First, the averaged reaction terms, $\langle G(\overline{p}+\tilde{p},\overline{s})\rangle$ and $\langle F(\overline{p}+\tilde{p},\overline{s})\rangle$, can differ from those specified originally. Second, the same is true regarding the averaged diffusivity, $\overline{\mu }$, and sensitivity, $\overline{\chi }$. Third, last but not least, the wind term, W, can enter the flux Equation (19). The usage of the name “wind” relies on an analogy to some features in the dynamics of compressible fluid [6].

Now let the following conditions hold

$$\exists \phi =\phi (p,s):d\phi =\mu (p,s)dp-\kappa p\chi (p,s)ds.$$

(22)

In analogy with fluid mechanics, we will call this function pressure (if any). Indeed, given an existing pressure, Equation (2) represents a conservation law:

$$q_t+\nu q=-{\partial }_x\phi (p,s).$$

Note in passing that the occurrence of shock waves is natural in such a case, see [5].

When the pressure exists, Equation (12) is written as

$${\tilde{q}}_{1\tau }=-{\partial }_{\xi }\left({\left.d\phi \right|}_{p=p_0,s=\overline{s}}(p_1,s_1)\right)-{\partial }_x\phi (p_0,\overline{s})-\nu q_0-q_{0t}.$$

Averaging the last equation again leads to the conservation law as follows

$${\overline{q}}_t+\nu \overline{q}=-{\partial }_x\overline{\phi }(\overline{p},\overline{s}),\overline{\phi }(\overline{p},\overline{s})=\langle \phi (\overline{p}+\tilde{p},\overline{s})\rangle .$$

(23)

Comparing Equation (23) to Equations (19)–(21) (where one can put $M=\phi$ now) brings the wind into a simpler form

$$W=\langle \mu (\overline{p}+\tilde{p},\overline{s}){\tilde{p}}_x-\kappa {\overline{s}}_x\chi (\overline{p}+\tilde{p},\overline{s})\tilde{p}\rangle ,$$

(24)

and formulae (20) define the averaged transport coefficients, $\overline{\mu }$ and $\overline{\chi }$, again. Generally, they satisfy condition (22) neither with pressure $\overline{\phi }$ nor with any others even if it has held originally.

The most tractable averaged system arises from the signal-dependent transport coefficients—that is, for

$$\mu =\mu \left(s\right),\chi =\chi \left(s\right).$$

At the same time, there are examples of considering this assumption for some real-life species in the biological literature. We allude to articles by Tytyunov et al. [4,7], and Menolascina et al. [8] for examples and references. In such a case system (9) reduces to two uncoupled linear wave equations where the sound speed is $c=c\left(\overline{s}(x,t)\right)=\sqrt{\mu \left(\overline{s}(x,t)\right)}$—that is, it is a constant relative to the fast variables. Then the solution reads as

$$\tilde{p}=\overline{p}(x,t)({\tilde{\psi }}^{+}(\xi +c\tau )+{\tilde{\psi }}^{-}(\xi -c\tau )),c=c\left(\overline{s}(x,t)\right)=\sqrt{\mu \left(\overline{s}(x,t)\right)}$$

(25)

for every pair of function ${\tilde{\psi }}^{\pm }$ (as the sound speed is a constant relative to the fast variables). To make the solution positive, choose functions ${\psi }^{\pm }$ that are $2\pi -$periodic and positive, and set

$$\frac{1}{2\pi }\underset{-\pi }{\overset{\pi }{\int }}{\psi }^{\pm }\left(\eta \right)d\eta ={\overline{\psi }}^{\pm },{\overline{\psi }}^{+}+{\overline{\psi }}^{-}=1,{\tilde{\psi }}^{\pm }={\psi }^{\pm }-{\overline{\psi }}^{\pm }.$$

Then the solution defined by equality (25) vanishes on average, while the leading term, $p_0=\overline{p}+\tilde{p}$, is positive for every positive function $\overline{p}$. The corresponding flux, $\tilde{q}$, reads as

$$\tilde{q}=c\overline{p}(x,t)({\tilde{\psi }}^{-}(\xi -c\tau )-{\tilde{\psi }}^{+}(\xi +c\tau )),$$

(26)

The envelope of this shortwave parcel, $\overline{p}$, the averaged flux, $\overline{q}$, and the prey density, $\overline{s}$, are a solution to the homogenized system. The last gets simplified substantially, namely,

$$\overline{\mu }=\mu ,\overline{\chi }=\chi ,W=({\mu }_s\left(\overline{s}\right)+\kappa \chi \left(\overline{s}\right))\langle {\tilde{s}}_{1\xi }\tilde{p}\rangle .$$

(27)

The wind vanishes when

$${\mu }_s\left(s\right)+\kappa \chi \left(s\right)=0.$$

(28)

This is the particular form of integrability condition (22) so that expression (24) has to determine the wind too and it is equal to zero as well.

Let us proceed with simplifying the model. Assume the reaction terms to be linear in the predator density (as in the Lotka–Volterra case)—that is, set

$$G(p,s)=G_1\left(s\right)p+G_0\left(s\right),F(p,s)=F_1\left(s\right)p+F_0\left(s\right).$$

(29)

Then the reaction terms in Equations (17) and (18) remain the same as we have defined by formulae (29). Equation (13) now reads as

$${\tilde{s}}_{1\tau }-{\tilde{s}}_{1\xi \xi }=G_1\left(\overline{s}\right)\tilde{p}.$$

Hence,

$${\tilde{s}}_1=\overline{p}{G}_1\left(\overline{s}\right)\sum _{k\ne 0}\left(\frac{{\psi }_k^{+}{\mathrm{e}}^{ik(\xi +c\tau )}}{k^2+ick}+\frac{{\psi }_k^{-}{\mathrm{e}}^{ik(\xi -c\tau )}}{k^2-ick}\right).$$

where numbers ${\psi }_k^{\pm }$ are the Fourier coefficients of functions ${\psi }^{\pm }$. Finally,

$$\begin{array}{cc}& W=\theta ({\psi }^{\pm },c)G_1\left(\overline{s}\right)({\mu }_s\left(\overline{s}\right)+\kappa \chi \left(\overline{s}\right)){\overline{p}}^2,c=c\left(\overline{s}\right)=\sqrt{\mu \left(\overline{s}\right)},\end{array}$$

(30)

$$\begin{array}{ccc}& \theta ({\psi }^{\pm },c)=2c\sum _{k=1}^{\infty }\frac{|{\psi }_k^{+}{|}^2-{\left|{\psi }_k^{-}\right|}^2}{c^2+k^2}.& \hfill \end{array}$$

(31)

Thus, in the case under consideration, the homogenized system takes the form as follows

$$\begin{array}{ccc}& {\overline{p}}_t+{\overline{q}}_x=F_1\left(\overline{s}\right)\overline{p}+F_0\left(\overline{s}\right),& \hfill \end{array}$$

(32)

$$\begin{array}{ccc}& {\overline{q}}_t+\nu \overline{q}+\mu \left(\overline{s}\right){\overline{p}}_x=\kappa \chi \left(\overline{s}\right)\overline{p}{\overline{s}}_x+W,& \hfill \end{array}$$

(33)

$$\begin{array}{ccc}& {\overline{s}}_t=G_1\left(\overline{s}\right)\overline{p}+G_0\left(\overline{s}\right),& \hfill \end{array}$$

(34)

where the wind term, W, is determined by formula (30).

Let us consider the maximization of factor $\theta$ in the absolute value for a fixed $c>0$. There are restrictions for functions ${\psi }^{\pm }$ to obey; namely, they must be positive and their sum (normalized by $2\pi$) must be a probabilistic density on the unit circumference. Hence, $|{\psi }_k^{\pm }|<1\forall k\in \mathbb{N}$ and the maximizing (minimizing) pair of sequences is $|{\psi }_k^{+}|=1, {\psi }_k^{-}=0$ ($|{\psi }_k^{-}|=1, {\psi }_k^{+}=0$). So, the supremum and infimum of factor $\theta$ read as

$$\pm 2c\sum _{k=1}^{\infty }\frac{1}{c^2+k^2}=\pm \frac{\pi coth\left(\pi c\right)c-1}{c}$$

(35)

It is worth noting that functional (31) probably does not attain these extremal values at some smooth densities ${\psi }^{\pm }$. On the other hand, it tends to one of these values when one concentrates density ${\psi }^{+}$ or ${\psi }^{-}$ near some singular point on the unit circumference. The similar limiting object on the whole real axis is known as the Dirac measure, or $\delta$, or the atomic measure. To build the same on the circumference, we construct the $2\pi$-periodic grid of the Dirac atoms while setting one atom on the period.

5. The Case of Constant Diffusivity and Sensitivity

For reducing the wind term (30) into its simplest form, let us set

$$\mu =const,\chi =const,\kappa =1,G(p,s)=s(\alpha -p-s),F=p(\gamma s-\beta ),$$

(36)

where $\alpha$, $\beta$ and $\gamma$ are positive constants. Then the homogenized system reads as

$$\begin{array}{ccc}& {\overline{p}}_t+{\overline{q}}_x=\overline{p}(\gamma \overline{s}-\beta ),& \hfill \end{array}$$

(37)

$$\begin{array}{ccc}& {\overline{q}}_t+\nu \overline{q}+\mu {\overline{p}}_x=\chi \overline{p}({\overline{s}}_x-\theta \overline{p}\overline{s}),& \hfill \end{array}$$

(38)

$$\begin{array}{ccc}& {\overline{s}}_t=\overline{s}(\alpha -\overline{p}-\overline{s}),& \hfill \end{array}$$

(39)

where $\mu =const>0$, $\chi =const>0$, while formula (31) determines the magnitude of factor $\theta$ for every $c=\sqrt{\mu }=const$. Since it depends on functions ${\psi }^{\pm }$ as well, we will be considering this factor as an additional parameter, which varies within bounds (35). Actually, we can restrict our considerations within positive values of parameter $\theta$ due to the system invariance with respect to the following transform

$$(x,t,\overline{p},\overline{q},\overline{s},\theta )\to (-x,t,\overline{p},-\overline{q},\overline{s},-\theta ).$$

For $\theta =0$, system (37)–(39) is nothing else than a version of the original system (1)–(3), which possesses the mirror symmetry $(x,t,p,q,s)\to (-x,t,p,-q,s)$. Note in passing that solutions (25) and (26), over which we do the homogenization, cannot possess the above invariance unless ${\psi }^{+}(-\eta )={\psi }^{-}\left(\eta \right)\forall \eta$. The last, however, entails $|{\psi }_k^{+}|=|{\psi }_k^{-}|\forall k$ and, hence, $\theta =0$ by (31), so that the homogenized system gets the mirror invariance.

Recall now that we aim at studying the slow evolution of the amplitudes of the short-wavelength parcels. In this respect, understanding the effect of the wind seems to be the keynote as the wind is the only difference between the homogenized system (37)–(39) and the original system. Fortunately, it is particularly easy under hypotheses (36), as the wind term, $-\theta \chi \overline{s}{\overline{p}}^2$, is proportional to only one scalar parameter, $\theta$. Hence, estimating the effect of wind reduces to estimating the effect of the magnitude of parameter $\theta$, which is moreover confined within the known finite interval. It is worth noting that the same is true when the transport coefficients depends on the prey density only, and the reaction terms are linear in the predator density.

6. Linear Stability

The homogenized system (37)–(39) allows an equilibrium

$$p_e=\alpha -\beta /\gamma ,s_e=\beta /\gamma ,q_e=-\theta s_e{p}_e^2/\nu .$$

(40)

For $\theta =0$, this is the equilibrium of the original system. Here by the original system we mean the formulation of system (1)–(3) arising from assumptions (36) in the case of non-diffusive (sedentary) prey (i.e., for $\delta =0$).

In general, from formula (21), it follows that $W(\overline{p},\overline{s})=const$ provided that $\overline{p}=const$, $\overline{s}=const$ and system (9) has a two-parametric family of solutions $\tilde{p}=P(\overline{p},\overline{s}),\tilde{q}=Q(\overline{p},\overline{s})$. Hence, the equilibria of the general homogenized system (if any) are to read as

$$\overline{p}\equiv {\overline{p}}_e=const,\overline{s}\equiv {\overline{s}}_e=const{\overline{q}}_e=W({\overline{p}}_e,{\overline{s}}_e)/\nu ,$$

where the equilibrium densities are to set to zero the averaged reaction terms—that is, the following equations are to hold

$$\langle F({\overline{p}}_e+\tilde{p},{\overline{s}}_e)\rangle =0,\langle G({\overline{p}}_e+\tilde{p},{\overline{s}}_e)\rangle =0,\tilde{p}=P({\overline{p}}_e,{\overline{s}}_e).$$

Every equilibrium of the homogenized system generates the unmodulated short-wave pattern by formula $\tilde{p}=P({\overline{p}}_e,{\overline{s}}_e),\tilde{q}=Q({\overline{p}}_e,{\overline{s}}_e)$. In a simpler situation, for $\mu =\mu \left(s\right)$, $\chi =\chi \left(s\right)$, such a quasi-equilibrium pattern results from formulae (25) and (26) where $\overline{p}={\overline{p}}_e$, $\overline{q}={\overline{q}}_e$.

Applying the rationales above to the case of the linear reactions in the predator density leads to the conclusion that the equilibrium densities of species are the same for both the homogenized and original systems, i.e.,

$${\overline{p}}_e=p_e,{\overline{s}}_e=s_e.$$

(where by original we mean the version of the governing system dealt with). This observation opens a relatively easy way towards unfolding the effect of the wind, at least partly. It consists of comparing the stability of the equilibrium with the wind or without it and deciding whether the wind exerts a stabilizing or destabilizing effect. Regarding system (37)–(39), the problem simplifies to estimating the effect of the parameter $\theta$ on the stability of the above equilibria.

Let us consider the linearization of system (37)–(39) near the equilibrium (40) and look for a special solution that reads as

$$(\widehat{p},\widehat{q},\widehat{s})exp(ikx+\lambda t),\widehat{p},\widehat{q},\widehat{s},\lambda \in \mathbb{C},k\in \mathbb{R},$$

(41)

It is common to name solutions (41) as the normal modes (of the corresponding equilibrium) and say that such a mode is stable (unstable, neutral) if the real part of decrement $\lambda$ is negative (positive, equal to zero). At that, the so-called spectral stability problem arises for which $\lambda$ is the eigenvalue and the corresponding normal mode is the eigenmode. It is also common to say that an equilibrium (or more general steady regime) is stable (unstable, neutral) provided its normal mode is stable (there exists an unstable (neutral) mode). Note that genuine stability in the Lyapunov sense often needs a much more accurate definition and does not always relate to this spectral test straightforwardly. This notice concerns the instability too.

Typically, the regimes we examine for stability depend on several parameters as, for example, equilibrium (40) does. Note that the wavelength of the eigenmode, k, stands as a parameter too. In the parametric space, every point belongs to either the area of instability or the area of stability or the critical submanifold that separates them. In the former (latter) area, there exist unstable normal modes (every normal mode is stable).

A smooth path in the parametric space that parameterizes a family of equilibria parameterizes the corresponding family of stability problems. We call a submanifold in the parametric space neutral if the transversal intersection of it by such a path generically entails the transversal intersection of the imaginary axis by a branch of eigenvalues of the corresponding stability problems. The unstable eigenmodes disappearing or arising upon crossing the neutral submanifold is what we call an occurrence of instability. We say that the instability is oscillatory (monotonic) if the corresponding branch of the eigenvalues is not real (is real).

In the case dealt with here, the stability test for an individual normal mode reduces to counting the eigenvalues of some complex matrix $3\times 3$ that belong to the right semi-plane. We allude to widely recognized textbook [9] for more details, though the procedure is routine. Nevertheless, it is worth noting that the characteristic polynomials of the mentioned matrices depend not on the parameter $\theta$ itself but on ${\theta }^2$ only—that is, the direction of the wind is nothing to do with the stability. This feature also agrees with the invariance of the homogenized system.

Here is the summary of our outcomes concerning the linear stability of equilibrium (40) for $\theta =0$—that is, for the original system.

• The neutral submanifold is the graph of function $\chi ={\chi }_{nt}(\gamma ,p,s,k^2)$, where

  $${\chi }_{nt}=\frac{\mu \nu }{p_e{s}_e}+\frac{\gamma p_e{s_e}^2+s_e{\nu }^2+{s_e}^2\nu }{p_e{s}_e{k}^2}.$$

  (42)
  For a specific equilibrium and a specific wave number k, the normal mode is stable for $\chi <{\chi }_{nt}$ and becomes unstable otherwise.
• The critical submanifold is the graph of function $\chi ={\chi }_{cr}=\mu \nu /p_e{s}_e$. A specific equilibrium is stable provided that $\chi <{\chi }_{cr}$ and becomes unstable otherwise. This instability is always oscillatory and occurs to shorter waves; namely, for every $\chi >{\chi }_{cr}$ there exists $k_{*}>0$ such that every eigenmode is unstable (stable) provided that $\left|k\right|>k_{*}$ ($\left|k\right|<k_{*}$).
• The above conclusions mainly withstand allowing the prey diffusion, except for the shortwave instability; namely, the neutral submanifold is the graph of function

  $$\chi ={\check{\chi }}_{nt}={\delta }^2{k}^{2}A+\delta B+{\chi }_{nt},$$

  where ${\chi }_{nt}$ is function (42) while factors $A,B$ do not depend on the parameters that enter the right-hand side explicitly.
• At that, the critical submanifold is the graph of function

  $$\chi ={\check{\chi }}_{cr}=\underset{k}{min}{\check{\chi }}_{nt}=\frac{\delta ({\nu }^2+2\nu s_e+\gamma p_e{s}_e+\nu s_e\sqrt{{\nu }^2+\nu s_e+\gamma p_e{s}_e})+\mu \nu }{p_e{s}_e}>{\chi }_{cr}.$$

Again, for every $\delta >0$ a specific equilibrium is stable provided that $\chi <{\breve{\chi }}_{cr}$ and becomes unstable otherwise. This instability is always oscillatory and there exist real numbers $k_1>k_0>0$ such that every eigenmode’s mode is unstable provided that its wavenumber obeys the inequality $k_0<\left|k\right|<k_1$. Interestingly, function ${\breve{\chi }}_{cr}$ is linear in $\mu$ and in $\delta$ too.

Let us pass to the homogenized system (37)–(39) and test the linear stability for $\chi ={\chi }_{cr}$ and $\theta$ equal to its maximal absolute value (35). It turns out that the equilibrium (40) is unstable while it is stable for $\theta =0$ in accordance with what we have been saying above. Hence, the wind can make the equilibria unstable. In more detail, let $\chi$ and $\theta$ be set as above. Then for every admissible setting of the other parameters there exist numbers $k_2>k_1>0$ such that there are unstable normal modes for every wavelength $k:\left|k\right|>k_2$ and all the normal modes are stable for every $k:\left|k\right|<k_1$. Figure 2 (the right frame) is to illustrate this conclusion. Note that generically there are three different eigenvalues $\lambda$ for every k. Their expressions in the problem parameters are commonly called the eigenvalues branches. These branches remains smooth unless the eigenvalues collide, while their collisions lead to singularities. The right frame in Figure 2 shows the real part of the eigenvalue branch that turns positive when the squared wavenumber increases over some threshold. This threshold tends to zero rapidly upon increasing the diffusivity coefficient. Two more branch do not appear there since they always remain in the left complex semi-plane. It is worth mentioning that the increment, $\mathrm{Re}\lambda$, tends to a limit for $\left|k\right|\to \infty$ (which depends on other parameters), while $\mathrm{Im}\lambda =\left|k\right|\sqrt{\mu }+{O\left(\right|k|}^{-1})$, $\left|k\right|\to \infty$.

7. Numerical Results for the Homogenized System

Figure 3 shows typical patterns of the wind-induced pumping of a perturbation (upper row) or decaying it without the wind (lower row), cf. Figure 1. They both result from the numerical solutions of the Cauchy problems for system (37)–(39). Their settings are the same, except for the value of parameter $\theta$, which is equal either to zero (lower row) or to its negative absolutely maximal value (upper row). It is the negative quantity given by formula (35) where $c=\sqrt{\mu }$. The initial conditions and the values of the other parameters are as follows

$${p|}_{t=0}=p_e+0.05{\mathrm{e}}^{-\frac{x^2}{4}},s{{|}_{t=0}=s_e,q|}_{t=0}=q_e=-\theta \chi s_e{p}_e^2/\nu ,s_e=\beta /\gamma <\alpha ,p_e=\alpha -s_e>0,$$

(43)

$$\chi ={\chi }_{cr}=\mu \nu /\left(p_e{s}_e\right),\nu =1,\mu =1/2,\beta =0.4\alpha \gamma ,\alpha =3/4,\gamma =2/3.$$

(44)

From the numerical experiments, it follows that the behaviour of small finite perturbation qualitatively agrees with the linear theory of Section 6. Namely, the linear theory states the stability provided that $\theta =0$ while the other parameters take the values specified in (44), and we have observed the decay of the small finite perturbations. At the same time, we see amplifying them once the magnitude of $\theta$ has attained its maximal value with no changes in the other parameters. This observation also agrees with the linear theory, which predicts the onset of instability upon following such a path in the parametric space.

We emphasize the typicality of the outlined picture. We have tried not a few parameter combinations with stronger or weaker resistance, faster growth of prey ($\alpha >1$), weaker prey diffusivity, a moderate overabundance of prey ($\alpha =0.6\beta \gamma$), a greater predation efficiency ($\gamma >1$) and smaller values of parameter $\theta$ (e.g., produced by a single harmonic). Nevertheless, we have not observed any patterns of behaviour, linear or nonlinear, which differ qualitatively from those discussed above.

For numerical solving of the Cauchy problems for system (37)–(39) with data (43) and (44), we have eliminated the unknown q from Equations (37) and (38) and then have arrived at a system that consists of the reaction–diffusion Equation (39) and a nonlinear wave equation for the predator density, p. For completing the reduction, it is necessary to set up an ad hoc initial condition of form $p_t=a$ when $t=0$, where function $a=a\left(x\right)$ is given. By Equation (37), $a=F(p,s)-q_x$, where F is the predators’ sources density, and functions p, q, s are specified by the initial condition set up for the first-order system (37)–(39). Given the initial data (43) and (44) and Equation (37), we set $p_t{|}_{t=0}=0.$ This is not the case for the Cauchy problem of Section 3, where we have addressed in fact Equations (37) and (38) for $\theta =0$ and initial data (5). Given this, we arrive at $p_t+q_x{|}_{t=0}=0.$

In fact, instead of the Cauchy problem, we had been solving an initial-boundary value problem with Neumann boundary conditions on a sufficiently wide spatial interval with the built-in solver of Maple. We set the interval of $(-25,25)$ for computing the results, which we have been discussing in this section, and $(-35,35)$ for those of Section 3.

8. Conclusions

We have studied a model for the predator–prey community that couples purely diffusive prey with a prey-sensitive predator. For the latter, we employ the Cattaneo model, which is a notable member of the hyperbolic models family. Motivated by the preliminary numerical results presented in Section 3, we have been seeking a mechanism for pumping the slow amplitude modulation of short-wave patterns. For this purpose, we have addressed the high-frequency limit of the model for vanishing prey diffusivity. It has turned out that there exists an asymptotic solution that represents a short-wave pattern with slowly varying amplitude, for which there is an explicit closed system of equations to obey. This homogenized system resembles the original model in many respects, except for the source-like term, which enters the flux equation. We have called it the acoustic wind.

The linear stability analysis has discovered that the wind always destabilizes the equilibria of the homogenized system and, hence, the corresponding quasi-equilibria short-wave patterns. In agreement with the linear theory, the numerical solutions to the homogenized system presented in Section 7 reveal the growth of the finite perturbations of the equilibria when the wind is on and their decay when it is off. Thus, we conclude that the acoustic wind is solely responsible for the excitement and growth of the amplitude modulation of the short-wave patterns.

Thus, we have found an intrinsic mechanism for destabilizing and ruining the short-wave patterns in the Cattaneo model. It stops working only for quite a narrow class of models that obey condition (28). It is worth emphasizing that Tytyunov et al. [4,7] had earlier proposed this condition using probabilistic reasoning for modelling copepods with special regard to their ability to form small-scale mosaics upon spreading. The mentioned articles invoked Patlak–Keller–Segel modelling with transport coefficients obeying condition (28). Given the results of the present article, the hyperbolic models under condition (28) are not less suitable for studying the formation of small-scale patterns in population dynamics, and it is a challenging topic for continuing investigations.