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

HAL
open science

Global stability in the Ricker model with delay and
stocking

Ziyad Alsharawi, Sadok Kallel

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

To cite this version:

Ziyad Alsharawi, Sadok Kallel. Global stability in the Ricker model with delay and stocking. 2023. [hal-04435644]

HAL Id: hal-04435644

https://hal.science/hal-04435644v1

Preprint submitted on 2 Feb 2024

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

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

Public Domain
# Global stability in the Ricker model with delay and stocking

Ziyad AlSharawi* and Sadok Kallel

Department of Mathematics and Statistics
American University of Sharjah, P. O. Box 26666
University City, Sharjah, UAE

January 31, 2024

## Abstract

We consider the Ricker model with delay and constant or periodic stocking. We found that the high stocking density tends to neutralize the delay effect on stability. Conditions are established on the parameters in order to ensure the global stability of the equilibrium solution in the case of constant stocking, as well as the global stability of the 2-periodic solution in the case of 2-periodic stocking. Our approach extensively relies on the utilization of the embedding technique. Whether constant stocking or periodic stocking, the mode has the potential to undergo a Neimark-Sacker bifurcation in both cases. However, the Neimark-Sacker bifurcation in the 2-periodic case results in the emergence of two invariant curves that collectively function as a single attractor. Finally, we pose open questions in the form of conjectures about global stability for certain choices of the parameters.

**AMS Subject Classification:** 39A10, 39A30, 92D25

**Keywords:** Ricker model, global stability, embedding, periodic solutions, stocking.

## 1 Introduction

As an early attempt to model a population of a single species with non-overlapping generations, Moran (1950) [1] proposed using a one-hump map that increases to a maximum and then decreases asymptotically to a non-negative value. Along this line of thought and based on experimental data, Ricker introduced a density-dependent model that became known in the scientific literature as the Ricker model [2]

$$x_{n+1} = x_n f(x_n) = x_n e^{r(1 - \frac{x_n}{K})}, \quad K, r > 0 \text{ and } x_0 \in \mathbb{R}_+ = [0, +\infty). \quad (1.1)$$

In this model, the parameters K and r represent the carrying capacity and the intrinsic growth rate, respectively. To learn more about the origin and significance of Ricker's model, and Ricker's pivotal contribution in the field of quantitative fishery science, we refer the reader to [3]. Rescaling can reduce the model to $y_{n+1} = y_n \exp(r - y_n)$. This model enjoyed a boost of success after the publication of the two studies by May [4, 5], and a third study by May and Oster [6], in which the model's intriguing dynamics were investigated, and the global stability of its positive equilibrium

<footnote>*Corresponding author: zsharawi@aus.edu</footnote>
for $0 < r < 2$ was observed. A rigorous proof of the global stability for $0 < r \le 2$ was developed later using various techniques such as Lyapunov functions [7], Singer's theorem [8] and enveloping [9].

While implicit time lags are embedded in discrete systems, explicit time lags caused by significant recruitment delays must be accounted for in the density-dependent function $f$. Time delays are necessary to accommodate the recruitment to the adult stage, which may vary based on the type of species [10–12]. This mechanism requires that a mathematical model be considered with a certain time lag $N$. As a result, Eq. (1.1) becomes

$$y_{n+1} = y_n f(y_{n-N}) = y_n e^{r-y_{n-N}}. \quad (1.2)$$

Due to implicit delays, solutions of discrete models are observed to overshoot and undershoot their equilibrium level and generate oscillations, periodic solutions, or chaos. This tendency is typically aggravated when explicit delays are incorporated into the system [10, 11, 13]. When a population exhibits convergence to its equilibrium regardless of its initial density, the phenomenon is called global stability. Proving mathematically that the equilibrium solution is globally stable becomes a challenging task when time delays are involved. Under general settings of the density-dependent map $f$, Liz et al. [14] found sufficient conditions to obtain global stability when $0 < r < \frac{3}{2(N+1)}$. This implies $0 < r < \frac{3}{4}$ when $N = 1$. In a lengthy and technical paper that included computer-assisted proofs, Bartha et al. [15] established global stability when $N = 1$ and $0 < r < 1$.

Stocking or harvesting can be used as a management tool to achieve various goals, including chaos reversal and re-stabilization of an equilibrium level [16, 17]. According to [18], constant rate harvesting is observed to transform contest competition into scramble competition. On the other hand, constant stocking slows the fluctuation and has a stabilizing effect on equilibrium solutions [19]. When constant stocking is applied to the Ricker model with delay, Eq. (1.2) becomes

$$y_{n+1} = y_n f(y_{n-1}) + h = y_n e^{r-y_{n-1}} + h = F(y_n, y_{n-1}), \quad h > 0. \quad (1.3)$$

It is also worthwhile to consider this equation in terms of perturbation, as was done in [20] in the absence of time delay. However, Eq. (1.3) becomes intriguing because it includes two significant factors, namely delay and stocking. The individual impact of delay and stocking on stability is known to reflect opposing effect, which motivates investigating the interaction between these factors and their influence on overall stability. Equations of the form $y_{n+1} = y_n f(y_{n-1}) + h$ have been considered in [18] under the assumption that $tf(t)$ is increasing; however, Eq. (1.3) does not belong to this category. If the stocking is done through seasonal quotas, the constant $h$ is replaced by a sequence $\{h_n\}$, which we assume to be $p$-periodic. Therefore, Eq. (1.3) becomes

$$y_{n+1} = y_n f(y_{n-1}) + h_n = y_n e^{r-y_{n-1}} + h_n = F_n(y_n, y_{n-1}), \quad (1.4)$$

where $h_n \ge 0$ is a $p$-periodic sequence and $F$ is defined on the positive orthant $\mathbb{R}_+^2$. We continue to write $f(y)$ instead of $e^{r-y}$ whenever we find it convenient. Unlike Eq. (1.1), which has two equilibrium points $x = 0$ and $x = r$, Eq. (1.3) has a unique non-negative equilibrium, which bifurcates into a $p$-periodic solution when periodicity is introduced in Eq. (1.4).

Our primary objective is to investigate the global stability of the equilibrium solution in Eq.
(1.3) and the p-periodic solution in Eq. (1.4). Our method, when it applies, provides a clear geometric understanding for the reason why stability (local or global) holds. This paper is structured as follows: In Section Two, preliminary findings pertaining to the utilization of the embedding technique in both autonomous and periodic cases are presented. In the third section, we leverage the outcomes from Section Two to apply them to Equation (1.3). This section delineates sufficient conditions on the variable h that ensure global stability in the autonomous case. We determine “almost completely” the regions of local and global stability in the (h, r) plane, modulo a well-identified restricted region of uncertainty. In Section Four, we turn our attention to Equation (1.4) and employ the embedding technique to investigate the 2-periodic case, ultimately establishing the global stability of the 2-periodic solution under certain sufficient conditions. This paper closes with a conclusion section that not only summarizes our key findings but also raises pertinent open questions for further exploration.

## 2 The embedding technique

Embedding a dynamical system in general into a higher dimensional dynamical system that can be utilized to classify certain characteristics of the original system is a known approach [21–23]. However, the challenge arises when it comes to the technical details of the approach, particularly with the partial order in the embedded system. In this section, we build the basic machinery needed in the sequel. Let V denote a partially ordered metric space, specifically the positive orthant, as defined for the purposes of this discussion, $\mathbb{R}_+^n$ or $[a, b]^n$ for some $n$. A recursive sequence in V, with delay k, is any sequence defined by

$$ \alpha_F : x_{n+1} = F(x_n, \dots, x_{n-k+1}), \quad k \ge -1, \qquad (2.1) $$

where $F: V^k \to V$ is a continuous function and the initial terms $x_0, x_{-1}, \dots, x_{-k+1}$ are given in V. These initial terms, together with F, determine the sequence uniquely. We write $S(V)$, the set of all recursive sequences in V. This can be topologized as a subspace of $V^w$, the infinite direct product of V. Notice that different functions F can give rise to the same recursive sequence (2.1), so only the sequence uniquely determines the system.

**Definition 2.1.** A continuous injection $\Psi: S(V_1) \hookrightarrow S(V_2)$, which sends convergent sequences to convergent sequences, is called an “embedding”. If $V_1$ is a subspace of $V_2$, then there is a canonical inclusion $S(V_1) \subset S(V_2)$.

It is pertinent to note that, generally, the convergence of $\Psi(\alpha)$ does not entail the convergence of $\alpha$ for a given embedding $\Psi: S(V_1) \to S(V_2)$.

### 2.1 Embedding in the autonomous case

Consider the two-dimensional difference equation $\alpha_F: x_{n+1} = F(x_n, x_{n-1})$, where $F: \mathbb{R}_+^2 \to \mathbb{R}$ is any continuous function. We define $T(x, y) = (F(x, y), x)$, so that $T(x_n, x_{n-1}) = (x_{n+1}, x_n)$ and its iterations produce a sequence in the plane

$$ (X_n) : (x_0, x_1), (x_1, x_0), \dots, X_n := (x_n, x_{n-1}), $$

which entirely describes the dynamics of the system. We refer to T as the “vector form” of the system $\alpha_F$. The system T has a global attractor means of course that the sequence $(X_n)$ converges, independently of the choice of the initial value within the given domain. The assignment
$(x_n) \mapsto (X_n)$ gives an embedding $\Psi : S(\mathbb{R}) \hookrightarrow S(\mathbb{R}^2)$. This is a strong embedding in the sense that $\Psi(\alpha)$ is convergent if and only if $\alpha$ is convergent.

Next, we write $V = \mathbb{R}_+^2$. For every sequence $\alpha_F \in S(V)$, we define a 4-dimensional sequence $\zeta_F \in S(V^4)$ with general term $\zeta_n = (X_n, X_n) = (x_n, x_{n-1}, x_n, x_{n-1})$, and initial term $\zeta_0 = (x_0, x_{-1}, x_0, x_{-1})$. This is the image of $(X_n)$ under the diagonal embedding $V^2 \to V^2 \times V^2$, $(x, y) \mapsto (x, y, x, y)$. This is a strong embedding as well $S(V) \hookrightarrow S(V^4)$, which we refer to as the “diagonal embedding”. Because this is a strong embedding, a good way to establish convergence of $\zeta_F$ will help us obtain convergence in our initial system.

