Efficient qudit based scheme for photonic quantum computing

Márton Karácsony1$\star$,
László Oroszlány1,2 and
Zoltán Zimborás3,4

1 Department of Physics of Complex Systems, Eötvös Loránd University, 1117 Budapest, Hungary

2 MTA-BME Lendület Topology and Correlation Research Group, Budapest University of Technology and Economics, 1521 Budapest, Hungary

3 Wigner Research Centre for Physics, H-1525, P.O.Box 49, Budapest, Hungary

4 Algorithmiq Ltd, Kanavakatu 3C, FI-00160 Helsinki, Finland

${}^\star$<EMAIL_ADDRESS>
August 27, 2024

§ ABSTRACT

Linear optics is a promising alternative for the realization of quantum computation protocols due to the recent advancements in integrated photonic technology. In this context usually qubit based quantum circuits are considered, however, photonic systems naturally allow also for $d$-ary, i.e., qudit based, algorithms. This work investigates qudits defined by the possible photon number states of a single photon in $d > 2$ optical modes. We demonstrate how to construct locally optimal non-deterministic many-qudit gates using linear optics and photon number resolving detectors, and explore the use of qudit cluster states in the context of a $d$-ary optimization problem. We find that the qudit cluster states require less optical modes and are encoded by a fewer number of entangled photons than the qubit cluster states with similar computational capabilities. We illustrate the benefit of our qudit scheme by applying it to the k-coloring problem.

§ INTRODUCTION

In recent years, many advancements have been made in the area of integrated photonic quantum technologies.
Integrated quantum photonic processors (IQP) with on-chip linear optics <cit.>, single photon sources <cit.> and photon number resolving detectors (PNRD) <cit.> have been demonstrated.
These integrated components and waveguides are usually printed onto silica, which makes it possible to miniaturize the chips and increase the number of on-chip components. Since photons couple very weakly to their environment, IQPs typically do not require millikelvin temperatures to operate as opposed to other quantum information processing architectures like superconducting devices <cit.>.
However, there are still many hurdles to overcome in order to achieve fully scalable and fault-tolerant quantum computation. Until then even IQPs will remain intermediate scaled and thus the amount of available quantum resources will be limited.
A common technique for reducing computational resources is to encode quantum information efficiently. Typically quantum information is encoded into two level systems, qubits which became standard due to their simplicity and analogy to classical computation. However, there are other ways to deal with quantum information. An alternative is to use qudits, $d > 2$ level quantum systems, which are the basis for high-dimensional quantum computation <cit.>.
Qudits have been considered in many areas of quantum computation, for example in topological quantum computation the braiding of $\mathbb{Z}_d$ parafermions <cit.> provide a natural way to implement qudits <cit.>.
There also have been attempts to realize superconducting qudits <cit.>, but because of the flexibility in how photon states can be interpreted, qudits are most prominent in photonics <cit.>. The use of qudits increases the number of dimensions per computation unit. This property of qudits in theory can provide some benefits over qubit systems by reducing circuit complexity of quantum algorithms or by increasing the channel capacity and noise tolerance of communication protocols <cit.>.

There has been considerable experimental progress on the way towards photonic quantum supremacy via boson sampling (see, for example, Refs. <cit.>), however, the realization of universal photonic quantum computation still seems far from reality.

This is because linear optics alone cannot be used to produce certain multi-photon states which are required in some form by all universal photonic schemes. The most well-known universal photonic scheme is the KLM scheme <cit.> which solves the problem of creating entanglement between photonic qubits by applying postselection based on ancilla measurements. This makes it possible to prepare entangled two qubit states non-deterministically using only linear optics and PNRDs. This measurement based method combined with the quantum teleportation of qubits was shown to be able to create near-deterministic two qubit gates. However, the complete implementation of the KLM scheme proved to be out of reach with current technology, and only some of the its techniques has been demonstrated experimentally. Another approach for realizing universal quantum computation is to use measurement based quantum computation (MBQC). In MBQC the quantum algorithms are performed by making single qubit measurements on a large entangled state, usually a cluster state <cit.>. For photonics MBQC provides a resource efficient and deterministic way to implement universal quantum computation.
Since using qudits instead of qubits are often beneficial in quantum computation, as discussed above, it is natural to ask whether photonic MBQC can also be improved by utilizing high-dimensional computational units.
In this work, we show that qudit cluster states can indeed be more resource efficient than qubit cluster states.

The present paper is structured as follows: In Sec. <ref> we go over a high-dimensional generalization of the KLM encoding and explain how to use passive linear optics to create single- and many-qudit gates. We generalize the post selection based approach of the KLM scheme to implement non-deterministic many-qudit gates. Then we use an optimization method to find locally optimal interferometer configurations for some two-qudit gates. In Sec. <ref> we introduce the theory of high-dimensional clusters states, and then in Sec. <ref> we show how to use high-dimensional cluster states to solve the $k$-coloring problem.

§ MULTI-RAIL QUDIT ENCODING

In this section, we describe a simple method for encoding d-level quantum system, i.e., qudits, into optical modes and show how to implement non-deterministic many-qudit gates using linear optics and PNRDs. Before  describing the qudit encoding, we present a few basic qudit gates and discuss some useful identities between them. Then, after the specific description of the multi-rail encoding, we also present an optimization method for increasing the success rate of non-deterministic photonic many-qudit gates.

§.§ Qudit gates

There are several different conventions in use for the generalization of the well known Pauli $X$ and $Z$ gates for qudit systems. We will use the following definitions

\begin{align}
X = \sum_{n = 0}^{d - 1} \ket{n \oplus 1} \bra{n}, && Z = \sum_{n = 0}^{d - 1} \omega^n \ket{n} \bra{n},
\end{align}

where $\oplus$ means addition modulo $d$ and $\omega = \exp\left(i 2\pi / d\right)$. One can immediately see that $X$ and $Z$ are unitary and  $X^d = Z^d = \mathbbm{1}$. However, it is important to note that they are not Hermitian for $d > 2$. Since the qudit gates $X$ and $Z$ defined above have the same eigenvalues, they satisfy a similarity relation,

\begin{equation}
\label{eq:XZsimilarity}
H X H^\dagger = Z,
\end{equation}

where the corresponding basis change operator defines the qudit Hadamard gate

\begin{equation}
H = \frac{1}{\sqrt{d}} \sum_{n = 0}^{d - 1} \sum_{m = 0}^{d - 1} \omega^{nm} \ket{n} \bra{m}.
\end{equation}

The Hadamard gate is often used to prepare the qudit state

\begin{equation}
\ket{+} \equiv H \ket{0} = \frac{1}{\sqrt{d}} \sum_{n = 0}^{d - 1} \ket{n},
\end{equation}

which plays a key role in many quantum algorithms.
The controlled $X$ and $Z$ gates are defined as

\begin{align}
\label{eq:cxcz}
CX = \sum_{n = 0}^{d - 1}  \ket{n} \bra{n} \otimes X^n, && CZ = \sum_{n = 0}^{d - 1}  \ket{n} \bra{n} \otimes Z^n.
\end{align}

From Eq. (<ref>) it follows that

\begin{equation}
CX = \left(\mathbbm{1} \otimes H^\dagger\right) CZ \left(\mathbbm{1} \otimes H\right).
\end{equation}

Equation (<ref>) means that $CX \ket{n} \otimes \ket{m} = \ket{n} \otimes \ket{n \oplus m}$, thus one can write

\begin{align}
CX^\dagger (\mathbbm{1} \otimes Z) CX \ket{n} \otimes \ket{m} &= CX^\dagger (\mathbbm{1} \otimes Z) \ket{n} \otimes \ket{n \oplus m}  \nonumber \\
&= CX^\dagger \omega^{n + m} \ket{n} \otimes \ket{n \oplus m} \nonumber \\
&= \omega^{n + m} \ket{n} \otimes \ket{m} = Z \otimes Z \ket{n} \otimes \ket{m},
\end{align}

where $n, m \in \{0, 1, 2, \ldots d - 1\}$, therefore

\begin{equation}
\label{eq:pZpZ}
CX^\dagger (\mathbbm{1} \otimes Z^m) CX = Z^m \otimes Z^m,
\end{equation}

and by a very similar derivation we obtain

\begin{equation}
\label{eq:pZmZ}
CX (\mathbbm{1} \otimes Z^{-m}) CX^\dagger = Z^m \otimes Z^{-m}.
\end{equation}

Equations (<ref>) and (<ref>) will prove to be useful later when one needs to decompose unitaries of the form $\exp\left(i\alpha Z \otimes Z\right)$. Another useful qudit gate is the one-level controlled-Z gate  $\overline{\mathrm{CZ}}$, which is defined by its action on the basis states $\ket{k} \otimes \ket{l}$ as

\begin{equation} \label{eq:single-level-cz}
\overline{\mathrm{CZ}} \ket{k} \otimes \ket{l} = \exp \left(i 2\pi l \delta_{k, d-1}/d\right) \ket{k} \otimes \ket{l},
\end{equation}

where $\delta_{k, l}$ denotes the Kronecker delta. In other words, the gate $\overline{\mathrm{CZ}}$ is another  generalization of the two-qubit controlled sign flip operation, it applies a phase factor $\exp(i 2\pi l /d)$ on the target (second) qudit if the control (first) qudit is in the state $\ket{d-1}$.

§.§ Optical multi-rail qudit protocol

One of the simplest way to encode a qudit in an optical system is to use $d$ different optical modes, labeled by the numbers $0, 1, \ldots, d-1$,

and map the qudit levels to the occupation of different modes via the mapping $\ket{i} \rightarrow \ket{0}_{0} \ket{0}_{1} \cdots \ket{1}_{i} \cdots \ket{0}_{d-1}$, where $\ket{i}$
denotes the $d$ different qudit levels. This encoding of qudit levels is a multi-rail encoding and has the advantage, that any single-qudit $U(d)$ transformation can be efficiently implemented using the Clements decomposition <cit.>, for illustration see Fig. <ref>.

Interferometer configuration for the implementation of a single qudit gate when $d = 5$. This configuration can implement any single qudit gate and minimizes the optical depth. Each crossing represents a general beamsplitter between neighboring optical modes. In general the configuration requires $d (d - 1) / 2$ number of beamsplitters. On the left the input state $\ket{i}$ encoded by single photon enters the interferometer and the output state $U \ket{i}$ appears on the right.

When encoding many qudit systems, one thus needs
$d$ optical modes for each qudit, and all the one-qubit operations can be performed by linear optical gates. To implement two-qudit gates, one can generalize the KLM construction, where a nonlinear sign flip operation is used to obtain non-deterministic two-qubit gates.

Moreover, in general not only the sign flip but any phase shift can be implemented this way. The nonlinear phase shift acts on a state of an optical mode that is in a linear combination of the vacuum, one- and two-photon states in the following way $NS_{\varphi} : \alpha_0 \ket{0} + \alpha_1 \ket{1} + \alpha_2 \ket{2} \rightarrow \alpha_0 \ket{0} + \alpha_1 \ket{1} + e^{i \varphi} \alpha_2 \ket{2}$.
It can be implemented in a probabilistic way with linear optics and PNRDs utilizing two ancilla modes.  The ancilla modes are prepared in the state $\ket{10}_{23}$, so that the initial state of the system looks like $\left(\alpha_0 \ket{0}_1 + \alpha_1 \ket{1}_1 + \alpha_2 \ket{2}_1 \right) \otimes \ket{10}_{23}$. After this preparation of the ancillas, a specific interferometer configuration $U$ is applied on the three modes and then finally the ancillas are measured. If the result of the measurement is $\ket{10}_{23}$ then the operation was successful otherwise it failed.

There are infinitely many choices for $U$ such that one obtains on the first mode the state  $ \alpha_0 \ket{0}_1 + \alpha_1 \ket{1}_1 + e^{i \varphi} \alpha_2 \ket{2}_1$ when measuring $\ket{10}_{23}$ on the ancilla modes, however, the optimal solution in terms of the probability of success is unique up to some phase factors of the rows and columns of the matrix. The set of solutions for $U$ is the following

\begin{align}
U_{11} &= 1 - \sqrt{1 - e^{i\varphi}}, && |U_{12}| \in \left(0, \sqrt{1 - |U_{11}|^2}\right), && U_{13} = \sqrt{1 - |U_{11}|^2 - |U_{12}|^2}, \nonumber \\
U_{21} &= \frac{\sqrt{P} (1 - U_{11})}{U_{12}}, &&  \phantom{|} U_{22} = \sqrt{P}, && U_{23} = \sqrt{P}\left(\frac{|U_{11}|^2 -U_{11}^* - |U_{12}|^2}{U_{12} U_{13}}\right), \nonumber\\
U_{3i} &= \sum_{j,k} \epsilon_{ijk} U_{1j}^* U_{2k}^*,
\end{align}

