*Article* **Synchronization of Epidemic Systems with Neumann Boundary Value under Delayed Impulse**

**Ruofeng Rao 1,\*, Zhi Lin 1,2, Xiaoquan Ai <sup>1</sup> and Jiarui Wu <sup>1</sup>**


**Abstract:** This paper reports the construction of synchronization criteria for the delayed impulsive epidemic models with reaction–diffusion under the Neumann boundary value. Different from the previous literature, the reaction–diffusion epidemic model with a delayed impulse brings mathematical difficulties to this paper. In fact, due to the existence of second-order partial derivatives in the reaction–diffusion model with a delayed impulse, the methods of first-order ordinary differential equations from the previous literature cannot be effectively applied in this paper. However, with the help of the variational method and an appropriate boundedness assumption, a new synchronization criterion is derived, and its effectiveness is illustrated by numerical examples.

**Keywords:** Neumann boundary value; delayed impulse; synchronization; reaction–diffusion epidemic models; variational methods

**MSC:** 34K24; 34K45

#### **1. Introduction**

The dynamics of epidemic models has always been a hot topic [1,2]. Ordinary differential equation epidemic dynamic models are the most common models, and fractional order models especially have been hot topics in recent research [3–6] whose ideas or methods have been applied to studying epidemic dynamic models. Moreover, the reaction–diffusion epidemic models have become one of the key topics because of the migration behavior of the population [7–10]. Usually, infectious diseases are controlled within a certain range, so we consider the Neumann boundary value, that is, there is no diffusion on the boundary of the infectious area because the disease area is usually isolated from the outside world by some measures, so the Neumann zero boundary value is considered in this paper. To prevent the spread of disease, the government or relevant departments often take impulse measures. This impulse management measure is not only aimed at the epidemic situation, but also considered impulse control measures for economic management, mechanical engineering and other issues [11–20]. Delayed impulse models have also been investigated by many researchers [11,12], for delayed impulse models better simulate the actual situation, that is, the impulse effect usually takes some time to appear. However, the models with a delayed impulse are usually ordinary differential systems, and reaction–diffusion systems with a delayed impulse are rarely seen in the existing literature. This inspired us to write this paper. In fact, due to the existence of second-order partial derivatives in the reaction–diffusion model with a delayed impulse, the methods of first-order ordinary differential equations in the existing literature cannot be effectively applied to partial differential equations. By means of the variational method, differential mean value theorem and convergence of sequence, the synchronization criterion of the delayed impulse reaction–diffusion epidemic model is obtained in this paper. Intuition tells us that the shorter the impulse interval, the greater the impulse frequency and the faster the impulse effect appears, and the greater

**Citation:** Rao, R.; Lin, Z.; Ai, X.; Wu, J. Synchronization of Epidemic Systems with Neumann Boundary Value under Delayed Impulse. *Mathematics* **2022**, *10*, 2064. https:// doi.org/10.3390/math10122064

Academic Editor: Quanxin Zhu

Received: 30 May 2022 Accepted: 13 June 2022 Published: 15 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the impulse intensity, the faster the synchronization between the models must be. These intuitive conclusions are affirmed in the synchronization criterion given in this paper.

The rest of this paper is organized as follows. In Section 2, we present some preliminaries about the reaction–diffusion epidemic model with a delayed impulse. In Section 3, we propose and derive the synchronization criterion for reaction–diffusion epidemic models under a delayed impulse. In Section 4, an illustrative numerical example is provided to show the effectiveness of the newly obtained criterion. Finally, some conclusions are written in Section 5.

The main contributions are as follows:


#### **2. System Description and Preliminaries**

In [1], the following epidemic system was studied:

$$\begin{cases} \frac{dS}{dt} = -S\beta(t)I, \\ \frac{dI}{dt} = S\beta(t)I - I\gamma(t), \\ \frac{dR}{dt} = I\gamma(t), \end{cases} \tag{1}$$

where the function *S*(*t*) is the fraction of the susceptible population, *I*(*t*) the infected fraction, *R*(*t*) the recovered fraction, and 0 < *S*(*t*) < 1, 0 < *I*(*t*) < 1, 0 < *R*(*t*) < 1. In addition, the disease transmission rate is denoted by *β*(*t*) and the recovery rate is *γ*(*t*). In 2020, the authors of [7] considered the epidemic system with inevitable diffusions:

$$\begin{cases} \frac{\partial S(t,\mathbf{x})}{\partial t} = \Delta[d\_1 S(t,\mathbf{x})] - I(t,\mathbf{x})\beta(t)S(t,\mathbf{x}),\\ \frac{\partial I(t,\mathbf{x})}{\partial t} = \Delta[d\_2 I(t,\mathbf{x})] + I(t,\mathbf{x})\beta(t)S(t,\mathbf{x}) - I(t,\mathbf{x})\gamma(t),\\ \frac{\partial R(t,\mathbf{x})}{\partial t} = \Delta[d\_3 R(t,\mathbf{x})] + I(t,\mathbf{x})\gamma(t), \end{cases} \tag{2}$$

where Δ is the Laplacian operator.

Generally speaking, <sup>Δ</sup>*ϕ*(*x*) = *<sup>n</sup>* ∑ *i*=1 *∂*2*ϕ ∂x*<sup>2</sup> *i* , for *<sup>x</sup>* = (*x*1, *<sup>x</sup>*2, ··· , *xn*)*<sup>T</sup>* <sup>∈</sup> <sup>R</sup>*n*.

Denote *X* = (*X*1, *X*2, *X*3)*<sup>T</sup>* with *X*<sup>1</sup> = *S*, *X*<sup>2</sup> = *I*, *X*<sup>3</sup> = *R*. The following impulsive epidemic model with a Neumann boundary value is investigated in this paper:

$$\begin{cases} \frac{\partial X(t, \mathbf{x})}{\partial t} = D \Delta X(t, \mathbf{x}) + A(t)X(t, \mathbf{x}) + f(t, X(t, \mathbf{x})), \quad \mathbf{x} \in \Omega, \ t \gg t\_0, \\ X(t\_k^+, \mathbf{x}) - X(t\_k^-, \mathbf{x}) = M\_k X(t\_k - \tau\_k, \mathbf{x}), \quad k \in \mathbb{N}, \ x \in \Omega, \\ \frac{\partial X(t, \mathbf{x})}{\partial \nu} = 0, \quad \mathbf{x} \in \partial \Omega, \ t \gg 0, \\ X(0, \mathbf{x}) = \boldsymbol{\varrho}(\mathbf{x}), \quad \mathbf{x} \in \Omega, \end{cases} \tag{3}$$

where <sup>N</sup> <sup>=</sup> {1, 2, 3, ···}, *<sup>t</sup>*<sup>0</sup> <sup>=</sup> 0 and *tk* is the impulse moment for *<sup>k</sup>* <sup>=</sup> 1, 2, ··· , satisfying 0 < *t*<sup>1</sup> < *t*<sup>2</sup> < ··· < *tk* < ··· and lim *<sup>k</sup>*→<sup>∞</sup> *tk* = +∞. For any impulse moment *tk*, *Mk* is a constant parameters matrix that quantifies the impulse strength at the moment *tk*. The time delay *<sup>τ</sup><sup>k</sup>* <sup>∈</sup> [0, *<sup>τ</sup>*] with (*tk* <sup>−</sup> *<sup>τ</sup>k*, *tk*) <sup>⊂</sup> (*tk*−1, *tk*) for each *<sup>k</sup>* <sup>∈</sup> <sup>N</sup>, *<sup>τ</sup>* <sup>=</sup> sup *<sup>k</sup>*∈<sup>N</sup> *τk*, and so *t*<sup>0</sup> *t*<sup>1</sup> − *τ*1.

<sup>Ω</sup> <sup>⊂</sup> <sup>R</sup>*N*(*<sup>N</sup>* <sup>2</sup>) is a bounded smooth domain with smooth boundary *<sup>∂</sup>*Ω.

Denote *ν* the external normal direction of *∂*Ω.

$$D = \begin{pmatrix} d\_1 & 0 & 0 \\ 0 & d\_2 & 0 \\ 0 & 0 & d\_3 \end{pmatrix}, \quad A(t) = \begin{pmatrix} 0 & 0 & 0 \\ 0 & -\gamma(t) & 0 \\ 0 & \gamma(t) & 0 \end{pmatrix}, \quad f(t, X) = \begin{pmatrix} -\beta(t)X\_1X\_2 \\ \beta(t)X\_1X\_2 \\ 0 \end{pmatrix}. \tag{4}$$

System (3) is the drive system, and its response system can be considered as follows,