Let us write $\zeta_F$ in recursive form as follows. Define the self-map of $V^4$ as

$$G(x, y, u, v) = (F(x, y), u, F(u, v), x). \qquad (2.2)$$

Starting with $\zeta_0$, we see that $\zeta_n = G^n(\zeta_0)$. The advantage of introducing the map $G$ is to obtain a form of monotonicity [22, 23]. Introduce the “southeast partial ordering” on $V \times V$ by $(x_1, y_1) \leq_{se} (x_2, y_2)$ if and only if $x_1 \leq x_2$ and $y_1 \geq y_2$. Using the southeast partial ordering on $V^2 \times V^2$ means using $\leq_{se}$ on $W \times W$ with $W = V^2$. A key observation now is that if $F$ is non-decreasing in its first component and non-increasing in its second one, i.e., $F(\uparrow, \downarrow)$, then the sequence $G$ is monotonic with respect to the $\leq_{se}$ partial ordering on $V^2 \times V^2$. More precisely, rewrite $G : V^4 \to V^4$ as a map $G : V^2 \times V^2 \to V^2 \times V^2$, $(X, U) \mapsto G(X, U)$. Then

$$(X_1, U_1) \leq_{se} (X_2, U_2) \implies G(X_1, U_1) \leq_{se} G(X_2, U_2).$$

This monotonicity gives a handy way of studying the convergence of the diagonal sequence $\zeta_F : \zeta_n = G^n(\zeta_0)$, and thus ultimately that of $\alpha_F$. We develop this approach next.

**Definition 2.2.** Write $A = (a, b)$ and $B = (b, a)$. We say that $X = (x, y)$ is inside a box with vertices $A$ and $B$ in $V^2$ if $a \leq b$ and $A \leq_{se} X \leq_{se} B$. Equivalently, if in $V^4$ we have

$$(A, B) \leq_{se} (X, X) \leq_{se} (B, A)$$

Suppose that for a choice of the map $F$ and the points $A$ and $B$, the map $G$ in Eq. (2.2) has the property that $(A, B) \leq_{se} G(A, B)$ and $G(B, A) \leq_{se} (B, A)$. This happens if $(a, b) \leq_{se} (F(a, b), F(b, a))$. Then, if $X_0$ (resp. some $X_k = T^k(X_0)$ for some $k \geq 0$) is inside a box with vertices $A, B$, we can see immediately that for positive $n$,

$$(A, B) \leq_{se} G^n(A, B) \leq_{se} G^n(X_0, X_0) \leq_{se} G^n(B, A) \leq_{se} (B, A), \qquad (2.3)$$

which means that our sequence $\zeta_F$ is eventually caught between the orbits of $G$ through $(A, B)$ and $(B, A)$. The sequences $G^n(A, B)$ and $G^n(B, A)$, being bounded and monotonic in the southeast ordering, converge to fixed points of $G$. A fixed point of $G$ must be of the form

$$\bar{x} = (x, y, y, x), \text{ where } x = F(x, y) \text{ and } y = F(y, x).$$

Such fixed points come in pairs if $x \neq y$, and so if $(F(x, y), F(y, x)) = (x, y)$ has a unique solution, then $x = y = \bar{x}$ is a fixed point of $F$. If $x \neq y$, then $(x, y)$ and $(y, x)$ are dubbed as pseudo fixed points of $F$ (these were called “artificial fixed points” in [24]). The above discussion leads us to the following result [22, 23].
**Proposition 2.3.** Let $F: V^2 \to V$, where $F(\uparrow, \downarrow)$ and consider the system $x_{n+1} = F(x_n, x_{n-1})$ with an initial condition $X_0 = (x_0, x_{-1}) \in V^2$. Suppose there is $(a, b) \in V^2$ such that $(a, b) \le_{se} (F(a, b), F(b, a))$ and $(a, b) \le_{se} (x_k, x_{k-1}) \le_{se} (b, a)$ for some $k \ge 0$. If $(F(x, y), F(y, x)) = (x, y)$ has a unique solution in $V^2$, then $\zeta_F$ converges to a global attractor in $V^4$ and $\alpha_F$ converges to a global attractor in $V^2$. If $(F(x, y), F(y, x)) = (x, y)$ has three solutions $(x^*, y^*)$, $(\tilde{y}, \tilde{y})$ and $(y^*, x^*)$, $x^* < y^*$, then

$$x^* \le \liminf \alpha_F \le \limsup \alpha_F \le y^*.$$

A convenient way of stating this result is as follows: Recall that $X_n = T^n(x_0, x_{-1}) = (x_n, x_{n-1})$. Suppose $F(\uparrow, \downarrow)$ and $(a, b) \le_{se} (F(a, b), F(b, a))$. Suppose that the sequence $(X_n)$ is eventually in the box with vertices $A = (a, b)$ and $B = (b, a)$. If $G: V^4 \to V^4$ has a unique fixed point, then it must be of the form $(x, x, x, x)$ and $\alpha_F$ must converge to $x$. Alternatively, the system $T$ has the global attractor $(x, x)$ in that box.

A $q$-cycle $C_q := \{\xi_0, \xi_1, \dots, \xi_{q-1}\}$ of the map $G$ in Eq. (2.2) is a periodic solution of the system $\xi_{n+1} = G(\xi_n)$ with minimal period $q$. Note that $C_q$ could be driven by a $q$-cycle $\{x_0, x_1, \dots, x_{q-1}\}$ of $F$, where

$$\xi_0 = (x_0, x_{q-1}, x_0, x_{q-1}) \quad \text{and} \quad \xi_j = (x_j, x_{j-1}, x_j, x_{j-1}), \quad j = 1, \dots, q-1.$$

Otherwise, we say that $F$ has a pseudo $q$-cycle. This notion has been introduced and classified in [24].

## 2.2 Periodic systems

We consider in this section dynamical systems of the form

$$\alpha: x_{n+1} = F_n(x_n, x_{n-1}), \quad (x_0, x_{-1}) \in V^2, \qquad (2.4)$$

where we assume $F_{n+p} = F_n$, $\forall n \ge 0$, and $p \ge 1$ is the minimal positive integer with such a property. We refer to such a sequence as being “$p$-periodic”. These sequences are treated in [24], and we give a general overview next. We stress that the starting time $n = n_0$ is crucial in non-autonomous systems because it dictates certain synchronization between time and space. In this paper, we consider $n_0 = 0$ throughout. The sequence $\alpha := \{(x_n)\}_{n \ge 0}$ in $V$ gives rise to the sequence $(X_n) := \{(x_n, x_{n-1})\}_{n \ge 0}$ in $V^2$ as before. To stress the role of the individual maps in the system (2.4), we denote it by $[F_0, F_1, \dots, F_{p-1}]$ and its vector form by $[T_0, T_1, \dots, T_{p-1}]$, where $T_n(x_n, x_{n-1}) := (F_n(x, y), x)$, $n \ge 0$. Considering first the period two case, let us write $T_{01} = T_0 \circ T_1$ and $T_{10} = T_1 \circ T_0$ the compositions whose associated sequences in $V^2$ are given by

$$(X_{2n+1}) : (x_{2n+1}, x_{2n}) = T_{01}(x_{2n-1}, x_{2n-2}) \quad \text{and} \quad (X_{2n}) : (x_{2n}, x_{2n-1}) = T_{10}(x_{2n-2}, x_{2n-3})$$

with initial values $(x_0, x_{-1})$ and $(x_1, x_0)$ respectively. These two sequences are interlocked as in Figure 1.

Decomposing the 2-periodic recursive sequence $(X_n)$ into these two subsequences give an embedding $S(V) \to S(V^2) \times S(V^2)$ whose image is $(X_{2n+1}, X_{2n})$. Concatenating these sequences side-by-side produces a recursive sequence $\zeta = \Psi(\alpha)$ in $S(V^4)$ given by the general term

$$\zeta_n := (x_{2n+1}, x_{2n}, x_{2n}, x_{2n-1}) \text{ where } \zeta_0 = (x_1, x_0, x_0, x_{-1}).$$

The above is not a strong embedding in general since a convergence of a pair of subsequences
<figure><img src="image_3.png" /><figcaption>Figure 1</figcaption></figure>

implies convergence of the original sequence unless they both converge to the same limit. In all cases, this construction extends in an obvious manner to p-periodic recursive sequences of finite delay and gives an embedding of the subset of recursive p-periodic sequences of delay N into S(V$^N$).

Starting with Eq. (2.4), p = 2, the system in vector form is written [T₀, T₁], and it breaks into the two “folded” systems T₁₀ and T₀₁ as already indicated. We will relate the k-cycles of T₁₀ and T₀₁ to the 2k-cycles of [T₀, T₁]. Note that this process is known as folding and unfolding, and it has been characterized for one-dimensional maps [25].

**Lemma 2.4.** Define α : x$_{n+1}$ = F$_{n}$ (x$_{n}$, x$_{n}$$_{−}$$_{1}$), (x₀, x$_{−}$₁) ∈ V² which is 2-periodic. Let C$_{k}$ (T₀₁) ⊂ (V²)k (resp. C$_{k}$ (T₁₀)) describe the set of k-cycles of T₀₁ (resp. T₁₀). Then T₀ maps C$_{k}$ (T₀₁) bijectively onto C$_{k}$ (T₁₀), with inverse T₁. If these sets are disjointed, their union is a 2k-cycle of [T₀, T₁]. In particular, [T₀, T₁] has a unique 2-cycle if and only if T₀₁ and T₁₀ have unique fixed points which are distinct.

*Proof.* This is an immediate consequence of the identities $T_1 \circ T_{01} = T_{10} \circ T_1$, and $T_0 \circ T_{10} = T_{01} \circ T_0$. The final claim is a consequence of the fact that the fixed points of $\alpha$ correspond to its one cycle. $\square$

In the next two illustrative examples, we ignore monotonicity, and focus on the structure of cycles.

**Example 2.5.** Let p = k = 2. Here's an example where the 2-cycles for T₀₁ and T₁₀ are different. Define α as in Eq. (2.4) with $F_0(x, y) = y + x^2 - 1$ and $F_1(x, y) = -y$. Then