where $\varphi$ is the desired phase shift and

\begin{equation} \label{eq:nonlinearprob}
P = \left[\frac{\left|\left(\frac{U_{11}^*(1 - U_{11})}{|U_{12}|^2} + 1\right)U_{12}\right|^2}{1 - |U_{11}|^2 - |U_{12}|^2} + \frac{|1 - U_{11}|^2}{|U_{12}|^2} + 1\right]^{-1}
\end{equation}

is the probability of success. The optimal solution can be found by simply finding the maximum of $P$ as a function of $|U_{12}|$.

Since all single-qudit gates are already established, it is enough to implement the two-qudit transformation $\overline{CZ}$ defined in Eq. (<ref>) to create a universal gate set <cit.>. The qudit gate $\overline{\mathrm{CZ}}$ can be constructed using beamsplitters and nonlinear phase shift operations, this is shown by Fig. <ref>. As can be seen the $\overline{CZ}$ gate is made up by $2(d {-} 1)$ number of non-linear phase shifts and the same number of beamsplitters. The nonlinear phase shifts have a success rate less than $1/4$ which is the global maximum of Eq. (<ref>), thus the probability of the $\overline{CZ}$ gate is bound by $P \leq 1/16^{d-1}$.
This construction proves that a universal qudit protocol is possible
in the same way as in the teleportation-based KLM scheme.
However, the proposed $\overline{\mathrm{CZ}}$ gate has a very low success rate, which implies that it is an unfeasible option for actual applications (especially when $d$ is high). This means that one needs an alternative method, which we will present in the next subsection.

Realization of the $\overline{\mathrm{CZ}}$ gate entangling two qudits defined by the multi-rail encoding. The red and blue boxes indicate symmetric beamsplitters with angles $\theta = \pm \pi/4$ and $\phi = 0$ (red is plus and blue is minus), and the white boxes denote nonlinear phase shifts with angles $\varphi_k = 2\pi k/d$.

§.§ Locally optimal many-qudit gates

We now present, as an alternative to the scheme of Sec. <ref>, an optimized method for obtaining photonic many-qudit gates, e.g., $\mathrm{CZ}$ or $\overline{\mathrm{CZ}}$ operations, with a considerably improved success rate. Similarly to the scheme presented in the previous subsection, the desired gate is executed by implementing a non-trivial unitary transformation and doing postselection on the measurement results, as shown by Fig. <ref>. The parameters of the unitary are optimized to enhance the success rate. Such a scheme was already employed for generic two-qubits gates   successfully <cit.>. Now we apply the method for generic many-qudit gates.

Configuration of the interferometer for a general nonlinear operation. The many-qudit state is encoded in the modes at the top and the ancilla modes are placed below them. The ancilla modes are measured using PNRDs and the measurement result is used to determine whether the operation was successful or not. In the event of success, the computational output is encoded in the many-qudit state $\ket{\psi_\text{out}}_c$.

In the multi-rail encoding the many-qudit basis states $\ket{i}$ are mapped to photon number states $\ket{n_1, n_2, \ldots, n_M} \equiv \ket{\mathbf{n}}$ with fixed photon number $N_c$ matching the number of qudits and a fixed number of optical modes $M = N_c d$.

The general problem we aim to solve is how to implement a many-qudit gate $T$ using linear optics and PNRDs. Given some computational input state $\ket{\psi_\text{in}}_c = \sum_\mathbf{n} \alpha_\mathbf{n} \ket{\mathbf{n}}$ encoded into the $M$ computational optical modes, the many-qudit gate $T$ transforms it into

\begin{equation}
\ket{\psi_\text{out}}_c = T \ket{\psi_\text{in}}_c = \sum_{\mathbf{m}} \sum_{\mathbf{n}} T_{\mathbf{m}\mathbf{n}} \alpha_\mathbf{n} \ket{\mathbf{m}},
\end{equation}

where $\mathbf{n}$ and $\mathbf{m}$ indices are referring to the photon number states $\ket{\mathbf{n}} = \ket{n_1, n_2, \ldots, n_M}$ and $\ket{\mathbf{m}} = \ket{m_1, m_2, \ldots, m_M}$ with $\sum_i n_i = \sum_i m_i = N_c$. Unfortunately for a general transformation $T$ there exists no linear interferometer which could execute such a transformation, i.e., most of the time some nonlinear effects are required. Here we will provide nonlinearities via photon number measurements and postselection.

A linear optical interferometer can be described by a unitary $U$ which transforms the creation operators via ${a'_i}^\dagger = \sum_{j} u_{ij} a^\dagger_j$ and the photon number states via

\begin{equation}
U\ket{\mathbf{n}} = \prod_i \left(\sum_j u_{ij} a^\dagger_j\right)^{n_i} \ket{0}.
\end{equation}

Since for most transformations $T$ there exist no $U$ such that $U\ket{\psi_\text{in}}_c = \ket{\psi_\text{out}}_c$, the next best thing we can achieve is

\begin{equation}
U \ket{\psi_\text{in}}_c \ket{1}^{\otimes N_a} \ket{0}^{\otimes N_v} = e^{i\phi}\sqrt{P} \ket{\psi_\text{out}}_c  \ket{1}^{\otimes N_a} \ket{0}^{\otimes N_v} + \text{terms with other ancilla values},
\end{equation}

where we introduced $N_a$ number of ancilla modes each occupied by a single photon and $N_v$ number of vacuum ancilla modes.
If we were to prepare the state $U \ket{\psi_\text{in}}_c \ket{1}^{\otimes N_a} \ket{0}^{\otimes N_v}$, we would find the ancillas in the state $\ket{1}^{\otimes N_a} \ket{0}^{\otimes N_v}$ with probability $P$. This means that by preparing $U \ket{\psi_\text{in}}_c \ket{1}^{\otimes N_a} \ket{0}^{\otimes N_v}$ and measuring the ancilla modes $T$ is implemented with probability $P$. The relation between $U$ and $T$ is

\begin{equation}
\label{eq:optmization}
\bra{\tilde{\mathbf{m}}} U \ket{\tilde{\mathbf{n}}} = \tilde{U}_{\mathbf{m}\mathbf{n}} = e^{i\phi} \sqrt{P} T_{\mathbf{m}\mathbf{n}},
\end{equation}

where $\ket{\tilde{\mathbf{n}}} = \ket{\mathbf{n}}  \ket{1}^{\otimes N_a} \ket{0}^{\otimes N_v}$. The Fock space matrix elements $\tilde{U}_{\mathbf{m}\mathbf{n}}$ can be expressed in terms of $u_{ij}$ using permanents<cit.>. Equation (<ref>) is a nonlinear polynomial equation for $u_{ij}$ and it can be solved for $u_{ij}$ using numerical optimization methods while maximizing the success rate $P$. Since $u$ is unitary it also has to satisfy $u u^\dagger = \mathbb{1}$. To solve this problem via optimization we can define the cost function

\begin{equation}
L = F + \lambda P + \sigma C,
\end{equation}

\begin{align}
F = \frac{\left| \mathrm{Tr}\left(\tilde{U}^\dagger T \right)\right|^2}{ \mathrm{Tr}\left(\tilde{U}^\dagger \tilde{U}\right) \mathrm{Tr} \left(T^\dagger T\right)}, && P = \frac{ \mathrm{Tr}\left(\tilde{U}^\dagger \tilde{U}\right)}{\mathrm{Tr} \left(T^\dagger T\right)}, && C = -\mathrm{Tr}\left\{\left(u u^\dagger - \mathbb{1}\right)^2\right\},
\end{align}

and $\lambda$ and $\sigma$ are positive real constants. The fidelity $F$ reaches its maximum if and only if Eq. (<ref>) is satisfied for some probability $P$, this is guaranteed by the Cauchy–Schwarz inequality. The unitary constraint $C$ reaches its maximum if and only if $U$ is a unitary matrix.

The optimization problem defined by $L$ was solved using a trust-region method <cit.>. The optimization process went as follows. First the problem was solved using $\lambda = 0$. Once a solution with $F = 1$ and $C = 0$ was found $\lambda$ was gradually increased to the largest possible value for which the solution still converged to $F = 1$.
At the beginning $N_a$ and $N_v$ were chosen such that the system of equations in Eq. (<ref>) would be under-defined. After finding a solution $N_a$ and $N_v$ would be decreased until no solution with $F = 1$ and $C = 0$ could be found. This process was carried out for the qutrit $\overline{CZ}$ and $CZ$ gates (see Table <ref>).
There is a significant improvement in the success rate when the gates are optimized. These optimized gates could be used to non-deterministically prepare entangled photon states which are necessary for many measurement based architectures like optical cluster state computing <cit.>. In the following sections we will look at how high-dimensional cluster states might be used in photonic qudit schemes in general and in particular for the $k$-coloring problem.

$T$                                   $P_\mathrm{naive}$       $P$                $N_a$     $N_v$
CZ (qubit) <cit.>     0.0625     0.0740 …    2         0

$\overline{\mathrm{CZ}}$ (qutrit)     $\approx 0.0016$         0.011              3         3

CZ (qutrit)                           $\approx 0.00000256$     0.000507           5         4

Optimized two-qutrit gates. $P_\text{naive}$ is the probability we get when the gates are decomposed into many nonlinear phase shifts. $P$ is the best probability found after doing the optimization process for many different randomly selected initial conditions.

§ HIGH-DIMENSIONAL CLUSTER STATES

The most common computational model is the quantum circuit model which resembles the logic gate based description of classical computation.
Measurement based quantum computation is an alternative to the circuit model, where instead of applying unitary gates one performs adaptive single qubit or qudit measurements on a large entangled state, called the cluster state.

At first this may seem completely different from the circuit model
but their equivalence have been proved

<cit.>. MBQC is in particular a promising paradigm for the realization of universal optical quantum computation, because if cluster states can be prepared efficiently there is no longer need for further non-deterministic gates to perform the computation. In other words, MBQC shifts the problem of implementing nonlinearities to the problem of state preparation.

Several methods have been developed
to generate cluster states efficiently,

MBQC is
one of the most viable options for realizing universal quantum computation using photonics <cit.>.

Cluster states have the form

\begin{equation}
\prod_{\{i, j\} \in C}CZ_{i,j} \ket{+}^{\otimes N},
\end{equation}

where $CZ$ is the controlled qudit $Z$ gate and $C$ is a graph whose edges represent the $CZ$ gates between qudits. The idea is to perform single-qudit measurements sequentially on the qudits found in the cluster, where the basis of each measurement depends on the results of the previous measurements. When a qudit is measured it is effectively removed from the cluster reducing its size. After the measurements, the remains of the cluster state will encode the computational state up to some known single-qudit Pauli errors which can be easily corrected for at the end of the computation.

To connect the cluster state picture with the quantum circuit model, we can use the teleportation identity

\begin{equation}
\label{eq:dittele}
\begin{quantikz}
\lstick{$\ket{\psi}$} & \ctrl{1} & \gate{P(\boldsymbol{\theta})} & \gate{H^\dagger} & \meter{} & = m\\
\lstick{$\ket{+}$} & \gate{Z}  & \qw & \qw  & \qw \rstick{$ X^m H P(\boldsymbol{\theta}) \ket{\psi}$}
\end{quantikz},
\end{equation}

where $\ket{\psi}$ is a general qudit state, $m$ is the result of the measurement performed on the top qudit, and $P(\boldsymbol{\theta})$ is a phase gate with $\boldsymbol{\theta} \in \mathbb{R}^d$ which acts like $P(\boldsymbol{\theta}) \ket{n} = \exp(i\theta_n) \ket{n}$ <cit.>. This circuit represents a single measurement step in the cluster state picture. It applies the single-qudit unitary $H P(\boldsymbol{\theta})$ and propagates the state to the qudit that connects to it. Any unitary can be decomposed into a product that only contains gates of the form $HP(\boldsymbol{\theta})$, thus chaining together many of these circuits one can perform any single-qudit unitary <cit.>. These chains directly correspond to measurements on a linear cluster state. An important detail is that after each measurement we introduce Pauli errors of the form $X^a Z^b$ where $a, b$ depend on the results of the previous measurements. In order to account for the accumulated Pauli errors we have to change the measurement angles to $\boldsymbol{\theta}'$ such that $HP(\boldsymbol{\theta}') X^aZ^b = X^{a'}Z^{b'} HP(\boldsymbol{\theta})$.