$$\begin{cases} \frac{\partial Y(t, \mathbf{x})}{\partial t} = f(t, Y(t, \mathbf{x})) + D\Delta Y(t, \mathbf{x}) + A(t)Y(t, \mathbf{x}), \quad \mathbf{x} \in \Omega, t \gg t\_{0, \mathbf{x}} \\ Y(t\_k^+, \mathbf{x}) = M\_k Y(t\_k - \tau\_k, \mathbf{x}) + Y(t\_k^-, \mathbf{x}), \quad k \in \mathbb{N}, \mathbf{x} \in \Omega, \\ \frac{\partial Y(t, \mathbf{x})}{\partial \nu} = 0, \quad \mathbf{x} \in \partial\Omega, t \gg 0, \\ \phi(\mathbf{x}) = Y(s, \mathbf{x}), \quad \mathbf{x} \in \Omega, s \in [-\tau, 0], \end{cases} (5)$$

and then the error system is proposed as follows,

$$\begin{cases} \frac{\partial \varepsilon(t, \mathbf{x})}{\partial t} = D \Delta \varepsilon(t, \mathbf{x}) + A(t)\varepsilon(t, \mathbf{x}) + F(t, \varepsilon(t, \mathbf{x})), \quad \mathbf{x} \in \Omega, \ t \geqslant t\_{0},\\ \varepsilon(t\_{k}^{+}, \mathbf{x}) - \varepsilon(t\_{k}^{-}, \mathbf{x}) = M\_{k}\varepsilon(t\_{k} - \tau\_{k}, \mathbf{x}), \quad k \in \mathbb{N}, \mathbf{x} \in \Omega, \\ \frac{\partial \varepsilon(t, \mathbf{x})}{\partial \nu} = 0, \quad t \geqslant 0, \mathbf{x} \in \partial \Omega, \\ \varepsilon(0, \mathbf{x}) = \varrho(\mathbf{x}) - \phi(\mathbf{x}), \quad \mathbf{x} \in \Omega, \end{cases} (6)$$

where *e* = *X* − *Y*. Moreover,

$$F(\mathbf{e}(t, \mathbf{x})) = f(t, X(t, \mathbf{x})) - f(t, Y(t, \mathbf{x})) = \begin{pmatrix} -\beta(t)X\_1X\_2 + \beta(t)Y\_1Y\_2 \\ \beta(t)X\_1X\_2 - \beta(t)Y\_1Y\_2 \\ 0 \end{pmatrix} \tag{7}$$

We assume in this paper that variables are left continuous at impulse moment *tk*, for example, *e*(*tk*, *x*) = *e*(*t* − *<sup>k</sup>* , *x*).

Obviously, −1 < *ei* < 1. Moreover, in consideration of the fact that population resources are limited, we can assume throughout this paper that their regional change rate is also limited, and so the change rate of the change rate is even limited:


**Lemma 1** (See, e.g., [21])**.** <sup>Ω</sup> <sup>⊂</sup> <sup>R</sup>*<sup>m</sup> is a bounded domain with its smooth boundary <sup>∂</sup>*<sup>Ω</sup> *that is of class* C2*. <sup>ξ</sup>*(*x*) ∈ *<sup>H</sup>*<sup>1</sup> <sup>0</sup> (Ω) *is a real-valued function and ∂ξ*(*x*) *∂ν* |*∂*<sup>Ω</sup> = 0*. Then,*

$$\int\_{\Omega} |\nabla \xi(\mathbf{x})|^2 d\mathbf{x} \geqslant \lambda\_1 \int\_{\Omega} |\xi(\mathbf{x})|^2 d\mathbf{x} \mathcal{A}$$

*where λ*<sup>1</sup> *is defined by the least positive eigenvalue of the problem:*

$$\begin{cases} \lambda \xi - \Delta \xi = 0, & \mathbf{x} \in \Omega, \\\frac{\partial \xi(\mathbf{x})}{\partial \nu} = 0, & \mathbf{x} \in \partial \Omega. \end{cases}$$

#### **3. Main Result**

**Theorem 1.** *Assume there exists a positive definite diagonal matrix Q and a constant q*<sup>0</sup> > 0 *such that*

$$Q \lessapprox q\_0 \mathcal{Z} \tag{8}$$

*and*

$$\sup\_{k \in \mathbb{N}} \left[ \| \mathcal{Z} + M\_k \| + \tau\_k \| M\_k \| \cdot \left( \| D \| \sqrt{\lambda\_{\text{max}} \mathbb{C}^2} + \| \tilde{A} \| + 2\beta \right) \right] \sqrt{\frac{\lambda\_{\text{max}} Q\_\rho}{\lambda\_{\text{min}} Q}} e^{\lambda \rho} \leqslant \rho\_0 < 1,\tag{9}$$

*then system (3) and system (5) are synchronized, where* I *is the identity matrix, C* = *diag*(*c*1, *c*2, *c*3) > 0 *with ci* > 0 *defined in (H1), ρ* = sup *<sup>k</sup>*∈<sup>N</sup> (−*tk*−<sup>1</sup> + *tk*)*, <sup>ζ</sup>* = inf *<sup>k</sup>*∈<sup>N</sup> (*tk* − *tk*−1) > 0,

$$
\tilde{A} = \begin{pmatrix} 0 & 0 & 0 \\ 0 & \gamma & 0 \\ 0 & \gamma & 0 \end{pmatrix},
\tag{10}
$$

$$\lambda = \left[ \frac{1}{\lambda\_{\text{min}} Q} \lambda\_{\text{max}} \left( -\lambda\_1 (QD + DQ) + Q\tilde{A} + \tilde{A}^T Q + Q + 4q\_0 \beta^2 \mathcal{I} \right) \right]. \tag{11}$$

Here, inequality (8) indicates that (*q*0I − *Q*) is a non-negative definite matrix. For any symmetric matrix *B*, the real numbers *λ*min*B* and *λ*max*B* represent the minimum and maximum eigenvalue of *<sup>B</sup>*, respectively. For a matrix *<sup>B</sup>*, *<sup>B</sup>* is its norm with *<sup>B</sup>* <sup>=</sup> <sup>9</sup>*λ*max(*BTB*).

**Proof.** Consider the following Lyapunov function:

$$V(t) = \int\_{\Omega} e^T(t, \mathbf{x}) Q e(t, \mathbf{x}) d\mathbf{x},$$

where *<sup>Q</sup>* is a positive definite symmetric matrix. Denote *<sup>η</sup>* <sup>2</sup> <sup>=</sup> <sup>3</sup> ∑ *i*=1 <sup>Ω</sup> |*ηi*(*x*)| <sup>2</sup>*dx* for any vector Lebesgue square-integrable function *η*(*x*)=(*η*1(*x*), *η*2(*x*), *η*3(*x*))*T*. It follows from 0 < *Xi* < 1 and 0 < *Yi* < 1 (*i* = 1, 2, 3) that

$$\mathbb{E}F^T(t,\varepsilon)QF(t,\varepsilon) \leqslant 2q\_0\beta^2 \cdot 2[(X\_2 - Y\_2)^2 + (X\_1 - Y\_1)^2] \leqslant 4q\_0\beta^2\varepsilon^T\varepsilon$$

So,

$$\begin{split} \mathcal{D}^{+}V(t) &\leqslant -\lambda\_{1} \int\_{\Omega} \varepsilon^{T}(QD+DQ) \mathrm{e} \mathrm{d}x + \int\_{\Omega} \varepsilon^{T} \Big( Q\tilde{A} + \tilde{A}^{T}Q \Big) \mathrm{e} \mathrm{d}x + \int\_{\Omega} \left( \varepsilon^{T} QF(t,\varepsilon) + F^{T}(t,\varepsilon)Q\varepsilon \right) \mathrm{d}x \\ &\leqslant \frac{1}{\lambda\_{\mathrm{min}}Q} \lambda\_{\mathrm{max}} \Big( -\lambda\_{1}(QD+DQ) + Q\tilde{A} + \tilde{A}^{T}Q + Q + 4q\_{0}\delta^{2}\mathcal{L} \Big) \int\_{\Omega} \varepsilon^{T}(t,\mathbf{x}) Q\varepsilon(t,\mathbf{x}) \mathrm{d}x, \; t \in (t\_{k-1}, t\_{k}]. \end{split}$$

which means

$$||\boldsymbol{\varepsilon}(t)||^{2} \prec\_{\boldsymbol{\varepsilon}} \frac{\lambda\_{\max} \mathbb{Q}}{\lambda\_{\min} \mathbb{Q}} e^{\lambda(t - t\_{k-1})} ||\boldsymbol{\varepsilon}(t\_{k-1}^{+})||^{2}, \quad \boldsymbol{t} \in (t\_{k-1}, t\_{k}].$$

Particularly,

$$\left\|\left\|\boldsymbol{\varepsilon}(t\_k)\right\|\right\|^2 = \left\|\boldsymbol{\varepsilon}(t\_k^-)\right\|^2 \leqslant \frac{\lambda\_{\max} Q}{\lambda\_{\min} Q} e^{\lambda(t - t\_{k-1})} \left\|\boldsymbol{\varepsilon}(t\_{k-1}^+)\right\|^2, \quad k \in \mathbb{N}.$$