$$T_0(x, y) = (y + x^2 - 1, x) \text{ and } T_1(x, y) = (-y, x).$$

We can check that $\{(1, y), (-1, y)\}$ is the unique 2-cycle of $T_{10} = T_1 T_0$, while $\{(y, 1), (y, -1)\}$ is the unique 2-cycle of $T_{01} = T_0 T_1$. Clearly, $T_0(1, y) = (y, 1)$ and $T_1(y, 1) = (-1, y)$. The system $[T_0, T_1]$ has a 4-cycle.

**Example 2.6.** Define α as in (2.4) with $F_0(x, y) = xy$ and $F_1(x, y) = \frac{x}{y}$. This is 2-periodic, with subsequences in vector form induced from $T_0(x, y) = (xy, x)$ and $T_1(x, y) = (\frac{x}{y}, x)$. A calculation shows that both systems $T_{10}$ and $T_{01}$ have each a unique three cycle that is common to both

$$\{(-1, -1), (1, -1), (-1, 1)\} = C_3(T_{10}) = C_3(T_{01}).$$

This is also the same 3-cycle for each individual map $T_j$ and for the 2-periodic system $[T_0, T_1]$. This is due to the fact that the periodicity of the system 2 and the periodicity of the cycle 3 are coprime numbers.

Similarly, treating the p-periodic case for any $p \ge 2$ is possible. In that case, the sequence $\alpha$ in (2.4) gives rise to p sequences in $V^2$, indexed by the action of the cyclic group $\mathbb{Z}_p$ with generator $\sigma$ acting on the ordered tuple $(p-1, p-2, \dots, 1, 0)$ by cyclic permutations of coordinates. Here
$\sigma^i(p-1, p-2, \dots, 1, 0) = (i-1, i-2, \dots, 0, \dots, i)$. To this corresponds the sequence $X_{\sigma^i}$ with general recursive term

$$ (x_{n+1}, x_n) = T_{\sigma^i}(x_n, x_{n-1}) \quad \text{where} \quad T_{\sigma^i} = T_{i-1}T_{i-2}\cdots T_0\cdots T_i, \qquad (2.5) $$

$T_i(x, y) = (F_i(x, y), x)$, with initial term $T_{i-1}\cdots T_0(x_0, x_{-1})$. The sequences $X_{\sigma^i}$ partition $(X_n)$. They are mirror copies of each other in the sense that either all converge simultaneously or diverge. If some sequence $X_{\sigma^i}$ has a $k$-cycle, then all other sequences have $k$-cycles. If some sequence has a global attractor, so do all other sequences (see [25], Lemma 2.1). It is of potential interest to understand the way in which $k$-cycles of the $T_{\sigma^i}$s' produce $q$-cycles of $[T_0, T_1, \dots, T_{p-1}]$. The $k$-cycles in $C_k(T_{\sigma^i})$, for different $i$, may overlap or even be equal as illustrated in Example 2.5 and Example 2.6. Note for example that if the $F_i$'s are injective, then a fixed point of $T_{\sigma^i}$ gives rise to a $p$-cycle of $[T_0, T_1, \dots, T_{p-1}]$. This notion is a straightforward generalization of the one-dimensional case in [25].

## 2.3 Embedding in the periodic case

In the autonomous case (Subsection 2.1), we extended a one dimensional system $\alpha_F$ as given in Eq. (2.1) into a four dimensional recursive system $\zeta_F : \zeta_{n+1} = G(\zeta_n)$, where $G$ is a self map of $V^4$. When working in a box (see Definition 2.2), if $G$ had a unique fixed point, then $\alpha_F$ converged to a unique fixed point as well. We seek an analogous result in the periodic case, i.e., Eq. (2.4). It turns out in this case that a unique fixed point of some corresponding embedded system in $V^4$ leads to the existence of a globally attracting cycle of $[T_0, T_1]$ (and $\alpha_F$). We will restrict below to the 2-periodic case, i.e., $p=2$.

Starting with (2.4), $F_0, F_1, F_1 \neq F_0$ and $F_i(\uparrow, \downarrow)$, define similarly

$$ G_0(x, y, u, v) = (F_0(x, y), u, F_0(u, v), x), \quad G_1(x, y, u, v) = (F_1(x, y), u, F_1(u, v), x). $$

We write $G_{01} = G_0 \circ G_1$ and $G_{10} = G_1 \circ G_0$. The four dimensional embedded system $[G_0, G_1]$ is defined by $\zeta_{n+1} = G_n \mod 2(\zeta_n)$. Again, the restriction of $[G_0, G_1]$ to the diagonal subset in $V^2 \times V^2$ is the “diagonal system” $([T_0, T_1], [T_0, T_1])$. The following is a direct check.

**Lemma 2.7.** Assume that $F_i(\uparrow, \downarrow)$ are strictly monotonic in each component. Then both $G_0$ and $G_1$ are one-to-one.

A fixed point of the system $[G_0, G_1]$ must be a fixed point of all the $G_i$'s. In other words, if $\bar{\zeta} \in V^4$ is a fixed point of $[G_0, G_1]$, then

$$ \bar{\zeta} = (\bar{x}, \bar{y}, \bar{y}, \bar{x}), \text{ where } F_j(\bar{x}, \bar{y}) = \bar{x} \text{ and } F_j(\bar{y}, \bar{x}) = \bar{y}, j = 0, 1. \qquad (2.6) $$

Note that a fixed point of $[G_0, G_1]$ is necessarily a fixed point of both $G_{01}$ and $G_{10}$. The converse is not always true. It could happen that $G_{01}$ and $G_{01}$ have different fixed points that form a 2 cycle of $[G_0, G_1]$. The relationship between the cycles of $[F_0, F_1]$ and $[G_0, G_1]$ turns out to be interesting. It has been characterized in [24], but for the sake of clarity and completeness, we illustrate it here. A 2-cycle $\{x_0, x_1\}$ of $[F_0, F_1]$ must satisfy

$$ F_0(x_1, x_0) = x_0 \text{ and } F_1(x_0, x_1) = x_1, \quad x_1 \neq x_0. \qquad (2.7) $$

A 2-cycle of $[G_0, G_1]$ is of the form $\{\zeta_0, \zeta_1\} \in V^4 \times V^4$ such that $G_0(\zeta_0) = \zeta_1$ and $G_1(\zeta_1) = \zeta_0$,
$\zeta_0 \neq \zeta_1$. Let $\zeta = (\bar{x}, \bar{y}, \bar{u}, \bar{v})$ be a fixed point of $G_{10}$, that is $G_1(G_0(\zeta)) = \zeta$. Then we must have

$$ (\bar{u}, \bar{v}) = (F_1(\bar{y}, \bar{x}), F_0(\bar{x}, \bar{y})) \quad \text{and} \quad (\bar{x}, \bar{y}) = (F_1(\bar{v}, \bar{u}), F_0(\bar{u}, \bar{v})). \qquad (2.8) $$

Note that if $\zeta = (\bar{x}, \bar{y}, \bar{u}, \bar{v})$ is a fixed point of $G_{10}$, then $\zeta' = (\bar{u}, \bar{v}, \bar{x}, \bar{y})$ is also a fixed point. A unique fixed point of $G_{10}$ is therefore of the form $(\bar{x}, \bar{x}, \bar{x}, \bar{x})$ or $(\bar{x}, \bar{y}, \bar{x}, \bar{y})$, $\bar{x} \neq \bar{y}$ and satisfying the conditions (2.8). We discuss next the various implications of the existence of the fixed points of $G_{10}$ on the cycles of $[G_0, G_1]$ and $[F_0, F_1]$. Assume $\xi_0 = (\bar{x}, \bar{y}, \bar{u}, \bar{v})$ satisfies the system of equations in (2.8).

(i) Let $\bar{x} = \bar{y}$. We have the following scenarios:

* $\bar{\xi}_0 = (\bar{x}, \bar{x}, \bar{x}, \bar{x})$ and $\bar{x}$ is a common fixed point for each individual map $F_j$. Obviously, $\bar{x}$ will be a fixed point of $[F_0, F_1]$.

* $\bar{\xi}_0 = (\bar{x}, \bar{x}, \bar{u}, \bar{u})$, and $G_0(\bar{\xi}_0) = \bar{\xi}_1 = (\bar{u}, \bar{u}, \bar{x}, \bar{x})$, where $\bar{x} \neq \bar{u}$. In this case, $\{\bar{\xi}_0, \bar{\xi}_1\}$ is a 2-cycle of $[G_0, G_1]$ and $\{\bar{x}, \bar{u}\}$ is a 2-cycle of the one dimensional 2-periodic system $[f_0, f_1]$, where $f_0(t) = F_0(t, t)$ and $f_1(t) = F_1(t, t)$. $\{(\bar{x}, \bar{x}), (\bar{u}, \bar{u})\}$ is dubbed as a pseudo 2-cycle of $[F_0, F_1]$.

(ii) Let $\bar{x} \neq \bar{y}$. We have the following scenarios:

* $\bar{\xi}_0 = (\bar{x}, \bar{y}, \bar{y}, \bar{x})$, where $\bar{\xi}_0$ is a common fixed point of each individual map $G_j$. In this case, $(\bar{x}, \bar{y})$ and $(\bar{y}, \bar{x})$ are dubbed as two pseudo common fixed points of each individual map $F_j$.

* $\bar{\xi}_0 = (\bar{x}, \bar{y}, \bar{x}, \bar{y})$ and $G_0(\bar{\xi}_0) = \bar{\xi}_1 = (\bar{y}, \bar{x}, \bar{y}, \bar{x})$. In this case, $\{\bar{\xi}_0, \bar{\xi}_1\}$ is a 2-cycle of $[G_0, G_1]$ and $\{\bar{y}, \bar{x}\}$ is a 2-cycle of $[F_0, F_1]$.

