Threshold dynamics scenario of a plants-pollinators cooperative system with impulsive effect on a periodically evolving domain

Abstract

Flowering plants depend on some animals for pollination and contribute to nourish the animals in natural environments. We call these animals pollinators and build a plants-pollinators cooperative model with impulsive effect on a periodically evolving domain. Next, we define the ecological reproduction index for single plant model and plants-pollinators system, respectively, whose threshold dynamics, including the extinction, persistence and coexistence, is established by the method of upper and lower solutions. Theoretical analysis shows that a large domain evolution rate has a positive influence on the survival of pollinators whether or not the impulsive effect occurs, and the pulse eliminates the pollinators even when the evolution rate is high. Moreover, some selective numerical simulations are still performed to explain our theoretical results.

1. Introduction and model formulation

In recent decades, interspecies interactions have received considerable attention since it is rare for a single species to survive independently in the natural realm. Based on the positive or negative influence of each species on the others, interactions between populations of different species are used to being classified as competition, predator-prey and cooperation. Animals resonate with plants in many forms during the species succession, forming a cross-species communication. A typical example is the cooperative relationships between the flowering plants and insects, such as bees and butterflies. Insects obtain nutrients and energy from pollens of plants, and they, as a pollinators, spread pollen to the plants (see Figure 1).

Figure 1. The schematic diagram of the plants-pollinators system.

The spread of species is widely present in nature. Interestingly, plant populations exhibit both the short-distance dispersal directly determined by the plants itself and the long-distance dispersal that is indirectly dependent on other media such as wind and animals [29]. Local and non-local diffusion operators have always been one of the theories to describe short-distance and long-distance dispersal. The results of [7, 13] suggest that although some current works focus on the non-local dispersal patterns of plants and their ecological significance, local Laplacian dispersal in the stochastic order sense, a basis for investigating more general dispersal phenomena, is an indispensable and effective tool for theoretical study of the diffusion phenomenon in plant populations.

The habitat of a species plays a crucial role in its dispersal. For example, the spatial structure of habitat can enhance the persistence of species survival [12]. When the boundary of varying habitat is unknown, differential equation with free boundaries is considered to characterise population dynamics. For instance, the spreading of invasive species was discussed in [11, 15, 38, 40, 41] and investigated the transmission of infected disease in [2, 16, 21, 37]. On the other hand, the range of habitats can present periodic variations. For example, the areas and depth of rivers and lakes change regularly due to the alternation of seasons. In summer, the water area becomes larger, while in winter, the water level drops, causing the habitat of aquatic species to extension and contraction. Although seasonal change is the direct source of periodic factors, the indirect coupling of habitat cannot be ignored. Like the model in [12], we would analyse the potential impacts of domain changes on long-term asymptotic behaviours such as population persistence and coexistence by partially introducing the intrinsic structure of habitats into the model, so as to achieve the purpose of proceeding with a comprehensive investigation of such ecological issues.

Therefore, habitat is called a growing or evolving domain. In reality, there are many recent advances concerning the dynamical behaviours of good mathematical models, as seen from [19, 28, 31, 36, 49]. Particularly, Montano et al. [28] discussed a diffusive two predators-one prey model with Holling-type II functional response, which showed that suitable conditions, depending on the domain evolution function and the space dimension, were introduced leading to the extinction of one predator and the stable coexistence of the surviving predator and its prey. References [31, 36] indicated that the small evolving rate had positive impact on prevention and control of disease. Nevertheless, [49] researched that the increase in domain evolution ratio would boost the spread of dengue fever.

Clearly, the size of a habitat’s population affects their survival and reproduction. For example, mammals give birth at a specific period, so the population experiences a birth pulse growth. At the same time, they will suffer the depletion of numbers mainly caused by capture. We suspect that part of the power from the harvesting pulse cancels out the power from the birth pulse, which may be one reason why species are kept at a certain number. The assumption of harvesting pulse only is also reflected in our model in the following text, which can be used to describe the artificial transient disturbance, including the release of natural enemies and the spraying of pesticides at certain fixed time points. In contrast to the evolution of the population, the impact of such disturbance of the pulse on the population system is transient and can be regarded as temporary, whereas the disturbance has a great impact on the population density or the number of individuals. Classical differential equations are not suitable for describing such phenomena, in which the important drivers are non-continuous process. Hence, impulsive differential equations are employed to describe the evolution of population under transient perturbations. Mil’man and Myshkis’ work in the 1960s [27] launched the theoretical study of impulsive differential equations, which has progressed since the 1980s.

Recently, based on the complicated dynamics induced by pulses, impulsive differential systems have been deeply explored by many scholars. Since drugs are frequently given into the body as pulses via oral or injection in the treatment of diseases, impulsive ordinary differential equations can be utilised to examine the dynamics of infectious diseases [17, 35]. In the dynamics of population ecology, the impulse equation model is often used to describe the occurrence of population numbers in a short period of time. Such as many species such as fish or large mammal populations will experience a birth pulse growth. Pulse models have been used in many population ecosystems, such as predator-prey systems, pest management systems, and systems with control strategies. With the change of seasons, the pulse phenomenon and population spread simultaneously affect the survival of the population. In particular, Lewis and Li [22] discussed how a seasonal birth pulse influences population dynamics, incorporating spreading speeds, travelling wave speeds, minimal domain size, as well as complicated bifurcations, in a response diffusion model with a seasonal birth pulse. Later, [46] considered the non-local dispersal stage of the system into account, established the threshold-type dynamics of the system with bounded domain, and proved the existence of a spreading speed in unbounded domain. Interestingly, impulsive harvesting was introduced into the free boundary problem and periodically evolving domain problem [24, 25], which aim to study the impacts of their combinations on persistence and extinction of species.

Flowering plants depend on pollinators (usually insects) for pollination and contribute to nourish them in natural environments. In this paper, we consider the cooperative relationship between plants and pollinators, and build a plants-pollinators model with impulsive effect on a periodically evolving domain. This model depicts the case in which impulsive effect, described by a function $g$ , occurs at every time $nT$ $(n=0,1 \ldots )$ with impulsive period $T\gt 0$ throughout the continuous growth and dispersal process of a population. During the dispersal stage, species $P$ and $H$ diffuse by the coefficients $d_1$ and $d_2$ $(\gt 0)$ , respectively. Inspired by [25, 42], the model is introduced as follows