On the other hand,

$$\| |D\Delta \mathbf{e}(\varrho\_k, \mathbf{x})| \| \lesssim \| |D| \| \cdot \| \Delta \mathbf{e}(\varrho\_k, \mathbf{x}) \| \lesssim \| |D| \| \sqrt{\lambda\_{\max} C^2} \cdot \| \| \mathbf{e}(\varrho\_k) \| \| \rho$$

and

$$\|\|F(t,e(t))\|\| \leqslant \sqrt{2}\beta \cdot \sqrt{2} \sqrt{\int\_{\Omega} [(X\_2 - Y\_2)^2 + (X\_1 - Y\_1)^2] dx} \leqslant \|e(t)\| \cdot 2\beta.$$

Now we can see it from the differential mean value theorem and definition of *A*(*t*) that there exists *<sup>ς</sup><sup>k</sup>* ∈ (*tk* − *<sup>τ</sup>k*, *tk*) ⊂ (*tk*−1, *tk*) such that

$$\begin{split} \|\boldsymbol{\varepsilon}(t\_{k}^{+})\| &= \|\boldsymbol{\varepsilon}(t\_{k}^{-},\boldsymbol{x}) + M\_{k}\boldsymbol{\varepsilon}(t\_{k} - \boldsymbol{\tau}\_{k},\boldsymbol{x})\| \\ &\leqslant \|\mathbb{Z} + M\_{k}\| \cdot \|\boldsymbol{\varepsilon}(t\_{k},\boldsymbol{x})\| + \|M\_{k}\| \cdot \|\boldsymbol{\varepsilon}(t\_{k},\boldsymbol{x}) - \boldsymbol{\varepsilon}(t\_{k} - \boldsymbol{\tau}\_{k},\boldsymbol{x})\| \\ &\leqslant \|\mathbb{Z} + M\_{k}\| \cdot \|\boldsymbol{\varepsilon}(t\_{k},\boldsymbol{x})\| + \tau\_{k} \|M\_{k}\| \cdot \left( \|\boldsymbol{D}\| \sqrt{\lambda\_{\max}\mathbb{C}^{2}} \cdot \|\boldsymbol{\varepsilon}(\boldsymbol{\varepsilon}\_{k})\| + \|\tilde{A}\| \cdot \|\boldsymbol{\varepsilon}(\boldsymbol{\varepsilon}\_{k})\| + 2\beta\|\boldsymbol{\varepsilon}(\boldsymbol{\varepsilon}\_{k})\| \right) \\ &\leqslant \Big{[} \|\mathbb{Z} + M\_{k}\| + \tau\_{k} \|M\_{k}\| \cdot \left( \|\boldsymbol{D}\| \sqrt{\lambda\_{\max}\mathbb{C}^{2}} + \|\tilde{A}\| + 2\beta\right) \Big{]} \sqrt{\frac{\lambda\_{\max}\mathbb{Q}}{\lambda\_{\min}\mathbb{Q}}} e^{\lambda(t\_{k} - t\_{k-1})} \Big{]} \|\boldsymbol{\varepsilon}(t\_{k-1}^{+})\| \\ &\leqslant \rho\_{0} \|\boldsymbol{\varepsilon}(t\_{k-1}^{+})\|, \quad \forall k = 1,2,3,\cdots \ \vdots \\ &\qquad \text{which means} \end{split}$$

$$\|e(t\_k^+) \prec\_\varepsilon \rho\_0^k\|e(0)\|\_\prime \quad \forall k = 1, 2, 3, \cdots \; . \tag{12}$$

Finally, for *<sup>t</sup>* ∈ (*tk*−1, *tk*],

$$\|\|\boldsymbol{\varepsilon}(t)\|\|^2 \lesssim \frac{\lambda\_{\max} Q}{\lambda\_{\min} Q} e^{\lambda(t - t\_{k-1})} \|\|\boldsymbol{\varepsilon}(t\_{k-1}^+)\|\|^2 \lesssim \frac{\lambda\_{\max} Q}{\lambda\_{\min} Q} e^{\lambda \rho} \rho\_0^{2(k-1)} \|\|\boldsymbol{\varepsilon}(0)\|\|^2 \lesssim \frac{\lambda\_{\max} Q}{\lambda\_{\min} Q} e^{\lambda \rho} \|\|\boldsymbol{\varepsilon}(0)\|\|^2 \lesssim \frac{\lambda\_{\max} Q}{\lambda\_{\max} Q} e^{\lambda \rho} \|\|\boldsymbol{\varepsilon}(0)\|\|^2 \lesssim \frac{\lambda\_{\max} Q}{\lambda\_{\max} Q} e^{\lambda \rho} \|\|\boldsymbol{\varepsilon}(0)\|\|^2$$

which, together with *<sup>t</sup>* > *tk*−1, implies

$$\|\|e(t)\|\| \leqslant \sqrt{\frac{\lambda\_{\text{max}} Q}{\lambda\_{\text{min}} Q}} e^{\lambda \rho} \|\|e(0)\|\| e^{-\lambda\_0(t-t\_0)}.\tag{13}$$

where *<sup>λ</sup>*<sup>0</sup> <sup>=</sup> <sup>−</sup><sup>1</sup> *<sup>ζ</sup>* ln *ρ*<sup>0</sup> > 0 . This completes the proof.

**Remark 1.** *Contrary to the existing literature related to impulsive reaction–diffusion epidemic models (see, e.g., [13,14]), the delayed impulse is firstly considered in the reaction–diffusion epidemic system in this paper. Indeed, although time delays were introduced in [13], the impulse was not delayed. However, in real life, the effectiveness of many defensive measures usually takes place after a period of time. Therefore, the delayed impulse epidemic model studied in this paper clearly has practical significance.*

**Remark 2.** *Introducing the delayed impulse into reaction–diffusion epidemic models means bringing new mathematical difficulties to this paper. Therefore, this paper adopts a method different from [13,14] to overcome the mathematical difficulties, and a new synchronization criterion is derived.*

#### **4. Numerical Example**

**Example 1.** *Let* <sup>Ω</sup> = (0, 1) <sup>×</sup> (0, 1) <sup>⊂</sup> <sup>R</sup>2*, then <sup>λ</sup>*<sup>1</sup> <sup>=</sup> *<sup>π</sup>*2. *Set <sup>C</sup>* <sup>=</sup> <sup>I</sup> <sup>=</sup> *diag*(1, 1, 1)*, <sup>γ</sup>*(*t*) = 0.1 sin2 *t*, *β*(*t*) = 0.1 cos2 *t, and β* = 0.1 = *γ*,

$$A(t) = \begin{pmatrix} 0 & 0 & 0 \\ 0 & -0.1\sin^2 t & 0 \\ 0 & 0.1\sin^2 t & 0 \end{pmatrix}, f(t,X) = \begin{pmatrix} -0.1\cos^2 tX\_1X\_2 \\ 0.1\cos^2 tX\_1X\_2 \\ 0 \end{pmatrix}, \\ \tilde{A} = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0.1 & 0 \\ 0 & 0.1 & 0 \end{pmatrix}.$$

**Case 1**: Let *ρ* = 0.1, *τ<sup>k</sup>* ≡ *τ* = 0.01, *Mk* ≡ −0.7I, and

$$D = \begin{pmatrix} 0.045 & 0 & 0 \\ 0 & 0.035 & 0 \\ 0 & 0 & 0.055 \end{pmatrix}.$$

Using a computer with Matlab's LMI toolbox results in the following feasibility datum:

$$Q = \begin{pmatrix} 0.0155 & 0 & 0\\ 0 & 0.0135 & 0\\ 0 & 0 & 0.0161 \end{pmatrix}$$

then, *q*<sup>0</sup> = *λ*max*Q* = 0.0161, *λ*min*Q* = 0.0135, and

$$\sup\_{k \in \mathbb{N}} \left[ \| \mathcal{Z} + M\_k \| + \pi\_k \| M\_k \| \cdot \left( \| D \| \sqrt{\lambda\_{\max} C^2} + \| \tilde{A} \| + 2 \beta \right) \right] \sqrt{\frac{\lambda\_{\max} Q}{\lambda\_{\min} Q}} \epsilon^{\lambda \rho} = 0.8196 \lesssim \rho\_0 < 1/\rho$$

where *ρ*<sup>0</sup> = 0.8196. According to Theorem 1, system (3) and system (5) are synchronized. **Case 2**: Let *ρ* = 0.05, *τ<sup>k</sup>* ≡ *τ* = 0.01, *Mk* ≡ −0.7I, and

$$D = \begin{pmatrix} 0.045 & 0 & 0 \\ 0 & 0.035 & 0 \\ 0 & 0 & 0.055 \end{pmatrix}.$$