* $\bar{\xi}_0 = (\bar{x}, \bar{y}, \bar{u}, \bar{v})$ and $\bar{\xi}_1 = (\bar{v}, \bar{u}, \bar{y}, \bar{x})$, where $(\bar{u}, \bar{v})$ is neither $(\bar{x}, \bar{y})$ nor $(\bar{y}, \bar{x})$. In this case, $(\bar{u}, \bar{v}) = (F_1(\bar{y}, \bar{x}), F_0(\bar{x}, \bar{y}))$ and $C_2 := \{\bar{\xi}_0, \bar{\xi}_1\}$ is a 2-cycle of $[G_1, G_0]$, while $\{(\bar{x}, \bar{y}), (\bar{v}, \bar{u})\}$, $\{(\bar{u}, \bar{v}), (\bar{y}, \bar{x})\}$ are dubbed as pseudo 2-cycles of $[F_0, F_1]$.

Next, we need to focus on stability. Given a unique fixed point of $G_{01}$, and provided we are working “in a box”, we deduce that $T_{01}$ also has a unique fixed point. Similarly for $T_{10}$. By Lemma 2.4, $[T_0, T_1]$ has a 2-cycle if the fixed point is not common. This must be a globally attracting cycle [24], and we extract the needed result below.

**Proposition 2.8.** Let $F_j : \mathbb{R}_+^2 \to \mathbb{R}_+$, where $F_j(\uparrow, \downarrow)$ for each $j = 0, 1$. Consider the system $[F_0, F_1]$ with initial conditions $X_0 = (x_0, x_{-1}) \in \mathbb{R}_+^2$. Suppose there exists $(a, b) \in \mathbb{R}_+^2$ such that $a < b$, $(a, b) \leq_{se} (x_k, x_{k-1}) \leq_{se} (b, a)$ for some $k \geq 0$ and

$$ a < \min\{F_0(a, b), F_1(F_0(a, b), b)\}, \quad b > \max\{F_0(b, a), F_1(F_0(b, a), a)\}. \qquad (2.9) $$

If $G_{10} = G_1 \circ G_0$ has a unique fixed point in $\mathbb{R}_+^4$, then the 2-periodic system $[F_0, F_1]$ has a global attractor, which can be an equilibrium solution or a 2-cycle.

*Proof.* The condition $(a, b) \leq_{se} (x_k, x_{k-1}) \leq_{se} (b, a)$ means that the sequence $(X_n)$ is eventually trapped in the box $[a, b]^2$. The condition (2.9) gives us $(A, B) < G_{10}(A, B)$, where $A = (a, b)$ and $B = (b, a)$. Since $G_{10}$ is monotonic with respect to $\leq_{se}$, we can use our formalism directly (Section 2.1, inequalities (2.3)) to infer that the sequence $G_{10}^n(X_0, X_0)$ converges to the unique fixed point of $G_{10}$ in $[a, b]^4$. If the unique fixed point of $G_{10}$ is of the form $\bar{\xi} = (\bar{x}, \bar{x}, \bar{x}, \bar{x})$, then $\bar{x}$ is an equilibrium
point of $[F_0, F_1]$ and both $\{x_{2n}\}, \{x_{2n-1}\}$ converge to $\bar{x}$. On the other hand, if the unique fixed point of $G_{10}$ is of the form $\bar{\xi} = (\bar{x}, \bar{y}, \bar{x}, \bar{y})$, $\bar{x} \neq \bar{y}$ then $\{\bar{x}, \bar{y}\}$ is a 2-cycle of $[F_0, F_1]$, and the even terms of the orbit $\{x_n\}$ converge to one branch of the 2-cycle, while the odd terms converge to the other branch. Hence, whether the system $[F_0, F_1]$ has an equilibrium or a 2-cycle, it serves as a global attractor. $\square$

We remark that the inequalities (2.9) can be replaced by the more simple inequalities $(a, b) <_{se} (F_j(a, b), F_j(b, a))$, $j = 0, 1$. However, simplicity here comes at the expense of generalization.

### 3 Stability under constant stocking

Since $h > 0$, all solutions of Eq. (1.3) must satisfy $y_n > h$ for all $n$. Furthermore, an equilibrium solution $\bar{y}$ satisfies $(1 - \frac{h}{\bar{y}}) = e^{r-\bar{y}}$, and it is obvious that we obtain a unique equilibrium solution. To stress the role of the parameter $h$ in our analysis, we denote the positive equilibrium solution by $\bar{y}_h$. The Jacobian matrix at $\bar{y}_h$ is given by

$$J(\bar{y}_h, \bar{y}_h) = \begin{bmatrix} e^{r-\bar{y}_h} & -\bar{y}_h e^{r-\bar{y}_h} \\ 1 & 0 \end{bmatrix} = \begin{bmatrix} 1 - \frac{h}{\bar{y}_h} & h - \bar{y}_h \\ 1 & 0 \end{bmatrix}.$$

Let the trace of the Jacobian matrix be $T$ and the determinant be $D$, then $0 < T = 1 - \frac{h}{\bar{y}_h} < 1$ and $D = \bar{y}_h - h > 0$. Furthermore, we have

$$D - T = \bar{y}_h - h + \frac{h}{\bar{y}_h} - 1 > -1.$$

Based on the Jury conditions for stability, the eigenvalues of $J$ are within the unit disk when $|T| < 1 + D < 2$, and consequently, the equilibrium $\bar{y}_h$ is locally asymptotically stable (LAS) when $D < 1$, i.e., $\bar{y}_h < 1 + h$. Since

$$0 < e^{r-\bar{y}_h} = 1 - \frac{h}{\bar{y}_h} < 1,$$

we obtain $\bar{y}_h > h$ and $\bar{y}_h > r$. In fact, $\bar{y}_0 = r$. The next lemma is relevant in the sequel and describes the behavior of the fixed point in terms of $r$ and $h$.

**Lemma 3.1.** The fixed point $\bar{y}_h$ of Eq. (1.3) is increasing in both $h$ and $r$.

*Proof.* This is an immediate consequence of implicitly differentiating $(1 - \frac{h}{\bar{y}_h}) = \frac{d}{dh}e^{r-\bar{y}_h}$ with respect to both $h$ and $r$. In the first case, we obtain

$$(\bar{y}_h - h + \frac{h}{\bar{y}_h}) \frac{d\bar{y}_h}{dh} = 1.$$

So that $\frac{d\bar{y}_h}{dh} > 0$. Similarly, $\frac{d\bar{y}_h}{dr} > 0$ and $\bar{y}_h$ is increasing in $r$. $\square$

Now, the question is: can $\bar{y}_h$ reach $h + 1$ and lose its stability? If yes, then we must have $h + 1 = e^{h+1-r}$. This is possible when

$$r = r_1 := h + 1 - \ln(h + 1). \qquad (3.1)$$

**Lemma 3.2.** Consider Eq. (1.3) with $h > 0$, and let $r_1$ be as defined in Eq. (3.1). Each of the following holds true:
(i) Eq. (1.3) has a unique non-negative equilibrium solution $\bar{y}_h$ and $\bar{y}_h > \max\{r, h\}$.

(ii) If $0 < r < r_1$, then $\bar{y}_h$ is LAS, while if $r > r_1$, then $\bar{y}_h$ is unstable.

*Proof.* Part (i) is clear. It remains to clarify Part (ii). Recall that $\bar{y}_h$ is LAS if and only if $\bar{y}_h < h+1$; in turn, we claim this happens iff $r < r_1$. But $h+1$ is the value of $\bar{y}_h$ when $r = r_1$, and so the claim is a consequence of the fact that $\bar{y}_h$ is increasing in $r$. When $r = r_1$, the claim, and its proof are standard. □

When $r = r_1$, Eq. (1.3) has the potential to go through a Neimark-Sacker bifurcation (see Fig. 2). Note that since $r_1$ is increasing in $h$, $0 < r < r_0 = 1$ gives LAS regardless of the value of $h$. Our next example illustrates the validity of Part (ii) of Lemma (3.2). Also, Fig. 5 shows the stability region in the $(h, r)$-plane.

**Example 3.3.** Let $h = e-1$ and $r = 2$. Then $r_1 = e-1 < 2 = r$. In this case, we obtain $\bar{y}_h \approx 2.898$, and the eigenvalues of the Jacobian matrix are $\lambda \approx 0.204 \pm 1.067i$. They are out of the unit disk, and consequently, the equilibrium is unstable. On the other hand, if we consider $h = e-1$ and $r = \frac{3}{2}$, then $r_1 = e-1 > r$, $\bar{y}_h \approx 2.589$ and, the eigenvalues of the Jacobian matrix are $\lambda \approx 0.168 \pm 0.918i$. The eigenvalues are located within the unit disk, and consequently, the equilibrium is LAS.

<figure><img src="image_4.png" /><figcaption>Figure 2: Let $r = r_1$. This figure shows the equilibrium solution $y_h$ (solid green circle) undergoes a Neimark-Sacker bifurcation at $h = 1.0$. The red color represents the even terms in the orbits, while the blue represents the odd terms. The horizontal axis is for $x_n$, and the vertical axis is for $x_{n-1}$.</figcaption></figure>

Next, we define the 4-dimensional map

$$G: \mathbb{R}_+^4 \to \mathbb{R}_+^4, \quad \text{where} \quad G(x, y, u, v) = (xf(y) + h, u, uf(v) + h, x).$$

Again here, $f(y) = e^{r-y}$. We solve $G(\xi) = \xi$ to find the fixed points of $G$. We obtain $\xi = (x^*, y^*, y^*, x^*)$, where $x^* = x^*f(y^*) + h$ and $y^* = y^*f(x^*) + h$, with $x^* \le y^*$. The case $x^* = y^*$ is when both are equal $\bar{y}_h$, and we write the symmetric solution $\tilde{\xi} = \tilde{\xi}_h = (\bar{y}_h, \bar{y}_h, \bar{y}_h, \bar{y}_h)$.
Let $L_1$ and $L_2$ denote the curves of $y = yf(x) + h$ and $x = xf(y) + h$, respectively. We discuss how these curves intersect since the intersection points (i.e. the “solutions”) correspond to the fixed points of $G$. The fixed point $(\bar{y}_h, \bar{y}_h)$ is always a solution. Clearly, $L_1$ has a vertical asymptote at $x = r$ and a horizontal asymptote at $y = h$, whereas $L_2$ has a vertical asymptote at $x = h$ and a horizontal asymptote at $y = r$. A variation of $h$ makes asymmetric solutions bifurcate from the symmetric ones. The slope $dy/dx$ of $L_1$ at $\bar{y}_h$ is $\frac{1}{h}\bar{y}_h(h - \bar{y}_h)$, and the reciprocal of this expression is the slope of $L_2$. Since the slope is negative and $\bar{y}_h > h$, we consider the solution of $\bar{y}_h(\bar{y}_h - h) = h$ to be