In general, we can translate any quantum circuit to a cluster state with measurement instructions using the  teleportation identity defined above. This way each qudit in the quantum circuit model will correspond to a linear subcluster made up of several physical qudits and the $CZ$ gates between qudits in the quantum circuit will become connections across the corresponding linear clusters. Although this construction is enough to translate any quantum circuit, it makes the translation easier if we allow the use of $CZ^\dagger$ connections between physical qudits in the cluster. This is justified by another teleportation identity slightly different from (<ref>)

\begin{equation}
\label{eq:dittele2}
\begin{quantikz}
\lstick{$\ket{\psi}$} & \ctrl{1} & \gate{P(\boldsymbol{\theta})} & \gate{H} & \meter{} & = m\\
\lstick{$\ket{+}$} & \gate{Z^\dagger}  & \qw & \qw  & \qw \rstick{$ X^m H^\dagger P(\boldsymbol{\theta}) \ket{\psi}$}
\end{quantikz}.
\end{equation}

These cluster states can be illustrated by depicting the linear clusters as rows with the $CZ$ or $CZ^\dagger$ gates between them shown as vertical lines. An example of this is shown by Fig. <ref>.

Qudit cluster state with the measurement instructions to implement the two qudit unitary $U = CX P(\boldsymbol{\theta}) CX^\dagger$.
The black and red edges denote the $CZ$ and $CZ^\dagger$ connections between the physical qudits. The initial state with its Pauli errors is encoded in the green qudits. The measurement order is indicated by the numbers. After the execution of all the measurements, the final state is encoded in the red qubits.
The arrows show how the qudit states are teleported from one qudit to the next, and the unitary of each arrow indicates what transformation should the teleportation perform.

§ THE K-COLORING PROBLEM

Graph coloring is one of the most studied topics of graph theory because of its applications in solving scheduling problems and compiler theory <cit.>.
Coloring refers to the assignment of labels (colors) from a label set to the vertices of a graph $G$. This assignment defines a map $C: V \rightarrow L$, where $V$ is the set containing the vertices of $G$ and $L$ is the label set. Arbitrarily ordering the vertices of $G$ allows us to represent a coloring as string of labels: If we call the $i$-th vertex $v_i \in V$ and its assigned color $C(v_i) = l_i \in L$, then the string $\mathbf{l} = l_1 l_2 l_3 \ldots l_{|V|}$ completely characterizes the coloring $C$. The $k$-coloring problem of a graph is to decide whether the vertices of the graph can be colored using $k$ colors, such that no adjacent vertices have the same color. The smallest $k$ for which the $k$-coloring is possible for a given graph is called the chromatic number of the graph. The $k$-coloring problem of an arbitrary graph is known to be NP-complete <cit.>. There has been rigorous proof showing the speed up of certain quantum algorithms compared to the best known classical, <cit.>, and even certain heuristic variational quantum algorithms for a similar problem have been  hinted to outperform classical methods <cit.>.

§.§ Qubit algorithm

The $k$-coloring problem can be reformulated as a combinatorial optimization problem where the goal is to minimize the number of edges connecting vertices with identical colors by searching through the possible $|V|$ length  strings of the label set $L = \{0, 1, 2, \ldots, k - 1\}$. On qubit-based architectures, binary optimization problems can be solved using the quantum approximate optimization algorithm (QAOA) <cit.>. Encoding the strings of $L$ as bit-strings means that we can use the QAOA to solve the $k$-coloring problem as well;
such a QAOA type of approach for the $k$-coloring (and the related $k$-SAT problem) has been studied in Refs. <cit.>.

The encoding of colorings into bit-strings is most easily done via one-hot encoding, this means that the strings of $L$ are encoded into $N = |V| \cdot k$ bits, where each set of $k$ consecutive bits encodes the color of a vertex. The color of the $i$-th vertex $l_i$ is indicated by the $i$-th set of $k$ bits of the form $00\ldots010\ldots00$, where the $l_i$-th bit is flipped and the rest of the $k$ bits is zero. We will denote the $j$-th bit of the $i$-th set of $k$ bits by $x_{i, j}$.

Minimizing the cost function $f : {\{0, 1\}}^{N} \rightarrow \mathbb{R}$

\begin{equation}
\label{eq:binary_cost}
f(\mathbf{x}) = C \sum_{n = 1}^{|V|} \left(1 - \sum_{i = 0}^{k-1} x_{n, i}\right)^2 + D \sum_{n = 1}^{|V|} \sum_{m = 1}^{|V|} \sum_{i = 0}^{k-1} A_{nm} x_{n, i} x_{m,i},
\end{equation}

is equivalent to finding a solution to the $k$-coloring problem,
where $\mathbf{x}$ is a bit-string of length $N$, $A_{nm}$ are the matrix elements of the adjacency matrix of $G$ and $C$ and $D$ are arbitrary positive real constants. The first term in Eq. (<ref>) weighted by $C$ is a penalty term, since not every bit-string of length $N$ corresponds to a coloring. For every $n \in \{1, 2, \ldots, |V|\}$ exactly one of the $k$ bits $x_{n, j}$, where $j \in \{0, 1, \ldots, k - 1\}$, has to equal one in order for $\mathbf{x}$ to encode a coloring. The remaining term simply counts the number of edges that connect vertices with the same color. To construct the QAOA cost Hamiltonian we can use the recipe

\begin{equation}
\label{eq:cost_ham}
H_C = \sum_{\mathbf{x} \in {\{0, 1\}}^N} f(\mathbf{x}) \ket{\mathbf{x}}\bra{\mathbf{x}},
\end{equation}

where $\ket{\mathbf{x}}$ is the computational state

\begin{equation*}
\bigotimes_{n = 1}^{|V|} \left(\bigotimes_{i = 0}^{k - 1} \ket{x_{n, i}}\right).
\end{equation*}

Substituting Eq. (<ref>) into Eq. (<ref>), the cost Hamiltonian becomes

\begin{equation} \label{eq:cost_ham2}
H_C = \frac{C}{4} \sum_{n = 1}^{|V|} \left(2\mathbbm{1} - \sum_{i = 0}^{k - 1} \left(\mathbbm{1} - Z_{n, i}\right)\right)^2 + \frac{D}{4} \sum_{n = 1}^{|V|} \sum_{m = 1}^{|V|} \sum_{i = 0}^{k - 1} A_{nm}\left(\mathbbm{1} - Z_{n, i}\right) \left(\mathbbm{1} - Z_{m, i}\right),
\end{equation}

where $Z_{n, i}$ is the Pauli Z matrix acting on the $(k \cdot (n - 1) + i)$-th qubit labeled by the tuple $(n, i)$. Obtaining the ground state of $H_C$ is the same as finding a bit-string that minimizes the cost function $f$ <cit.>.

QAOA minimizes the expectation value $\bra{\psi\left(\boldsymbol{\alpha}, \boldsymbol{\beta}\right)} H_C \ket{\psi\left(\boldsymbol{\alpha}, \boldsymbol{\beta}\right)}$, where

\begin{equation}
\ket{\psi\left(\boldsymbol{\alpha}, \boldsymbol{\beta}\right)} = \left(\prod_{n = 1}^p e^{i \alpha_n H_M} e^{i \beta_n H_C}\right) \ket{+}^{\otimes N},
\end{equation}

$\boldsymbol{\alpha}, \boldsymbol{\beta} \in \mathbb{R}^p$ and

\begin{equation*}
H_M = \sum_{n = 1}^{|V|} \sum_{i = 0}^{k - 1} X_{n, i}
\end{equation*}

is the mixing Hamiltonian, where $p$ denotes the number of QAOA layers.
Preparation of $\ket{\psi\left(\boldsymbol{\alpha}, \boldsymbol{\beta}\right)}$ requires the implementation of two parameterized unitary $U_M(\alpha) = e^{i\alpha H_M}$ and $U_C(\beta) = e^{i\beta H_C}$. $U_M$ can be realized using single qubit rotations and only $U_C$ requires entangling gates. $U_C$ can be further separated into the product of two unitaries, where one of the unitary requires only single qubit phase gates and the other unitary requires only controlled rotations. This can be achieved by separating the cost Hamiltonian into a non-interacting and an interacting part, $H_C = H_0 + H_1$, where

\begin{equation}
\label{eq:cost_int}
H_1 = \frac{C}{2} \sum_{n = 1}^{|V|} \sum_{i <j = 0}^{k - 1} Z_{n, i} Z_{n, j} + \frac{D}{2} \sum_{\{n, m\} \in E} \sum_{i = 0}^{k - 1} Z_{n, i} Z_{m, i}.
\end{equation}

In Eq. (<ref>) $E$ denotes the set containing the edges of the graph. The Hamiltonians $H_0$ and $H_1$ commute since they only contain Pauli Z matrices, thus

\begin{equation}
U_C(\alpha) = e^{i \alpha H_0} e^{i \alpha H_1}.
\end{equation}

The unitary $e^{i\alpha H_1}$ can be decomposed into the product of $|V| \binom{k}{2} + k |E|$ number of controlled rotations

\begin{equation}
\label{eq:qubit_decomp}
e^{i\alpha H_1} = \prod_{n = 1}^{|V|} \left(\prod_{i < j = 1}^k CX_{n,i;n,j} P_{n,j}\left(\alpha C/2\right) CX_{n,i;n,j}\right)
\times \prod_{\{n, m\} \in E}\left(\prod_{i = 0}^{k - 1} CX_{n, i; m, i} P_{m,i}\left(\alpha D / 2\right)CX_{n, i; m, i}\right),
\end{equation}

where $P_{n, j}(\alpha) = e^{i\alpha Z_{n, j}}$ is the phase gate acting on the qubit labeled by the pair of indices $(n, j)$ and $CX_{n, i;n,j}$ is a controlled NOT gate between the control qubit $(n, i)$ and the target qubit $(n, j)$. Eq. (<ref>) can be derived using the identity $e^{i\alpha Z\otimes Z} = CX \left(\mathbbm{1} \otimes P(\alpha)\right) CX$.

The physical realization of the unitary $e^{i\alpha H_1}$ requires by far the most amount of physical resources compared to the other parts of the algorithm, since only this part  needs entangling gates. Therefore, focusing only on the implementation details of $e^{i\alpha H_1}$ can give us a good lower bound on the resources necessary to carry out the QAOA.

§.§ Graph coloring using qudits

When using qudits with dimension $d = k$, there is a one-to-one correspondence between the computational basis and the possible graph colorings with $k$ colors via the mapping

\begin{equation}
l_1 l_2 l_3 \ldots l_{|V|} \rightarrow \ket{l_1} \otimes \ket{l_2} \otimes \ket{l_3} \otimes \cdots \otimes  \ket{l_{|V|}},
\end{equation}

where $l_i \in \{0, 1, 2, \ldots, k - 1\}$. This means that the coloring problem can be formulated as an unconstrained $d$-ary optimization problem which can be solved by the QAOA generalized to qudits <cit.>. Because the optimization problem is unconstrained, the decomposition of the QAOA layers will turn out to be simpler due to lack of penalty terms.

The QAOA for qudits is very similar to the qubit version only the form of the mixing and cost Hamiltonian is different. The mixing Hamiltonian is still a simple sum of single-qudit terms, e.g., the $r$-nearby-values single-qudit mixer $H_M = \sum_{n = 1}^{|V|}\sum_{i = 1}^r \left(X^i_n + {X^\dagger_n}^i\right)$ as in Ref. <cit.>, and the cost Hamiltonian is

\begin{equation}
\label{eq:qudit_cost_ham}
H_C = \sum_{n = 1}^{|V|} \sum_{m = 1}^{|V|} \sum_{i = 0}^{k - 1} A_{nm} Z_n^i Z_m^{k - i} = \sum_{\left\{n, m\right\} \in E} \sum_{i = 0}^{k-1} Z_n^i Z_m^{k - i},
\end{equation}

where $Z_n$ is the qudit Pauli $Z$ gate acting on the $n$-th qudit and $A_{nm}$ are the matrix elements of the adjacency matrix of the graph $G$. The ground states of the Hamiltonian $H_C$ correspond to optimal $k$-colorings of $G$. Similarly to the qubit case, we can decompose the unitary $U_C(\alpha) = e^{i\alpha H_C}$ into $|E|$ number of controlled qudit rotations

\begin{equation}
\label{eq:qudit_decomp}
e^{i\alpha H_C} = \prod_{\left\{n, m\right\} \in E} CX_{n;m} P_m(\alpha) CX^\dagger_{n;m},
\end{equation}

where $P_m(\alpha)$ is a single-qudit phase gate which acts on the $m$-th qudit the following way

\begin{equation}
P_m(\alpha) \ket{\mathbf{x}} = e^{i \alpha \sum_{n = 0}^{k-1} Z^n_m} \ket{\mathbf{x}} =
\begin{cases}
e^{i\alpha k} \ket{\mathbf{x}} & \text{if} \ x_m = 0,\\
\ket{\mathbf{x}} & \text{otherwise.}
\end{cases}
\end{equation}