Using Matlab's LMI toolbox results in the following feasibility datum:

$$Q = \begin{pmatrix} 0.0111 & 0 & 0 \\ 0 & 0.0131 & 0 \\ 0 & 0 & 0.0112 \end{pmatrix}$$

then, *q*<sup>0</sup> = *λ*max*Q* = 0.0131, *λ*min*Q* = 0.0111, and

$$\sup\_{k \in \mathbb{N}} \left[ \| \mathcal{Z} + M\_k \| + \tau\_k \| M\_k \| \cdot \left( \| D \| \sqrt{\lambda\_{\max} C^2} + \| \tilde{A} \| + 2 \beta \right) \right] \sqrt{\frac{\lambda\_{\max} Q}{\lambda\_{\min} Q}} e^{\lambda \rho} = 0.7988 \le \rho\_0 < 1.5$$

where *ρ*<sup>0</sup> = 0.7988. According to Theorem 1, system (3) and system (5) are synchronized.

**Remark 3.** *Table 1 reveals that the bigger the impulse frequency, the faster thesynchronization speed.*

**Table 1.** Comparison of the influence from different impulse frequencies when other data are unchanged.


**Case 3**: Let *ρ* = 0.1, *τ<sup>k</sup>* ≡ *τ* = 0.01, *Mk* ≡ −0.8I, and

$$D = \begin{pmatrix} 0.045 & 0 & 0 \\ 0 & 0.035 & 0 \\ 0 & 0 & 0.055 \end{pmatrix}.$$

Using Matlab's LMI toolbox results in the following feasibility datum:

$$Q = \begin{pmatrix} 0.0166 & 0 & 0\\ 0 & 0.0163 & 0\\ 0 & 0 & 0.0165 \end{pmatrix}$$

then, *q*<sup>0</sup> = *λ*max*Q* = 0.0163, *λ*min*Q* = 0.0166, and

$$\sup\_{k \in \mathbb{N}} \left[ \| \mathcal{Z} + M\_k \| + \pi\_k \| M\_k \| \cdot \left( \| D \| \sqrt{\lambda\_{\text{max}} \mathbf{C}^2} + \| \tilde{A} \| + 2\beta \right) \right] \sqrt{\frac{\lambda\_{\text{max}} Q}{\lambda\_{\text{min}} Q}} \varepsilon^{\lambda \rho} = 0.5466 \lesssim \rho\_0 < 1, \rho\_0 \in \mathbb{R}$$

where *ρ*<sup>0</sup> = 0.5466. According to Theorem 1, system (3) and system (5) are synchronized.

**Remark 4.** *Table 2 implies that the bigger the impulse intensity, the faster the synchronization speed.*


**Table 2.** Comparison of the influence from different impulse intensities when other data are unchanged.

**Case 4**: Let *ρ* = 0.1, *τ<sup>k</sup>* ≡ *τ* = 0.001, *Mk* ≡ −0.7I, and

$$D = \begin{pmatrix} 0.045 & 0 & 0 \\ 0 & 0.035 & 0 \\ 0 & 0 & 0.055 \end{pmatrix}.$$

Using Matlab's LMI toolbox results in the following feasibility datum:

$$Q = \begin{pmatrix} 0.0155 & 0 & 0\\ 0 & 0.0135 & 0\\ 0 & 0 & 0.0161 \end{pmatrix}$$

then, *q*<sup>0</sup> = *λ*max*Q* = 0.0161, *λ*min*Q* = 0.0135, and

$$\sup\_{k \in \mathbb{N}} \left[ \| \mathcal{Z} + M\_k \| + \mathfrak{r}\_k \| \| M\_k \| \cdot \left( \| D \| \sqrt{\lambda\_{\max} \mathbf{C}^2} + \| \tilde{A} \| + 2 \mathfrak{k} \right) \right] \sqrt{\frac{\lambda\_{\max} \mathbf{Q}}{\lambda\_{\min} \mathbf{Q}}} \mathbf{e}^{\lambda \rho} = 0.7188 \le \rho\_0 < 1, \rho\_0 \in \mathbb{R}$$

where *ρ*<sup>0</sup> = 0.7188. According to Theorem 1, system (3) and system (5) are synchronized.

**Remark 5.** *Table 3 means that the smaller the time delays of the impulse effect, the faster the synchronization speed.*

**Table 3.** Comparison of the influence of different time delays of the impulse effect when other data are unchanged.


Below, another numerical example is presented to show the validity of Theorem 1 via very simple computations.

**Example 2.** *Set C* = 15I = *diag*(15, 15, 15)*, Mk* ≡ −0.9I*, D* = 0.1I*, Q* = I, *τ<sup>k</sup>* ≡ *τ* = 0.01*, β* = 0.1*, and then direct computations lead to*

$$\sqrt{\lambda\_{\text{max}} C^2} = 15, \; \|\mathbb{Z} + M\_{\text{k}}\| \equiv 0.1,\\ \|M\_{\text{k}}\| \equiv 0.9,\\ \|D\| = 0.1, \lambda\_{\text{max}} Q = \lambda\_{\text{min}} Q = 1 = q\_0. \tag{14}$$

$$\text{Let } \Omega = (0,1) \times (0,1) \subset \mathbb{R}^2; \text{ then, } \lambda\_1 = \pi^2. \text{ Set } \gamma = 0.1; \text{ then,}$$

$$
\check{A} = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0.1 & 0 \\ 0 & 0.1 & 0 \end{pmatrix} \text{ and } \|\check{A}\| = \sqrt{0.0200} = 0.1414. \tag{15}
$$

Hence, it follows from (14) and (15) that

$$\begin{split} \lambda &= \left[ \frac{1}{\lambda\_{\text{min}} Q} \lambda\_{\text{max}} \left( -\lambda\_1 (QD + DQ) + Q\tilde{A} + \tilde{A}^T Q + Q + 4q\rho\theta^2 \mathcal{Z} \right) \right] \\ &= \lambda\_{\text{max}} \left( -\lambda\_1 (QD + DQ) + Q\tilde{A} + \tilde{A}^T Q + Q + 4q\rho\theta^2 \mathcal{Z} \right) \\ &= \lambda\_{\text{max}} \left( -\pi^2 (0.1 \mathcal{Z} + 0.1 \mathcal{Z}) + \tilde{A} + \tilde{A}^T + \mathcal{Z} + 0.04 \mathcal{Z} \right) \\ &= -0.6925, \end{split} \tag{16}$$

where

$$\left(-\pi^2(0.1\mathcal{Z}+0.1\mathcal{Z}) + \tilde{A} + \tilde{A}^T + \mathcal{Z} + 0.04\mathcal{Z}\right) = \begin{pmatrix} -0.9339 & 0 & 0\\ 0 & -0.7339 & 0.1000\\ 0 & 0.1000 & -0.9339 \end{pmatrix}.$$

Now, letting *ρ* = 0.3 and *ρ*<sup>0</sup> = 0.9, we can get by (14)–(16) that

$$\begin{split} & \sup\_{k \in \mathbb{N}} \left[ \| \mathcal{Z} + M\_{k} \| + \tau\_{k} \| \| M\_{k} \| \cdot \left( \| D \| \sqrt{\lambda\_{\text{max}} \mathbf{C}^{2}} + \| \tilde{A} \| \| + 2 \xi \right) \right] \sqrt{\frac{\lambda\_{\text{max}} Q}{\lambda\_{\text{min}} \mathbf{Q}}} e^{\lambda \rho} \\ & \equiv \left[ \| \mathcal{Z} + M\_{k} \| + \tau\_{k} \| M\_{k} \| \cdot \left( \| D \| \sqrt{\lambda\_{\text{max}} \mathbf{C}^{2}} + \| \tilde{A} \| + 2 \xi \right) \right] \sqrt{\frac{\lambda\_{\text{max}} Q}{\lambda\_{\text{min}} \mathbf{Q}}} e^{\lambda \rho} \\ & = \left[ 0.1 + 0.01 \times 0.9 \times \left( 0.1 \times 15 + 0.1414 + 0.2 \right) \right] \sqrt{e^{-0.6925 \times 0.3}} \\ & < \left[ 0.1 + 0.01 \times 0.9 \times \left( 0.1 \times 15 + 0.1414 + 0.2 \right) \right] \times 1 \\ & = 0.1166 < \rho\_{0} < 1, \end{split}$$

which implies (9) holds.

According to Theorem 1, system (3) and system (5) are synchronized (see Figures 1–3).

**Figure 1.** Computer simulation of *X*<sup>1</sup> in (3) and *Y*<sup>1</sup> in (5).

**Figure 2.** Computer simulation of *X*<sup>2</sup> in (3) and *Y*<sup>2</sup> in (5).

**Figure 3.** Computer simulation of *X*<sup>3</sup> in (3) and *Y*<sup>3</sup> in (5).