$$H^* := \frac{1}{2} \left( h + \sqrt{h^2 + 4h} \right). \qquad (3.2)$$

In this case,

$$r = r_2 := H^* + \ln(H^* - h) - \ln(H^*), \quad (\text{and } \bar{y}_h = H^*). \qquad (3.3)$$

Notice that we always have (we skip the computation) that

$$r_2 < r_1 \quad \text{and} \quad r_1 < h \quad \text{if} \quad h > e - 1. \qquad (3.4)$$

To recapitulate, for a given $h$, and for $r = r_2$, the fixed point of Eq. (1.3) is $\bar{y}_h = H^*$, and at the solution $x = y = H^*$, both curves have slope $-1$. The number of solutions of $L_1 \cap L_2$ depends on whether $\bar{y}_h$ is larger or smaller than $H^*$, but this also depends on whether $h < r$ or $h > r$. From a combinatorial standpoint, there are four options; however, $(h, H^*) <_{se} (r, \bar{y}_h)$ is not viable according to the following result.

**Lemma 3.4.** If $h \le r$, then $H^* \le \bar{y}_h$.

*Proof.* Because $1 + t \le e^t$ for all $t \in \mathbb{R}$, then $1 + h - \bar{y}_h \le e^{h-\bar{y}_h}$. Since $h \le r$, we obtain that $1 + h - \bar{y}_h \le e^{r-\bar{y}_h}$, and consequently,

$$1 - e^{r-\bar{y}_h} \le \bar{y}_h - h \iff \bar{y}_h(1 - e^{r-\bar{y}_h}) \le \bar{y}_h(\bar{y}_h - h).$$

From the fact that $\bar{y}_h$ is an equilibrium point, we obtain $h = \bar{y}_h(1 - e^{r-\bar{y}_h})$, and therefore, $h \le \bar{y}_h(\bar{y}_h - h)$, which gives us $\bar{y}_h \ge H^*$. $\square$

The next lemma summarizes and settles all cases.

**Lemma 3.5.** Consider $r_2$ as defined in Eq. (3.3). Let $L_1$ and $L_2$ denote the curves of $y = F(y, x) = yf(x) + h$ and $x = F(x, y) = xf(y) + h$, respectively. Each of the following holds true:

(i) If $h < r$, then $H^* < \bar{y}_h$ and $L_1$ and $L_2$ intersect at the unique point $(\bar{y}_h, \bar{y}_h)$.

(ii) If $r < r_2 < h$, then $\bar{y}_h < H^*$ and $L_1$ and $L_2$ intersect at the unique point $(\bar{y}_h, \bar{y}_h)$.

(iii) If $r_2 < r < h$, then $H^* < \bar{y}_h$ and $L_1$ and $L_2$ intersect at three points denoted by $(\bar{y}_h, \bar{y}_h)$, $(x^*, y^*)$ and $(y^*, x^*)$, where $x^* \ne y^*$.

*Proof.* (i) The two curves $L_1$ and $L_2$ always intersect at the point $(\bar{y}_h, \bar{y}_h)$ along the diagonal. Since the curves are symmetric images of one another with respect to the $y = x$ axis, we can confine ourselves to the case $x > \bar{y}_h > \max\{r, h\}$ and investigate other intersection points. Express $y$ explicitly for both $L_1$ and $L_2$ to find, respectively, that

$$y_1 = \frac{h}{1 - e^{r-x}} \quad \text{and} \quad y_2 = r - \ln \left( 1 - \frac{h}{x} \right).$$
We show $y_2 - y_1 > 0$ for all $x > \bar{y}_h$ when $h < r$. For a fixed value of $x$, we show that $y_2 - y_1$ is positive and minimum at $h = r$. We have

$$ \frac{d}{dh}(y_2 - y_1) = \frac{-1}{1 - e^{r-x}} + \frac{1}{x - h} $$

and $e^{r-x} > 1 + r - x$, for all $x > r$. Therefore, $\frac{1}{1 - e^{r-x}} > \frac{1}{x - r}$, and since $h < r < \bar{y}_2 < x$, we obtain

$$ \frac{d}{dh}(y_2 - y_1) < \frac{-1}{x - r} + \frac{1}{x - h} < 0. $$

So $y_2 - y_1$ is decreasing and takes its minimum at $h = r$. At $h = r$,

$$ y_2 - y_1 = \frac{-r}{1 - e^{r-x}} + r - \ln\left(1 - \frac{r}{x}\right) > 0. $$

Hence, the minimum of $y_2 - y_1$ is strictly positive, and $y_2 \neq y_1$. This shows that $L_1$ and $L_2$ cannot intersect again for $x > \bar{y}_h > \max\{r, h\}$, and also by symmetry, for $x < \bar{y}_h$. This completes the proof of Part (i).

To show (ii) and (iii), we have shown in Lemma 3.1 that the fixed point $\bar{y}_h$ is an increasing function of $r$. Therefore if $r < r_2 < h$, then necessarily $\bar{y}_h < H^*$, since $H^*$ is the fixed point when $r = r_2$. This establishes the first part of the claims in (ii) and (iii). The second part in each case follows along the same outline as in (i). □

Figures 3 and 4 illustrate all three possible scenarios in Lemma 3.5.

<figure><img src="image_5.png" /><figcaption>Figure 3: Part (i) of this figure shows the curves of $L_1$ and $L_2$ when $H^* \le \bar{y}_h$ and $h < r$. The unique intersection point is $(\bar{y}_h, \bar{y}_h)$. Part (ii) shows the unique intersection point $(\bar{y}_h, \bar{y}_h)$ but when $\bar{y}_h \le H^*$ and $h > r$. The shaded region in Part (ii) represents the solution of $(x, y) \le_{se} (F(x, y), F(y, x))$ in which $x \le y$.</figcaption></figure>

**Remark 3.1.** The inequalities $(x, y) \le_{se} (F(x, y), F(y, x))$ and $x \le y$ have a feasible solution in case (ii) and case (iii) of Lemma 3.5, that is when $h > r$. This is essential for applying the embedding technique. Concerning the intersections between the curves $L_1$ and $L_2$, let $g_1(t) = \frac{h_0}{1 - f(t)}$, $t > r$. The function $g_1$ is decreasing in $t$. Since an intersection point $(a, b)$ must satisfy $(a, b) = (g_1(b), g_1(a))$, the fixed point of $g_1$ is the fixed point of $F$, and the 2-cycles of $g_1$ form the pseudo fixed points of $F$. Since a period doubling bifurcation occurs when the fixed point of $g_1$ loses its stability, the existence of pseudo fixed points can be investigated through the local stability of the fixed point of $g_1$.
<figure><img src="image_6.png" /><figcaption>Figure 4: This figure shows the curves of $L_1$ and $L_2$ when $H^* \le \bar{y}_h$ and $r < h$. The three intersection points are emphasized, and for the given choice of our parameters, they are (2.741, 4.969), (3.424, 3.424) and (4.969, 2.741). The shaded region represents the solution of $(x, y) \le_{se} (F(x, y), F(y, x))$ in which $x \le y$.</figcaption></figure>

After establishing the needed machinery and lemmas, we can present one of our first main results. Recall that the fixed points of the embedded map $G$ are associated with the intersection points of $L_1$ and $L_2$. In particular, $G$ has a unique fixed point given by $\xi = \bar{\xi}_h$ when cases (i) and (ii) of Lemma 3.5 are satisfied, and $G$ has three fixed points in case (iii) of that Lemma. As done earlier, we denote the asymmetric fixed points of $G$ by $\bar{\xi}_1 = (x^*, y^*, y^*, x^*)$ and $\bar{\xi}_2 = (y^*, x^*, x^*, y^*)$ where $x^* < y^*$. In other words, $(x^*, y^*)$ and $(y^*, x^*)$ are the pseudo fixed points of $F$.

**Theorem 3.6.** Consider $r_2$ as given in Eq. (3.3). Each of the following holds true for Eq. (1.3):

(i) If $r \le r_2$, then $\bar{y}_h$ is globally asymptotically stable.

(ii) If $r_2 < r < h$, then the compact box $[x^*, y^*]^2$ forms an absorbing region. In particular,

$$x^* \le \liminf\{x_n\} \le \limsup\{x_n\} \le y^*.$$

*Proof.* (i) Since $r \le r_2$, Lemma 3.5 guarantees a unique intersection between $L_1$ and $L_2$, and consequently, the embedding $G$ has a unique fixed point, which we denoted by $\bar{\xi}_h = (\bar{y}_h, \bar{y}_h, \bar{y}_h, \bar{y}_h)$. Since $r_2 < r_1$, $\bar{y}_h$ is LAS by Lemma 3.2. To show global stability, notice that the inequalities $(a, b) \le_{se} (F(a, b), F(b, a))$ and $a < b$ have a feasible region (as illustrated in Part (ii) of Fig. 3 and Remark 3.1). For any initial condition $(x_0, x_{-1})$ of Eq. (1.3), there is $(a, b)$ in the feasible region such that Proposition 2.3 is applicable. This implies that for all $\xi_0 = (y_0, y_{-1}, y_0, y_{-1})$ such that

$$\xi = (a, b, b, a) \le_{se} \xi_0 \le_{se} (b, a, a, b) = \xi^t, \quad (3.5)$$

$G^n(\xi_0)$ converges to $\bar{\xi}_h$ and Eq. (1.3) with the initial condition $(y_0, y_{-1})$ converges to $\bar{y}_h$. This verifies Part (i). Part (ii) is similar and refers to Part (iii) of Lemma 3.5 to obtain the three fixed points of $G$ denoted by $\bar{\xi}_1, \bar{\xi}_2$ and $\bar{\xi}_h$. We have the same scenario as in Condition (3.5) for initial conditions $(x_0, x_{-1})$ and the existence of $(a, b)$ in the feasible region as above so that $G^n(\xi_0)$ must converge to a fixed point of $G$, and by the embedding result, the orbit of Eq. (1.3) that starts
at $(y_0, y_{-1})$ will be eventually squeezed between $x^*$ and $y^*$ (Proposition 2.3). This completes the proof. □