§.§ Cluster states for the qubit and qudit implementations

Having established all the necessary techniques to implement the QAOA solution of the $k$-coloring problem both using qubits and qudits, we can now compare the physical resources used by the two methods. A single layer of the qudit QAOA algorithm without the mixing part can be executed on the qudit cluster state $C_d$ constructed from many copies of $R_d$, where $R_d$ is a cluster with the topology depicted by Fig. <ref> made up by $d$-dimensional qudits. When the number of colors is $k$, one has to use $d{=}k$-dimensional qudits for the  algorithm. The construction of $C_d$ is shown by Fig. <ref>. The construction of $C_2$, the qubit equivalent of $C_d$, works differently. In the qubit case, one has to add $k$ copies of $R_2$ for each edge of $G$ and also $k(k - 1)/2$ copies for each node of $G$. These construction rules follow from the decompositions in Eqs. (<ref>) and (<ref>).

The number of qudits in $C_d$ is approximately

\begin{equation}
\label{eq:cluster_size_approx_qdit}
|C_d| \approx |R_d| \cdot |E| = 8 |E|,
\end{equation}

where $|E|$ is the number of edges in $G$, which is equal to the number of copies of $R_d$ used to build $C_d$, and $|R_d|$ is the number of qudits in $R_d$. The approximation can  be made exact by carefully considering the precise sequence of node additions and removals in the cluster growing process, however, for large enough graphs Eq. (<ref>) provides a highly accurate estimate.
Similarly, the number of qubits in $C_2$ is

\begin{equation}
|C_2| \approx |R_2| \cdot \left(\binom{k}{2}|V| + k |E|\right) = 8\binom{k}{2} |V| + 8k |E|.
\end{equation}

In general, the ratio $|C_2| / |C_d|$ is strictly larger than $k$, and in the large dense graph limit, when $|E| \gg |V|$, this ratio approaches $k$. This is illustrated by Fig. <ref> for Erdős-Rényi random graphs.  If one uses the multi-rail encoding described in Sec. <ref>, the photon numbers of the cluster states is equal to their size. The number of optical modes required to encode the cluster state $C_d$ is $d |C_d|$ including the qubit case $d = 2$. Thus, the number of photons decreases by a $k$-fold and the number of optical modes halves, when using multi-rail qudits instead of KLM qubits to encode the $k$-coloring problem.

Qudit cluster state for a single layer of the qudit QAOA. Each edge of $G$ corresponds to a controlled qudit rotation in the decomposition of $\exp(i\alpha H_C)$ described by Eq. (<ref>). Thus one can chain many copies of  $R_d$ to form the complete qudit cluster. The two input nodes of $R_d$ are merged with two output nodes of the unfinished cluster state based on the topology of $G$. Each row of physical qudits represents a logical qudit labeled by the numbers. Each logical qudit encodes the color of a node in $G$.

(a) Number of physical qubits and qudits versus the number of edges in $G$, the graph to be colored. $|C_2|$ is the number of physical qubits in the qubit cluster state implementation of $\exp(i\alpha H_1)$, where $H_1$ is defined by Eq. (<ref>). $|C_d|$ is the number of physical qudits in the cluster state implementation of $\exp(i\alpha H_C)$, where $H_C$ is defined in Eq. (<ref>). (b) Ratio of the qubit and qudit cluster sizes when $|V| = 500$ and $k = 3$, where $|V|$ is the number of vertices in $G$ and $k$ is the number of colors used for the coloring.

§ CONCLUSION AND OUTLOOK

In this work, we investigated a KLM-like qudit encoding for universal photonic quantum computation. First,  the optimization of non-deterministic many-qudit gates was discussed, then the solution of the $k$-coloring problem as an application of high-dimensional cluster states was examined.

The physical implementation of the single-qudit gates by the described  protocol should pose no problem in experiments, since they require only the use of passive linear optics. However, the low probabilities of the entangling many-qudit gates make this scheme applied with postselection challenging from the  universal quantum information processing point of view. Nevertheless, the approach may still be used to non-deterministically prepare highly entangled photon states, which are required by many other universal photonic schemes, most notably photonic cluster state computation. We present a numerical optimization method to improve the success rate of a general many-qudit gate. The method is capable of finding locally optimal solutions with perfect fidelity, improving the success rates typically by orders of magnitude compared to the naive implementations when the gates are decomposed into nonlinear phase shift operations. We only considered qudit gates where photon number of the input states is fixed, but in principle it could also be used when this restriction is lifted.
For example, it can be applied to numerically find an optimal solution for the nonlinear phase shift. Furthermore, it can be employed for more general hybrid qudit gates <cit.>, when the dimensions of the qudits are not the same.

To amend the problems of the non-deterministic postselection based approach, we used cluster states encoded in the multi-rail encoding to give an implementation for the high-dimensional QAOA.
We demonstrated our method on the $k$-coloring problem which can be solved by either using qubits or qudits with $d{=}k$.
We show that the high-dimensional cluster states proposed in this work to solve the $k$-coloring problem are made up of fewer photons and require less optical modes to encode than the KLM qubit cluster states needed to perform the same task.

The reduced number of photons can help with quantum storage since usually the fidelity of quantum memories drop quickly with number of photons due to photon loss <cit.>.
The real difficulty of cluster state computation is in the generation of the cluster states, and all results involving cluster states rely on the assumption that they can be generated efficiently.
In a recent work it has been shown that in theory one can produce these type of high-dimensional cluster states deterministically using quantum emitters <cit.>. The methods shown in Ref. <cit.> combined with the result of this work give a convincing argument for the usefulness of multi-rail encoded high-dimensional cluster states.

§ ACKNOWLEDGEMENTS

Funding information
MK and LO was supported by the Ministry of Culture and Innovation and the National Research, Development and Innovation Office within the Quantum Information National Laboratory of Hungary (Grant No. 2022-2.1.1-NL-2022-00004) and grants K131938, K142179, FK 135220. ZZ was supported supported by TKP Project no. TKP2021-NVA-04 financed under the TKP2021-NVA funding scheme.
LO also acknowledges support from of the NRDI Office of Hungary and the Hungarian Academy of Sciences through the Bolyai and Bolyai+ scholarships.


# Generic Poincaré-Bendixson Theorem for singularly perturbed monotone systems
with respect to cones of rank-$2$††thanks: Supported by NSF of China
No.11825106, 11871231, 11771414 and 11971232, CAS Wu Wen-Tsun Key Laboratory
of Mathematics, University of Science and Technology of China.

Lin Niua and Xizhuang Xiea,b

aSchool of Mathematical Sciences
University of Science and Technology of China
Hefei, Anhui, 230026, P. R. China
bSchool of Mathematical Sciences
Huaqiao University
Quanzhou, Fujian, 362021, P. R. China
Corresponding author.

###### Abstract

We investigate the singularly perturbed monotone systems with respect to cones
of rank $2$ and obtain the so called Generic Poincaré-Bendixson theorem for
such perturbed systems, that is, for a bounded positively invariant set, there
exists an open and dense subset $\mathcal{P}$ such that for each
$z\in\mathcal{P}$, the $\omega$-limit set $\omega(z)$ that contains no
equilibrium points is a closed orbit.

Keywords: Monotone systems with respect to $2$-cones, Regular perturbation,
Singular perturbation, Generic Poincaré-Bendixson Theorem

AMS Subject Classification (2020): 34C12, 34D15, 37C65

##  1 Introduction

In this paper, we are concerned with the dynamics of singularly perturbed
monotone systems with respect to a cone $C$ of rank $2$. Roughly speaking, a
cone $C$ of rank $2$ is a closed subset of $\mathbb{R}^{n}$ that contains a
linear subspace of dimension $2$ and no linear subspaces of higher dimension.
The concept of cones of rank $k$ was introduced by Fusco and Oliva [2] in a
finite-dimensional space, and Kransosel’skii et al. [15] in a Banach space.

Monotone dynamical systems whose flows respect an order structure have been
studied for a long time. It was Hirsch [6, 7, 8, 9, 10, 11] who started the
full description of the main dynamical properties in the setting of
cooperative and competitive systems, whose flows preserve the partial order
induced by a proper (i.e., convex, closed, solid, pointed) cone (see [12, 13,
25, 26]). Among others, the celebrated result of a classical monotone flow
with strong order preserving property is the so called Hirsch’s Generic
Convergence Theorem, which concludes that the set of all $x\in X$, for which
the omega-limit set $\omega(x)$ belongs to the set of equilibria, is generic
(open-dense, residual) in $X$ (see e.g., [7, 25, 13]). Furthermore, for a
classical smooth strongly monotone systems, precompact semi-orbits are
generically convergent to equilibria in the continuous-time case [21, 20, 27]
or to cycles in the discrete-time case [5, 22, 31, 33].

The classical monotone flow can be viewed as a monotone system with respect to
a cone of rank-$1$. And the long time behaviors of classical monotone flows
are determined by a one dimensional space in the cone of rank-$1$. For the
smooth monotone flows with respect to cones of rank $k$ $(k\geq 2)$ in
$\mathbb{R}^{n}$, Sánchez [24] first studied the structure of an omega limit
set. By using cones of rank $2$, Sánchez [24] projected the dynamics into
planes and proved a Poincaré-Bendixson type theorem: Any omega limit set of a
pseudo-ordered orbit that contains no equilibrium is a closed orbit. Here, an
orbit is called pseudo-ordered if it possesses one pair of distinct ordered
points. Besides, Sánchez’s work is strongly related to the theory of R. A.
Smith in [28, 29, 30], which obtained a Poincaré-Bendixson theorem for systems
in a higher dimensional space.

Later, Feng, Wang and Wu [3] established the Poincaré-Bendixson theorem of
pseudo-ordered orbits for strongly monotone semiflows (with respect to cones
of rank $2$) on a Banach space without any smoothness assumption. Very
recently, Weiss and Margaliot [34] described an important class of systems
that are monotone with respect to a cone of rank $k$ and presented the
Poincaré-Bendixson property for any bounded trajectory in the case $k=2$.
Feng, Wang and Wu [4] obtained the “generic Poincaré-Bendixson Theorem” for
strongly monotone flows with respect to cones of rank $2$, that is, for
generic (belonging to an open and dense set) points, the $\omega$-limit set
containing no equilibrium is a single closed orbit.

It deserves to point out that $C^{1}$-Closing Lemma was used as an effective
tool in perturbation theory both for Hirsch’s work [7] and Sánchez’s dynamical
result [24]. As a consequence, the perturbation of a monotone system began to
attract considerable attention. Hirsch [7] obtained that the regular
perturbation of a cooperative irreducible vector field is at most eventually
cooperative rather than cooperative. While, for singular perturbations of a
cooperative system, Sontag and Wang [32] showed that such perturbed system is
also eventually cooperative. Later, Niu [17] showed the $C^{1}$-regular
perturbation of a competitive irreducible vector field is at most eventually
competitive. For the high-rank cones, Sánchez [24] considered the regular
perturbation in the sense that a sequence $F_{n}$ of $C^{1}$-vector field
converges to $F$ in the $C^{1}$-topology.

Based on these results, one can see that the perturbed systems (eventual
cases) inherit sort of dynamical properties of unperturbed cooperative or
competitive systems (see [7, 18]). Thus, it is interesting to consider the
(both regular and singular) perturbations of a monotone system with respect to
high-rank cones and investigate the dynamics of the perturbed systems .

In the present paper, we try to carry out this task by focusing on a
singularly perturbed monotone system (with respect to cones of rank $2$)
having the form:

$\begin{cases}\frac{dx}{dt}&=f_{0}(x,y,\epsilon),\\\
\epsilon\frac{dy}{dt}&=g_{0}(x,y,\epsilon),\end{cases}$ (1.1)

with a positive parameter $\epsilon$ near zero and $(x,y)\in U\times V$, where
$U\subset\mathbb{R}^{n}$, $V\subset\mathbb{R}^{m}$ are open and bounded sets
such that $U\times V$ contains a compact set $D_{\epsilon}$. Under a series of
hypotheses (A1)-(A6) (see the details in Section $2$), including that the
limiting system ($\epsilon=0$) is monotone with respect to a cone of rank $2$
in $\mathbb{R}^{n}$, we are able to deduce the dynamics of perturbed system
from the unperturbed monotone system ($\epsilon=0$) with respect to a cone of
rank $2$. More precisely, the Poincaré-Bendixson type theorem is inherited for
the perturbed system:

Theorem A For singular perturbed system (1.1) satisfying (A1)-(A6), there
exists an open and dense subset $\mathcal{P}_{\epsilon}\subset D_{\epsilon}$
such that for each $z\in\mathcal{P}_{\epsilon}$, the $\omega$-limit set
$\omega(z)$ containing no equilibrium is a closed orbit.