(1.1) \begin{equation} \left \{\begin{aligned} &\frac{\partial P}{\partial t}-d_1\frac{\partial ^2P}{\partial x^2}=P\left (\gamma _1-cP+\frac{\alpha _1H}{aP+bH}\right ),&&x\in (0,l(t)), t\gt 0,\\[5pt] &\frac{\partial H}{\partial t}-d_2 \frac{\partial ^2H}{\partial x^2}=H\left (-\gamma _2+\frac{\alpha _2P}{aP+bH}\right ),&&x\in (0,l(t)), t\in ((nT)^+, (n+1)T], \\[5pt] &P(t,0)=P(t,l(t))=0, H(t,0)=H(t,l(t))=0, &&t\gt 0,\\[5pt] &P(0,x)=P_0(x), H(0,x)=H_0(x), &&x\in [0,l(0)],\\[5pt] &H((nT)^+, x)=g(H(nT,x)), &&x\in (0,l(0)), \end{aligned}\right. \end{equation}

where $P$ and $H$ denote the population densities of plants and pollinators, respectively. $\gamma _1$ represents the intrinsic growth rate of the plants, $\gamma _2$ is the pollinators’ per capita mortality rate, $a$ denotes the half-saturation constant of plants, $b$ is called the saturation effect of pollinators, $c$ is the carrying capacity of plants. $\alpha _1$ represents the plants’ efficiency in translating plants-pollinators interactions into fitness, while $\alpha _2$ shows the corresponding value for the pollinators. Moreover, it suffices to guarantee that this solution makes sense in the case where $\gamma _2$ cannot be too large. In fact, we need only take $\gamma _2\lt \frac{\alpha _2}{a}$ . The term $H((nT)^+, x)=g(H(nT,x))$ shows that the density of pollinators at the end of the pulse is the function $g$ of the density of pollinators at the start of the pulse. We always take $n=0,1,2,\ldots$ unless otherwise stated. All coefficients are positive.

In the current paper, we make the following assumptions about the impulsive function $g$ :

1. (A ₁) $g(H)$ is the first order continuously differentiable for $H\geq 0$ , $g(0)=0$ , $g'(0)\gt 0$ , and for $H\gt 0$ , $g(H)\gt 0$ , $g(H)/H$ is nonincreasing with respect to $H$ and $0\lt g(H)/H\leq 1$ .

2. (A ₂) $g(H)$ is nondecreasing with respect to $H\geq 0$ .

3. (A ₃) There are positive constants $D, \nu \gt 1$ and small $\sigma$ such that $g(H)\geq g'(0)H-DH^\nu$ for $0\leq H\leq \sigma$ .

The impulsive function $g$ satisfying assumptions ( $\boldsymbol{\mathrm{A}_1}$ ), ( $\boldsymbol{\mathrm{A}_2}$ ) and ( $\boldsymbol{\mathrm{A}_3}$ ) usually take the form of linear function $g(H)=H$ and the Beverton-Holt function:

(1.2) \begin{equation} g(H)=\frac{n_1H}{n_2+H} \end{equation}

with $n_1\gt 0$ and $n_2\gt 0$ as in [4].

Similar to [10, 32], we let $(0,l(t))$ be a periodically evolving domain with the moving boundary $l(t)$ , and assume that any point $x(t)\in (0,l(t))$ still satisfies $x(t+T)=x(t)$ . Due to the principle of mass conservation and Reynolds transport theorem [1], the evolution of domain $l(t)$ generates the spacial flow velocity $\textbf{a}$ . Inspired by [3], the evolution of domain introduces two types of extra terms, one is the dilution terms $P(\nabla \cdot \textbf{a})$ and $H(\nabla \cdot \textbf{a})$ in terms of local volume expansion (see details for [5]), another is the advection terms $\textbf{a}\cdot \nabla P$ and $\textbf{a}\cdot \nabla H$ , which represent the transport of material around $(0,l(t))$ at a rate determined by the flow $\textbf{a}$ . Therefore, problem (1.1) can be converted to the following problem

(1.3) \begin{equation} \left \{\begin{aligned} &\frac{\partial P}{\partial t}-d_1 \frac{\partial ^2P}{\partial x^2}+\textbf{a}\cdot \nabla P+P(\nabla \cdot \textbf{a})\\[5pt] &=P\left (\gamma _1-cP+\frac{\alpha _1H}{aP+bH}\right ),&&x\in (0,l(t)), t\gt 0,\\[5pt] &\frac{\partial H}{\partial t}-d_2 \frac{\partial ^2H}{\partial x^2}+\textbf{a}\cdot \nabla H+H(\nabla \cdot \textbf{a})\\[5pt] &=H\left (-\gamma _2+\frac{\alpha _2P}{aP+bH}\right ),&&x\in (0,l(t)), t\in ((nT)^+, (n+1)T], \\[5pt] &P(t,0)=P(t,l(t))=0,H(t,0)=H(t,l(t))=0, &&t\gt 0,\\[5pt] &P(0,x)=P_0(x), H(0,x)=H_0(x), &&x\in [0,l(0)],\\[5pt] &H((nT)^+, x)=g(H(nT,x)), &&x\in (0,l(0)). \end{aligned}\right. \end{equation}

To circumvent the complexities caused by the advection and dilution terms, we modify problem (1.3) from the evolving domain into the fixed domain by employing Lagrangian transformations [3, 26]. Hence, we assume that the evolution of domain is uniform and isotropic; that is, the domain evolves by the same ratio in all directions as time flies. One possibility can be denoted as

(1.4) \begin{equation} x(t)=\rho (t)y,\, y\geq 0, \end{equation}

where the positive continuous function $\rho (t)$ represents the evolution rate of domain, and $\rho (t)$ is $T$ -periodic in time $t$ , that is,

\begin{equation*} \rho (t)=\rho (t+T) \end{equation*}

for some $T\gt 0$ and $\rho (0)\equiv 1$ . Assuming $l(0)=l_0$ , we rewrite the evolving domain as $(0,l(t))=(0,\rho (t)l_0)$ . Thus, $P$ and $H$ can be mapped as a new function with the definition

(1.5) \begin{equation} P(x,t)=p(y,t),\quad H(x,t)=h(y,t), \end{equation}

then problems (1.4) and (1.5) yield that

\begin{equation*} p_t=\frac {\partial P}{\partial t}+\textbf {a}\cdot \nabla P,\quad h_t=\frac {\partial H}{\partial t}+\textbf {a}\cdot \nabla H, \end{equation*}

\begin{equation*} \nabla \textbf {a}=\frac {\dot \rho (t)}{\rho (t)},\quad P_{xx}=\frac {1}{\rho ^2(t)} p_{yy},\quad H_{xx}=\frac {1}{\rho ^2(t)} h_{yy}. \end{equation*}

As a result, problem (1.3) is changed into the following problem on a fixed domain

(1.6) \begin{equation} \left \{ \begin{aligned} &p_t-\frac{d_1}{\rho ^2(t)} p_{yy}=-\frac{\dot \rho (t)}{\rho (t)}p+p\left (\gamma _1-cp+\frac{\alpha _1h}{ap+bh}\right ), &&y\in (0,l_0), t\gt 0,\\[5pt] &h_t-\frac{d_2}{\rho ^2(t)} h_{yy}=-\frac{\dot \rho (t)}{\rho (t)}h+h\left (-\gamma _2+\frac{\alpha _2p}{ap+bh}\right ), &&y\in (0,l_0), t\in ((nT)^+, (n+1)T],\\[5pt] &p(t,0)=p(t,l_0)=0, h(t,0)=h(t,l_0)=0, &&t\gt 0,\\[5pt] &p(0,y)=p_0(y), h(0,y)=h_0(y), &&y\in [0,l_0],\\[5pt] &h((nT)^+, y)=g(h(nT,y)), &&y\in (0,l_0). \end{aligned}\right. \end{equation}

The structure of the current paper is organised as follows. In the next Section, we are devoted to investigating the threshold dynamics scenario of the plants-pollinators system, which is the core work of the current paper. Firstly, we discuss the threshold dynamics of plant model without pollinators, and present extinction-persistence phenomenon by the ecological reproduction index $R_0^1$ . Then we are concerned with investigating the threshold dynamics of plants-pollinators system, the ecological reproduction index $R_0^2$ of pulse problem is introduced by an explicit formula. Finally, we consider threshold-type results for the asymptotic behaviour of the solution to problem (1.6). Furthermore, numerical simulations are performed to understand the impacts of the domain evolution rate and impulsive effect on the dynamics of the population in Section 3. In Section 4, we end our investigation with a brief discussion.

2. Threshold scenario of plants-pollinators system

In this section, we focus on investigating the threshold dynamics scenario of problem (1.6). In plants world without pollinators, we discuss the dynamical behaviours of (1.6) with $h\equiv 0$ and present the sharp persistence-extinction dichotomy of plants. Besides, we are devoted to studying the dynamical behaviours of plants-pollinators system, in which we first establish the ecological reproduction index $R_0^2$ , then obtain the extinction-coexistence dichotomy of (1.6) by $R_0^2$ .

2.1. The threshold dynamics in plants world

In this subsection, as the starting point of the further investigation, we first discuss the dynamical behaviours in plants world in the absence of pollinators and present a considerably widespread persistence-extinction phenomenon of plants by the ecological reproduction index $R_0^1$ , which is analogous to the basic reproduction number in epidemiology [45].

In reality, we consider the corresponding periodic problem for system (1.6) with $h\equiv 0$ as follows

(2.1) \begin{equation} \left \{ \begin{aligned} &p_t-\frac{d_1}{\rho ^2(t)} p_{yy}=-\frac{\dot \rho (t)}{\rho (t)}p+p(\gamma _1-cp), &&y\in (0,l_0),t\gt 0,\\[5pt] &p(t,0)=p(t,l_0)=0,&&t\gt 0,\\[5pt] &p(0,y)=p(T,y), &&y\in [0,l_0]. \end{aligned}\right. \end{equation}

Furthermore, we focus on discussing the dynamical behaviours to problem (2.1). To address this, linearising problem (2.1) at $p=0$ , we obtain the following periodic parabolic eigenvalue problem

(2.2) \begin{equation} \left \{ \begin{aligned} &\psi _t-\frac{d_1}{\rho ^2(t)} \psi _{yy}=\frac{\gamma _1}{\mu }\psi -\frac{\dot \rho (t)}{\rho (t)}\psi, &&y\in (0,l_0),t\gt 0,\\[5pt] &\psi (t,0)=\psi (t,l_0)=0,&&t\gt 0,\\[5pt] &\psi (y,0)=\psi (y,T), &&y\in [0,l_0]. \end{aligned}\right. \end{equation}

Denote $R_0^1\;:\!=\;\mu _0^1$ , where $\mu _0^1$ is the principal eigenvalue of problem (2.2). We have the following statements and can refer to [44, 49] for more details.

Lemma 2.1. sign $(1-R_0^1)$ =sign $(\lambda _0^1)$ , where $\lambda _0^1$ is the principal eigenvalue of the following eigenvalue problem

(2.3) \begin{equation} \left \{ \begin{aligned} &\psi _t-\frac{d_1}{\rho ^2(t)} \psi _{yy}=\gamma _1\psi -\frac{\dot \rho (t)}{\rho (t)}\psi +\lambda _0\psi, &&y\in (0,l_0),t\gt 0,\\[5pt] &\psi (t,0)=\psi (t,l_0)=0,&&t\gt 0,\\[5pt] &\psi (y,0)=\psi (y,T), &&y\in [0,l_0]. \end{aligned}\right. \end{equation}

Proof. For any fixed $\mu \gt 0$ , we first consider the eigenvalue problem

\begin{equation*} \left \{ \begin {aligned} &\psi _t-\frac {d_1}{\rho ^2(t)} \psi _{yy}=\frac {\gamma _1}{\mu }\psi -\frac {\dot \rho (t)}{\rho (t)}\psi +\lambda _0\psi, &&y\in (0,l_0),t\gt 0,\\[5pt] &\psi (t,0)=\psi (t,l_0)=0,&&t\gt 0,\\[5pt] &\psi (y,0)=\psi (y,T), &&y\in [0,l_0], \end {aligned}\right. \end{equation*}

and regard $\lambda _0^*$ as its principal eigenvalue. According to [9], we can deduce that $\lambda _0^*(\mu )$ is continuous and strictly increasing for $\mu$ . Furthermore, it follows from the uniqueness of principal eigenvalue that $\lambda _0^1=\lambda _0^*(1)$ and $\lambda _0^*(\mu _0^1)=0$ .

As noted by [9] that $\lambda _0^*(\mu )$ satisfies

\begin{equation*} \lim _{\mu \rightarrow 0^+}\lambda _0^*(\mu )\lt 0,\quad \lim _{\mu \rightarrow \infty }\lambda _0^*(\mu )\gt 0, \end{equation*}

and the monotonicity of $\lambda _0^*(\mu )$ can provide the conclusion that $R_0^1=\mu _0^1$ is the unique positive solution for the equation $\lambda _0^*(\mu )=0$ . Due to

\begin{equation*} \lambda _0^1=\lambda _0^*(1)-\lambda _0^*(\mu _0^1)=\lambda _0^*(1)-\lambda _0^*(R_0^1), \end{equation*}

together with the monotonicity implies sign $(1-R_0^1)$ =sign $(\lambda _0^1)$ .

Next, we focus on the existence and attractivity of the periodic solutions of problem (2.1), and similar results can be found in [19].

Theorem 2.2.

1. (i) If $R_0^1\gt 1,$ problem (2.1) admits a unique positive steady-state solution $p^*,$ which is globally asymptotically stable;

2. (ii) If $R_0^1\leq 1,$ the trivial solution $0$ of problem (2.1) is globally asymptotically stable.

Proof. (i) If $R_0^1\gt 1$ , that is, $\lambda _0^1\lt 0$ , we let $\underline{p}=\varepsilon ^*\psi$ , where $\psi$ is the positive principal eigenfunction corresponding to the principal eigenvalue $\lambda _0^1$ of periodic parabolic eigenvalue problem (2.3), we normalise $\psi$ such that $\Vert \psi \Vert _{C([0,l_0]\times [0,T])}=1$ . Then, $\underline{p}=\varepsilon ^*\psi$ is a lower solution of problem (2.1) for small enough $\varepsilon ^*$ with $0\lt \varepsilon ^*\leq \frac{-\lambda _0^1}{c}$ . Assume that $\overline{p}$ is a constant satisfying $\overline{p}\gt \max \limits _{t\in [0,T]}\left \{\frac{\gamma _1-\frac{\dot \rho (t)}{\rho (t)}}{c} \right \}$ , thus $\overline{p}$ is a upper solution of problem (2.1). As a result, a positive steady-state solution $p^*$ of problem (2.1) exists and satisfies the following equation

\begin{equation*} \left \{ \begin {aligned} &p^*_t-\frac {d_1}{\rho ^2(t)} p^*_{yy}=(\gamma _1-cp^*)p^*-\frac {\dot \rho (t)}{\rho (t)}p^*, &&y\in (0,l_0), t\gt 0,\\[5pt] &p^*(t,0)=p^*(t,l_0)=0,&&t\gt 0,\\[5pt] &p^*(0,y)=p^*(T,y), &&y\in [0,l_0]. \end {aligned}\right. \end{equation*}

Motivated by [18, Theorem 27.1], employing the concavity of the nonlinearity term $-\frac{\dot \rho (t)}{\rho (t)}p+p(\gamma _1-cp)$ in problem (2.1), we can obtain that the positive steady-state solution $p^*$ of problem (2.1) is unique and attracts all positive solutions of problem (2.1).

(ii) When $R_0^1\leq 1$ , Lemma 2.1 implies $\lambda _0^1\geq 0$ . Then, one can follow the statement in the proof of [18, Theorem 28.1] to deduce that trivial solution 0 is globally asymptotically stable. This completes the proof.

What’s more, the following result can also be regarded as a sufficient condition to ensure $R_0^1\gt 1$ .

Remark 2.3. Suppose that $p(y,t)$ is a positive $T$ -periodic solution of problem (2.1). Then

\begin{equation*} p_t-\frac {d_1}{\rho ^2(t)} p_{yy}-\left (\gamma _1-\frac {\dot \rho (t)}{\rho (t)}\right )p=-cp^2\lt 0. \end{equation*}

According to [18, Theorem 16.6 and Remark 16.7], we can obtain $\lambda _0^1\lt 0$ , which together with Lemma 2.1 gives $R_0^1\gt 1$ .

2.2. The threshold dynamics in plants-pollinators world

Since the case that $R_0^1\leq 1$ only leads to the extinction of plants, and in view of the pollinators’ reliance on plants, henceforth, we focus on the comprehensive investigations on the threshold dynamics of plants-pollinators system given that $R_0^1\gt 1$ .

2.2.1. The ecological reproduction index

We are going to define the ecological reproduction index $R_0^2$ associated with the impulsive problem (1.6), which plays a more critical role than $R_0^1$ . Linearising the problem (1.6) around $(p, h)=(p^*, 0)$ , then there exists

(2.4) \begin{equation} \left \{ \begin{aligned} &p_t-\frac{d_1}{\rho ^2(t)} p_{yy}=\frac{\alpha _1}{a}h+(\gamma _1-2cp^*)p-\frac{\dot \rho (t)}{\rho (t)}p, &&y\in (0,l_0), t\gt 0,\\[5pt] &h_t-\frac{d_2}{\rho ^2(t)} h_{yy}=\frac{\alpha _2}{a}h-\gamma _2h-\frac{\dot \rho (t)}{\rho (t)}h, &&y\in (0,l_0), t\in ((nT)^+, (n+1)T],\\[5pt] &p(t,0)=p(t,l_0)=0, h(t,0)=h(t,l_0)=0, &&t\gt 0,\\[5pt] &p(0,y)=p_0(y), h(0,y)=h_0(y), &&y\in [0,l_0],\\[5pt] &h((nT)^+, y)=g'(0)h(nT,y), &&y\in (0,l_0). \end{aligned}\right. \end{equation}

Furthermore, the second equation of problem (2.4) can be rewritten as

(2.5) \begin{equation} \left \{ \begin{aligned} &h_t-\frac{d_2}{\rho ^2(t)} h_{yy}=\left (\frac{\alpha _2}{a}-\gamma _2-\frac{\dot \rho (t)}{\rho (t)}\right )h, &&y\in (0,l_0), t\in ((nT)^+, (n+1)T],\\[5pt] &h(t,0)=h(t,l_0)=0, &&t\gt 0,\\[5pt] &h((nT)^+, y)=g'(0)h(nT,y), &&y\in (0,l_0). \end{aligned}\right. \end{equation}

To begin with, we consider the following auxiliary problem

(2.6) \begin{equation} \left \{ \begin{aligned} &h_t=\left (\frac{\alpha _2}{a}-\gamma _2-\frac{\dot \rho (t)}{\rho (t)}\right )h, &&t\in ((nT)^+, (n+1)T],\\[5pt] &h((nT)^+)=g'(0)h(nT). \end{aligned}\right. \end{equation}

Motivated by [8], denote $B(t,\omega )$ the evolution operator of problem

(2.7) \begin{equation} \left \{ \begin{aligned} &h_t=-\gamma _2-\frac{\dot \rho (t)}{\rho (t)}h, &&t\in ((nT)^+, (n+1)T], \\[5pt] &h((nT)^+)=g'(0)h(nT). \end{aligned}\right. \end{equation}

Based on the theory of linear impulsive equations in [6], the evolution operator $B(t,\omega )$ , $t\geq \omega$ , associated with problem (2.7), can be written as

\begin{equation*} B(t,\omega )=e^{-\int _\omega ^t[\gamma _2+\frac {\dot \rho (\tau )}{\rho (\tau )}]\mathrm {d}\tau }(g'(0))^k, \end{equation*}

where $k$ denotes the number of impulsive points on $[\omega,t)$ . Due to the boundedness of $\int _\omega ^t\frac{\dot \rho (\tau )}{\rho (\tau )}\mathrm{d}\tau$ , there exists a positive constant $K$ such that

\begin{equation*} \Vert B(t,\omega )\Vert \leq K, \,t\geq \omega, \,\omega \in \mathbb {R}. \end{equation*}

Furthermore, let $C_T$ be the Banach space given by

\begin{equation*} C_T=\left \{\varsigma |\varsigma \in C\left ((nT,(n+1)T], \varsigma (t+T)=\varsigma (t) \,\text {for}\, t\in \mathbb {R}, \varsigma ((nT)^+)=\varsigma (((n+1)T)^+), n\in \mathbb {Z}\right ) \right \}, \end{equation*}

which is equipped with the maximum norm $\Vert \varsigma \Vert =\sup _{t\in [0,T]}|\varsigma (t)|$ and the positive cone $C_T^+\;:\!=\;\left \{\varsigma \in C_T|\varsigma (t)\geq 0,\forall t\in \mathbb{R}\right \}$ .

As in [23, 34], the linear operator $\mathfrak{L}\;:\;C_T\rightarrow C_T$ can be introduced by

\begin{equation*} [\mathfrak {L}\varsigma ](t)=\int _0^{+\infty }\frac {\alpha _2}{a} B(t,t-\omega )\varsigma (t-\omega )\mathrm {d}\omega, \end{equation*}

which is called as the next-generation operation. It is easily seen that $\mathfrak{L}$ is continuous, compact on $C_T\times C_T$ and positive (namely, $\mathfrak{L}(C_T^+\times C_T^+)\subset (C_T^+\times C_T^+)$ ). Therefore, we define the spectral radius of $\mathfrak{L}$ as the basic reproduction number of problem (2.6), that is,

\begin{equation*} R_0\;:\!=\;r(\mathfrak {L}). \end{equation*}

Besides, we have the following significant conclusions.

Lemma 2.4. $R_0^2=\mu _0^2$ , where $\mu _0^2$ is the principal eigenvalue of the following periodic parabolic eigenvalue problem

(2.8) \begin{equation} \left \{ \begin{aligned} &\phi _t-\frac{d_2}{\rho ^2(t)} \phi _{yy}=\frac{1}{R_0^2}\frac{\alpha _2}{a}\phi -\gamma _2\phi -\frac{\dot \rho (t)}{\rho (t)}\phi, &&y\in (0,l_0), t\in ((nT)^+, (n+1)T], \\[5pt] &\phi (t,0)=\phi (t,l_0)=0, &&t\gt 0,\\[5pt] &\phi (0,y)=\phi (T,y), &&y\in [0,l_0],\\[5pt] &\phi ((nT)^+,y)=g'(0)\phi (nT,y), &&y\in (0,l_0). \end{aligned}\right. \end{equation}

For impulsive problem, its basic reproduction number theory is not established completely. However, motivated by the above, we can present the ecological reproduction index $R_0^2$ of impulsive problem by solving problem (2.8), which can provide an explicit formula.

Theorem 2.5. The ecological reproduction index of problem (2.8) can be specifically represented as

(2.9) \begin{equation} R_0^2=\frac{\frac{\alpha _2}{a}}{\gamma _2+\frac{d_2\lambda ^*}{T}\int _0^T \frac{1}{\rho ^2(t)}\mathrm{d}t-\frac{1}{T}\ln g'(0)}, \end{equation}

where $\lambda ^*$ $(\gt 0)$ is the principal eigenvalue of $-\partial _{yy}$ in $(0, l_0)$ under the Dirichlet boundary condition.

Proof. Set

\begin{equation*} \phi (y,t)=q(t)\varphi (y), \end{equation*}

where $\varphi (y)$ is the eigenfunction corresponding to $\lambda ^*$ in the Cauchy problem

(2.10) \begin{equation} \left \{ \begin{aligned} &-\varphi _{yy}=\lambda ^*\varphi, &&y\in (0,l_0),\\[5pt] &\varphi (0)=\varphi (l_0)=0. \end{aligned}\right. \end{equation}

Hence, problem (2.8) becomes

(2.11) \begin{equation} \left \{ \begin{aligned} &\dot{q}(t)\varphi (y)-\frac{d_2}{\rho ^2(t)}q(t) \varphi _{yy}\\[5pt] &=\left (\frac{1}{R_0^2}\frac{\alpha _2}{a}-\gamma _2-\frac{\dot \rho (t)}{\rho (t)}\right )q(t)\varphi (y), &&y\in (0,l_0), t\in ((nT)^+, (n+1)T],\\[5pt] &\varphi (0)=\varphi (l_0)=0, \\[5pt] &q(0)=q(T), \\[5pt] &q((nT)^+)=g'(0)q(nT). \\[5pt] \end{aligned}\right. \end{equation}

Substituting into problem (2.10), then the first equation of (2.11) becomes

\begin{equation*} \frac {\dot {q}(t)}{q(t)}=\frac {1}{R_0^2}\frac {\alpha _2}{a}-\gamma _2-\frac {\dot \rho (t)}{\rho (t)}-\frac {d_2\lambda ^*}{\rho ^2(t)}, \end{equation*}

by solving the above equation, we derive

\begin{equation*} q(t)=Ce^{\int _0^t[\frac {1}{R_0^2}\frac {\alpha _2}{a}-\gamma _2-\frac {\dot \rho (\tau )}{\rho (\tau )}-\frac {d_2\lambda ^*}{\rho ^2(\tau )}]\mathrm {d}\tau }, \end{equation*}

where the initial value $C$ satisfies $C=q(0^+)=g'(0)q(0)$ . We rewrite

\begin{equation*} q(T)=g'(0)q(0)e^{\int _0^T[\frac {1}{R_0^2}\frac {\alpha _2}{a}-\gamma _2-\frac {\dot \rho (\tau )}{\rho (\tau )}-\frac {d_2\lambda ^*}{\rho ^2(\tau )}]\mathrm {d}\tau }, \end{equation*}

due to the periodicities of $q$ and $\rho$ , which yields

\begin{equation*} \frac {1}{g'(0)}=e^{\int _0^T[\frac {1}{R_0^2}\frac {\alpha _2}{a}-\gamma _2-\frac {d_2\lambda ^*}{\rho ^2(\tau )}]\mathrm {d}\tau }. \end{equation*}

Thus, we have

\begin{equation*} R_0^2=\frac {\frac {\alpha _2}{a}}{\gamma _2+\frac {d_2\lambda ^*}{T}\int _0^T \frac {1}{\rho ^2(t)}\mathrm {d}t-\frac {1}{T}\ln g'(0)}. \end{equation*}

We notice that in the trivial case where $g(h)=h$ , one yields

(2.12) \begin{equation} R_0^2=\frac{\frac{\alpha _2}{a}}{\gamma _2+\frac{d_2\lambda ^*}{T}\int _0^T \frac{1}{\rho ^2(t)}\mathrm{d}t}, \end{equation}

which naturally implies that the impulsive effect does not occur.

Remark 2.6. It suffices to emphasise $g'(0)\leq 1$ to guarantee $R_0^2\gt 0$ . If $g'(0)\gt 1$ , the ecological reproduction index $R_0^2$ is meaningless without positivity. For the second equation of problem (1.6), we have the following problem

\begin{equation*} \left \{ \begin {aligned} &h_t-\frac {d_2}{\rho ^2(t)} h_{yy}\\[5pt] &=\left (\frac {\alpha _2p}{ap+bh}+M^*\right )h-\frac {\dot \rho (t)}{\rho (t)}h-\gamma _2h-M^*h, &&y\in (0,l_0), t\in ((nT)^+, (n+1)T], \\[5pt] &h(t,0)=h(t,l_0)=0, &&t\gt 0,\\[5pt] &h(0,y)=h_0(y), &&y\in [0,l_0],\\[5pt] &h((nT)^+,y)=g(h(nT,y)), &&y\in (0,l_0).\\[5pt] \end {aligned}\right. \end{equation*}

Similarly, we also consider the following periodic eigenvalue problem

\begin{equation*} \left \{ \begin {aligned} &\phi _t-\frac {d_2}{\rho ^2(t)} \phi _{yy}\\[5pt] &=\frac {1}{R_0^*}\left (\frac {\alpha _2}{a}+M^*\right )\phi -\gamma _2\phi -\frac {\dot \rho (t)}{\rho (t)}\phi -M^*\phi, &&y\in (0,l_0), t\in ((nT)^+, (n+1)T], \\[5pt] &\phi (t,0)=\phi (t,l_0)=0, &&t\gt 0,\\[5pt] &\phi (0,y)=\phi (T,y), &&y\in [0,l_0],\\[5pt] &\phi ((nT)^+,y)=g'(0)\phi (nT,y), &&y\in (0,l_0).\\[5pt] \end {aligned}\right. \end{equation*}

From the above analysis, the ecological reproduction index can be re-obtained by

\begin{equation*} R_0^*=\frac {\frac {\alpha _2}{a}+M^*}{\gamma _2+\frac {d_2\lambda ^*}{T}\int _0^T \frac {1}{\rho ^2(t)}\mathrm {d}t-\frac {1}{T}\ln g'(0)+M^*}, \end{equation*}

where $M^*=\frac{1}{T}|\ln g'(0)|$ can guarantee $R_0^*\gt 0$ .

The following conclusion is well-known, which can be found in [23, 47]. Furthermore, it follows from [25] that the result also holds for our impulsive problem.

Lemma 2.7. $sign(1-R_0^2)=sign(\lambda _0^2)$ , where $\lambda _0^2$ is the principal eigenvalue of the following eigenvalue problem

(2.13) \begin{equation} \left \{ \begin{aligned} &\phi _t-\frac{d_2}{\rho ^2(t)} \phi _{yy}=\frac{\alpha _2}{a}\phi -\gamma _2\phi -\frac{\dot \rho (t)}{\rho (t)}\phi +\lambda _0^2\phi, &&y\in (0,l_0), t\in ((nT)^+, (n+1)T],\\[5pt] &\phi (t,0)=\phi (t,l_0)=0, &&t\gt 0,\\[5pt] &\phi (0,y)=\phi (T,y), &&y\in [0,l_0],\\[5pt] &\phi ((nT)^+,y)=g'(0)\phi (nT,y), &&y\in (0,l_0).\\[5pt] \end{aligned}\right. \end{equation}

Similarly, we can also obtain $\lambda _0^2=\gamma _2+\frac{d_2\lambda ^*}{T}\int _0^T\frac{1}{\rho ^2(t)}\mathrm{d}t-\frac{\alpha _2}{a}-\frac{1}{T}\ln g'(0)$ .

To present our coexistence results, we provide the following critical lemma.

Lemma 2.8. The principal eigenvalue $\lambda _0^2$ of problem (2.13) is also an eigenvalue for the following eigenvalue problem with some strict positive eigenfunctions $(\Psi _0,\Phi _0)$ ,

(2.14) \begin{equation} \left \{ \begin{aligned} &\Psi _t-\frac{d_1}{\rho ^2(t)} \Psi _{yy}=\frac{\alpha _1}{a}\Phi +(\gamma _1-2cp^*)\Psi -\frac{\dot \rho (t)}{\rho (t)}\Psi +\Lambda \Psi, &&y\in (0,l_0), t\gt 0,\\[5pt] &\Phi _t-\frac{d_2}{\rho ^2(t)} \Phi _{yy}=\frac{\alpha _2}{a}\Phi -\gamma _2\Phi -\frac{\dot \rho (t)}{\rho (t)}\Phi +\Lambda \Phi, &&y\in (0,l_0), t\in ((nT)^+, (n+1)T], \\[5pt] &\Psi (t,0)=\Psi (t,l_0)=0, \Phi (t,0)=\Phi (t,l_0)=0, &&t\gt 0,\\[5pt] &\Psi (0,y)=\Psi (T,y),\Phi (0,y)=\Phi (T,y), &&y\in [0,l_0],\\[5pt] &\Phi ((nT)^+, y)=g'(0)\Phi (nT,y), &&y\in (0,l_0) \end{aligned}\right. \end{equation}

provided that $\lambda _0^2\lt 0$ .

Proof. Assume that $(\lambda _0^2,\varphi _0)$ is the eigenpair of problem (2.13) with $\lambda _0^2\lt 0$ and $\varphi _0\gt 0$ . Thus, $(\lambda _0^2,\varphi _0)$ satisfies

\begin{equation*} \left \{ \begin {aligned} &\Phi _t-\frac {d_2}{\rho ^2(t)} \Phi _{yy}=\frac {\alpha _2}{a}\Phi -\gamma _2\Phi -\frac {\dot \rho (t)}{\rho (t)}\Phi +\Lambda \Phi, &&y\in (0,l_0), t\in ((nT)^+, (n+1)T],\\[5pt] &\Phi (t,0)=\Phi (t,l_0)=0, &&t\gt 0,\\[5pt] &\Phi (0,y)=\Phi (T,y), &&y\in [0,l_0],\\[5pt] &\Phi ((nT)^+, y)=g'(0)\Phi (nT,y), &&y\in (0,l_0). \end {aligned}\right. \end{equation*}

We then take into account the following problem

(2.15) \begin{equation} \left \{ \begin{aligned} &\Psi _t-\frac{d_1}{\rho ^2(t)} \Psi _{yy}=\frac{\alpha _1}{a}\varphi _0+(\gamma _1-2cp^*)\Psi -\frac{\dot \rho (t)}{\rho (t)}\Psi +\lambda _0^2\Psi, &&y\in (0,l_0), t\gt 0,\\[5pt] &\Psi (t,0)=\Psi (t,l_0)=0,&&t\gt 0,\\[5pt] &\Psi (0,y)=\Psi (T,y), &&y\in [0,l_0]. \end{aligned}\right. \end{equation}

Since $p^*$ solves

\begin{equation*} \left \{ \begin {aligned} &p^*_t-\frac {d_1}{\rho ^2(t)} p^*_{yy}=(\gamma _1-cp^*)p^*-\frac {\dot \rho (t)}{\rho (t)}p^*, &&y\in (0,l_0), t\gt 0,\\[5pt] &p^*(t,0)=p^*(t,l_0)=0,&&t\gt 0,\\[5pt] &p^*(0,y)=p^*(T,y), &&y\in [0,l_0], \end {aligned}\right. \end{equation*}

as mentioned in [47] the monotonicity of the principal eigenvalue, we can obtain the following problem

\begin{equation*} \left \{ \begin {aligned} &\Psi _t-\frac {d_1}{\rho ^2(t)} \Psi _{yy}=(\gamma _1-2cp^*)\Psi -\frac {\dot \rho (t)}{\rho (t)}\Psi +\Lambda \Psi, &&y\in (0,l_0), t\gt 0,\\[5pt] &\Psi (t,0)=\Psi (t,l_0)=0,&&t\gt 0,\\[5pt] &\Psi (0,y)=\Psi (T,y), &&y\in [0,l_0] \end {aligned}\right. \end{equation*}

has a positive principal eigenvalue $\Lambda _0\gt 0$ .

Therefore, recalling the positivity of $\frac{\alpha _1}{a}$ and $\varphi _0$ together with [18, Theorem 16.6], we deduce that problem (2.15) admits a unique solution $\Psi _0(t, y)$ satisfying $\Psi _0(t, y)\gt 0$ for all $(t, y)\in [0,T]\times [0,l_0]$ . In conclusion, if the principal eigenvalue $\lambda _0^2\lt 0$ , then it is still an eigenvalue of the eigenvalue problem (2.14) with strict positive eigenfunctions $(\Psi _0,\Phi _0)$ = $(\Psi _0,\varphi _0)$ . This completes the proof of lemma.

2.2.2. The extinction dynamics of plants-pollinators populations

For the further investigations, we first present the global existence, uniqueness and some estimates of the solution $(p,h)(t,y)$ .

Hereafter, we always assume that the mortality $\gamma _2$ is large sufficiently to compensate for the positive effect due to the domain shrinking, i.e. $\gamma _2\gt A_\rho \;:\!=\;\min \limits _{t\in [0,T]}\left \{\frac{\dot \rho (t)}{\rho (t)} \right \}$ .

Lemma 2.9. Problem (1.6) admits a unique global solution $(p,h)(t,y)$ , and

\begin{equation*} \begin {aligned} &(p,h)(t,y)\in C^{1,2}((0,+\infty )\times (0,l_0))\times PC^{1,2}((0,+\infty )\times (0,l_0))\\[5pt] &\;:\!=\;\{(p,h)(t,y)\mid (p,h)(t,y)\in C^{1,2}((0,+\infty )\times (0,l_0))\times C^{1,2}((nT,(n+1)T]\times (0,l_0))\}. \end {aligned} \end{equation*}

In addition, there are positive constants $S_p$ and $S_h$ such that the solution $(p, h)(t,y)$ satisfying

\begin{equation*} (0, 0)\lt (p, h)(t, y)\leq (S_p, S_h) \end{equation*}

holds for all $t\geq 0$ , $y\in [0,l_0]$ , provided that $(0,0)\not \equiv,\leq (p_{0}, h_{0})(y)\leq (S_p, S_h)$ for $y\in [0,l_0]$ .

Proof. First, we notice that the reaction function $\frac{ph}{ap+bh}$ is locally Lipschitz continuous in the whole first quadrant by extending the definition to be zero when either $p=0$ or $h=0$ . Therefore, by employing the methods used in [48], we obtain the local existence, uniqueness and regularity of the solution of problem (1.6) defined for some $T\gt 0$ in $C^{1,2}([0,T]\times [0,l_0])\times C^{1,2}([0,T]\times [0,l_0])$ .

Next, we exhibit the estimates of solution $(p,h)(t,y)$ . Recalling the strong maximum principle, it suffices to ensure the strict positivity of solution. Moreover, it follows from the first equation of problem (1.6), we obtain that

\begin{equation*} \begin {aligned} p_t-\frac {d_1}{\rho ^2(t)} p_{yy} &=-\frac {\dot \rho (t)}{\rho (t)}p+p\left (\gamma _1-cp+\frac {\alpha _1h}{ap+bh}\right ) \\[5pt] &\leq p\left ( \gamma _1+\frac {\alpha _1}{b}-\frac {\dot \rho (t)}{\rho (t)}-cp \right )\\[5pt] &\leq p\left ( \gamma _1+\frac {\alpha _1}{b}-A_\rho -cp \right ) \end {aligned} \end{equation*}

for all $t\geq 0$ , $y\in [0,l_0]$ , which yields that

\begin{equation*} p(t,y)\leq \max \left \{\frac {1}{c}\left ( \gamma _1+\frac {\alpha _1}{b}-A_\rho \right ), \Vert p_0 \Vert _\infty \right \}{\triangleq S_p} \end{equation*}

for all $t\geq 0$ , $y\in [0,l_0]$ , where $A_\rho$ is defined above.

And similarly, from the second equation of problem (1.6), for $t\in ((nT)^+, (n+1)T],y\in (0,l_0)$ , we obtain that

\begin{equation*} \begin {aligned} h_t-\frac {d_2}{\rho ^2(t)} h_{yy} &=-\frac {\dot \rho (t)}{\rho (t)}h+h\left (-\gamma _2+\frac {\alpha _2p}{ap+bh}\right ) \\[5pt] &{\leq h\left ( -\gamma _2+\frac {\alpha _2S_p}{bh}-A_\rho \right )}. \end {aligned} \end{equation*}

Particularly, for $t\in (0^+,T]$ , $y\in (0,l_0)$ , consider initial function $h(0^+,y)=g(h_0(y))$ and use the comparison principle to deduce that $h(t,y)\leq \widehat{H}(t)$ , where $\widehat{H}(t)$ satisfies the following problem

\begin{equation*} \left \{ \begin {aligned} &\frac {\mathrm {d}\widehat {H}(t)}{\mathrm {d}t}=h\left (-A_{\rho }-\gamma _2+\frac {\alpha _2S_p}{bh}\right ), &&t\in (0^+, T],\\[5pt] &\widehat {H}(0^+)=\Vert g(h_0(y)) \Vert _\infty. \end {aligned}\right. \end{equation*}

Hence, we have

\begin{equation*}h(t,y)\leq \sup \limits _{0\lt t\leq T}\widehat {H}(t)=\max \left \{\frac {\alpha _2S_p}{b(\gamma _2+A_\rho )}, \Vert g(h_0(y)) \Vert _\infty \right \}\end{equation*}

for all $t\in [0^+,T]$ , $y\in [0,l_0]$ . Due to the inequality $0\lt g(H)/H\leq 1$ in $(\textrm{A}_1)$ , it is easily shown that

\begin{equation*}h(t,y)\leq \max \left \{\frac {\alpha _2S_p}{b(\gamma _2+A_\rho )}, \Vert h_0 \Vert _\infty \right \}\triangleq S_h\end{equation*}

for $t\in [0,T]$ , $y\in [0,l_0]$ . Taking $n=1,2,3,\ldots$ and using the same procedures, we get that $h(t,y)\leq S_h$ for $t\geq 0$ , $y\in [0,l_0]$ . The proof is completed.

The above-mentioned preliminaries allow us to investigate the asymptotic behaviours of the solution to problem (1.6).

Theorem 2.10. If $R_0^2\leq 1$ , then the solution $(p, h)(t,y)$ of problem (1.6) satisfies

\begin{equation*} \lim _{t\rightarrow \infty }(p,h)(t,y)=(p^*, 0) \end{equation*}

uniformly for $y\in [0,l_0]$ .

Proof. Let

\begin{equation*} h^\diamond (t,y)=Me^{{\frac {\alpha _2}{a}\left (1-\frac {1}{R_0^2}\right ) t}}\phi (t,y), \end{equation*}

where $\phi (t,y)$ satisfying $\Vert \phi \Vert _\infty =1$ is the normalised eigenfunction associated with $R_0^2$ , $M$ is a sufficient large positive constant to be chosen later. It deduces that

\begin{equation*} \begin {aligned} &h^\diamond _t-\frac {d_2}{\rho ^2(t)} h^\diamond _{yy}-\frac {\alpha _2p}{ap+bh^\diamond }h^\diamond +\gamma _2h^\diamond +\frac {\dot \rho (t)}{\rho (t)}h^\diamond \\[5pt] \geq &h^\diamond _t-\frac {d_2}{\rho ^2(t)} h^\diamond _{yy}-\frac {\alpha _2}{a}h^\diamond +\gamma _2h^\diamond +\frac {\dot \rho (t)}{\rho (t)}h^\diamond \\[5pt] =\;&{\frac {\alpha _2}{a}(1-\frac {1}{R_0^2}) Me^{\frac {\alpha _2}{a}(1-\frac {1}{R_0^2}) t}\phi +Me^{\frac {\alpha _2}{a}(1-\frac {1}{R_0^2}) t}\phi _t-Me^{\frac {\alpha _2}{a}(1-\frac {1}{R_0^2}) t}\frac {d_2}{\rho ^2(t)} \phi _{yy}}\\[5pt] &{-\frac {\alpha _2}{a}Me^{\frac {\alpha _2}{a}(1-\frac {1}{R_0^2}) t}\phi +\gamma _2Me^{\frac {\alpha _2}{a}(1-\frac {1}{R_0^2}) t}\phi +\frac {\dot \rho (t)}{\rho (t)}Me^{\frac {\alpha _2}{a}(1-\frac {1}{R_0^2}) t}\phi }\\[5pt] =\;&{\frac {\alpha _2}{a}(1-\frac {1}{R_0^2}) Me^{\frac {\alpha _2}{a}(1-\frac {1}{R_0^2}) t}\phi +Me^{\frac {\alpha _2}{a}(1-\frac {1}{R_0^2}) t}\left [\frac {d_2}{\rho ^2(t)} \phi _{yy}+\frac {1}{R_0^2}\frac {\alpha _2}{a}\phi -\gamma _2\phi -\frac {\dot \rho (t)}{\rho (t)}\phi \right ]}\\[5pt] &{-Me^{\frac {\alpha _2}{a}(1-\frac {1}{R_0^2}) t}\frac {d_2}{\rho ^2(t)} \phi _{yy} -\frac {\alpha _2}{a}Me^{\frac {\alpha _2}{a}(1-\frac {1}{R_0^2}) t}\phi +\gamma _2Me^{\frac {\alpha _2}{a}(1-\frac {1}{R_0^2}) t}\phi +\frac {\dot \rho (t)}{\rho (t)}Me^{\frac {\alpha _2}{a}(1-\frac {1}{R_0^2}) t}\phi }\\[5pt] =\;&{h^\diamond \left [ \frac {\alpha _2}{a}\left (1-\frac {1}{R_0^2}\right )+\frac {1}{R_0^2}\frac {\alpha _2}{a}-\frac {\alpha _2}{a} \right ]}\\[5pt] =\;&0. \end {aligned} \end{equation*}

Now we can easily verify that $h^\diamond$ is a solution of the following problem

(2.16) \begin{equation} \left \{ \begin{aligned} &h_t-\frac{d_2}{\rho ^2(t)} h_{yy}=\frac{\alpha _2}{a}h-\gamma _2h-\frac{\dot \rho (t)}{\rho (t)}h, &&y\in (0,l_0), t\in ((nT)^+, (n+1)T],\\[5pt] &h(t,0)=h(t,l_0)=0, &&t\gt 0,\\[5pt] &h(0,y)=M\phi (0,y), &&y\in [0,l_0],\\[5pt] &h((nT)^+, y)=g'(0)h(nT,y), &&y\in (0,l_0). \end{aligned}\right. \end{equation}

Recalling the assumption ( $\boldsymbol{\mathrm{A}_1}$ ), we obtain that

\begin{equation*} \frac {g(h)}{h}\leq \lim _{\varepsilon \rightarrow 0}\frac {g(\varepsilon )-g(0)}{\varepsilon }=g'(0) \end{equation*}

holds for small enough $\varepsilon \gt 0$ , which deduces that

\begin{equation*} h^\diamond ((nT)^+,y)=g'(0)h^\diamond (nT,y)\geq g(h^\diamond (nT,y)). \end{equation*}

What’s more, for any initial function $h(0,y)$ , we select a large enough constant $M$ such that $h^\diamond (0,y)\geq h(0,y)$ . Since the reaction term in problem (2.16) is larger than that in problem (1.6), it yields that $h^\diamond (t,y)$ is a upper solution of problem (1.6). Due to the comparison principle, we have

\begin{equation*} h(t,y)\leq h^\diamond (t,y), \,\,\,t\geq 0, y\in [0, l_0]. \end{equation*}

If $R_0^2\leq 1$ , one can obtain that $\lim \limits _{t\rightarrow \infty }h^\diamond (t,y)=0$ for all $y\in [0,l_0]$ . Hence, we have $\lim \limits _{t\rightarrow \infty }h(t,y)=0$ uniformly for $y\in [0,l_0]$ .

In addition, by the nearly parallel approach adopted in [31], we can also provide that $\lim \limits _{t\rightarrow \infty }p(t,y)=p^*(t,y)$ holds uniformly for all $t\geq 0, y\in [0,l_0]$ .

2.2.3. The coexistence dynamics of plants-pollinators populations

To explore the periodic steady-state coexistence solutions of problem (1.6) and their attractivity, we first consider the following periodic problem

(2.17) \begin{equation} \left \{ \begin{aligned} &p_t-\frac{d_1}{\rho ^2(t)} p_{yy}=-\frac{\dot \rho (t)}{\rho (t)}p+p\left (\gamma _1-cp+\frac{\alpha _1h}{ap+bh}\right ), &&y\in (0,l_0), t\gt 0,\\[5pt] &h_t-\frac{d_2}{\rho ^2(t)} h_{yy}=-\frac{\dot \rho (t)}{\rho (t)}h+h\left (-\gamma _2+\frac{\alpha _2p}{ap+bh}\right ), &&y\in (0,l_0), t\in ((nT)^+, (n+1)T],\\[5pt] &p(t,0)=p(t,l_0)=0, h(t,0)=h(t,l_0)=0, &&t\gt 0,\\[5pt] &p(0,y)=p(T,y), h(0,y)=h(T,y), &&y\in [0,l_0],\\[5pt] &h((nT)^+, y)=g(h(nT,y)), &&y\in (0,l_0). \end{aligned}\right. \end{equation}

The definition of upper and lower solutions to (2.17) with pulses are presented as follows.

Definition 2.11. Let $(\widetilde{p},\widetilde{h})$ and $(\widehat{p},\widehat{h})$ be a pair of ordered upper and lower solutions of problem (2.17), if $(0,0)\leq (\widehat{p},\widehat{h})\leq (\widetilde{p},\widetilde{h})\leq (S_p,S_h)$ and

(2.18) \begin{equation} \left \{ \begin{aligned} &\widetilde{p}_t-\frac{d_1}{\rho ^2(t)} \widetilde{p}_{yy}\geq -\frac{\dot \rho (t)}{\rho (t)}\widetilde{p}+\widetilde{p}\left (\gamma _1-c\widetilde{p}+\frac{\alpha _1\widetilde{h}}{a\widetilde{p}+b\widetilde{h}}\right ), &&y\in (0,l_0), t\gt 0,\\[5pt] &\widetilde{h}_t-\frac{d_2}{\rho ^2(t)} \widetilde{h}_{yy}\geq -\frac{\dot \rho (t)}{\rho (t)}\widetilde{h}+\widetilde{h}\left (-\gamma _2+\frac{\alpha _2\widetilde{p}}{a\widetilde{p}+b\widetilde{h}}\right ), &&y\in (0,l_0), t\in ((nT)^+, (n+1)T], \\[5pt] &\widehat{p}_t-\frac{d_1}{\rho ^2(t)} \widehat{p}_{yy}\leq -\frac{\dot \rho (t)}{\rho (t)}\widehat{p}+\widehat{p}\left (\gamma _1-c\widehat{p}+\frac{\alpha _1\widehat{h}}{a\widehat{p}+b\widehat{h}}\right ), &&y\in (0,l_0), t\gt 0,\\[5pt] &\widehat{h}_t-\frac{d_2}{\rho ^2(t)} \widehat{h}_{yy}\leq -\frac{\dot \rho (t)}{\rho (t)}\widehat{h}+\widehat{h}\left (-\gamma _2+\frac{\alpha _2\widehat{p}}{a\widehat{p}+b\widehat{h}}\right ), &&y\in (0,l_0), t\in ((nT)^+, (n+1)T], \\[5pt] &\widehat{p}(t,y)=0\leq \widetilde{p}(t,y), \widehat{h}(t,y)=0\leq \widetilde{h}(t,y), &&y\in [0,l_0], t\gt 0,\\[5pt] &\widetilde{p}(0,y)\geq \widetilde{p}(T,y), \widetilde{h}(0,y)\geq \widetilde{h}(T,y), &&y\in [0,l_0],\\[5pt] &\widehat{p}(0,y)\leq \widehat{p}(T,y), \widehat{h}(0,y)\leq \widehat{h}(T,y), &&y\in [0,l_0],\\[5pt] &\widetilde{h}((nT)^+, y)\geq g(\widetilde{h}(nT,y)),&&y\in (0,l_0), \\[5pt] &\widehat{h}((nT)^+, y)\leq g(\widehat{h}(nT,y)), &&y\in (0,l_0). \end{aligned}\right. \end{equation}

For further analysis, we let $f_1=-\frac{\dot \rho (t)}{\rho (t)}p+p\left (\gamma _1-cp+\frac{\alpha _1h}{ap+bh}\right )$ , $f_2=-\frac{\dot \rho (t)}{\rho (t)}h+h\left (-\gamma _2+\frac{\alpha _2p}{ap+bh}\right )$ and choose

\begin{equation*} k_1=\max _{t\in [0,T]}\left \{\frac {\dot {\rho }(t)}{\rho (t)}\right \}+S_p, \quad k_2=\max _{t\in [0,T]}\left \{\frac {\dot {\rho }(t)}{\rho (t)}\right \}+\gamma _2 \end{equation*}

such that

\begin{equation*} F_1(p,h)=k_1p+f_1(p,h), \quad F_2(p,h)=k_2h+f_2(p,h). \end{equation*}

Then, problem (2.17) can be rewritten as

\begin{equation*} \left \{ \begin {aligned} &p_t-\frac {d_1}{\rho ^2(t)} p_{yy}+k_1p=F_1(p,h), &&y\in (0,l_0), t\gt 0,\\[5pt] &h_t-\frac {d_2}{\rho ^2(t)} h_{yy}+k_2h=F_2(p,h), &&y\in (0,l_0), t\in ((nT)^+, (n+1)T],\\[5pt] &p(t,0)=p(t,l_0)=0, h(t,0)=h(t,l_0)=0, &&t\gt 0,\\[5pt] &p(0,y)=p(T,y), h(0,y)=h(T,y), &&y\in [0,l_0],\\[5pt] &h((nT)^+, y)=g(h(nT,y)), &&y\in (0,l_0). \end {aligned}\right. \end{equation*}

It is easy to verify that both $F_1(p,h)$ and $F_2(p,h)$ are nondecreasing with respect to $p$ and $h$ .

Furthermore, we consider the following iteration process associated with the initial values $(\overline{p}^{(0)},\overline{h}^{(0)})$ = $(\widetilde{p},\widetilde{h})$ , $(\underline{p}^{(0)},\underline{h}^{(0)})$ = $(\widehat{p},\widehat{h})$ and the iteration sequences $\big \{(\overline p^{(m)},\overline h^{(m)})\big \}$ and $\big \{(\underline p^{(m)},\underline h^{(m)})\big \}$ by the following process

(2.19) \begin{equation} \left \{ \begin{aligned} &\overline{p}_t^{(m)}-\frac{d_1}{\rho ^2(t)} \overline{p}_{yy}^{(m)}+k_1\overline{p}^{(m)}=F_1(\overline{p}^{(m-1)},\overline{h}^{(m-1)}), &&y\in (0,l_0),t\gt 0,\\[5pt] &\overline{h}_t^{(m)}-\frac{d_2}{\rho ^2(t)} \overline{h}_{yy}^{(m)}+k_2\overline{h}^{(m)}=F_2(\overline{p}^{(m-1)},\overline{h}^{(m-1)}), &&y\in (0,l_0), t\in ((nT)^+, (n+1)T], \\[5pt] &\underline{p}_t^{(m)}-\frac{d_1}{\rho ^2(t)} \underline{p}_{yy}^{(m)}+k_1\underline{p}^{(m)}=F_1(\underline{p}^{(m-1)},\underline{h}^{(m-1)}), &&y\in (0,l_0),t\gt 0,\\[5pt] &\underline{h}_t^{(m)}-\frac{d_2}{\rho ^2(t)} \underline{h}_{yy}^{(m)}+k_2\underline{h}^{(m)}=F_2(\underline{p}^{(m-1)},\underline{h}^{(m-1)}), &&y\in (0,l_0), t\in ((nT)^+, (n+1)T], \\[5pt] &\overline{p}^{(m)}(t,y)=\overline{h}^{(m)}(t,y)=\underline{p}^{(m)}(t,y)=\underline{h}^{(m)}(t,y)=0,&&y\in [0,l_0],t\gt 0,\\[5pt] &\overline{p}^{(m)}(y,0)=\overline{p}^{(m-1)}(y,T),\,\overline{h}^{(m)}(y,0)=\overline{h}^{(m-1)}(y,T),&&y\in [0,l_0],\\[5pt] &\underline{p}^{(m)}(y,0)=\underline{p}^{(m-1)}(y,T),\,\underline{h}^{(m)}(y,0)=\underline{h}^{(m-1)}(y,T),&&y\in [0,l_0],\\[5pt] &\overline{h}^{(m)}((nT)^+,y)=g(\overline{h}^{(m-1)}((n+1)T,y)), &&y\in (0,l_0),\\[5pt] &\underline{h}^{(m)}((nT)^+,y)=g(\underline{h}^{(m-1)}((n+1)T,y)), &&y\in (0,l_0).\\[5pt] \end{aligned}\right. \end{equation}

Motivated by [30], we present the following lemma to expound the monotone property of the iteration sequences.

Lemma 2.12. Let $(\widetilde{p},\widetilde{h})$ and $(\widehat{p},\widehat{h})$ be a pair of ordered upper and lower solutions of problem (2.17) $,$ respectively. Then, the sequence $\big \{(\overline p^{(m)},\overline h^{(m)})\big \}$ decreases and converges monotonically to $(\overline p,\overline h)$ which is a maximal T-periodic solution of problem (2.17) $,$ while the sequence $\big \{(\underline p^{(m)},\underline h^{(m)})\big \}$ increases and converges monotonically to $(\underline p,\underline h)$ which is a minimal $T$ -periodic solution of problem (2.17) $,$ that is $,$

(2.20) \begin{equation} \begin{aligned} (\widehat{p},\widehat{h})&\leq (\underline{p}^{(m)},\underline{h}^{(m)})\leq (\underline{p}^{(m+1)},\underline{h}^{(m+1)})\leq (\underline{p},\underline{h})\\[5pt] &\leq (\overline{p},\overline{h})\leq (\overline{p}^{(m+1)},\overline{h}^{(m+1)})\leq (\overline{p}^{(m)},\overline{h}^{(m)})\leq (\widetilde{p},\widetilde{h}). \end{aligned} \end{equation}

Theorem 2.13. If $R_0^2\gt 1,$ then there exists the following statements hold:

1. (i) There are a pair of minimal and maximal positive $T$ -periodic solutions $(\underline{p},\underline{h})\leq (\overline{p},\overline{h})$ of problem (2.17) over $(p^*,0),$ besides, if $(\underline{p},\underline{h})(0,y)=(\overline{p},\overline{h})(0,y)$ , then $(\underline{p},\underline{h})=(\overline{p},\overline{h})\;:\!=\;(p^{**},h^{**})$ is the unique positive $T$ -periodic solution to problem (2.17);

2. (ii) Let $(p,h)(t,y;p_0,h_0)$ be the solution of problem (1.6) with bounded and continuous initial conditions $(0,0)\not \equiv,\leq (p_0,h_0)(y)\leq (S_p,S_h)$ on $[0,l_0].$ Then $(\underline{p},\underline{h})\leq (\overline{p},\overline{h})$ is attractive in the sense that

   (2.21) \begin{equation} \begin{aligned} (\underline{p},\underline{h})(t,y)&=\liminf _{m\rightarrow \infty }(p,h)(t+mT,y;\, p_0,h_0)\\[5pt] &\leq \limsup _{m\rightarrow \infty }(p,h)(t+mT,y;\, p_0,h_0)=(\overline{p},\overline{h})(t,y) \end{aligned} \end{equation}
   holds on $[0,\infty )\times [0,l_0],$ that is $,$
   \begin{equation*} \lim _{m\rightarrow \infty }(p,h)(t+mT,y;\, p_0,h_0)\rightarrow (p^{**},h^{**})(t,y). \end{equation*}

Proof. (i) The proof of this part is not particularly difficult but is too long; thus, we divide it into two steps to help the reader understand.

Step 1 The existence of the positive $T$ -periodic solution to problem (2.17).

We first construct the upper solution of problem (2.17). Let $(\widetilde{p},\widetilde{h})$ = $(M_1U(t),M_2V(t))$ , $M_1, M_2\gt 1$ , where $(U(t),V(t))$ satisfies

(2.22) \begin{equation} \left \{ \begin{aligned} &U_t(t)=-\frac{\dot{\rho }(t)}{\rho (t)}U(t)+\gamma _1U(t)+\frac{\alpha _1}{b}U(t)-cU^2(t), &&t\gt 0,\\[5pt] &V_t(t)=-\frac{\dot{\rho }(t)}{\rho (t)}V(t)-\gamma _2V(t)+{\frac{\alpha _2U(t)}{aU(t)+bV(t)}V(t)}, &&t\in ((nT)^+, (n+1)T],\\[5pt] &U(t)=U(t+T), V(t)=V(t+T),&&t\geq 0,\\[5pt] &V((nT)^+)=g'(0)V(nT)\geq g(V(nT)).\\[5pt] \end{aligned}\right. \end{equation}

Now, it is necessary to explain the existence of solutions $U(t)$ and $V(t)$ to problem (2.22). For the equation of $U$ , by direct calculation, we can conclude that

\begin{equation*} U(t)=\frac {e^{(\gamma _1+\frac {\alpha _1}{b})t}\rho ^{-1}(t)(e^{(\gamma _1+\frac {\alpha _1}{b})T}-1)} {\int _{nT}^t \frac {ce^{(\gamma _1+\frac {\alpha _1}{b})\tau }}{\rho (\tau )}\mathrm {d}\tau (e^{(\gamma _1+\frac {\alpha _1}{b})T}-1)+e^{(\gamma _1+\frac {\alpha _1}{b})nT}\int _0^T \frac {ce^{(\gamma _1+\frac {\alpha _1}{b})\tau }}{\rho (\tau )}\mathrm {d}\tau }. \end{equation*}

And for the equation of $V$ , since $R_0^2\gt 1$ , $V(t)$ can also be obtained by using the upper and lower solution method similar to the logistic equation.

Furthermore, we have

\begin{equation*} \begin {aligned} &\frac {\partial \widetilde {p}}{\partial t}-\left [\frac {d_1}{\rho ^2(t)}\widetilde {p}_{yy}-\frac {\dot \rho (t)}{\rho (t)}\widetilde {p}+\gamma _1\widetilde {p}-c\widetilde {p}^2 +\frac {\alpha _1\widetilde {h}\widetilde {p}}{a\widetilde {p}+b\widetilde {h}}\right ]\\[5pt] \gt &\frac {\partial \widetilde {p}}{\partial t}-\left [\frac {d_1}{\rho ^2(t)}\widetilde {p}_{yy}-\frac {\dot \rho (t)}{\rho (t)}\widetilde {p}+\gamma _1\widetilde {p}-c\widetilde {p}^2 +\frac {\alpha _1}{b}\widetilde {p}\right ]\\[5pt] =\;&M_1\left [-\frac {\dot {\rho }(t)}{\rho (t)}U(t)+\gamma _1U(t)+\frac {\alpha _1}{b}U(t)-cU^2(t)\right ]-\left [\frac {d_1}{\rho ^2(t)}\widetilde {p}_{yy}-\frac {\dot \rho (t)}{\rho (t)}\widetilde {p}+\gamma _1\widetilde {p}-c\widetilde {p}^2 +\frac {\alpha _1}{b}\widetilde {p}\right ]\\[5pt] =\;&M_1\left [-\frac {\dot {\rho }(t)}{\rho (t)}U(t)+\gamma _1U(t)+\frac {\alpha _1}{b}U(t)-cU^2(t)\right ]-\left [-\frac {\dot \rho (t)}{\rho (t)}\widetilde {p}+\gamma _1\widetilde {p}-c\widetilde {p}^2 +\frac {\alpha _1}{b}\widetilde {p}\right ]\\[5pt] =\;&0, \end {aligned} \end{equation*}

which implies that $\widetilde{p}$ is the upper solution of the first equation to problem (2.17). Since the nonlinear term $\frac{\alpha _2p}{ap+bh}$ is strictly increasing in $p$ , and considering the cooperative relationship, we can deduce that $\widetilde{h}$ is the upper solution of the second equation to problem (2.17). Hence, we obtain that $(\widetilde{p},\widetilde{h})$ = $(M_1U(t),M_2V(t))$ with $M_1$ and $M_2\gt 1$ is an upper solution of problem (2.17).

In the following, we aim to consider the lower solution and define $(\widehat{p},\widehat{h})=(p^*+\delta \Psi _{\tilde{\varepsilon }},\delta \Phi _{\tilde{\varepsilon }})$ with $\widehat{h}$ satisfying

(2.23) \begin{equation} \widehat{h}(t,y)= \left \{ \begin{aligned} &\delta \Phi _{\tilde{\varepsilon }}(nT,y), &&t=nT,\\[5pt] &\delta \frac{\rho _1}{g'(0)}\Phi _{\tilde{\varepsilon }}((nT)^+,y), &&t=(nT)^+,\\[5pt] &\delta \frac{\rho _1}{g'(0)}e^{\left [-(\lambda _0^2)_{\tilde{\varepsilon }}-\xi \right ](t-nT)}\Phi _{\tilde{\varepsilon }}(t,y), &&t\in ((nT)^+, (n+1)T], \end{aligned}\right. \end{equation}

where $\xi +(\lambda _0^2)_{\tilde{\varepsilon }}\lt 0$ with a positive constant $\xi$ and $\rho _1=e^{((\lambda _0^2)_{\tilde{\varepsilon }}+\xi )T}g'(0)$ such that $\widehat{h}(nT,y)=\widehat{h}((n+1)T,y)$ . $\delta$ is a small enough positive constant to be chosen later, and the positive eigenfunctions $(\Psi _{\tilde{\varepsilon }},\Phi _{\tilde{\varepsilon }})$ satisfy

(2.24) \begin{equation} \left \{ \begin{aligned} &(\Psi _{\tilde{\varepsilon }})_t-\frac{d_1}{\rho ^2(t)} (\Psi _{\tilde{\varepsilon }})_{yy}\\[5pt] &=\frac{\alpha _1p^*}{\tilde{\varepsilon }+ap^*}\Phi _{\tilde{\varepsilon }}+(\gamma _1-2cp^*)\Psi _{\tilde{\varepsilon }} -\frac{\dot \rho (t)}{\rho (t)}\Psi _{\tilde{\varepsilon }}+\Lambda \Psi _{\tilde{\varepsilon }}, &&y\in (0,l_0), t\gt 0,\\[5pt] &(\Phi _{\tilde{\varepsilon }})_t-\frac{d_2}{\rho ^2(t)} (\Phi _{\tilde{\varepsilon }})_{yy}\\[5pt] &=\frac{\alpha _2p^*}{\tilde{\varepsilon }+ap^*}\Phi _{\tilde{\varepsilon }}-\gamma _2\Phi _{\tilde{\varepsilon }} -\frac{\dot \rho (t)}{\rho (t)}\Phi _{\tilde{\varepsilon }}+\Lambda \Phi _{\tilde{\varepsilon }}, &&y\in (0,l_0), t\in ((nT)^+, (n+1)T], \\[5pt] &\Psi _{\tilde{\varepsilon }}(t,0)=\Psi _{\tilde{\varepsilon }}(t,l_0)=0, \Phi _{\tilde{\varepsilon }}(t,0)=\Phi _{\tilde{\varepsilon }}(t,l_0)=0, &&t\gt 0,\\[5pt] &\Psi _{\tilde{\varepsilon }}(0,y)=\Psi _{\tilde{\varepsilon }}(T,y),\Phi _{\tilde{\varepsilon }}(0,y)=\Phi _{\tilde{\varepsilon }}(T,y), &&y\in [0,l_0],\\[5pt] &\Phi _{\tilde{\varepsilon }}((nT)^+, y)=g'(0)\Phi _{\tilde{\varepsilon }}(nT,y), &&y\in (0,l_0), \end{aligned}\right. \end{equation}

which is similar to problem (2.14) and can be obtained by perturbation theory for sufficiently small $\tilde{\varepsilon }$ , and one could refer to Kato [20] for more details about this point. By Lemma 2.8, we can also obtain that $(\lambda _0^2)_{\tilde{\varepsilon }}\lt 0$ if $\tilde{\varepsilon }$ is sufficiently small. For $t\in ((nT)^+, (n+1)T]$ and $y\in (0,l_0)$ , if $\delta \lt \delta _1$ , we can obtain that

\begin{equation*} \begin {aligned} &\frac {\partial \widehat {h}}{\partial t}-\left [ \frac {d_2}{\rho ^2(t)}\widehat {h}_{yy}-\frac {\dot \rho (t)}{\rho (t)}\widehat {h}+\widehat {h}\left (-\gamma _2+\frac {\alpha _2\widehat {p}}{a\widehat {p}+b\widehat {h}}\right )\right ]\\[5pt] =\;&\left [-(\lambda _0^2)_{\tilde {\varepsilon }}-\xi \right ]\times \delta \frac {\rho _1}{g'(0)}e^{[-(\lambda _0^2)_{\tilde {\varepsilon }} -\xi ](t-nT)}\Phi _{\tilde {\varepsilon }}\\[5pt] &+\delta \frac {\rho _1}{g'(0)}e^{\left [-(\lambda _0^2)_{\tilde {\varepsilon }}-\xi \right ](t-nT)}\times \left [\frac {d_2}{\rho ^2(t)} (\Phi _{\tilde {\varepsilon }})_{yy}+\frac {\alpha _2p^*}{\tilde {\varepsilon }+ap^*}\Phi _{\tilde {\varepsilon }}-\gamma _2\Phi _{\tilde {\varepsilon }}-\frac {\dot \rho (t)}{\rho (t)}\Phi _{\tilde {\varepsilon }} +(\lambda _0^2)_{\tilde {\varepsilon }}\Phi _{\tilde {\varepsilon }}\right ]\\[5pt] &-\delta \frac {\rho _1}{g'(0)}e^{\left [-(\lambda _0^2)_{\tilde {\varepsilon }}-\xi \right ](t-nT)}\times \left [\frac {d_2}{\rho ^2(t)}(\Phi _{\tilde {\varepsilon }})_{yy} -\left (\frac {\dot \rho (t)}{\rho (t)}+\gamma _2\right )\Phi _{\tilde {\varepsilon }}\right ]\\[5pt] &-\delta \frac {\rho _1}{g'(0)}e^{\left [-(\lambda _0^2)_{\tilde {\varepsilon }}-\xi \right ](t-nT)}\Phi _{\tilde {\varepsilon }}\times \left [\frac {\alpha _2(p^* +\delta \Psi _{\tilde {\varepsilon }})}{a(p^*+\delta \Psi _{\tilde {\varepsilon }})+b\delta \frac {\rho _1}{g'(0)}e^{\left [-(\lambda _0^2)_{\tilde {\varepsilon }} -\xi \right ](t-nT)}\Phi _{\tilde {\varepsilon }}} \right ]\\[5pt] =\;&\delta \frac {\rho _1}{g'(0)}e^{\left [-(\lambda _0^2)_{\tilde {\varepsilon }}-\xi \right ](t-nT)}\Phi _{\tilde {\varepsilon }} \\[5pt] &\times \left \{\left [-(\lambda _0^2)_{\tilde {\varepsilon }}-\xi \right ]+(\lambda _0^2)_{\tilde {\varepsilon }}+\frac {\alpha _2p^*}{\tilde {\varepsilon }+ap^*} -\frac {\alpha _2(p^*+\delta \Psi _{\tilde {\varepsilon }})}{a(p^*+\delta \Psi _{\tilde {\varepsilon }})+b\delta \frac {\rho _1}{g'(0)}e^{\left [-(\lambda _0^2)_{\tilde {\varepsilon }} -\xi \right ](t-nT)}\Phi _{\tilde {\varepsilon }}} \right \}\\[5pt] =\;&\delta \frac {\rho _1}{g'(0)}e^{\left [-(\lambda _0^2)_{\tilde {\varepsilon }}-\xi \right ](t-nT)}\Phi _{\tilde {\varepsilon }} \\[5pt] &\times \left \{ -\xi +\frac {-\alpha _2\tilde {\varepsilon }p^*+\delta \left [a\alpha _2p^*\Psi _{\tilde {\varepsilon }}+b\alpha _2p^*\frac {\rho _1}{g'(0)}e^{\left [-(\lambda _0^2)_{\tilde {\varepsilon }} -\xi \right ](t-nT)}\Phi _{\tilde {\varepsilon }}-\alpha _2(\tilde {\varepsilon }+ap^*)\Psi _{\tilde {\varepsilon }} \right ] }{(\tilde {\varepsilon }+ap^*)\left [a(p^*+\delta \Psi _{\tilde {\varepsilon }})+b\delta \frac {\rho _1}{g'(0)}e^{\left [-(\lambda _0^2)_{\tilde {\varepsilon }} -\xi \right ](t-nT)}\Phi _{\tilde {\varepsilon }}\right ]} \right \}. \end {aligned} \end{equation*}

If taking $\delta =0$ in the last term of the equation above, we have $\frac{-\alpha _2\tilde{\varepsilon }}{a(\tilde{\varepsilon }+ap^*)}\lt 0$ . One thus could choose $\delta _1$ small sufficiently such that the last term of the equation above is negative uniformly.

Thus, this yields that

\begin{equation*} \begin {aligned} &\frac {\partial \widehat {h}}{\partial t}-\left [ \frac {d_2}{\rho ^2(t)}\widehat {h}_{yy}-\frac {\dot \rho (t)}{\rho (t)}\widehat {h}+\widehat {h}\left (-\gamma _2+\frac {\alpha _2\widehat {p}}{a\widehat {p}+b\widehat {h}}\right )\right ]\lt 0. \end {aligned} \end{equation*}

Besides, if $\delta \lt \delta _2\;:\!=\;\left ( \frac{g'(0)-\rho _1}{D}\right )^{\frac{1}{\nu -1}}$ , it follows from the assumption $\boldsymbol{\mathrm{A}_3}$ that

\begin{equation*} \begin {aligned} g(\widehat {h}(nT,y))-\widehat {h}((nT)^+,y)&=g(\widehat {h}(nT,y))-\delta \frac {\rho _1}{g'(0)}\Phi _{\tilde {\varepsilon }}((nT)^+,y)\\[5pt] &=g(\widehat {h}(nT,y))-\rho _1\widehat {h}(nT,y)\\[5pt] &\geq (g'(0)-\rho _1)\widehat {h}(nT,y)-D(\delta \Phi _{\tilde {\varepsilon }}(nT,y))^\nu \\[5pt] &=[(g'(0)-\rho _1)-D(\delta \Phi _{\tilde {\varepsilon }}(nT,y))^{\nu -1}]\delta \Phi _{\tilde {\varepsilon }}(nT,y)\\[5pt] &\geq 0. \end {aligned} \end{equation*}

Similarly, we can also verify that

\begin{equation*} \begin {aligned} &\frac {\partial \widehat {p}}{\partial t}-\left [ \frac {d_1}{\rho ^2(t)}\widehat {p}_{yy}-\frac {\dot \rho (t)}{\rho (t)}\widehat {p}+\widehat {p}\left (\gamma _1-c\widehat {p}+\frac {\alpha _1\widehat {h}}{a\widehat {p}+b\widehat {h}}\right )\right ]\lt 0. \end {aligned} \end{equation*}

Consequently, we infer that $(\widehat{p},\widehat{h})=(p^*+\delta \Psi _{\tilde{\varepsilon }},\delta \Phi _{\tilde{\varepsilon }})$ is the lower solution of problem (2.17).

Next, we select the $(\overline{p}^{(0)},\overline{h}^{(0)})=(\widetilde{p},\widetilde{h})$ and $(\underline{p}^{(0)},\underline{h}^{(0)})=(\widehat{p},\widehat{h})$ as initial iteration, the sequences $\big \{(\overline p^{(m)},\overline h^{(m)})\big \}$ and $\big \{(\underline p^{(m)},\underline h^{(m)})\big \}$ are defined by (2.19). It follows from Lemma 2.12 that we have

\begin{equation*} (\widehat {p},\widehat {h})\leq (\underline {p}^{(m)},\underline {h}^{(m)})\leq (\underline {p}^{(m+1)},\underline {h}^{(m+1)})\leq (\overline {p}^{(m+1)},\overline {h}^{(m+1)})\leq (\overline {p}^{(m)},\overline {h}^{(m)})\leq (\widetilde {p},\widetilde {h}). \end{equation*}

Based on the monotone convergence theorem, we obtain that the limits of the sequences $\big \{(\overline p^{(m)},\overline h^{(m)})\big \}$ and $\big \{(\underline p^{(m)},\underline h^{(m)})\big \}$ exist and

\begin{equation*} \lim _{m\rightarrow \infty }(\overline {p}^{(m)},\overline {h}^{(m)})=(\overline {p},\overline {h}), \lim _{m\rightarrow \infty }(\underline {p}^{(m)},\underline {h}^{(m)})=(\underline {p},\underline {h}), \end{equation*}

where $(\overline{p},\overline{h})$ and $(\underline{p},\underline{h})$ are $T$ -periodic solutions of problem (2.17) satisfying $(p^*,0)\leq (\underline{p},\underline{h})\leq (\overline{p},\overline{h})$ . Moreover,

\begin{equation*} \begin {aligned} (\widehat {p},\widehat {h})&\leq (\underline {p}^{(m)},\underline {h}^{(m)})\leq (\underline {p}^{(m+1)},\underline {h}^{(m+1)})\leq (\underline {p},\underline {h})\\[5pt] &\leq (\overline {p},\overline {h})\leq (\overline {p}^{(m+1)},\overline {h}^{(m+1)})\leq (\overline {p}^{(m)},\overline {h}^{(m)})\leq (\widetilde {p},\widetilde {h}). \end {aligned} \end{equation*}

Now we claim that $(\overline{p},\overline{h})$ and $(\underline{p},\underline{h})$ are the maximal and minimal positive $T$ -periodic solutions of problem (2.17). In fact, for any positive periodic solution $(p^{**},h^{**})$ of problem (2.17) over $(p^*,0)$ satisfies $(\widehat{p},\widehat{h})\leq (p^{**},h^{**})\leq (\widetilde{p},\widetilde{h})$ . Employing the same iteration as problem (2.19), we choose $(\widetilde{p},\widetilde{h})$ and $(p^{**},h^{**})$ as the initial iteration with $(\overline{p}^{(0)},\overline{h}^{(0)})=(\widetilde{p},\widetilde{h})$ and $(\underline{p}^{(0)},\underline{h}^{(0)})=(p^{**},h^{**})$ , it follows that

\begin{equation*} (p^{**},h^{**})\leq (\overline {p},\overline {h}),\quad t\geq 0, y\in [0,l_0], \end{equation*}

thus, $(\overline{p},\overline{h})$ is the maximal positive $T$ -periodic solution of problem (2.17). Similarly, $(\underline{p},\underline{h})$ is the minimal positive $T$ -periodic solution of problem (2.17).

Step 2 we now present the uniqueness of the positive $T$ -periodic solution of problem (2.17).

Indeed, since the two components of the system are weakly coupled, we can refer to [39] to prove the uniqueness.

We remark that $(\overline{p},\overline{h})=(\underline{p},\underline{h})\;:\!=\;(p^{**},h^{**})$ provided with $(\underline{p},\underline{h})(0,y)=(\overline{p},\overline{h})(0,y)$ . In fact, choosing the initial condition $p(y,0)=p_0(y)$ , one can regard problem (2.17) as an initial boundary value problem and then acquire its uniqueness condition through the standard existence-uniqueness theorem on the initial boundary value parabolic problem. On the other hand, assume that $h_1$ and $h_2$ are the two solutions and define set

\begin{equation*} \eta =\left \{\varsigma \in [0,1],\varsigma h_1\leq h_2, t=0, t=0^+, t\in (0^+,T], y\in [0,l_0]\right \}, \end{equation*}

which can be shown that $\eta$ possesses a right neighbourhood around 0. We say that $1\in \eta$ . Suppose not, then we have that $\varsigma _0=\sup \eta \lt 1$ . We note that $F(p,h,t)=f(p,h,t)+k_2h$ is nondecreasing and $\frac{f(p,h,t)}{h}$ is decreasing in $h$ on $[0,\max \limits _{[0,l_0]\times [0,T]} h_2]$ , it yields that

\begin{equation*} \begin {aligned} (h_2-\varsigma _0h_1)_t-\frac {d_2}{\rho ^2(t)}(h_2-\varsigma _0h_1)_{yy}+k_2(h_2-\varsigma _0h_1) &\geq \left (-\frac {\dot \rho (t)}{\rho (t)}-\gamma _2\right )(h_2-\varsigma _0h_1)+k_2(h_2-\varsigma _0h_1)\\[5pt] &=f_2(h_2,t)+k_2h_2-\varsigma _0(f_2(h_1,t)+k_2h_1)\\[5pt] &\geq f_2(\varsigma _0h_1,t)+k_2\varsigma _0h_1-\varsigma _0(f_2(h_1,t)+k_2h_1)\\[5pt] &\geq 0 \end {aligned} \end{equation*}

for $t\in (0^+,T]$ and $y\in (0,l_0)$ . By assumptions $\boldsymbol{\mathrm{A}_1}$ and $\boldsymbol{\mathrm{A}_2}$ , we deduce that

\begin{equation*} \begin {aligned} h_2(0^+,y)-\varsigma _0h_1(0^+,y)&=g(h_2(0,y))-\varsigma _0g(h_1(0,y))\\[5pt] &\geq g(\varsigma _0h_1(0,y))-\varsigma _0g(h_1(0,y))\geq 0 \end {aligned} \end{equation*}

for $y\in (0,l_0)$ . However, for $t\gt 0$ ,

\begin{equation*} \begin {aligned} h_2(t,0)-\varsigma _0h_1(t,0)=h_2(t,l_0)-\varsigma _0h_1(t,l_0)=0. \end {aligned} \end{equation*}

Due to the strong maximum principle [33], we have significant statements as follows:

1. (a) $h_2-\varsigma _0h_1\gt 0$ holds for $t=0^+, t\in (0^+,T]$ and $y\in (0,l_0)$ . Since $h_1$ and $h_2$ are $T$ -periodic solutions, that is, $h_1(0,y)=h_1(T,y)$ and $h_2(0,y)=h_2(T,y)$ for $y\in (0,l_0)$ , and utilising the strong maximum principle implies $h_2-\varsigma _0h_1\gt 0$ for $t\in (0,T]$ and $y\in (0,l_0)$ . Based on the Hopf’s boundary lemma, we deduce that $\frac{\partial }{\partial \textbf{n}}\mid _{y=0}(h_2-\varsigma _0h_1)\gt 0$ and $\frac{\partial }{\partial \textbf{n}}\mid _{y=l_0}(h_2-\varsigma _0h_1)\lt 0$ , where $\textbf{n}$ is the outward unit normal vector. Then, there is a constant $\epsilon \gt 0$ such that $h_2-\varsigma _0h_1\gt 0\geq \epsilon h_1$ , which leads to $\varsigma _0+\epsilon \in \eta$ . This contradicts the maximality of $\varsigma _0$ .

2. (b) $h_2-\varsigma _0h_1\equiv 0$ for $t=0^+, t\in (0^+,T]$ and $y\in (0,l_0)$ . In this case, we have $f_2(h_2,t)=\varsigma _0f_2(h_1,t)$ . However, recalling $\varsigma _0\lt 1$ , $f_2(h_2,t)=f_2(\varsigma _0h_1,t)\gt \varsigma _0f_2(h_1,t)$ ; thus, it is also impossible.

To sum up, problem (2.17) admits a unique positive $T$ -periodic solution $(p^{**},h^{**})$ .

(ii) Due to the Hopf’s boundary lemma, we obtain that $\phi _y(0,0)\gt 0$ and $\phi _y(0,l_0)\lt 0$ , and we select a small enough $\delta$ to make sure $\delta \phi (0,y)\leq h(0,y)$ . Meanwhile, a large enough $M_2$ can be chosen such that $h(0,y)\leq M_2U(0)$ . For given $\delta$ , $M_1$ and $M_2$ , the function $(\widetilde{p},\widetilde{h})$ = $(M_1U(t),M_2V(t))$ with $(U(t),V(t))$ defined in (2.22) and $(\widehat{p},\widehat{h})=(p^*+\delta \psi,\delta \phi )$ with $\widehat{h}$ defined in (2.23), satisfies

\begin{equation*} {(0,\widehat {h})(0,y)\leq (p,h)(0,y)\leq (\widetilde {p},\widetilde {h})(0,y),\quad y\in [0,l_0].} \end{equation*}

It follows from ( $\boldsymbol{\mathrm{A}_2}$ ) that $g$ is nondecreasing with respect to $h$ , we obtain that

\begin{equation*} \widehat {h}(0^+,y)\leq g(\widehat {h}(0,y))\leq g(h(0,y))=h(0^+,y)\leq g(\widetilde {h}(0,y))=\widetilde {h}(0^+,y). \end{equation*}

The classical comparison principal yields $\widehat{h}(t,y)\leq h(t,y)\leq \widetilde{h}(t,y)$ , $t\in (0^+,T]$ , $y\in [0,l_0]$ . Induction reveals that $\widehat{h}(t,y)\leq h(t,y)\leq \widetilde{h}(t,y)$ , $t=nT, (nT)^+, t\in ((nT)^+,(n+1)T]$ , $y\in [0,l_0]$ .

Therefore,

\begin{equation*} \underline {h}^{(0)}(t,y)\leq h(t,y)\leq \overline {h}^{(0)}(t,y), \,\,t=nT, (nT)^+, t\in ((nT)^+,(n+1)T], y\in [0,l_0]. \end{equation*}

Moreover,

(2.25) \begin{equation} \underline{h}^{(0)}(T,y)\leq h(T,y)\leq \overline{h}^{(0)}(T,y), \,\, y\in [0,l_0], \end{equation}

which together with $\underline{h}^{(1)}(0,y)=\underline{h}^{(0)}(T,y)$ and $\overline{h}^{(1)}(0,y)=\overline{h}^{(0)}(T,y)$ yields

\begin{equation*} \underline {h}^{(1)}(0,y)\leq h(T,y)\leq \overline {h}^{(1)}(0,y), \,\, y\in [0,l_0]. \end{equation*}

By the assumption ( $\boldsymbol{\mathrm{A}_2}$ ) and (2.25), we obtain that

\begin{equation*} g(\underline {h}^{(0)}(T,y))\leq g(h(T,y))\leq g(\overline {h}^{(0)}(T,y)), \,\, y\in [0,l_0]. \end{equation*}

Recalling problems (2.19) and (1.6), ones deduce that

\begin{equation*} \begin {aligned} \underline {h}^{(1)}(0^+,y)&=g(\underline {h}^{(0)}(T,y))\leq g(h(T,y))\\[5pt] &=h(T^+,y)\leq g(\overline {h}^{(0)}(T,y))=\overline {h}^{(1)}(0^+,y), \,\, y\in [0,l_0], \end {aligned} \end{equation*}

that is,

\begin{equation*} \underline {h}^{(1)}(0^+,y)\leq h(T^+,y)\leq \overline {h}^{(1)}(0^+,y), \,\, y\in [0,l_0]. \end{equation*}

The comparison principle shows that $\underline{h}^{(1)}(t,y)\leq h(t+T,y)\leq \overline{h}^{(1)}(t,y)$ holds for $t\in (0^+,T]$ and $y\in [0,l_0]$ . Utilising induction again, we have

\begin{equation*} \underline {h}^{(1)}(t,y)\leq h(t+T,y)\leq \overline {h}^{(1)}(t,y),\,\,t=nT, (nT)^+, t\in ((nT)^+,(n+1)T], y\in [0,l_0], \end{equation*}

and combining with the last two equations in problem (2.19) that

\begin{equation*} \underline {h}^{(m)}(t,y)\leq h(t+mT,y)\leq \overline {h}^{(m)}(t,y),\,\,t\geq 0, y\in [0,l_0], \end{equation*}

since the above inequality holds for $m=0$ and $m=1$ . Similarly, we do some same work for $p$ .

Recalling the uniqueness of the periodic solution of problem (2.17) provided with $\lim \limits _{m\rightarrow \infty }(\underline{p}^{(m)},\underline{h}^{(m)})(t,y)=\lim \limits _{m\rightarrow \infty }(\overline{p}^{(m)},\overline{h}^{(m)})(t,y)=(p^{**},h^{**})$ in (i), we have

\begin{equation*} \lim _{m\rightarrow \infty }(p,h)(t+mT,y;p_0,h_0)\rightarrow (p^{**},h^{**})(t,y). \end{equation*}

Remark 2.14. Considering the cooperative characteristics of problem (2.17), one can still use the subhomogeneity of the nonlinear term to obtain the uniqueness of the positive periodic equilibrium solution of problem (2.17). However, in the current paper, we prefer to apply more detailed differential equation analysis techniques based on the maximum principle in order to reveal more details of problem (2.17).

3. the impacts of the evolution rate and impulsive effect

In this section, we will perform some numerical simulations to verify the theoretical results obtained in the previous section. We aim to investigate how the domain evolution and impulse affect these dynamical behaviours in plants-pollinators world and set $R_0^2(\rho )$ to emphasise such dependence. In all simulations, we always consider the interval $[0,l(t)]=[0,\rho (t)l_0]$ , where $l_0=\pi$ , and we first fix

\begin{equation*} d_1=0.3, d_2=0.005, \alpha _1=0.1, \alpha _2=0.01, \gamma _1=0.1, \gamma _2=0.09, a=0.047, b=0.01, c=0.01 \end{equation*}

and then provide $\lambda ^*=\left (\frac{\pi }{l_0}\right )^2=1$ . We select

\begin{equation*}p_0(y)=5\sin (y), \quad h_0(y)=0.5\sin (y)+0.2\sin (3y)\end{equation*}

as initial function.

In fact, we must note that, if $\rho (t)\equiv 1$ and $l(t)=l(0)$ , then we can easily obtain the corresponding problem on fixed domain for problem (1.6). The ecological reproduction index is represented as $\Re _0^2=R_0^2(1)=\frac{\frac{\alpha _2}{a}}{\gamma _2+d_2\lambda ^*-\frac{1}{T}\ln g'(0)}$ by the identical arguments with Lemma (2.4), and we have the threshold scenario completely similar to Theorems (2.10) and (2.13). Here we omit it and refer to [43] for more details.

The following assertions can help to explain more graphically the influence of the periodically evolving domain on the persistence of pollinators, which is obtained directly from Theorem (2.5). Specifically, if the evolution rate of domain is small, then we have $R_0^2\lt \Re _0^2$ , which means the evolving domain is not conducive to pollinators survival, that is, the evolving domain has a negative influence on the persistence of pollinators. While the evolution rate of domain is large, then we have $R_0^2\gt \Re _0^2$ , which implies the evolving domain can promote diffusion of pollinators such that pollinators have more space for transmission. If $\rho (t)=1$ , then we have $R_0^2=\Re _0^2$ , which gives that pollinators can persist on the evolving domain at the same scope of diffusive rate on the fixed domain, that is, the evolving domain has no influence on the persistence of pollinators.

3.1. The impact of the evolution rate

We select different $\rho (t)$ to emphasise the impact of the domain evolution rate on the dynamical behaviours of pollinators when the impulsive effect occurs every time $T=2$ . We first fix $n_1=8$ and $n_2=10$ in (1.2).

Example 3.1. Take $g(H)=\frac{8H}{10+H}$ and then $g'(0)=0.8$ . Let $\rho _1(t)=e^{-0.1(1-\cos (\pi t))}$ . Then from (2.9), we calculate that

\begin{equation*} R_0^2(\rho _1)=\frac {\frac {\alpha _2}{a}}{\gamma _2+\frac {d_2\lambda ^*}{T}\int _0^T \frac {1}{\rho ^2(t)}\mathrm {d}t-\frac {1}{T}\ln g'(0)}\approx 0.9627\lt 1. \end{equation*}

In Figure 2 , there is a clear trend of $h$ decreases $0$ , which is consistent with Theorem 2.10 that pollinators suffer extinction eventually.

Next, let $\rho _2(t)=e^{0.1(1-\cos (\pi t))}$ in (2.9). Then we have

\begin{equation*} R_0^2(\rho _2)=\frac {\frac {\alpha _2}{a}}{\gamma _2+\frac {d_2\lambda ^*}{T}\int _0^T \frac {1}{\rho ^2(t)}\mathrm {d}t-\frac {1}{T}\ln g'(0)}\approx 1.0343\gt 1. \end{equation*}

It is shown in Figure 3 that $h$ approaches a positive periodic steady state. And it agrees with Theorem 2.13 that pollinators can coexist with plants.

The results, as shown in Example 3.1, indicate that the evolution rate of the domain is crucial to the persistence and extinction of the population. Specifically, the larger the evolution rate is conducive for pollinators‘ survival, that is, a large domain evolution rate has a positive influence on the survival of pollinators when impulsive effect takes place. Nevertheless, that pollinators eventually extinct in a periodically evolving domain with a small evolution rate.

Figure 2. $\rho _1(t)=e^{-0.1(1-\cos (\pi t))}$ , $n_1=8$ and $n_2=10$ . The domain is periodically evolving with $\rho _1$ and $R_0^2\lt 1$ . Graphs $(a)$ - $(c)$ show that the population $H(t,x)$ decays to 0. Graphs $(b)$ and $(c)$ are the cross-sectional view and projection of $H$ on the $t-H-$ plane, respectively. The colour bar in graph $(b)$ shows the density of $H(t,x)$ .

Figure 3. $\rho _2(t)=e^{0.1(1-\cos (\pi t))}$ , $n_1=8$ , $n_2=10$ , and $R_0^2\gt 1$ . Graph $(a)$ shows the dynamics of pollinators $H(t,x)$ , which implies that pollinators tend to a positive periodic steady state, it also shows pollinators can coexist with plants on the periodically evolving domain. Graph $(b)$ is the cross-sectional view and indicates the periodic evolution of the domain. The appearance of impulsive effect every time $T=2$ can be seen in graph $(c)$ , which is the projection of $H$ on the $t-H-$ plane.

3.2. The impact of impulsive effect

In order to study how pulse affects the dynamical behaviours of pollinators in a periodically evolving domain, we first employ numerical simulations to illustrate the case when impulsive effect does not occur.

Example 3.2. Fix $\rho _1(t)=e^{-0.1(1-\cos (\pi t))}$ . From (2.12), it yields that $R_0^2\lt 1$ and without impulsive effect. Observing Figure 4 , we find that pollinators suffer extinction eventually. Comparing Figure 2 , we obtain that pollinators suffer extinction at a faster speed when impulsive effect occurs.

Next, we fix $\rho _2(t)=e^{0.1(1-\cos (\pi t))}$ , and assume that impulsive effect does not occur. By (2.12), we have $R_0^2\gt 1$ . It is easily seen from Figure 5 that pollinators stabilise to a positive periodic steady state.

Finally, we select the impulsive function with $n_1=5$ and $n_2=10$ , that is, $g(H)=5H/(10+H)$ . One can see from Figure 6 that pollinators now decay to extinction. Figures 5 and 6 show that pollinators survive in an evolving domain with a large evolution rate, but vanishes when the impulsive effect takes place.

Figure 4. $\rho _1(t)=e^{-0.1(1-\cos (\pi t))}$ and without pulse. In this case, $R_0^2\lt 1$ . Graphs $(a)$ - $(c)$ imply that pollinators $H(t,x)$ decreases to 0. Graphs $(b)$ and $(c)$ are the cross-sectional view and projection of $H$ on the plane $t-H$ , respectively.

Figure 5. $\rho _2(t)=e^{0.1(1-\cos (\pi t))}$ and without pulse. In this situation, $R_0^2\gt 1$ . Graphs $(a)$ - $(c)$ show that population $H(t,x)$ approaches to a positive periodic steady state.

Figure 6. In the case with impulsive effect, $\rho _2(t)=e^{0.1(1-\cos (\pi t))}$ and $g(H)$ is chosen with $n_1=5$ and $n_2=10$ , this implies $R_0^2\lt 1$ . Graph $(a)$ suggests that the population will go to extinction eventually. Graph $(b)$ is the case where the domain is periodically evolving. We can also clearly observe the impact of impulsive effect every time $T=2$ from graph $(c)$ , in which the population suffers extinction.

Example 3.2 reveals that when pollinators live in a periodically evolving domain with a small evolution rate, impulsive effect can speed up the extinction of pollinators. Taken together, impulsive effect has a negative influence on the survival of pollinators and, eventually, leads to the extinction of pollinators.

4. Discussion

In this paper, we have combined the periodic evolution of domain and impulsive effect into the plants-pollinators cooperative system, which makes it more reasonable for describing the persistence and extinction of species. Based on the interdependence of pollinators and plants, through the current paper, we focus on discussing the case when plants are survival given that $R_0^1\gt 1$ . The main purpose is to examine the threshold dynamics scenario of pollinators under the influence of plants.

Firstly, we define the ecological reproduction index $R_0^2$ of pulse problem and provide an explicit formula. Then, utilising the monotone iteration technique with the proper upper and lower solutions, we establish dynamical behaviours of the solution to problem (1.6) when the impulsive function is monotone. We conclude that in the case of $R_0^2\leq 1$ , the solution $(p,h)(t,y)=(p^*, 0)$ , which sees details from Theorem 2.10. On the contrary, when $R_0^2\gt 1$ , Theorem 2.13 implies that the solution converges to a positive periodic steady state, indicating that pollinators can coexist with plants. In addition, our numerical simulations further illustrate that pollinators suffer extinction with a small evolution rate of domain (see Figure 2), but survive in one with a large evolution rate of domain (see Figure 3). Meanwhile, a large evolution rate of domain is beneficial to the survival of the pollinators. Another notable result is that the impulsive effect can speed up the pollinators‘ extinction (see Figures 2 and 4) and has a negative influence on the pollinators’ survival (see Figures 5 and 6), and, eventually, leads to the extinction of pollinators (see Figure 6). Which is consistent with our theoretical results.

We still note that the current analysis is based on the one-dimensional case for plants-pollinators cooperative system. Fazly et al. [14] have recently expanded the aforementioned findings of seasonal birth pulse in [22] to a higher dimensional for logistic model without domain evolution. They provided domain parameters and discussed species extinction and persistence as a function of domain shape and size. However, it remains to be further investigated whether our approach can be applied to higher dimensional evolving domains with impulsive effects in plant-pollinator systems. In addition, it is also worth considering that the combination of free boundary and impulsive effect in plant-pollinator systems.