It is worth stressing that $r_2 < r_1$ always, and $0 < r < r_1$ guarantees the local stability of $\bar{y}_h$. However, when $r_2 < r < \min\{h, r_1\}$, we obtain the pseudo fixed points, which cripple the embedding technique. Nevertheless, an absorbing region was obtained as given in Part (ii) of Theorem 3.6. Finally, when $h < r < r_1$, $\bar{y}_h$ is LAS. However, the absence of a feasible solution for the system $(x, y) \le_{se} (F(x, y), F(y, x))$ and $x \le y$ led to the unsuccessful outcome of our embedding approach for addressing global stability. Figure 5 illustrates and summarizes our stability results in the parameter space.

<figure><img src="image_7.png" /><figcaption>Figure 5: This figure shows the stability regions in the $(h, r)$-plane. The red curve is $r = r_1$, which represents the solution of $\bar{y}_h = 1 + h$, while the blue curve is $r = r_2$, which represents the solution of $H^* = H^*f(H^*) + h$. We have local stability in the region below $r = r_1$, but the embedding technique fails to address global stability when $r_2 < r < h$. However, we conjecture that we also obtain global stability when $r_2 < r < r_1$. Observe that global stability in the green part of the $r$-axis, i.e., $h = 0$ and $0 < r < 1$ has been addressed in [15].</figcaption></figure>

# 4 Stability under periodic stocking

In this section, we consider the Ricker model with delay and $p$-periodic stocking as given in Eq. (1.4), i.e., the minimal period of the stocking sequence is $p$. Recall that we give ourselves the liberty to use $f(x)$ instead of $\exp(r - x)$ for convenience of our writing. Define the sequence of 2-dimensional maps $\{T_j\}$,

$$T_j : \mathbb{R}_+^2 \to \mathbb{R}_+^2 \quad \text{as} \quad T_j(x, y) = (xf(y) + h_j, x).$$

The equation $X_{n+1} = T_{n \mod p}(X_n)$ is just a vector form of Eq. (1.4). Since

$$\begin{aligned} y_{n+2} &= y_n f(y_n) f(y_{n-1}) + h_0 f(y_n) + h_1 \\ &\le (e^r + h_0)e^r + h_1, \end{aligned}$$

the composition operator $\tilde{T} := T_{p-1} \circ T_{p-2} \circ \cdots \circ T_0$ maps a nonempty convex compact set into itself. Therefore, $\tilde{T}$ has a fixed point in that domain. Because $T_i(x, y) = T_j(x, y)$ if and only if $i = j$, then a fixed point of $\tilde{T}$ must form a $p$-periodic solution of $X_{n+1} = T_{n \mod p}(X_n)$, which reflects a $p$-periodic solution of Eq. (1.4). We are interested in the global stability of the obtained $p$-periodic solution, and we focus on the case $p = 2$. The existing 2-periodic solution is in fact the 2-cycle of
the system $[F_0, F_1]$, and we denote it by $\{\bar{z}_0 = \bar{z}_0(h_0, h_1), \bar{z}_1 = \bar{z}_1(h_0, h_1)\}$, where

$$ \bar{z}_0 = \bar{z}_1 f(\bar{z}_0) + h_1 \quad \text{and} \quad \bar{z}_1 = \bar{z}_0 f(\bar{z}_1) + h_0. \qquad (4.1) $$

Observe that $\bar{z}_0 > h_1$ and $\bar{z}_1 > h_0$. Furthermore, $\bar{z}_0$ and $\bar{z}_1$ contribute to the formation of two fixed points of $\tilde{T} = T_1 \circ T_0$, namely $\bar{X}_0 := (\bar{z}_0, \bar{z}_1)$ and $\bar{X}_1 := (\bar{z}_1, \bar{z}_0)$. Now, we give the following fact:

**Proposition 4.1.** $h_0 > h_1$ if and only if $\bar{z}_1 > \bar{z}_0$.

*Proof.* From Eqs. (4.1), we obtain

$$ \bar{z}_1 - \bar{z}_0 = \bar{z}_0(f(\bar{z}_1)) - \bar{z}_1 f(\bar{z}_0) + h_0 - h_1. $$

This gives us

$$ \frac{h_0 - h_1}{\bar{z}_1 - \bar{z}_0} = 1 - \bar{z}_0 \frac{f(\bar{z}_1) - f(\bar{z}_0)}{\bar{z}_1 - \bar{z}_0} + f(\bar{z}_0) > 0. $$

Hence, $sign(h_0 - h_1) = sign(\bar{z}_1 - \bar{z}_0)$, which completes the proof. $\square$

**Remark 4.1.** Since $h_0 < h_1$ and $h_1 < h_0$ give similar behaviour, we focus on the case $h_1 < h_0$. In this case, we have $h_1 < \bar{z}_0 < \bar{z}_1$.

To investigate the stability of the 2-cycle, it is enough to investigate the stability of the fixed point $\bar{X}_0 = (\bar{z}_0, \bar{z}_1)$ under $\tilde{T} = T_1 \circ T_0$. Since

$$ \tilde{T}(x, y) = \begin{bmatrix} (x f(y) + h_0) f(x) + h_1 \\ x f(y) + h_0 \end{bmatrix}, $$

the Jacobian matrix at $\bar{X}_0$ is

$$ J(\bar{X}_0) = \begin{bmatrix} \bar{z}_1 f'(\bar{z}_0) + f(\bar{z}_0) f(\bar{z}_1) & \bar{z}_0 f(\bar{z}_0) f'(\bar{z}_1) \\ f(\bar{z}_1) & \bar{z}_0 f'(\bar{z}_1) \end{bmatrix}. $$

The trace and determinant here are given by

$$ Tr = \bar{z}_1 f'(\bar{z}_0) + \bar{z}_0 f'(\bar{z}_1) + f(\bar{z}_0) f(\bar{z}_1) \quad \text{and} \quad Det = \bar{z}_0 \bar{z}_1 f'(\bar{z}_1) f'(\bar{z}_0). $$

We re-write the two expressions as

$$ Tr = (h_1 - \bar{z}_0) + (h_0 - \bar{z}_1) + \left(1 - \frac{h_0}{\bar{z}_1}\right) \left(1 - \frac{h_1}{\bar{z}_0}\right), \qquad (4.2) $$

$$ Det = (\bar{z}_1 - h_0)(\bar{z}_0 - h_1). \qquad (4.3) $$

Obviously, $Det > 0$ and, $Det < 1$ when $(\bar{z}_1 - h_0)(\bar{z}_0 - h_1) < 1$. It turns out that this condition is sufficient for the local stability of the 2-cycle as the following result shows.

**Lemma 4.2.** Let $h_0 \neq h_1$ and both are non-negative. The 2-cycle of $[F_0, F_1]$ (i.e., $\{\bar{z}_0, \bar{z}_1\}$) is LAS stable if $(\bar{z}_1 - h_0)(\bar{z}_0 - h_1) < 1$.

*Proof.* We check the Jury conditions here. We have $0 < Det < 1$. It remains to show that

$$ Det - Tr > -1 \quad \text{and} \quad Det + Tr > -1. $$
The first inequality is satisfied because

$$ Det - Tr + 1 = (\bar{z}_1 - h_0 + 1)(\bar{z}_0 - h_1 + 1) - \left(1 - \frac{h_0}{\bar{z}_1}\right) \left(1 - \frac{h_1}{\bar{z}_0}\right) > 0. \quad (4.4) $$

Now, we use this fact to write

$$ (\bar{z}_1 - h_0 + 1)(\bar{z}_0 - h_1 + 1) > \left(1 - \frac{h_0}{\bar{z}_1}\right) \left(1 - \frac{h_1}{\bar{z}_0}\right), $$

or equivalently,

$$ \bar{z}_1 \left(1 - \frac{1}{\bar{z}_1 - h_0}\right) \bar{z}_0 \left(1 - \frac{1}{\bar{z}_0 - h_1}\right) > 1. $$

Next, write

$$ Det + Tr + 1 = (\bar{z}_0 - h_1 - 1)(\bar{z}_1 - h_0 - 1) + \left(1 - \frac{h_1}{z_0}\right) \left(1 - \frac{h_0}{z_1}\right) \\ = \left(1 - \frac{h_1}{z_0}\right) \left(1 - \frac{h_0}{z_1}\right) \left[ \bar{z}_1 \left(1 - \frac{1}{\bar{z}_1 - h_0}\right) \bar{z}_0 \left(1 - \frac{1}{\bar{z}_0 - h_1}\right) + 1 \right], $$

which is positive. This completes the proof. □

Note that the condition $(\bar{z}_1 - h_0)(\bar{z}_0 - h_1) < 1$ can be replaced by the more simple (but unnecessary) condition $r \le 1$. Also, the simple conditions $\bar{z}_j \le h_{j+1} + 1$ for $j = 0, 1$ are sufficient to make $Det < 1$. For the reader's convenience, we summarize some special cases in the following corollary:

**Corollary 4.3.** Each of the following cases is valid:

(i) If $r \le 1$, then $Det < 1$.

(ii) If $\bar{z}_j \le h_{j+1} + 1$ for all $j = 0, 1$, then the 2-cycle is LAS.

(iii) If $\bar{z}_j > h_{j+1} + 1$ for all $j = 0, 1$, then the 2-cycle is unstable.

Now, we provide some examples that illustrate the validity of Lemma 4.2.

**Example 4.4.** Consider Eq. (1.4) in each one of the following cases:

(i) Let $h_0 = 2.0, h_1 = 6.444$ and $r = 3.0$. We obtain $\bar{z}_0 \approx 6.635$ and $\bar{z}_1 \approx 7.237$. In this case, the condition of Lemma 4.2 is not satisfied. Furthermore, $Det \approx 1.000$ and $Tr \approx -5.407$. This makes $Det + Tr + 1 < 0$, and consequently, the 2-cycle $\{\bar{z}_0, \bar{z}_1\}$ of the system $[F_0, F_1]$ is unstable. In fact, we obtain the 4-cycle