We call this theorem as “generic Poincaré-Bendixson Theorem” for singular
perturbed monotone system with respect to cones of rank $2$. By the geometric
singular perturbation theory (see [1, 14, 16]) under the assumptions
(A1)-(A6), there is a slow manifold which attracts all flows in
$D_{\epsilon}$. The flow on the slow manifold can be treated as a
$C^{1}$-regular perturbation of the limiting flow for $\epsilon=0$. Further,
all the dynamics in $D_{\epsilon}$ could be tracked by the flows on the slow
manifold.

As one may know, one of the typical examples of a monotone flow with respect
to a cone of rank $2$ is the three-dimensional competitive systems. If the
limiting flow is competitive, system (1.1) can be treated as a singular
perturbed competitive system, and the flow on the slow manifold is eventually
competitive (see [17]). For an eventually competitive system, the asymptotic
behaviors of the complete (full) orbits have been investigated, due to a
priori drawback restriction of the so called Non-oscillation Principle for the
full orbits (see Niu and Wang [18]). To the best of our knowledge, very
limited work has been done on the dynamics (more specifically, the structure
of $\omega$-limit sets) of the singular perturbed competitive system. However,
in the view of cones of rank $2$, the asymptotic behaviors of the positive
orbits can be investigated in this paper. By applying Theorem A, we can obtain
the Generic Poincaré-Bendixson theorem for a singularly perturbed three-
dimensional competitive systems. We will present a concrete example to
illustrate our main result of generic Poincaré-Bendixson Theorem.

This paper is organized as follows. In section $2$, we introduce some basic
definitions and assumptions for singularly perturbed monotone systems with
respect to cones of rank $2$, and present the main result (see Theorem 2.3).
Section $3$ is devoted to give the detailed proof of the main theorem.
Finally, in Section $4$, an example for singular perturbations is presented
and numerical examples are given to illustrate the effectiveness of our
theoretical results.

##  2 Basic definitions and the main theorem

A nonempty closed set $C\subset\mathbb{R}^{n}$ is called a cone of rank $k$
($k$-cone) if it satisfies $\alpha C\subset C$ for all $\alpha\in\mathbb{R}$
and $\max\\{\text{dim}W:C\supset W\ \ \text{linear}\ \ \text{subspace}\\}=k$.
A cone $C$ is solid if $\text{Int }C\neq\emptyset$; and $C$ is called
$k$-solid if there is a $k$-dimensional linear subspace $W$ such that
$W\backslash\\{0\\}\subset\text{Int }C$. We write

$\displaystyle x$ $\displaystyle\thicksim y\quad$ if $\displaystyle\ \ y-x\in
C,$ $\displaystyle x$ $\displaystyle\thickapprox y$ if $\displaystyle\ \
y-x\in\text{Int }C.$

A flow $\phi_{t}$ is called monotone with respect to a $k$-solid cone $C$ if
$\phi_{t}(x)\thicksim\phi_{t}(y)$ whenever $x\thicksim y$ and $t\geq 0$. And
$\phi_{t}$ is called strongly monotone with respect to $C$ if $\phi_{t}$ is
monotone with respect to $C$ and $\phi_{t}(x)\thickapprox\phi_{t}(y)$ whenever
$x\neq y$, $x\thicksim y$ and $t>0$.

Consider the system:

$\frac{dx}{dt}=F(x),$ (2.1)

for which $F:O\rightarrow\mathbb{R}^{n}$ is a $C^{1}$-vector field and
$O\subset\mathbb{R}^{n}$ is an open convex set. We denote by $\phi_{t}$ the
flow generated by (2.1).

Let $U_{F}^{pq}(t)$ be the solution of

$\begin{cases}\frac{d}{dt}U=A_{F}^{pq}(t)U\\\ U(0)=I,\end{cases}$ (2.2)

where $A_{F}^{pq}(t)=\int^{1}_{0}DF(s\phi_{t}(p)+(1-s)\phi_{t}(q))ds$ for any
$p,q\in O$.

######  Definition 2.1.

The system (2.1) is called $C$-cooperative if for any $p,q\in O$, the matrix
$U_{F}^{pq}(t)$ satisfies $U_{F}^{pq}(t)[C\backslash\\{0\\}]\subset\text{Int
}C$ for all $t>0$.

######  Remark 2.2.

By Sánchez [24, Proposition 1], if system (2.1) is $C$-cooperative, then the
flow of (2.1) is strongly monotone with respect to cone $C$.

In this paper, we are concerned with the study of the singular perturbed
differential equations of the form

$\begin{cases}\frac{dx}{dt}=f_{0}(x,y,\epsilon),\\\
\epsilon\frac{dy}{dt}=g_{0}(x,y,\epsilon).\end{cases}$ (2.3)

System (2.3) can be reformulated with a change of time scale as

$\begin{cases}\frac{dx}{d\tau}=\epsilon f_{0}(x,y,\epsilon),\\\
\frac{dy}{d\tau}=g_{0}(x,y,\epsilon),\end{cases}$ (2.4)

where $\tau=t/\epsilon$. The time scale given by $\tau$ is said to be fast
whereas that for $t$ is slow. When $\epsilon$ is small enough, we call (2.3)
the slow system and (2.4) the fast system. And the two systems are equivalent
as long as $\epsilon\neq 0$.

Let $C^{r}_{b}$ denote a class of functions such that if a function $f$ is in
$C^{r}$ and its derivatives up to order $r$ as well as $f$ are bounded, then
$f\in C^{r}_{b}$. Throughout this paper, let $Df(x_{0})$ denote the
derivatives of $f$ evaluated at $x_{0}$ with respect to the variable $x$.
Meanwhile, $D_{x}f(x_{0},y_{0})$ and $D_{y}f(x_{0},y_{0})$ denote the partial
derivatives of $f$ with respect to $x$ and $y$ evaluated at $(x_{0},y_{0})$,
respectively. A family of sets $D_{\epsilon}\subset\mathbb{R}^{n}$ is upper
semicontinuous at $\epsilon_{1}\in[0,\epsilon_{0}]$ if given a neighborhood
$U$ of $D_{\epsilon_{1}}$, there exists a $\delta>0$ such that
$D_{\epsilon}\subset U$ for all
$\epsilon\in(\epsilon_{1}-\delta,\epsilon_{1}+\delta)$. $D_{\epsilon}$ is
lower semicontinuous at $\epsilon_{1}\in[0,\epsilon_{0}]$ if given any open
set $N$ such that $N\cap D_{\epsilon_{1}}$ is non-empty, there exists a
$\delta>0$ such that $N\cap D_{\epsilon}$ is non-empty for all
$\epsilon\in(\epsilon_{1}-\delta,\epsilon_{1}+\delta)$. It is continuous at
$\epsilon$ if it is both lower and upper semicontinuous at $\epsilon$ and
continuous if it is continuous at every $\epsilon\in[0,\epsilon_{0}]$. Then,
we have the following assumptions, where the integer $r>1$ and the positive
number $\epsilon_{0}$ are fixed from now on.

  1. (A1)

Let $U\subset\mathbb{R}^{n}$ and $V\subset\mathbb{R}^{m}$ be open and bounded
sets. The functions

$f_{0}:U\times V\times[0,\epsilon_{0}]\rightarrow\mathbb{R}^{n}$

and

$g_{0}:U\times V\times[0,\epsilon_{0}]\rightarrow\mathbb{R}^{m}$

are both of class $C^{r}_{b}$.

  2. (A2)

There is a function

$h_{0}:U\rightarrow V$

in $C^{r}_{b}$ such that $g_{0}(x,h_{0}(x),0)=0$ for all $x$ in $U$.

  3. (A3)

All eigenvalues of the matrix $D_{y}g_{0}(x,h_{0}(x),0)$ have negative real
parts for every $x\in U$.

  4. (A4)

There exists a family of convex compact sets $D_{\epsilon}\subset U\times V$,
which depend continuously on $\epsilon\in[0,\epsilon_{0}]$, such that (2.3) is
positively invariant on $D_{\epsilon}$ for $\epsilon\in(0,\epsilon_{0}]$.

  5. (A5)

For each $x\in U$, the system

$\frac{dz}{d\tau}=g_{1}(x,z,0)\triangleq g_{0}(x,z+h_{0}(x),0)$ (2.5)

is defined on $\\{z\in\mathbb{R}^{m}:z+h_{0}(x)\in V\\}$. And the steady state
$z=0$ of (2.5) is globally asymptotically stable on $\\{z:z+h_{0}(x)\in V\\}$.

  6. (A6)

Let $K_{0}$ be the projection of $D_{0}\cap\\{(x,y):y=h_{0}(x),x\in U\\}$ onto
the $x$-axis. For $x\in K_{0}$, the system

$\frac{dx}{dt}=f_{0}(x,h_{0}(x),0)$ (2.6)

is $C$-cooperative, where $C\subset\mathbb{R}^{n}$ is a cone of rank $2$.

Now, we present the generic Poincaré-Bendixson Theorem:

######  Theorem 2.3.

Assume that (A1)-(A6) hold. Then there exists a positive constant
$\epsilon^{*}<\epsilon_{0}$ such that for each $\epsilon\in(0,\epsilon^{*})$,
system (2.3) has the following property: there exists an open and dense subset
$\mathcal{P}_{\epsilon}\subset D_{\epsilon}$ such that, for each
$z\in\mathcal{P}_{\epsilon}$, $\omega$-limit set $\omega(z)$ that contains no
equilibrium is a single closed orbit.

##  3 Details of the proof

Our approach to prove the main theorem consists of two steps. First, we focus
on the regular perturbations of system (2.6). And then, we utilize the
geometric construction (see, e.g., [1, 23]) of system (2.3) to show the
generic dynamics.

###  3.1 Regularly perturbed monotone systems with respect to high-rank cone

For a general case, we consider the regular perturbations of system (2.1) in
this subsection. We focus on the system of ODE’s

$\frac{dx}{dt}=G(x),$ (3.1)

for which $G:O\rightarrow\mathbb{R}^{n}$ is a $C^{1}$-vector field. Let
$W\subset O$ be a convex compact subset. Then we obtain the following lemma.

######  Lemma 3.1.

Let system (2.1) be a $C$-cooperative system. Then there exists $\delta>0$
with the following property: If $\parallel F(z)-G(z)\parallel+\parallel
DF(z)-DG(z)\parallel<\delta$ for all $z\in O$, and $W$ is positively invariant
under the flow $\psi_{t}$ generated by $G$, then there exists $t_{*}>0$ such
that the matrix $U_{G}^{pq}(t)$ satisfies
$U_{G}^{pq}(t)[C\backslash\\{0\\}]\subset\text{Int }C$ for any $p,q\in W$ and
$t\geq t_{*}$.

###### Proof.

Pick $t_{*}>0$, we consider the matrices

$A_{G}^{pq}(t)=\int^{1}_{0}DG(s\psi_{t}(p)+(1-s)\psi_{t}(q))ds$

and the solution $U_{G}^{pq}(t)$ of

$\begin{cases}\frac{d}{dt}U=A_{G}^{pq}(t)U\\\ U(0)=I,\end{cases}$ (3.2)

for any $p,q\in W$.

We first prove the positiveness of $U_{G}^{pq}(t)$ for $t\in[t_{*},2t_{*}]$,
i.e., $U_{G}^{pq}(t)[C\backslash\\{0\\}]\subset\text{Int }C$ for
$t\in[t_{*},2t_{*}]$. Since system (2.1) is $C$-cooperative, there exists a
$\delta_{1}>0$ such that $B_{\delta_{1}}(U_{F}^{pq}(t)v)\subset\text{Int }C$
for all $v\in C\backslash\\{0\\}$ with $|v|=1$. A positive $\delta$ can be
found such that if $\parallel F(z)-G(z)\parallel+\parallel
DF(z)-DG(z)\parallel<\delta$ then $\parallel
U_{F}^{pq}(t)-U_{G}^{pq}(t)\parallel\leq\frac{\delta_{1}}{2}$ hold for
$t\in[t_{*},2t_{*}]$.

In fact, by the definition of flow, we have

$\begin{cases}\frac{d}{dt}\phi_{t}(z)=F(\phi_{t}(z)),\\\
\frac{d}{dt}\psi_{t}(z)=G(\psi_{t}(z)).\end{cases}$

Then,

$\displaystyle\parallel\phi_{t}(z)-\psi_{t}(z)\parallel$
$\displaystyle\leq\int^{t}_{0}\parallel F(\phi_{s}(z))-G(\psi_{s}(z))\parallel
ds$ $\displaystyle\leq\int^{t}_{0}\parallel
F(\phi_{s}(z))-F(\psi_{s}(z))\parallel ds+\int^{t}_{0}\parallel
F(\psi_{s}(z))-G(\psi_{s}(z))\parallel ds.$