**Remark 6.** *Example 2 illustrates that the validity of Theorem 1 can be easily verified even without using Matlab's LMI toolbox.*

#### **5. Conclusions**

This paper reported the synchronization control of two epidemic systems with a Neumann boundary value under a delayed impulse. Different from the previous relevant literature in which the effect of the impulse control was immediate, our impulse effect was delayed, which is in line with the actual situation during an epidemic. At the same time, the newly obtained criterion and numerical examples illustrate that the shorter the time delay of the pulse effect, the faster the synchronization speed. In addition, the smaller

the pulse interval, the faster the synchronization. On the other hand, Remarks 1 and 2 illustrated the novelty of this paper by comparing the related literature with this paper.

**Author Contributions:** Writing—original draft and revising, R.R.; participating in the discussion of the literature and polishing English, Z.L., X.A. and J.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the 2019 provincial undergraduate innovation and entrepreneurship training program of Chengdu Normal University (S201914389037).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


### *Article* **Analysis of Equilibrium Points in Quantized Hill System**

**Abdullah A. Ansari 1, Sawsan Alhowaity 2, Elbaz I. Abouelmagd 3,\* and Shiv K. Sahdev <sup>4</sup>**


**Abstract:** In this work, the quantized Hill problem is considered in order for us to study the existence and stability of equilibrium points. In particular, we have studied three different cases which give all whole possible locations in which two points are emerging from the first case and four points from the second case, while the third case does not provide a realistic locations. Hence, we have obtained four new equilibrium points related to the quantum perturbations. Furthermore, the allowed and forbidden regions of motion of the first case are investigated numerically. We demonstrate that the obtained result in the first case is a generalization to the classical one and it can be reduced to the classical result in the absence of quantum perturbation, while the four new points will disappear. The regions of allowed motions decrease as the value of the *Jacobian constant* increases, and these regions will form three separate areas. Thus, the infinitesimal body can never move from one allowed region to another, and it will be trapped inside one of the possible regions of motion with the relative large values of the *Jacobian constant*.

**Keywords:** Hill problem; quantum correction; equilibrium points; stability

**MSC:** 70F05; 70F07; 70F10; 70F15; 70H14

#### **1. Introduction**

Hill's problem is a particularly limiting case for the restricted three-body problem (RTBP). Researchers can obtain the Hill problem by using some scales and transformations while taking limits, as mass parameter tends to zero. Hence, it is an interesting application based problem, and many scientists have studied different versions of this problem by considering different perturbation forces in the classical Hill problem. This means the primary bodies possess point masses and move in circular orbits around their common centre of mass or in elliptical trajectory, while the third body moves in space under the effect of gravitational forces of the primary bodies without affecting their motions [1–4].

In [5], the authors have studied the Hill stability of satellites by utilising the RTBP configuration. However, in [6–8], the authors have studied the same configuration with various perturbations as radiation pressure and oblateness of the primaries. Additionally, in [9–12], the authors have explored and analysed the Hill four bodies problem with its application to the Earth–Moon–Sun system and satellite motion about binary asteroids. In this context, Hill's problem, with oblate secondary in three dimensions, has been illustrated in [13], where the equilibrium points and their stability have been determined.

Further to the precedent work, the radiation pressure effect of the bigger primary and the secondary oblateness on the new version of Hill's problem are investigated in [14],

**Citation:** Ansari, A.A.; Alhowaity, S.; Abouelmagd, E.I.; Sahdev, S.K. Analysis of Equilibrium Points in Quantized Hill System. *Mathematics* **2022**, *10*, 2186. https://doi.org/ 10.3390/math10132186

Academic Editors: Quanxin Zhu

Received: 21 May 2022 Accepted: 21 June 2022 Published: 23 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

where the authors illustrated that their study is more appropriate for astronomical application. They also used iterative methods to identify the locations of equilibrium points and used the linear stability analysis method to examine their stability properties. They proved that all the equilibrium points are unstable for this model. In [15], the authors investigated Hill's problem because space missions required the knowledge of orbits with some properties, where periodic solutions are illustrated numerically due to the non-integrability of this problem.

With the continuous contributions analysing the Hill body problem, the existence of positions and stability of collinear equilibrium points in its generalized version under radiation pressure and oblateness effects are studied in [16]; the authors also performed the basins of attraction through the Newton–Raphson method for many values of used parameters. Furthermore, in [17] the author investigated the basins of convergence in the aforementioned problem; his numerical analysis revealed the extraordinary and beautiful formations on the complex plane. In [18], the authors have performed the Hill's problem by assuming the primaries as the source of radiation pressure; they have determined the asymptotic orbits at collinear points and the same to the lyapunov periodic orbits.

The spatial or planar restricted three-body problem (RTBP) under any kind of perturbation is called the perturbed model. Otherwise, it can be called the phot–gravitational, relativistic, or quantized problem in the case that the system is analysed under the effect of radiation pressure, relativistic, or quantum corrections perturbation, respectively [19–23]. The analysis of the spatial quantized RTBP (i.e., the spatial of RTBP under the effect of quantum corrections) is studied in [24], where the locations of equilibrium points and the allowed and forbidden regions of motions are examined. Furthermore, the quantized RTBP is developed to construct a new version of the Hill problem [25], where the equations of motion for the Hill problem are evaluated under the quantum corrections. Thus, the obtained system is called quantized Hill problem (QHP).

Recently, in [26], the authors investigated the Hill's problem by assuming that the infinitesimal body varies its mass according to Jeans law, they investigated numerically the location of equilibrium points, regions of motion, and basins of attraction and also examined the stability status of these points by using Meshcherskii's space–time transformation. Furthermore, in [27] the authors investigated the differences and similarities among the classical perturbation theory, Poincaré–Lindstedt technique, multiple scales method, the KB averaging method, and averaging theory, while the latter is used to find periodic orbits in the framework of the spatial QHP. They stated that this model can be utilized to develop a lunar theory and families of periodic orbits.

In the framework of RTBP, which can be reduced to the Hill model, some effective contributions are outlined in [28–30], where the effects of lack of sphericity body shape and radiation pressure on the primaries are studied. In addition, the effect of mass variation in the frame of RTBP is investigated in [31–34], where the authors have also studied the impact of these perturbations on the positions of equilibrium points, Poincaré surfaces of section, regions of possible and forbidden motion, and basins of attraction and examined the stability of these equilibrium points such that it is proven that, in most cases, these points are unstable.

In general, the Hill body problem has a great significance in both stellar and solar systems and in dynamical astronomy; it has received a considerable analysis in its own literature. Primarily, it is formulated as a model to analyse the Moon's motion around the Earth under the effect of Sun perturbation. Furthermore, its model, with simple modifications, can also serve as a model for the motion of a star in a star cluster under the created perturbations from the galaxy. The importance of this problem motivated us to study and analyse the Hill body problem under the perturbation of quantum corrections.

In this work, the QHP is considered to study the existence of equilibrium points alongside examining their stability. Under the effect of quantum corrections, the locations of equilibrium points have been analysed. In particular, we have studied three different cases which give all possible locations, where two points are emerging from the first case and they are considered a generalization for the classical two points, as well as four points from the second case, while the third case does not provide any realistic locations. Hence, we have obtained four new equilibrium points related to the quantum perturbations. Furthermore, We demonstrate that the obtained result in the first case can be reduced to the classical result, while the four new points will disappear in the absence of quantum perturbation.

The paper is organized in six sections as follows: The literature surrounding the problem is given in Section 1. The equations of motion are preformed in Section 2. In Section 3, we have determined the positions of equilibrium points. The stability of equilibrium points are studied in Section 4. Furthermore, the numerical results are estimated in Section 5. Finally, the conclusion of the work is presented in Section 6.

#### **2. Equations of Motion**

Following the same notations and procedure in [25], we can write the equations of motion of the quantized Hill problem in the synodic coordinates system as:

$$\begin{aligned} \ddot{x} - 2\dot{y} &= V\_{x'}\\ \ddot{y} + 2\dot{x} &= V\_{y'}\\ \ddot{z} &= V\_{z'} \end{aligned} \tag{1}$$

where

$$V = \frac{1}{2} \left[ 3\,\mathrm{x}^2 + 4(\mathfrak{a}\_1 - \mathfrak{a}\_{11})\mathrm{x} - z^2 \right] + \frac{1}{r} \left( 1 + \frac{\mathfrak{a}\_{21}}{r} + \frac{\mathfrak{a}\_{22}}{r^2} \right) ,\tag{2}$$

and

$$r^2 = x^2 + y^2 + z^2.\tag{3}$$

By integrating System (1), one can write the *Jacobian integral* as

$$
\dot{x}^2 + \dot{y}^2 + \dot{z}^2 = 2\,\text{V} - \text{J}\_{\text{C}} \tag{4}
$$

where *JC* is the *Jacobian constant*.

In System (1), the parameters *α*1, *α*11, and *α*<sup>21</sup> represent very small amounts with order of O(1/*c*2), but *<sup>α</sup>*<sup>22</sup> is of order O(1/*c*3) where *<sup>c</sup>* is the speed of light. Therefore, the value of *α*<sup>1</sup> − *α*<sup>11</sup> will tends to zero [27]. Hence *α*<sup>1</sup> − *α*<sup>11</sup> ∼= 0. In this context, System (1) can be rewritten as:

$$\begin{aligned} \ddot{x} - 2\dot{y} &= x[3 - q(r)], \\ \ddot{y} + 2\dot{x} &= -yq(r), \\ \ddot{z} &= -z[1 + q(r)], \end{aligned} \tag{5}$$

where

$$q(r) = \frac{1}{r^3} \left( 1 + \frac{2\alpha\_{21}}{r} + \frac{3\alpha\_{22}}{r^2} \right)^3$$

We would like to provide the reader with the following investigations about the aforementioned perturbation parameters. In fact, the parameters *α*1, *α*11, and *α*<sup>21</sup> identify the size of the relativistic effect, while *α*<sup>22</sup> estimates the quantum correction contribution. However, all of these effects tend to zero in the case of large distances [21,35].

#### **3. Analysis of Equilibrium Points**

The equilibrium points can be obtained by equating all the derivatives with respect to time by zero in system (5), hence

$$\begin{aligned} \,\_2\mathbf{x}\left[3-q(r)\right]&=0,\\ \,\_2\mathbf{y}\,(r)&=0,\\ \,\_2\mathbf{z}\left[1+q(r)\right]&=0. \end{aligned} \tag{6}$$

In the classical case, we mean that the quantum effect will be neglected, the parameter *α*<sup>21</sup> = *α*<sup>22</sup> = 0, then the equilibrium points are given by (*x*, 0, 0) where *re* = 1/ <sup>√</sup><sup>3</sup> <sup>3</sup> and *<sup>y</sup>* <sup>=</sup> *<sup>z</sup>* <sup>=</sup> 0, hence *<sup>x</sup>* <sup>=</sup> <sup>±</sup>1/√<sup>3</sup> 3.

To find the equilibrium points under the quantized effect (*α*<sup>21</sup> = 0 and *α*<sup>22</sup> = 0), we have to find the solutions of System (6); there are some cases which can be applied to analyse the solutions of this system.

**First case:** *q*(*r*) = 3, then *y* = *z* = 0, and *x* = 0.

In this case, one obtains

$$\frac{1}{r^3} \left( 1 + \frac{2\alpha\_{21}}{r} + \frac{3\alpha\_{22}}{r^2} \right) = 3 \tag{7}$$

Equation (7) gives a quintic equation in the following form

$$2\mathbf{3}^5 - r^2 - 2\mathbf{a}\_{21}\mathbf{r} - \mathbf{3}\mathbf{a}\_{22} = \mathbf{0} \tag{8}$$

The solution of the fifth degree equation is generally too complicated, however the equation has at least one real root. Instead, numerical approximations can be evaluated using a root-finding algorithm for polynomials.

In fact, it is not our aim to find a solution of a quintic equation, but we aim to find the quantum corrections' impact on the locations of equilibrium points. Thus, we impose that *r* = *rq* = *re* + *ε*, where *ε* is a very small quantity which embodies the effect of quantum correction on the locations of equilibrium points after substituting *r* = *rq* = *re* + *ε* into Equation (7) or Equation (8), keeping all terms with coefficients of *ε* and *ε*<sup>2</sup> only, and neglecting all terms with an order of O(*ε*3) or more. Hence, *<sup>ε</sup>* will satisfy two values, *<sup>ε</sup>*<sup>1</sup> and *ε*2, which are given by

$$\varepsilon = \frac{15r\_{\varepsilon}^4 - 2r\_{\varepsilon} - 2a\_{21} \pm \sqrt{4a\_{21}^2 - 12a\_{22} + 180a\_{21}r\_{\varepsilon}^4 + 360a\_{22}r\_{\varepsilon}^3 - 135r\_{\varepsilon}^8 + 72r\_{\varepsilon}^5}}{2\left(1 - 30r\_{\varepsilon}^3\right)},\tag{9}$$

Substituting *re* <sup>=</sup> 1/√<sup>3</sup> 3 in Equation (9), one obtains

$$\varepsilon = \frac{1}{18} \left( 2a\_{21} - 3^{2/3} \pm \sqrt{4a\_{21}^2 + 20\sqrt[3]{9}a\_{21} + 108a\_{22} + 3\sqrt[3]{3}} \right) \tag{10}$$

As *<sup>α</sup>*<sup>21</sup> and *<sup>α</sup>*<sup>22</sup> are very small quantities with order of O(*c*2) and O(*c*23), respectively, we keep only terms with order of O(*α*21) and O(*α*22) and neglect the remanning terms. Thereby, the approximated values of the perturbed parameter *ε* are governed by

$$\begin{aligned} \varepsilon\_{11} &= \frac{1}{3} \left( 2a\_{21} + 3\sqrt[3]{3}a\_{22} \right), \\ \varepsilon\_{12} &= -\frac{1}{3\sqrt[3]{3}} \left( 1 + \frac{4\sqrt[3]{3}}{3}a\_{21} + 3\sqrt[3]{9}a\_{22} \right) \end{aligned} \tag{11}$$

The parameter of *ε* embodies the effect of the quantum corrections and it must equal zero in the absence of these corrections, i.e., when *α*<sup>21</sup> = 0 and *α*<sup>22</sup> = 0. However, the obtained solution of *ε*<sup>12</sup> does not equal zero and gives an inconvenient solution; thus, the value of *ε*<sup>12</sup> is rejected. Hence, the proper approximated value of the parameter *ε* is given by

$$
\varepsilon\_{11} = \frac{1}{3} \left( 2a\_{21} + 3\sqrt[3]{3}a\_{22} \right) \tag{12}
$$

Utilizing Equation (12) with relation to *rq*<sup>1</sup> = *re* + *ε*11, then the distance *rq*<sup>1</sup> at the quantized equilibrium point is

$$r\_{\eta\_1} = \frac{1}{\sqrt[3]{3}} \left( 1 + \frac{2\sqrt[3]{3}}{3} \alpha\_{21} + \sqrt[3]{9} \alpha\_{22} \right) \tag{13}$$

As *x* = |*rq*<sup>1</sup> |, we have two possible values for *x* and *x*<sup>1</sup> = *rq*<sup>1</sup> , *x*<sup>2</sup> = −*rq*<sup>1</sup> . Thus, the quantized equilibrium points are (*x*1, 0, 0) and (*x*2, 0, 0), which is considered a generalization of the classical case and can be reduced to the classical one when *α*<sup>21</sup> = 0 and *α*<sup>22</sup> = 0.

**Second case:** *q*(*r*) = 0, then *x* = *z* = 0, and *y* = 0.

This case could occur when the parameters of quantum corrections are negative, i.e., the values of *α*<sup>21</sup> and *α*<sup>22</sup> are negative [24]. Hence, *q*(*r*) = 0 when the solutions of the following quadratic equation are possible

$$r^2 + 2a\_{21}r + 3a\_{22} = 0.\tag{14}$$

The possible solutions of Equation (14) are

$$\begin{aligned} r\_{\theta\_2} &= -\mathfrak{a}\_{21} - \sqrt{a\_{21}^2 - 3\mathfrak{a}\_{22}} \\ r\_{\theta3} &= -\mathfrak{a}\_{21} + \sqrt{a\_{21}^2 - 3\mathfrak{a}\_{22}} \end{aligned} \tag{15}$$

The solutions in Equation (15) are valid if the values of *rq*<sup>2</sup> and *rq*<sup>3</sup> are positive. To investigate this property, first we remark that *α*21, *α*22, and *α*22/*α*<sup>21</sup> have values with order of O(1/*c*2), O(1/*c*3), and O(1/*c*). Then, the approximated series solutions of Equation (15) can be written as

$$\begin{aligned} r\_{q\_2} &= -2a\_{21} \left[ 1 - \frac{3}{4} \left( \frac{a\_{22}}{a\_{21}} \right) - \frac{9}{16} \left( \frac{a\_{22}}{a\_{21}} \right)^2 \right] + O(\frac{1}{c^5}) \\ r\_{q\_3} &= -\frac{3}{2} a\_{22} \left[ 1 + \frac{3}{4} \left( \frac{a\_{22}}{a\_{21}} \right) \right] + O(\frac{1}{c^5}) \end{aligned} \tag{16}$$

It is clear that from Equation (16) the values of *rq*<sup>2</sup> and *rq*<sup>3</sup> are very small and positive when *α*<sup>21</sup> and *α*22, respectively, take negative values. Then, we have four new equilibrium points corresponding to the second case under the perturbation of quantum corrections, where *y*<sup>2</sup> = |*rq*<sup>2</sup> | and *y*<sup>3</sup> = |*rq*<sup>3</sup> |. The new four points are (0, *y*2, 0), (0, −*y*2, 0), (0, *y*3, 0), and (0, −*y*3, 0)

**Third case:** *q*(*r*) = −1, then *x* = *y* = 0, and *z* = 0.

In this case, one obtains

$$\frac{1}{r^3} \left( 1 + \frac{2\alpha\_{21}}{r} + \frac{3\alpha\_{22}}{r^2} \right) = -1 \tag{17}$$

Equation (17) gives also a quintic equation in the following form

$$r^5 + r^2 + 2a\_{21}r + 3a\_{22} = 0\tag{18}$$

To find the solution of Equation (18), we impose that *r* = *rq* = *re* + *δ*, where *δ* is very small quantity which embodies the effect of quantum correction on the locations of equilibrium points in the current case after substituting *r* = *rq* = *re* + *δ* into Equation (18) and keeping all terms with coefficients of *δ* and *δ*<sup>2</sup> only, neglecting all terms with order of O(*δ*3) or more. Hence, *<sup>δ</sup>* will satisfy two values, *<sup>ε</sup>*<sup>31</sup> and *<sup>ε</sup>*32, which are given by

$$\begin{aligned} \varepsilon\_{31} &= -\frac{3}{26} \left( \frac{11}{3\sqrt[3]{3}} + 2a\_{21} + \sqrt{4a\_{21}^2 - \frac{20a\_{21}}{\sqrt[3]{3}} - 52a\_{22} - \frac{29}{3\sqrt[3]{9}}} \right) \\ \varepsilon\_{32} &= -\frac{3}{26} \left( \frac{11}{3\sqrt[3]{3}} + 2a\_{21} - \sqrt{4a\_{21}^2 - \frac{20a\_{21}}{\sqrt[3]{3}} - 52a\_{22} - \frac{29}{3\sqrt[3]{9}}} \right) \end{aligned} \tag{19}$$

It is clear that from Equation (19) the obtained values of the perturbation parameter *δ* is complex, which mean that the assumption of the third case does not lead to realistic situations. Thus, this case does not give real equilibrium points and it is rejected.

#### **4. Stability Status of Equilibrium Points**

Next, to check the equilibrium points stability, we have to write the equations of motion in to phase space. Thus, System (5) can be rewritten in the following form

$$\begin{aligned} \ddot{x} - 2\dot{y} &= H\_{\text{x}}, \\ \ddot{y} + 2\dot{x} &= H\_{y'} \\ \ddot{z} &= H\_z. \end{aligned} \tag{20}$$

where

$$H = \frac{1}{2} \left[ 3 \,\mathrm{x}^2 - z^2 \right] + \frac{1}{r} \left( 1 + \frac{a\_{21}}{r} + \frac{a\_{22}}{r^2} \right),\tag{21}$$

Here, the *Jacobian integral* can be rewitten as

$$
\dot{x}^2 + \dot{y}^2 + \dot{z}^2 = 2\,\text{H} - \text{J}\_{\text{C}'} \tag{22}
$$

The motion in the proximity of any of the equilibrium points (*a*, *b*, and *c*) can be studied by putting *x* = *a* + *ξ*, *y* = *b* + *η*, and *z* = *c* + *ζ* in Equations (20) and (21). Then, we can rewrite the equations of motion in the phase space as

$$\begin{aligned} \dot{\xi} &= \xi\_1, \\ \dot{\eta} &= \eta\_1, \\ \dot{\xi} &= \zeta\_1, \\ \dot{\xi}\_1 &= 2\eta\_1 + H\_{xx}^0 \xi + H\_{xy}^0 \eta + H\_{xz}^0 \zeta, \\ \dot{\eta}\_1 &= -2\xi\_1 + H\_{yx}^0 \xi + H\_{yy}^0 \eta + H\_{yz}^0 \zeta, \\ \zeta\_1 &= H\_{zx}^0 \xi + H\_{zy}^0 \eta + H\_{zz}^0 \zeta. \end{aligned} \tag{23}$$

where the superscript zero means that the second derivatives of *H* are evaluated at the related equilibrium point.

The characteristic polynomial of Equation (23) will be

$$f(\lambda) = \lambda^6 + H\_5\lambda^5 + H\_4\lambda^4 + H\_3\lambda^3 + H\_2\lambda^2 + H\_1\lambda + H\_0. \tag{24}$$

where

*H*<sup>0</sup> = *H*<sup>0</sup> *x z H*<sup>0</sup> *y y H*<sup>0</sup> *z x* <sup>−</sup> *<sup>H</sup>*<sup>0</sup> *x y H*<sup>0</sup> *y z H*<sup>0</sup> *z x* <sup>−</sup> *<sup>H</sup>*<sup>0</sup> *x z H*<sup>0</sup> *y z H*<sup>0</sup> *z y* + *H*<sup>0</sup> *x x H*<sup>0</sup> *y z H*<sup>0</sup> *z y* + *H*<sup>0</sup> *x y H*<sup>0</sup> *y x H*<sup>0</sup> *z z* <sup>−</sup> *<sup>H</sup>*<sup>0</sup> *x x H*<sup>0</sup> *y y H*<sup>0</sup> *z z*, *H*<sup>1</sup> = 2 (*H*<sup>0</sup> *x z H*<sup>0</sup> *z y* <sup>−</sup> *<sup>H</sup>*<sup>0</sup> *y z H*<sup>0</sup> *z x* <sup>−</sup> *<sup>H</sup>*<sup>0</sup> *x y H*<sup>0</sup> *z z* + *H*<sup>0</sup> *y x H*<sup>0</sup> *z z*), *<sup>H</sup>*<sup>2</sup> <sup>=</sup> <sup>−</sup> <sup>4</sup> *<sup>H</sup>*<sup>0</sup> *z z* <sup>−</sup> *<sup>H</sup>*<sup>0</sup> *x y H*<sup>0</sup> *y x* + *H*<sup>0</sup> *x x H*<sup>0</sup> *y y* <sup>−</sup> *<sup>H</sup>*<sup>0</sup> *x z H*<sup>0</sup> *z x* <sup>−</sup> *<sup>H</sup>*<sup>0</sup> *y z H*<sup>0</sup> *z y* + *H*<sup>0</sup> *x x H*<sup>0</sup> *z z* + *H*<sup>0</sup> *y y H*<sup>0</sup> *z z*, *H*<sup>3</sup> = 2 (*H*<sup>0</sup> *x y* <sup>−</sup> *<sup>H</sup>*<sup>0</sup> *y x*), *<sup>H</sup>*<sup>4</sup> <sup>=</sup> <sup>4</sup> <sup>−</sup> *<sup>H</sup>*<sup>0</sup> *x x* <sup>−</sup> *<sup>H</sup>*<sup>0</sup> *y y* <sup>−</sup> *<sup>H</sup>*<sup>0</sup> *z z*, *H*<sup>5</sup> = 0. (25)

We will examine the stability of equilibrium points in two cases only because there are no equilibrium points in the third case.

#### *4.1. First Case*

In this case the equilibrium points (*x*1, 0, 0) and (*x*2, 0, 0) are in symmetry about the *Y*-axis, therefore it is enough to examine the stability of only one of these two points. In this context, we have to evaluate the values of *Hi*<sup>1</sup> corresponding to (*x*1, 0, 0), which are as follows:

$$\begin{aligned} H\_{01} &= -\frac{3}{\mathbf{x}\_1^3} - \frac{5}{\mathbf{x}\_1^6} - \frac{2}{\mathbf{x}\_1^4} \left( 3 + \frac{11}{\mathbf{x}\_1^3} + \frac{7}{\mathbf{x}\_1^6} \right) \\ &- \frac{3\alpha\_{22}}{\mathbf{x}\_1^5} \left( 3 + \frac{12}{\mathbf{x}\_1^4} + \frac{8}{\mathbf{x}\_1^6} \right) - \frac{2}{\mathbf{x}\_1^9}, \\ H\_{11} &= 0, \\ H\_{21} &= 1 - \frac{3}{\mathbf{x}\_1^3} - \frac{8\alpha\_{21}}{\mathbf{x}\_1^4} \left( 1 + \frac{2}{\mathbf{x}\_1^3} \right) \\ &- \frac{15\alpha\_{22}}{\mathbf{x}\_1^5} \left( 1 + \frac{2}{\mathbf{x}\_1^3} \right) - \frac{3}{\mathbf{x}\_1^6}, \\ H\_{31} &= 0, \\ H\_{41} &= 2 - \frac{2\alpha\_{21}}{\mathbf{x}\_1^4} - \frac{6\alpha\_{22}}{\mathbf{x}\_1^5}, \\ H\_{51} &= 0. \end{aligned} \tag{26}$$

From Equations (24) and (26), we find

$$f(\lambda) = \lambda^6 + H\_{41}\lambda^4 + H\_{21}\lambda^2 + H\_{01\prime} \tag{27}$$

Here, *H*<sup>41</sup> > 0, *H*<sup>21</sup> < 0, and *H*<sup>01</sup> < 0 show that the sign changes occur one at a time time, thus there exists at least one positive real root. Therefore, the equilibrium point will be unstable in this case.

#### *4.2. Second Case*

In this case, the equilibrium points (0, *y*2, 0) and (0, *y*3, 0) are symmetrical about the *X*-axis, hence it is sufficient to examine the stability of only two of these four points. Additionally, we have to evaluate the values of *Hi*<sup>2</sup> and *Hi*3, *i* = 0, 1, 2, 3, 4, and 5 corresponding to (0, *y*2, 0) and (0, *y*3, 0), which are as follows:

$$\begin{aligned} H\_{02} &= \frac{6}{y\_2^5} + \frac{4}{y\_2^6} + \frac{2}{y\_2^4} \left( 9 + \frac{10}{y\_2^5} - \frac{7}{y\_2^6} \right) \\ &+ \frac{12a\_{22}}{y\_2^5} \left( 3 + \frac{3}{y\_2^3} - \frac{2}{y\_2^6} \right) - \frac{2}{y\_2^9}, \\ H\_{12} &= 0, \\ H\_{22} &= 1 + \frac{6}{y\_2^2} + \frac{16a\_{21}}{y\_2^4} \left( 1 - \frac{1}{y\_2^3} \right) \\ &+ \frac{30a\_{22}}{y\_2^5} \left( 1 - \frac{1}{y\_2^5} \right) - \frac{3}{y\_2^6}, \\ H\_{32} &= 0, \\ H\_{42} &= 2 - \frac{2}{y\_2^4} - \frac{6}{y\_2^5}, \\ H\_{52} &= 0, \end{aligned} \tag{28}$$

and

$$\begin{aligned} H\_{03} &= \frac{6}{y\_3^5} + \frac{4}{y\_3^6} + \frac{2}{y\_3^4} \left( 9 + \frac{10}{y\_3^5} - \frac{7}{y\_3^6} \right) \\ &+ \frac{12 \, a\_{22}}{y\_3^5} \left( 3 + \frac{3}{y\_3^4} - \frac{2}{y\_3^6} \right) - \frac{2}{y\_3^9}, \\ H\_{13} &= 0, \\ H\_{23} &= 1 + \frac{6}{y\_3^5} + \frac{16 \, a\_{21}}{y\_3^4} \left( 1 - \frac{1}{y\_3^3} \right) \\ &+ \frac{30 \, a\_{22}}{y\_3^5} \left( 1 - \frac{1}{y\_3^7} \right) - \frac{3}{y\_3^6}, \\ H\_{33} &= 0, \\ H\_{43} &= 2 - \frac{2 \, a\_{21}}{y\_3^4} - \frac{6 \, a\_{22}}{y\_3^5}, \\ H\_{53} &= 0. \end{aligned} \tag{29}$$

From Equations (24) and (28), we find

$$f(\lambda) = H\_{6k}\lambda^6 + H\_{4k}\lambda^4 + H\_{2k}\lambda^2 + H\_{0k} \tag{30}$$

Here, *H*6*<sup>k</sup>* = 1 > 0, *H*4*<sup>k</sup>* < 0, *H*2*<sup>k</sup>* < 0, and *H*0*<sup>k</sup>* < 0, where *k* = 2, 3, show that the sign changes occur one at a time, thus there exists at least one positive real root. Therefore, the equilibrium point will be unstable in this case.

#### **5. Numerical Results**

In this section, we illustrate some dynamical properties numerically for the proposed system (i.e., the quantized Hill system) such as the equilibrium points and the allowed and forbidden regions of motion under the quantum corrections. In order to avoid the reparation, we will present the numerical analysis on the first case of equilibrium points; the same procedure can be carried out for the second case.

The locations of equilibrium points are shown in Figure 1, for which we have taken zero as the derivatives with respect to time in Equation (5). Then, with the help of the well-known *Mathematica Software*, the collinear equilibrium points *L*<sup>1</sup> and *L*<sup>2</sup> under the quantum corrections, as well as the unperturbed equilibrium points *L*¯ <sup>1</sup> and *L*¯ 2, are estimated numerically. Both points exist either side of the origin on the *X*-axis and are in symmetry about *Y*-axis. However, we mark that the distance between the perturbed points is more than the distance between the unperturbed points. Of course, this perturbation will affect the other dynamical properties.

**Figure 1.** Locations of equilibrium points.

One of the most dynamical properties which can be identified by the *Jacobian integral* is the possible and forbidden regions of infinitesimal body motions, which are restricted to the locations of *<sup>v</sup>*<sup>2</sup> = <sup>2</sup>*<sup>H</sup>* − *CJ* ≥ 0 where *<sup>v</sup>* is the velocity of the infinitesimal body. Hence, Equation (22) can be used to determine the allowed or forbidden regions of motions, as in Figure 2, where the coloured green areas identify the regions of possible motions, while the white determine forbidden regions.

It is clear from Figure 2a that when the *Jacobian constant* is relatively small there is one large area for possible region of motion, and the body could move from any region point to another (or from *L*<sup>1</sup> (*L*¯ 1) to *L*<sup>2</sup> (*L*¯ 2)). When *CJ* becomes larger, the forbidden region is extended, as in Figure 2b. With further increase in the value of *CJ*, the forbidden region becomes larger, while the possible region of motion forms three septate areas starting from the perturbed equilibrium points *L*<sup>1</sup> and *L*2, as in Figure 2c. In addition, the body cannot move from one to another, because the three areas are not connected. With further increase in the value of *CJ*, the inner and two outer regions decrease while the separate areas start from the unperturbed equilibrium points *L*¯ <sup>1</sup> and *L*¯ 2, as in Figure 2d. We remark that the infinitesimal can never move from one allowed region to another, and the body will be trapped inside one of the possible regions of motion with the relative large values of the *Jacobian constant*, as in the case of Figure 2c,d.

The condition of *<sup>v</sup>*<sup>2</sup> ≥ 0 does not provide information about the size or shape of the orbit or the trajectory of the body; it can only identify the region where the infinitesimal body could move.

**Figure 2.** Regions of allowed (green area) and forbidden (white area) motion.

#### **6. Conclusions**

In this work, the quantized Hill problem is considered to study the existence of equilibrium points alongside examining their stability. Under the effect of quantum corrections, the locations of equilibrium points have been analysed, we have studied three different cases which give all possible locations, where two points emerge from the first case, taking a place on the *X*-axis, and four points dos so from the second case and lie on *Y*-axis. The third case does not provide a realistic location. Hence, we have obtained four new equilibrium points related to the quantum perturbations.

In this context, we have tested the stability status of all of the equilibrium points and we have found that all points are unstable. Further, we have illustrated the locations of equilibrium points for the first case and the related allowed regions of motion numerically. Similarly, we can perform these illustrations for the second case. Here, we found two equilibrium points which are either side of the origin on the *X*-axis and in symmetry about the *Y*-axis, as in Figure 1. The regions of possible and forbidden motion are investigated for different values of *Jacobian constant*, as in Figure 2.

Finally, we demonstrate that the obtained result in the first case is a generalization of the classical one, and it can be reduced to the classical result, while the four new points will disappear in the absence of quantum perturbation. The regions of possible motions decrease with the increasing value of *Jacobian constant* and these regions will form three separate areas. Thus, the infinitesimal body can never move from one allowed region to another, and it will be trapped inside one of the possible regions of motion with the relative large values for the *Jacobian constant*.

**Author Contributions:** Formal analysis, A.A.A., S.A., E.I.A. and S.K.S.; investigation, A.A.A., S.A., E.I.A. and S.K.S.; methodology, A.A.A., S.A., E.I.A. and S.K.S.; project administration, E.I.A.; software, E.I.A.; Validation, A.A.A., S.A., E.I.A. and S.K.S.; visualization, E.I.A.; Writing—original draft, E.I.A.; and writing—review and editing, A.A.A., S.A., E.I.A. and S.K.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by the National Research Institute of Astronomy and Geophysics (NRIAG), Helwan 11421, Cairo, Egypt. The third author, therefore, acknowledges his gratitude for NRIAG's technical and financial support. Moreover, this paper was supported by the National Natural Science Foundation of China (NSFC), grant no. 12172322.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The study does not report any data.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Mathematics* Editorial Office E-mail: mathematics@mdpi.com www.mdpi.com/journal/mathematics

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34

www.mdpi.com ISBN 978-3-0365-6323-7