$$ \{7.049, 19.611, 6.479, 2.000\}. $$

(ii) Let $r = 2.0, h_0 = 2.156$ and $h_1 = 2.720$. We obtain $\bar{z}_0 \approx 3.462$ and $\bar{z}_1 \approx 3.199$. In this case, Lemma 4.2 is satisfied, and the 2-cycle $\{\bar{z}_0, \bar{z}_1\}$ of the system $[F_0, F_1]$ is LAS. Furthermore, $Det \approx 0.774, Tr \approx -1.715$ and $Det + Tr + 1 \approx 0.059 > 0$, and consequently, the eigenvalues are $\approx -0.858 \pm 0.196i$, which are within the unit disk.

(iii) Let $h_0 = 0.820$ and $h_1 = 1.800$. We let $r = 1.5$ to obtain $\bar{z}_0 \approx 2.552$ and $\bar{z}_1 \approx 2.151$. In this case, the condition of Lemma 4.2 is not satisfied. However, by considering $r$ as a bifurcation
parameter, a Neimark-Sacker bifurcation occurs near $r = 1.5$. Indeed, the 2-cycle bifurcates into two invariant curves that serve as one attractor. Fig. 6 illustrates this scenario.

Based on Lemma 4.2, The 2-cycle of Eq. (1.4) has the potential to lose its stability after going through a Neimark-Sacker bifurcation. The period-doubling observed in Part (i) of Example 4.4 comes after the 2-cycle exits the stability region through the Neimark-Sacker bifurcations. Figure 6 Shows an illustrative case of the Neimark-Sacker bifurcation that occurs.

<figure><img src="image_8.png" /><figcaption>Figure 6: Let $h_0 = 0.820$ and $h_1 = 1.800$. When $r = 1.5$, the 2-cycle bifurcates into two curves through a Neimark-Sacker bifurcation. The two curves serve as one attractor. The red represents the even terms in the orbit, while the blue represents the odd terms. The green solid circles represent the 2-cycle. The horizontal axis is for $x_n$, and the vertical axis is for $x_{n-1}$.</figcaption></figure>

Next, we proceed to find conditions under which we obtain global stability. Recall $F_j(x, y) = xf(y) + h_j$ for $j = 0, 1$. First, we investigate the existence of pseudo 2-cycles. Solving the equation $G_1(G_0(\xi)) = \xi$, where $\xi = (x, y, u, v)$ gives us (see Subsection 2.3)

$$ (x, y) = (F_1(v, u), F_0(u, v)) \quad \text{and} \quad (u, v) = (F_1(y, x), F_0(x, y)). $$

Since each value of $(x, y)$ determines a unique value of $(u, v)$, it is enough to focus on the equations

$$ \begin{cases} x - h_1 = (xf(y) + h_0)f(yf(x) + h_1) \\ y - h_0 = (yf(x) + h_1)f(xf(y) + h_0). \end{cases} \quad (4.5) $$

Observe that the 2-cycle $\{\bar{z}_0, \bar{z}_1\}$ (see Eq. (4.1)) satisfies System (4.5). For other solutions, our best bargain here is to investigate the various scenarios based on the curves of the two equations. Let $\ell_1$ and $\ell_2$ be the curves of the first and second equations in System 4.5, respectively. The $\ell_1$-curve has a horizontal asymptote at $y = 2r - h_1$ and a vertical asymptote at $x = h_1$. Similarly, the $\ell_2$-curve has a horizontal asymptote at $y = h_0$ and a vertical asymptote at $x = 2r - h_0$. Multiple intersections between $\ell_1$ and $\ell_2$ guarantee the existence of artificial cycles, and in this case, the
embedding technique fails to help us establish global stability. Thus, a unique intersection between $l_1$ and $l_2$ is crucial to our embedding strategy. Next, we appeal to the results of Subsection 2.3 and Proposition 2.8 in particular. Our primary objective is to identify values for $(a, b)$ such that $a < b$ and the condition (2.9) of Proposition 2.8 are satisfied. By Part (ii) and Part (iii) of Lemma 3.5, and the consequent Remark 3.1, the inequalities $(a, b) <_{se} (F_0(a, b), F_0(b, a))$ have a feasible region when $h_0 > r$. However, we need the feasible region to overlap with the feasible region of the inequalities

$$a \le F_1(F_0(a, b), b) \quad \text{and} \quad b \ge F_1(F_0(b, a), a). \qquad (4.6)$$

A quick observation here is that if $(a, b) <_{se} (F_j(a, b), F_j(b, a))$ for each $j$, then the inequalities in (4.6) have a feasible region; however, we can avoid this strong constraint. Before we proceed, we give some illustrative computational examples.

**Example 4.5.** (i) Consider $r = 1.0, h_0 = 2.0$ and $h_1 = 1.5$. The 2-cycle is $\{\bar{z}_0, \bar{z}_1\} \approx \{2.230, 2.498\}$. From Eqs. (4.2) and (4.3), we obtain $Tr \approx -1.163$, and $Det \approx 0.364$. The eigenvalues are $\approx -0.582 \pm 0.161i$, which are within the unit disk. Thus, the 2-cycle $\{\bar{z}_0, \bar{z}_1\}$ of the system $[F_0, F_1]$ is locally asymptotically stable. The curves $l_1$ and $l_2$ are given in Part (i) of Fig. 7.

(ii) When $r = 1.0, h_0 = 2.0$ and $h_1 = 1.0$, we obtain the 2-cycle $[\bar{z}_0, \bar{z}_1] \approx [2.455, 1.950]$ of the system $[F_0, F_1]$. The trace and determinant from Eqs. (4.2) and (4.3) are given by $Tr \approx -1.42$ and $Det \approx -0.0732$. The eigenvalues are $\approx -1.470$ and $0.050$. Obviously, the 2-cycle of the system $[F_0, F_1]$ is unstable. Next, we investigate the solution of the system $G_1(G_0(\xi)) = \xi$. As clarified in Subsection 2.3, we explore the following scenarios:

* $G_1(\xi) = G_0(\xi) = \xi$, where $\xi = (x, y, y, x)$. $x = y$ leads to an equilibrium solution of Eq. (1.4), while $x \ne y$ leads to an artificial equilibrium solution. Both scenarios are not possible since $h_0 \ne h_1$.

* $\xi_0 = (x, x, y, y), x \ne y, G_0(\xi_0) = \xi_1 = (y, y, x, x)$ and $G_1(\xi_1) = \xi_0$. This means $\{x, y\}$ is a 2-cycle for each one of the one-dimensional maps $f_j(x) = F_j(x, x), j = 0, 1$. Again, this is not possible since $h_0 \ne h_1$.

* $\xi_0 = (x, y, x, y), x \ne y, G_0(\xi_0) = \xi_1 = (y, x, y, x)$ and $G_1(\xi_1) = \xi_0$. This means $\{\xi_0, \xi_1\}$ is a 2-cycle of $[G_0, G_1]$ and $\{x, y\}$ is a 2-cycle of $[F_0, F_1]$. This is always possible in our system since we verified the existence of a p-cycle for the p-periodic system $[F_0, F_1, \dots, F_{p-1}]$. For the specific choice of our parameters, we already computed the 2-cycle above.

* $\bar{\xi}_0 = (\bar{x}, \bar{y}, \bar{u}, \bar{v}), \bar{x} \ne \bar{y}, (\bar{x}, \bar{y}) \ne (\bar{u}, \bar{v}), G_0(\bar{\xi}_0) = \bar{\xi}_1 = (\bar{v}, \bar{u}, \bar{y}, \bar{x})$ and $G_1(\bar{\xi}_1) = \bar{\xi}_0$. This scenario is possible under certain choices of the parameters, and it leads to what we called artificial cycles. Indeed, when $r = 1.0, h_0 = 2.0$, and $h_1 = 1.0$, we obtain the following: The fixed points of $G_{10} = G_1 \circ G_0$ are $\bar{\xi}_0 = (\bar{x}, \bar{y}, \bar{u}, \bar{v})$ and $\bar{\xi}_1 = (\bar{u}, \bar{v}, \bar{x}, \bar{y})$, where

$$\bar{x} \approx 1.109, \bar{y} \approx 3.306, \bar{u} \approx 3.966 \quad \text{and} \quad \bar{v} \approx 2.110.$$

On the other hand, the fixed points of $G_{01} = G_0 \circ G_1$ are $\bar{\eta}_0 = (\bar{v}, \bar{u}, \bar{y}, \bar{x})$ and $\bar{\eta}_1 = (\bar{y}, \bar{x}, \bar{v}, \bar{u})$. Therefore, $[G_0, G_1]$ has two 2-cycles that are not generated by 2-cycles of $[F_0, F_1]$, namely $\{\bar{\xi}_0, \bar{\eta}_0\}$ and $\{\bar{\xi}_1, \bar{\eta}_1\}$. Those cycles lead to what we called artificial cycles of $[F_0, F_1]$, namely $\{(\bar{x}, \bar{y}), (\bar{v}, \bar{u})\}$ and $\{(\bar{u}, \bar{v}), (\bar{y}, \bar{x})\}$. In this case, the embedding technique fails to tackle global stability in the system $[F_0, F_1]$. However, we can eventually squeeze the orbits of $G_{10}$ between $\bar{\xi}_0$ and $\bar{\xi}_1$ (with respect to $\le_{se}$). Similarly, the orbits of $G_{01}$ can be eventually squeezed between $\bar{\eta}_0$ and $\bar{\eta}_1$ (again with respect to $\le_{se}$).
When it comes to the system $[F_0, F_1]$, the even terms of the orbits are eventually squeezes between $\bar{x}$ and $\bar{u}$, while the odd terms of the orbits are eventually squeezed between $\bar{v}$ and $\bar{y}$.