Since $F$ and $G$ are $C^{1}$ and $t_{*}\leq t\leq 2t_{*}$, there exists some
$M>0$ such that $\parallel F(\phi_{s}(z))-F(\psi_{s}(z))\parallel\leq
M\parallel\phi_{s}(z)-\psi_{s}(z)\parallel$. Then
$\parallel\phi_{t}(z)-\psi_{t}(z)\parallel\leq\int^{t}_{0}M\parallel\phi_{s}(z)-\psi_{s}(z)\parallel
ds+\delta t$.

Using the Gronwall’s inequality, we obtain

$\parallel\phi_{t}(z)-\psi_{t}(z)\parallel\leq\delta 2t_{*}e^{M2t_{*}}.$

Hence, for $s\in[0,1]$,

$\displaystyle\ \ \ \
\parallel[s\phi_{t}(p)+(1-s)\phi_{t}(q)]-[s\psi_{t}(p)+(1-s)\psi_{t}(q)]\parallel$
$\displaystyle\leq
s\parallel\phi_{t}(p)-\psi_{t}(p)\parallel+(1-s)\parallel\phi_{t}(q)-\psi_{t}(q)\parallel$
$\displaystyle\leq s\delta 2t_{*}e^{M2t_{*}}+(1-s)\delta
2t_{*}e^{M2t_{*}}=\delta 2t_{*}e^{M2t_{*}}.$

Note that $U_{F}^{pq}(t)$ and $U_{G}^{pq}(t)$ are the solutions of (2.2) and
(3.2), respectively, then

$\displaystyle\parallel U_{F}^{pq}(t)-U_{G}^{pq}(t)\parallel\leq$
$\displaystyle\int^{t}_{0}\parallel
A_{F}^{pq}(s)U_{F}^{pq}(s)-A_{G}^{pq}(s)U_{G}^{pq}(s)\parallel ds$
$\displaystyle\leq$ $\displaystyle\int^{t}_{0}\parallel
A_{F}^{pq}(s)\parallel\parallel U_{F}^{pq}(s)-U_{G}^{pq}(s)\parallel ds$
$\displaystyle+\int^{t}_{0}\parallel
A_{F}^{pq}(s)-A_{G}^{pq}(s)\parallel\parallel U_{G}^{pq}(s)\parallel ds.$

Using the Gronwall’s inequality again, it follows that

$\parallel U_{F}^{pq}(t)-U_{G}^{pq}(t)\parallel\leq e^{\int^{t}_{0}\parallel
A_{F}^{pq}(s)\parallel ds}\int^{t}_{0}\parallel
A_{F}^{pq}(s)-A_{G}^{pq}(s)\parallel\parallel U_{G}^{pq}(s)\parallel ds.$

There exist $N>0$ and $L>0$ such that $\parallel U_{G}^{pq}(s)\parallel\leq N$
and $\parallel A_{F}^{pq}(s)\parallel\leq L$ for all $t\in[t_{*},2t_{*}]$.
Meanwhile,

$\displaystyle\parallel A_{F}^{pq}(t)-A_{G}^{pq}(t)\parallel\leq$
$\displaystyle\int^{1}_{0}\parallel
DF(s\phi_{t}(p)+(1-s)\phi_{t}(q))-DG(s\psi_{t}(p)+(1-s)\psi_{t}(q))\parallel
ds$ $\displaystyle\leq$ $\displaystyle\int^{1}_{0}\parallel
DF(s\phi_{t}(p)+(1-s)\phi_{t}(q))-DF(s\psi_{t}(p)+(1-s)\psi_{t}(q))\parallel$
$\displaystyle+\parallel
DF(s\psi_{t}(p)+(1-s)\psi_{t}(q))-DG(s\psi_{t}(p)+(1-s)\psi_{t}(q))\parallel
ds$ $\displaystyle\leq$ $\displaystyle\int^{1}_{0}\parallel
DF(s\phi_{t}(p)+(1-s)\phi_{t}(q))-DF(s\psi_{t}(p)+(1-s)\psi_{t}(q))\parallel
ds+\delta.$

Since $F$ is $C^{1}$ and
$\parallel[s\phi_{t}(p)+(1-s)\phi_{t}(q)]-[s\psi_{t}(p)+(1-s)\psi_{t}(q)]\parallel\leq\delta
2t_{*}e^{M2t_{*}}$, one can choose
$\delta<\frac{1}{2}\cdot\frac{\delta_{1}}{2}\cdot\frac{1}{2Nt_{*}e^{2t_{*}L}}$
small enough such that $\parallel
DF(s\phi_{t}(p)+(1-s)\phi_{t}(q))-DF(s\psi_{t}(p)+(1-s)\psi_{t}(q))\parallel<\frac{1}{2}\cdot\frac{\delta_{1}}{2}\cdot\frac{1}{2Nt_{*}e^{2t_{*}L}}$.
So, $\parallel
A_{F}^{pq}(t)-A_{G}^{pq}(t)\parallel\leq\frac{\delta_{1}}{2}\cdot\frac{1}{2Nt_{*}e^{2t_{*}L}}$,
which implies that $\parallel
U_{F}^{pq}(t)-U_{G}^{pq}(t)\parallel\leq\frac{\delta_{1}}{2}$. Thus, we obtain
that $U_{G}^{pq}(t)[C\backslash\\{0\\}]\subset\text{Int }C$ for
$t\in[t_{*},2t_{*}]$ and $p,q\in W$.

For $t>2t_{*}$, let us write $t=t_{0}+kt_{*}$ with $t_{0}\in[t_{*},2t_{*})$.
Define $t_{j}=jt_{*}$, $p_{j}=\psi_{t_{j}}(p)$, $q_{j}=\psi_{t_{j}}(q)$ with
$j=1,2,\cdots,k$. It is clear that $p_{j},q_{j}\in W$ if $W$ is positively
invariant. Since $\psi_{t}(p)-\psi_{t}(q)=U_{G}^{pq}(t)(p-q)$, we obtain that

$U_{G}^{pq}(t)=U_{G}^{p_{k}q_{k}}(t_{0})U_{G}^{p_{k-1}q_{k-1}}(t_{*})\cdots
U_{G}^{p_{1}q_{1}}(t_{*})U_{G}^{pq}(t_{*}).$

By the preceding proof, $U_{G}^{pq}(t)[C\backslash\\{0\\}]\subset\text{Int }C$
for $t>2t_{*}$.

Thus, we have proved that the system (3.1) is $C$-cooperative for $t\geq
t_{*}$. ∎

###  3.2 Proof of the main theorem

Before proving Theorem 2.3, we give the following lemma about the singularly
perturbed systems on
$\mathbb{R}^{n}\times\mathbb{R}^{m}\times[0,\epsilon_{0}]$, which is a
restatement of the Theorem 2.1 and Theorem 3.1 in Sakamoto [23].

######  Lemma 3.2.

Consider the system

${}\begin{cases}\frac{dx}{d\tau}=\epsilon f(x,y,\epsilon),\\\
\frac{dy}{d\tau}=g(x,y,\epsilon),\end{cases}$ (3.3)

where
$f:\mathbb{R}^{n}\times\mathbb{R}^{m}\times[0,\epsilon_{0}]\rightarrow\mathbb{R}^{n}$
is $C^{r}_{b}$ and
$g:\mathbb{R}^{n}\times\mathbb{R}^{m}\times[0,\epsilon_{0}]\rightarrow\mathbb{R}^{m}$
is $C^{r}_{b}$. For $x\in\mathbb{R}^{n}$, there is a $C^{r}_{b}$ function
$\overline{h}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ such that
$g(x,\overline{h}(x),0)=0$. And there exists a positive constant $\mu$ such
that all eigenvalues of the matrix $D_{y}g(x,\overline{h}(x),0)$ have negative
real parts less than $-\mu$ for every $x\in\mathbb{R}^{n}$.

Then, there exists a positive number $\epsilon_{1}<\epsilon_{0}$ such that for
every $\epsilon\in(0,\epsilon_{1}]$:

  1. $1)$

There exists a $C_{b}^{r-1}$ function
$h:\mathbb{R}^{n}\times[0,\epsilon_{1}]\rightarrow\mathbb{R}^{m}$ such that
the set $\mathbb{C}_{\epsilon}$ defined by

$\mathbb{C}_{\epsilon}=\\{(x,h(x,\epsilon)):x\in\mathbb{R}^{n}\\},\
\epsilon\in(0,\epsilon_{1}]$

is invariant under the flow generated by (3.3) and

$\underset{x\in\mathbb{R}^{n}}{\sup}\\{|h(x,\epsilon)-\overline{h}(x)|:x\in\mathbb{R}^{n}\\}\sim
O(\epsilon),as\ \epsilon\rightarrow 0.$

In particular, we have $h(x,0)=\overline{h}(x)$ for all $x\in\mathbb{R}^{n}$.

  2. $2)$

There is an $(n+m)$-dimensional $C^{r-1}$ submanifold
$W^{s}(\mathbb{C}_{\epsilon})$. It is characterized by

$W^{s}(\mathbb{C}_{\epsilon})=\\{(x_{0},y_{0}):\underset{\tau\geq
0}{\sup}|y(\tau;x_{0},y_{0})-h_{\epsilon}(x(\tau;x_{0},y_{0}))|e^{\frac{\mu\tau}{4}}<\infty\\},$

where $(x(\tau;x_{0},y_{0}),y(\tau;x_{0},y_{0}))$ is the solution of (3.3)
passing through $(x_{0},y_{0})$ and $h_{\epsilon}(x)=h(x,\epsilon)$ is the
function defining $\mathbb{C}_{\epsilon}$.

  3. $3)$

The manifold $W^{s}(\mathbb{C}_{\epsilon})$ is a disjoint union of the
$m$-dimensional $C^{r-1}$ manifold $W^{s}(\xi)$:

$W^{s}(\mathbb{C}_{\epsilon})=\underset{\xi\in\mathbb{R}^{n}}{\cup}W^{s}(\xi).$

Moreover, $W^{s}(\xi)$ is characterized as

$W^{s}(\xi)=\\{(x_{0},y_{0}):\underset{\tau\geq
0}{\sup}|\widetilde{x}(\tau)|e^{\frac{\mu\tau}{4}}<\infty,\ \underset{\tau\geq
0}{\sup}|\widetilde{y}(\tau)|e^{\frac{\mu\tau}{4}}<\infty\\},$

where $\widetilde{x}(\tau)=x(\tau;x_{0},y_{0})-H_{\epsilon}(\xi)(\tau)$,
$\widetilde{y}(\tau)=y(\tau;x_{0},y_{0})-h_{\epsilon}(H_{\epsilon}(\xi)(\tau))$,
where $H_{\epsilon}(\xi)(\tau)$ stands for a unique solution of
$\frac{dx}{d\tau}=\epsilon f(x,h_{\epsilon}(x),\epsilon)$,
$x(0)=\xi\in\mathbb{R}^{n}$.

  4. $4)$

There is a constant $\delta_{0}>0$ such that if a solution $(x(\tau),y(\tau))$
of (3.3) satisfies

$\underset{\tau\geq 0}{\sup}|y(\tau)-h_{\epsilon}(x(\tau))|<\delta_{0},$

then $(x(0),y(0))\in W^{s}(\mathbb{C}_{\epsilon})$.

  5. $5)$

The fibers are positively invariant in the sense that

$W^{s}(H_{\epsilon}(\xi)(\tau))=\\{(x(\tau;,x_{0},y_{0}),y(\tau;x_{0},y_{0})):(x_{0},y_{0})\in
W^{s}(\xi)\\},$

for each $\tau\geq 0$.

  6. $6)$

The fibers restricted to the $\delta_{0}$ neighborhood of
$\mathbb{C}_{\epsilon}$, denoted by $W^{s}_{\epsilon,\delta_{0}}$, can be
parameterized as follows. Let
$L_{\delta_{0}}=\\{z\in\mathbb{R}^{m}:|z|\leq\delta_{0}\\}$. There are two
$C^{r-1}_{b}$ functions

$P_{\epsilon,\delta_{0}}:\mathbb{R}^{n}\times
L_{\delta_{0}}\rightarrow\mathbb{R}^{n},$
$Q_{\epsilon,\delta_{0}}:\mathbb{R}^{n}\times
L_{\delta_{0}}\rightarrow\mathbb{R}^{m},$

and a map

$T_{\epsilon,\delta_{0}}:\mathbb{R}^{n}\times
L_{\delta_{0}}\rightarrow\mathbb{R}^{n}\times\mathbb{R}^{m}$