<figure><img src="image_9.png" /><figcaption>Figure 7: Part (i) of this figure shows the curves $\ell_1$ and $\ell_2$ of the first and second equations in System (4.5), respectively. The unique intersection is the 2-cycle $\{\bar{z}_0, \bar{z}_1\}$ of the system $[F_0, F_1]$, where $\bar{z}_0 \approx 2.498$ and $\bar{z}_1 \approx 2.230$. The shaded thin region represents the points $(a, b)$ that satisfy $(a, b, b, a) <_{se} G_1(G_0(a, b, b, a))$ and $a < b$. Part (ii) shows the same curves but when $h_1$ is changed to 1.0. The three intersections between the curves represent the 2-cycle $\{\bar{z}_0 \approx 2.455, \bar{z}_1 \approx 1.950\}$ of the system $[F_0, F_1]$ and two other artificial cycles generated from the intersections (1.109, 3.306) and (3.966, 2.110).</figcaption></figure>

Now, we revisit the system of inequalities in (4.6), and re-write them as

$$ \begin{cases} a(1 - (f(b)))^2 \leq h_0 f(b) + h_1 \\ b(1 - (f(a)))^2 \leq h_0 f(a) + h_1. \end{cases} $$

For $a, b > r$, the factors on the left hand side are positive, which motivate us to define the map $g_2(t) = \frac{h_0 f(t) + h_1}{1 - (f(t))^2}$, $t > r$. Also, recall from Remark 3.1 that we defined the map $g_1(t) = \frac{h_0}{1 - f(t)}$. Both $g_1$ and $g_2$ are decreasing in $t$. The inequalities in (4.6) become $a \leq g_2(b)$ and $b \geq g_2(a)$. It is straightforward to check the following result.

**Lemma 4.6.** Assume $h_0 > r$. There exists an unbounded set of points $(a, b)$, $a < b$ that satisfy

$$ a \leq g_1(b) < g_2(b) < g_2(a) \quad \text{whenever} \quad h_0 < h_1 $$

and

$$ b \geq g_1(a) > g_2(a) > g_2(b) \quad \text{whenever} \quad h_0 > h_1. $$

Finally, we arrive at the main result of this section, which provides a generalization of Part (i) of Theorem 3.6.

**Theorem 4.7.** Consider Eq. (1.4) with $h_0 \neq h_1$. If both $h_0, h_1 > r$ and the solution of System (4.5) is unique, then the 2-cycle $\{\bar{z}_0, \bar{z}_1\}$ is global asymptotically stable.

*Proof.* For each initial condition $(x_0, x_{-1})$ of Eq. (1.4), we need to find a point $(a, b)$ that makes Proposition 2.8 applicable. Without loss of generality, we can consider the initial conditions to be $(x_2, x_1)$. In this case, we guarantee that $x_1 > h_0 > r$ and $x_2 > h_1 > r$. If $h_0 > h_1$, we consider $(a, b)$
such that $a < b, b = g_1(a)$ and $a$ is sufficiently small so that $r < a < h_1$. In this case and based on Lemma 4.6, we obtain

$$b \geq g_1(a) > g_2(a) > g_2(b) > h_1 > a.$$

If $h_0 < h_1$, we consider $(a, b)$ such that $a < b, a = g_1(b)$ and $b$ is sufficiently large so that $g_2(h_0) < b$. In this case and based on Lemma 4.6, we obtain

$$a \leq g_1(b) < g_2(b) < g_2(a) = g_2(g_1(b)) < g_2(h_0) < b.$$

Observe that this choice of $(a, b)$ can be done to have $(x_2, x_1) \in [a, b]^2$. All the hypotheses of Proposition 2.8 are satisfied, and we invoke the Lemma to obtain the result. □

The uniqueness condition in Theorem 4.7 is to guarantee that the global attractor is the 2-cycle of $[F_0, F_1]$. However, if we just require $h_0, h_1 > r$, then the theorem gives us either a globally attracting 2-cycle or a trapping box in which its boundaries are determined by two artificial 2-cycles of $[F_0, F_1]$.

# 5 Conclusion

We investigated two systems in this paper: the Ricker model with a time delay and constant stocking $x_{n+1} = F(x_n, x_{n-1}) = x_n \exp(r - x_{n-1}) + h$, and the system involving periodic stocking $x_{n+1} = F_n(x_n, x_{n-1}) = x_n \exp(r - x_{n-1}) + h_n$. In both cases, we embedded the 2-dimensional system into a 4-dimensional system that is monotonic with respect to a partially ordered set. In the first system, a unique positive equilibrium exists, which we denote by $\bar{y}_h$. The equilibrium $\bar{y}_h$ transitions into a cycle in the second system, and the cycle has the same period as the system. We utilized the embedding technique to analyze global stability in both systems.

In the constant stocking case, we obtained local stability when $\bar{y}_h < 1 + h$, which leads to the condition $0 < r < r_1 = h + 1 - \ln(h + 1)$. On the other hand, we proved global stability under the constraint

$$\bar{y}_h < \frac{1}{2} \left( h + \sqrt{h^2 + 4h} \right),$$

which leads to the constraint

$$0 < r < r_2 = \frac{1}{2} \left( h + \sqrt{h^2 + 4h} \right) + \ln \left( 1 - \frac{2h}{h + \sqrt{h^2 + 4h}} \right).$$

It's worth noting that the global stability condition we established is more stringent than the local stability condition, i.e., $r_2 < r_1$. Therefore, the common assertion that “local stability implies global stability” is still an open problem for all values of $r$ between $r_1$ and $r_2$ as we illustrated in Fig. 5.

For the 2-periodic system $x_{n+1} = F_n \mod 2(x_n, x_{n-1})$, we established conditions that ensure both local and global stability of the existing 2-cycle. We illustrated a scenario in which the 2-cycle undergoes a Neimark-Sacker bifurcation, leading to the emergence of two invariant curves that act as a unified attractor. Finding the explicit form of the cycle is not possible, but we provided illustrative numerical examples that clarify various possible scenarios. As in the constant stocking case, sufficient conditions for the global stability of the 2-cycle is obtained in the 2-periodic case. Nonetheless, identifying necessary and sufficient conditions remains an open problem. We also conjecture here that local stability of the 2-cycle is sufficient to ensure its global stability.
References

[1] P.A.P. Moran. Some remarks on animal population dynamics. *Biometrics*, 6(3):250-258, 1950.

[2] W. E. Ricker. Stock and recruitment. *J. Fish. Res. Board Can.*, 11(5):559-623, 1954.

[3] J. T. Schnute. Curiosity, recruitment, and chaos: A tribute to Bill Ricker's inquiring mind. *Environ Biol Fish*, 75:95-110, 2006.

[4] Robert M. May. Biological populations with nonoverlapping generations: Stable points, stable cycles, and chaos. *Science*, 186(4164):645-647, 1974.

[5] Robert M. May. Simple mathematical models with very complicated dynamics. *Nature*, 261:459-467, 1976.

[6] Robert M. May and George F. Oster. Bifurcations and dynamic complexity in simple ecological models. *American Naturalist*, 110:573-599, 1976.

[7] M. E. Fisher, B. S. Goh, and T. L. Vincent. Some stability conditions for discrete-time single species models. *Bull. Math. Biol.*, 41(6):861-875, 1979.

[8] David Singer. Stable orbits and bifurcation of maps of the interval. *SIAM J. Appl. Math.*, 35(2):260-267, 1978.

[9] Paul Cull. Enveloping implies global stability. In *Difference equations and discrete dynamical systems*, pages 71-85. World Sci. Publ., Hackensack, NJ, 2005.

[10] E.C. Pielou. *Population and Community Ecology*. Gordon and Breach, New York, 1974.

[11] Simon A. Levin and Robert M. May. A note on difference-delay equations. *Theoret. Population Biol.*, 9(2):178-187, 1976.

[12] Fred Brauer and Carlos Castillo-Chavez. *Mathematical models in population biology and epidemiology*, volume 40 of *Texts in Applied Mathematics*. Springer, New York, second edition, 2012.

[13] Mark Kot. *Elements of mathematical ecology*. Cambridge University Press, Cambridge, 2001.

[14] Eduardo Liz, Victor Tkachenko, and Sergei Trof'mchuk. Global stability in discrete population models with delayed-density dependence. *Math. Biosci.*, 199(1):26-37, 2006.

[15] Ferenc A. Bartha, Ábel Garab, and Tibor Krisztin. Local stability implies global stability for the 2-dimensional Ricker map. *J. Difference Equ. Appl.*, 19(12):2043-2078, 2013.

[16] James F. Selgrade. Using stocking or harvesting to reverse period-doubling bifurcations in discrete population models. *J. Differ. Equations Appl.*, 4(2):163-183, 1998.

[17] James F. Selgrade and James H. Roberds. Reversing period-doubling bifurcations in models of population interactions using constant stocking or harvesting. volume 6, pages 207-231. 1998. Geoffrey J. Butler Memorial Conference in Differential Equations and Mathematical Biology (Edmonton, AB, 1996).

[18] Z. AlSharawi and A. M. Amleh. Harvesting and stocking in discrete-time contest competition models with open problems and conjectures. *Palest. J. Math.*, 5(Special Issue):238-249, 2016.
[19] Z. AlSharawi. A global attractor in some discrete contest competition models with delay under the effect of periodic stocking. *Abstr. Appl. Anal.*, pages 7, Art. ID 101649, 2013.

[20] Elena Braverman and Damir Kinzebulatov. On linear perturbations of the Ricker model. *Math. Biosci.*, 202(2):323-339, 2006.

[21] J.-L. Gouzé and K. P. Hadeler. Monotone flows and order intervals. *Nonlinear World*, 1(1):23-34, 1994.

[22] H. L. Smith. The discrete dynamics of monotonically decomposable maps. *J. Math. Biol.*, 53(4):747-758, 2006.

[23] H. L. Smith. Global stability for mixed monotone systems. *J. Difference Equ. Appl.*, 14(10-11):1159-1164, 2008.

[24] Ziyad AlSharawi. Embedding and global stability in periodic 2-dimensional maps of mixed monotonicity. *Chaos Solitons Fractals*, 157:Paper No. 111933, 10, 2022.

[25] Z. AlSharawi, J. Cánovas, and A. Linero. Folding and unfolding in periodic difference equations. *J. Math. Anal. Appl.*, 417(2):643-659, 2014.