mapping $(\xi,\eta)$ to $(x,y)$, where

$x=\xi+P_{\epsilon,\delta_{0}}(\xi,\eta),\
y=h(x,\epsilon)+Q_{\epsilon,\delta_{0}}(\xi,\eta),$

such that

$W^{s}_{\epsilon,\delta_{0}}(\xi)=T_{\epsilon,\delta_{0}}(\xi,L_{\delta_{0}}).$

######  Remark 3.3.

  1. (1)

The $\delta_{0}$ in property 4) of Lemma 3.2 can be chosen uniformly for
$\epsilon\in(0,\epsilon_{0}]$.

  2. (2)

The property 3) of Lemma 3.2 is often referred to as the asymptotic phase
property in the way

$|x(\tau;x_{0},y_{0})-H_{\epsilon}(\xi)(\tau)|\rightarrow 0,$
$|y(\tau;x_{0},y_{0})-h_{\epsilon}(H_{\epsilon}(\xi)(\tau)|\rightarrow 0,$

as $\tau\rightarrow\infty$.

In order to use Sakamoto’s results in Lemma 3.2, we firstly extend the vector
fields from $U\times V$ to $\mathbb{R}^{n}\times\mathbb{R}^{m}$ for
$\epsilon\in[0,\epsilon_{0}]$. The technique is standard by [16], which can
also be found in [32], such that the extended system:

$\begin{cases}\frac{dx}{d\tau}=\epsilon f(x,y,\epsilon),\\\
\frac{dy}{d\tau}=g(x,y,\epsilon),\end{cases}$

satisfies the assumptions $(A1)$-$(A6)$ and the assumptions for the geometric
singular perturbation in Lemma 3.2. Moreover, $\overline{h}(x)$ coincides with
$h_{0}(x)$ on $K$, $f$ and $g$ coincide with $f_{0}$, $g_{0}$ on
$\Omega_{d_{1}}$, respectively. Where $K$ is a compact set with $K_{0}\subset
K\subset U$, $\Omega_{d_{1}}\triangleq\\{(x,y):x\in K,y\in V,|y-h_{0}(x)|\leq
d_{1}\\}$ and $d_{1}>0$ is fixed such that $\delta_{0}$ in Lemma 3.2 is less
than $d_{1}$.

###### Proof of Theorem 2.3.

First, we consider the solutions on the invariant manifold
$\mathbb{C}_{\epsilon}$ satisfying

$\begin{cases}\frac{dx}{dt}=f(x,h_{\epsilon}(x),\epsilon),\\\
y(t)=h_{\epsilon}(x(t)).\end{cases}$ (3.4)

For brevity, we only focus on the $x$-direction on the invariant manifold
$\mathbb{C}_{\epsilon}$, since $y=h_{\epsilon}(x)$. Clearly, the limiting
equation of (3.4) is (2.6) as $\epsilon$ approaches zero. For system (2.6),
the matrix $U_{f_{0}}^{pq}(t)$ satisfies
$U_{f_{0}}^{pq}(t)[C\backslash\\{0\\}]\subset\text{Int }C$ for all $p,q\in
K_{0}$ and $t>0$. By the continuity of $D_{\epsilon}$ and $h_{\epsilon}(x)$ at
$\epsilon=0$, we can pick an $\epsilon_{2}<\epsilon_{1}$ small enough such
that $U_{f_{0}}^{pq}(t)$ satisfies
$U_{f_{0}}^{pq}(t)[C\backslash\\{0\\}]\subset\text{Int }C$ for $p,q\in
K_{\epsilon}$ and $\epsilon\in(0,\epsilon_{2})$, where $K_{\epsilon}$ is the
projection of $\mathbb{C}_{\epsilon}\cap D_{\epsilon}$ to the $x$-axis. Note
also that $D_{\epsilon}$ is positively invariant under (3.4) and
$\mathbb{C}_{\epsilon}$ is an invariant manifold. Then, $K_{\epsilon}$ is
positively invariant under the flow $\psi^{\epsilon}_{t}$ of (3.4). Applying
Lemma 3.1, we obtain that there exist an $\epsilon_{3}\leq\epsilon_{2}$ and
some $t_{*}>0$ such that for each $\epsilon\in(0,\epsilon_{3})$, system (3.4)
is $C$-cooperative for $t\geq t_{*}$.

By [24, Proposition 1], it is clear that the flow of system (3.4) satisfies
the assumption (FWW) (see [4, p.4]) for $t\geq t_{*}$. Then [4, Theorem 5.3]
implies that there exists an open and dense subset $M_{\epsilon}\subset
K_{\epsilon}$ such that for any $x\in M_{\epsilon}$, the $\omega$-limit set
$\omega(x)$ containing no equilibrium is a single closed orbit.

It follows from the Lemma $7$ in [32] that there exist $\epsilon_{4}>0$ and
$d\in(0,\delta_{0})$ such that, for each $\epsilon\in(0,\epsilon_{4})$,
$(x_{0},y_{0})\in W^{s}(\mathbb{C}_{\epsilon})$ whenever $(x_{0},y_{0})\in
D_{\epsilon}$ satisfies $|y_{0}-h_{\epsilon}(x_{0})|<d$. Moreover, Lemma 3.2
$3)$ guarantees that $(x_{0},y_{0})\in W^{s}_{\epsilon,d}(\xi)$, where $\xi\in
K_{\epsilon}$ and $W^{s}_{\epsilon,d}(\xi)$ is the stable fiber restricted to
the $d$ neighborhood of $\mathbb{C}_{\epsilon}$.

Let $N_{\epsilon}=\\{(x,y)\in D_{\epsilon}:|y-h_{\epsilon}(x)|<d\\}$. In the
following, we will show that the set $\underset{\xi\in
M_{\epsilon}}{\cup}W^{s}_{\epsilon,d}(\xi)$ is open and dense in
$N_{\epsilon}$. By Lemma 3.2 $6)$, the stable fiber can be characterized as
$W^{s}_{\epsilon,d}(\xi)=T_{\epsilon,d}(\xi,L_{d})$, where
$L_{d}=\\{z\in\mathbb{R}^{m}:|z|<d\\}$. In other words, $\underset{\xi\in
M_{\epsilon}}{\cup}W^{s}_{\epsilon,d}(\xi)=T_{\epsilon,d}(M_{\epsilon},L_{d})$.
For any $(x_{0},y_{0})\in N_{\epsilon}$, there exists $(\xi_{0},\eta_{0})\in
K_{\epsilon}\times L_{d}$ such that
$(x_{0},y_{0})=T_{\epsilon,d}(\xi_{0},\eta_{0})$. Since $M_{\epsilon}$ is
dense in $K_{\epsilon}$, there exists a sequence $\\{\xi_{n}\\}\subset
M_{\epsilon}$ such that $\xi_{n}\to\xi_{0}$. Let $\eta_{n}=\eta_{0}$ and
$(x_{n},y_{n})=T_{\epsilon,d}(\xi_{n},\eta_{n})$ for $n=1,2,\cdots$, then
$(\xi_{n},\eta_{n})\to(\xi_{0},\eta_{0})$ and $(x_{n},y_{n})\in
T_{\epsilon,d}(M_{\epsilon},L_{d})$. Since the mapping $T_{\epsilon,d}$ is
continuous, one has $(x_{n},y_{n})\to(x_{0},y_{0})$ as $n\to\infty$. Thus, the
set $T_{\epsilon,d}(M_{\epsilon},L_{d})$ is dense in $N_{\epsilon}$. In order
to prove the openess of $T_{\epsilon,d}(M_{\epsilon},L_{d})$, we only need to
show the inverse of the mapping $T_{\epsilon,d}$ is continuous, because both
$M_{\epsilon}$ and $L_{d}$ are open sets. Clearly, the mapping
$T_{\epsilon,d}$ is both injective and surjective. And by the proof of the
Claim $3.7$ in [23], the inverse of $T_{\epsilon,d}$ is continuous (In fact,
there exists a mapping $\phi^{*}$, see [23, Lemma 3.9], such that
$\xi_{0}=x_{0}-\phi^{*}(0)$ and $\eta_{0}=y_{0}-h_{\epsilon}(x_{0})$ hold for
$(x_{0},y_{0})\in N_{\epsilon}$). Thus, the set
$T_{\epsilon,d}(M_{\epsilon},L_{d})$ is open and dense in $N_{\epsilon}$.

Moreover, there exist some positive $\tau_{0}$ and an
$\epsilon_{5}<\epsilon_{4}$ such that
$|y(\tau_{0})-h_{\epsilon}(x(\tau_{0}))|<d$ for all
$\epsilon\in(0,\epsilon_{5})$, whenever $(x(\tau),y(\tau))$ is the solution to
(2.4) with $(x(0),y(0))=(x_{0},y_{0})\in D_{\epsilon}$. This implies that
$(x(\tau_{0}),y(\tau_{0}))\in N_{\epsilon}$. The corresponding proof is
similar as Lemma 8 in [32] with the positive invariance of $D_{\epsilon}$ (In
fact,
$|y(\tau)-h_{\epsilon}(x(\tau))|\leq|y(\tau)-h_{0}(x(\tau))|+|h_{0}(x(\tau))-h_{\epsilon}(x(\tau))|$,
the assumption $A5$ and Lemma 3.2 $1)$ imply that there exist a $\tau_{0}>0$
and $\epsilon_{5}<\epsilon_{4}$ such that
$|y(\tau_{0})-h_{0}(x(\tau_{0}))|<\dfrac{d}{2}$ and
$|h_{0}(x(\tau_{0}))-h_{\epsilon}(x(\tau_{0}))|<\dfrac{d}{2}$). Let
$\phi^{\epsilon}_{\tau}$ be the flow of (2.4), and
$F(D_{\epsilon})=\\{\phi^{\epsilon}_{\tau_{0}}((x,y)):(x,y)\in
D_{\epsilon}\\}\subset N_{\epsilon}$. Since $\phi^{\epsilon}_{\tau_{0}}$ is
diffeomorphism, the set $G_{\epsilon}=F(\text{Int }D_{\epsilon})\cap
T_{\epsilon,d}(M_{\epsilon},L_{d})$ is open and dense in $F(D_{\epsilon})$.
So, $F^{-1}(G_{\epsilon})$ is open and dense in $D_{\epsilon}$. Let
$\mathcal{P}_{\epsilon}=F^{-1}(G_{\epsilon})$. For
$(x,y)\in\mathcal{P}_{\epsilon}$, the fact that
$\phi^{\epsilon}_{\tau_{0}}((x,y))\in T_{\epsilon,d}(M_{\epsilon},L_{d})$
implies the $\omega$-limit set $\omega((x,y))$ containing no equilibrium is a
single closed orbit by the asymptotic phase property.

We have completed the proof of Theorem 2.3 by taking
$\epsilon^{*}=\min\\{\epsilon_{3},\epsilon_{5}\\}$. ∎

##  4 An Example

In this section, we consider the following singular perturbed ODE’s system:

$\begin{cases}\frac{dx}{dt}=x-y-\frac{3}{2}xz^{2}-\frac{x(x^{2}+y^{2})}{2}+\epsilon
xw,\\\ \frac{dy}{dt}=x+y-\frac{3}{2}yz^{2}-\frac{y(x^{2}+y^{2})}{2}+\epsilon
yw,\\\ \frac{dz}{dt}=-z-\frac{z^{3}}{2}-\frac{3}{2}z(x^{2}+y^{2})+\epsilon
zw,\\\ \epsilon\frac{dw}{dt}=-w+x+y+z,\end{cases}$ (4.1)

where the parameter $\epsilon<\epsilon_{0}=0.1$. Let

$\begin{split}D_{\epsilon}=D=\\{&(x,y,z,w):|x|\leq 4,|y|\leq 4,|z|\leq
4,|w|\leq 16\\}.\end{split}$

Assume that $U\subset\mathbb{R}^{3}$, $V\subset\mathbb{R}^{1}$ are open
bounded sets such that $D_{\epsilon}\subset U\times V$.

We consider the limiting system:

$\begin{cases}\frac{dx}{dt}=x-y-\frac{3}{2}xz^{2}-\frac{x(x^{2}+y^{2})}{2},\\\
\frac{dy}{dt}=x+y-\frac{3}{2}yz^{2}-\frac{y(x^{2}+y^{2})}{2},\\\
\frac{dz}{dt}=-z-\frac{z^{3}}{2}-\frac{3}{2}z(x^{2}+y^{2}).\end{cases}$ (4.2)

This system is inspired by Ortega and Sánchez [19], where they found a stable
closed orbit of so called $P$-competitive systems (see [19, P2914]) such that
every orbit tends to the closed orbit or the origin. It is easy to see that
for the limiting system (4.2), the origin is the unique equilibrium point and
the linearization in the origin has eigenvalues $-1$ and $1\pm i$.

Let $P$ be a diagonal matrix

$P=\begin{pmatrix}-1&0&0\\\ 0&-1&0\\\ 0&0&1\end{pmatrix}.$

Let $\lambda:\mathbb{R}^{3}\to\mathbb{R}$ be a continuous function (not
necessarily positive) such that the matrices

$PDF(\xi)+DF(\xi)^{*}P+\lambda(\xi)P$ (4.3)

are negatively definite, for each $\xi\in U$, where $F$ is the vector field of
(4.2) and $DF(\xi)^{*}$ stands for the transpose of $DF(\xi)$. In fact, a
straightforward calculation yeilds that for any function $\lambda$ with
$3x^{2}+3y^{2}+3z^{2}-2<\lambda(x,y,z)<3x^{2}+3y^{2}+3z^{2}+2$, the matrices
(4.3) are negatively definite.

Define

$C=\\{\xi\in\mathbb{R}^{3}:\langle P\xi,\xi\rangle\leq 0\\},$

where $\langle\cdot,\cdot\rangle$ is the inner product in $\mathbb{R}^{3}$. It
is clear that $C$ is a $2$-solid cone and system (4.2) is $C$-cooperative (see
[24] or [3, 4]).

One can check that the assumptions (A1)-(A6) hold for system (4.1), since the
fast variable $w$ has the simple case that the fourth equation of (4.1) is
linear for $w$. Thus, we are able to obtain the following generic Poincaré-
Bendixson Theorem:

######  Theorem 4.1.

There exists a positive constant $\epsilon^{*}<\epsilon_{0}$ such that for
each $\epsilon\in(0,\epsilon^{*})$, system (4.1) has the following property:
there exists an open and dense subset $\mathcal{P}_{\epsilon}\subset
D_{\epsilon}$ such that, for each $z\in\mathcal{P}_{\epsilon}$, $\omega$-limit
set $\omega(z)$ that contains no equilibrium is a single closed orbit.

By a numerical simulation with $\epsilon=0.05$, we illustrate that the flow
$\phi$ of system (4.1) with initial value
$(x_{0},y_{0},z_{0},w_{0})=(2,2,3,12)\in D_{\epsilon}$ tends to a periodic
orbit in Figure 1. It is worth pointing out that this initial value is not
specially selected, actually, the orbits initiated from generic points in
$D_{\epsilon}$ will be attracted to closed orbits.

Fig 1: The solution curves of system (4.1) with initial value
$(x_{0},y_{0},z_{0},w_{0})$

$=(2,2,3,12)$. The small graph is a detail view when t $\in[0,1.5]$.

The corresponding phase portraits are shown in Figure 2. They illustrate that
the projections of the flow $\phi_{t}(x_{0},y_{0},z_{0},w_{0})$ onto
$3$-dimensional spaces tend to periodic orbits, respectively. This implies
that the flow $\phi_{t}(x_{0},y_{0},z_{0},w_{0})$ converges to a closed orbit
in four dimensional space.

(a) The projection of the flow $\phi_{t}(x_{0},y_{0},z_{0},w_{0})$ in space
O-xyz and the flow $\psi_{t}(x_{0},y_{0},z_{0})$ of limiting system (4.2). (b)
The projection of the flow $\phi_{t}(x_{0},y_{0},z_{0},w_{0})$ in space O-xyw.
(c) The projection of the flow $\phi_{t}(x_{0},y_{0},z_{0},w_{0})$ in space
O-yzw. (d) The projection of the flow $\phi_{t}(x_{0},y_{0},z_{0},w_{0})$ in
space O-xzw.

Fig 2: The corresponding phase portraits of Fig 1.

######  Remark 4.2.

Define the set $K=\\{\xi\in\mathbb{R}^{3}:\langle P\xi,\xi\rangle\geq
0,\langle\xi,v_{+}\rangle\geq 0\\}$, where $v_{+}$ is an eigenvector of $P$
with respect to the positive eigenvalue $1$. It is clear that $K$ is a usual
convex cone with nonempty interior. Moreover, the set
$C=\\{\xi\in\mathbb{R}^{3}:\langle P\xi,\xi\rangle\leq 0\\}$ can be written as
$C=\overline{\mathbb{R}^{3}\backslash(K\cup(-K))}$. We say a flow $\psi_{t}$
is strongly competitive with respect to the partial ordering induced by $K$ if
$y-x\in\text{Int }K$ whenever $\psi_{t}(y)-\psi_{t}(x)\in K\backslash\\{0\\}$
and $t>0$.

We emphasize that a flow $\psi_{t}$ is strongly monotone with respect to $C$
if and only if $\psi_{t}$ is strongly competitive; and moreover,
$D\psi_{t}(x)[C\backslash\\{0\\}]\subset\text{Int }C$ for $x\in U$ and $t>0$
if and only if $D\psi_{-t}(\psi_{s}(x))[K\backslash\\{0\\}]\subset\text{Int
}K$ for $x\in U$, $0<t\leq s$.

As a consequence, the fact that system (4.2) is $C$-cooperative implies that
its flow $\psi_{t}^{0}$ satisfies that
$D\psi_{-t}^{0}(\xi)[K\backslash\\{0\\}]\subset\text{Int }K$ for $t>0$ and
$-t\in I(\xi)$, where $I(\xi)\subset\mathbb{R}$ denotes the maximal interval
of existence of the solution passing though $\xi$. By our previous work in
[18, 17], there exists a positive constant $\epsilon_{*}<\epsilon_{0}$ such
that for each $\epsilon\in(0,\epsilon_{*})$, the $\omega$-limit set of (4.1)
is equivalent to an $\omega$-limit set of an eventually competitive system.
Moreover, if the orbit passing though any point $\zeta=(x,y,z,w)\in
D_{\epsilon}$ is tracked by a full orbit in $K_{\epsilon}$, then the
$\omega$-limit set $\omega(\zeta)$ containing no equilibrium is a single
closed orbit. The main result in the present paper provides a different view
that a flow which is competitive can be viewed as a monotone system with
respect to a high-rank cone, by which we can ignore the priori drawback
restriction of the so called Non-oscillation Principle for the full orbits in
competitive systems.

## Acknowledgments

The authors are greatly indebted to Professor Yi Wang for valuable suggestions
which led to much improvement of this paper.

## References

  * [1] N. Fenichel, Geometric singular perturbation theory for ordinary differential equations, J. Diff. Equ., 31(1979), 53-98.
  * [2] G. Fusco and W. Oliva, A Perron theorem for the existence of invariant subspaces, Ann. Mat. Pura. Appl., 160(1991), 63-76.
  * [3] L. Feng, Y. Wang and J. Wu, Semiflows “monotone with respect to highrank cones” on a Banach space, SIAM J. Math. Anal., 49(2017), 142-161.
  * [4] L. Feng, Y. Wang and J. Wu, Generical behavior of flows strongly monotone with respect to high-rank cones, J. Differential Equations, 275 (2021), 858–881.
  * [5] P. Hess and P. Poláčik, Boundeness of prime periods of stable cycles and convergence to fixed points in discrete monotone dynamical systems, SIAM J. Math. Anal. 24(1993), 1312-1330.
  * [6] M. W. Hirsch, Systems of differential equations which are competitive or cooperative I: limit sets, SIAM J. Appl. Math., 13 (1982), 167–179.
  * [7] M. W. Hirsch, Systems of differential equations which are competitive or cooperative II: convergence almost everywhere, SIAM J. Math. Anal., 16(1985), 423-439.
  * [8] M. W. Hirsch, Systems of differential equations which are competitive or cooperative III: compting species, Nonlinearity, 1 (1988), 51–71.
  * [9] M. W. Hirsch, Systems of differential equations which are competitive or cooperative IV: Structural stability in three dimensional systems, SIAM J. Math. Anal., 21 (1990), 1225–1234.
  * [10] M. W. Hirsch, Systems of differential equations that are competitive or cooperative V: convergence in $3$-dimensional systems, J. Differential Equations, 80 (1989), 94–106.
  * [11] M. W. Hirsch, Systems of differential equations that are competitive or cooperative VI: a local $C^{r}$ closing lemma for $3$-dimensional systems, Ergodic Theory Dynam. Systems, 11 (1991), 443–454.
  * [12] M. W. Hirsch and H. L. Smith, Monotone Systems, A Mini-review, Proceedings of the First Multidisciplinary Symposium on Positive Systems (POSTA 2003), Luca Benvenuti, Alberto De Santis and Lorenzo Farina (Eds.) Lecture Notes on Control and Information Sciences Vol. 294, Springer-Verlag, Heidelberg, 2003.
  * [13] M. W. Hirsch and H. L. Smith, Monotone dynamical systems, Handbook of Differential Equations: Ordinary Differential Equations, Vol. 2, Elsevier, Amsterdam 2005.
  * [14] C.K.R.T. Jones, Geometric singular perturbation theory, Dynamical Systems (Montecatini Terme,1994). Lect. Notes in Math., 1609, Springer, Berlin, 1995.
  * [15] M. A. Krasnosel’skii, J. A. Lifshits and A. V. Sobolev, Positive Linear Systems, the Method of Positive Operators, Heldermann Verlag, Berlin, 1989.
  * [16] K. Nipp, Smooth attractive invariant manifolds of singularly perturbed ODE’s, Research Report, (1992), 92–13.
  * [17] Lin Niu, Eventually competitive systems generated by perturbations, Electronic Journal of Differential Equations, 121(2019), 1-12.
  * [18] Lin Niu and Yi Wang, Non-oscillation principle for eventually competitive and cooperative systems, Discrete Contin. Dyn. Syst. B, 24 (2019), 6481-6494.
  * [19] R. Ortega and L. Sánchez, Abstract competitive systems and orbital stability in $\mathbb{R}^{3}$, Proc. Amer. Math. Soc., 128(2000), 2911-2919.
  * [20] P. Poláčik, Convergence in smooth strongly monotone ows defined by semilinear parabolic equations, J. Differential Equations, 79(1989), 89-100.
  * [21] P. Poláčik, Generic properties of strongly monotone semiflows defined by ordinary and parabolic differential equations, Qualitative theory of differential equations (Szeged 1988), 519-530, Colloq. Math. Soc. János Bolyai, 53, North-Holland, Amsterdam, 1990.
  * [22] P. Poláčik and I. Tereščák, Convergence to cycles as a typical asymptotic behavior in smooth strongly monotone discrete-time dynamical systems, Arch. Ration. Mech. Anal. 116(1992), 339-360.
  * [23] K. Sakamoto, Invariant manifolds in singular perturbation problems for ordinary differential equations, Proc. R. Soc. Edinb. A, 116(1990), 45-78.
  * [24] L. A. Sánchez, Cones of rank $2$ and the Poincaré-Bendixson property for a new class of monotone systems, J. Differential Equations, 216(2009), 1170-1190.
  * [25] H. L. Smith, Monotone Dynamical Systems, an introduction to the theory of competitive and cooperative systems, Math. Surveys and Monographs, 41, Amer. Math. Soc., Providence, Rhode Island 1995.
  * [26] H. L. Smith, Monotone dynamical systems: Reflections on new advances and applications, Discrete Contin. Dyn. Syst., 37(2017), 485-504.
  * [27] H. Smith and H. Thieme, Quasi convergence and stability for strongly order-preserving semiflows, SIAM J. Math. Anal. 21(1990), 673-692.
  * [28] R. A. Smith, The Poincaré-Bendixson theorem for certain differential equations of higher order, Proc. of Royal Soc. of Edinburgh, 83A(1979), 63-79.
  * [29] R. A. Smith, Existence of periodic orbits of autonomous ordinary differential equations, Proc. of Royal Soc. of Edinburgh A, 85(1980), 153-172.
  * [30] R. A. Smith, Orbital stability for ordinary differential equations, J. Differential Equations, 69(1987), 265-287.
  * [31] I. Tereščák, Dynamics of $C^{1}$ smooth strongly monotone discrete-time dynamical systems, preprint, Comenius University, Bratislava, 1994.
  * [32] L. Wang and E. D. Sontag, Singularly perturbed monotone systems and an application to double phosphorylation cycles, J. Nonlinear Sci., 18(2008), 527-550.
  * [33] Y. Wang and J. Yao, Dynamics alternatives and generic convergence for $C^{1}$-smooth strongly monotone discrete dynamical systems, J. Differential Equations, 269(2020), 9804-9818.
  * [34] E. Weiss and M. Margaliot, A generalization of linear positive systems with applications to nonlinear systems: Invariant sets and the Poincaré-Bendixson property, Automatica J.IFAC, 123, 2021.