Grotunits=360

# Bifurcation of the traveling wave solutions in a perturbed
$(1+1)$-dimensional dispersive long wave equation via a geometric approach

Hang Zheng1,2 Yonghui Xia1111Corresponding author: Email<EMAIL_ADDRESS><EMAIL_ADDRESS>
1\. Department of Mathematics, Zhejiang Normal University, Jinhua, 321004,
China
<EMAIL_ADDRESS><EMAIL_ADDRESS>
2\. Department of Mathematics and Computer, Wuyi University, Wuyishan, 354300,
China
<EMAIL_ADDRESS><EMAIL_ADDRESS>

###### Abstract

Choosing ${\kappa}$ (horizontal ordinate of the saddle point associated to the
homoclinic orbit) as bifurcation parameter, bifurcations of the travelling
wave solutions is studied in a perturbed $(1+1)$-dimensional dispersive long
wave equation. The solitary wave solution exists at a suitable wave speed $c$
for the bifurcation parameter
${\kappa}\in(0,1-\frac{\sqrt{3}}{3})\cup(1+\frac{\sqrt{3}}{3},2)$, while the
kink and anti-kink wave solutions exist at a unique wave speed
$c^{*}=\sqrt{15}/3$ for $\kappa=0$ or $\kappa=2$. The methods are based on the
geometric singular perturbation (GSP, for short) approach, Melnikov method and
invariant manifolds theory. Interestingly, not only the explicit analytical
expression of the complicated homoclinic Melnikov integral is directly
obtained for the perturbed long wave equation, but also the explicit
analytical expression of the limit wave speed is directly given. Numerical
simulations are utilized to verify our mathematical results.

Keywords: Wave equation, solitary wave solutions, geometric singular
perturbation.

MSC2020: 35Q35; 35L05;74J30;34D15;

##  1 Introduction

###  1.1 Model formulation

The water wave equations (e.g., KdV equation, BBM equation and CH or rotation-
CH equation) have attracted a lot of scholars’ attentions. These nonlinear
wave equations have been extensively used to describe dynamical behaviour of
nonlinear waves in shallow water. Since the interaction of nonlinear and
dispersion factors, long wave in shallow water would admit many special
characteristics. One of the important property of long waves, is that, they
retain their shapes and forms after mutual interactions and collisions. In
1996, Wu and Zhang [37] derived an equation describing nonlinear dispersive
long gravity waves travelling in two horizontal directions on shallow waters
of uniform depth format formulated as

$\begin{cases}\theta_{t}+\theta\theta_{x}+u\theta_{y}+v_{x}=0,\\\ u_{t}+\theta
u_{x}+uu_{y}+v_{y}=0,\\\ v_{t}+(\theta
v)_{x}+(uv)_{y}+\frac{1}{3}(\theta_{xxx}+\theta_{xyy}+u_{xxy}+u_{yyy})=0.\end{cases}$
(1)

where $\theta$ (resp., $u$) is the surface velocity of water along the $x$
(resp., $y$) direction and $v$ is the elevation of the water wave. For
dispersive long wave equation (1), Chen, Tang and Lou [17] obtained a special
type of multisoliton solution by the Weiss-Tabor-Carnvale Painlevé truncation
expansion. Eq. (1) reduces to the following $(1+1)$-dimensional dispersive
long wave equation (for short, DLWE) by symmetry reduction and scale
transition:

$\begin{cases}u_{t}=-uu_{x}-v_{x},\\\
v_{t}=-vu_{x}-uv_{x}-\frac{1}{3}u_{xxx}.\end{cases}$ (2)

Due to it models the nonlinear water wave availably, $(1+1)$-dimensional DLWE
(2) is often applied to coastal design and harbor construction. So far, a lot
of articles have been concerned with the exact solutions of
$(1+1)$-dimensional DLWE (2) by using various methods, such as a modified
Conte’s invariant Painlevé expansion (Chen and Lou [6]), a new Jacobi elliptic
function rational expansion method (Wang, Chen and Zhang [36]), a new general
algebraic method with symbolic computation (Chen and Wang [10]), a generalized
extended rational expansion method (Zeng and Wang [47]) and extended tanh-
function method (Fan [17]).

However, in the real world application, particularly for nonlinear wave, there
are many factors of uncertainty and unpredictable perturbations. Thus, small
perturbation terms are usually added to describe such unpredictable
perturbations when modelling the nonlinear waves. For shallow water wave
equation, weak backward diffusion $u_{xx}$ and dissipation $u_{xxxx}$ are
called Kuramoto-Sivashinsky (KS, for short) perturbation. However, to the best
of our knowledge, there are no papers considering the perturbed
$(1+1)$-dimensional DLWE. In this paper, we study a perturbed
$(1+1)$-dimensional DLWE described by

$\begin{cases}u_{t}=-uu_{x}-v_{x},\\\
v_{t}=-vu_{x}-uv_{x}-\frac{1}{3}u_{xxx}+\varepsilon(u_{xx}+u_{xxxx}),\end{cases}$
(3)

where $0<\varepsilon\ll 1$ (a sufficiently small parameter), $u_{xx}$ is
backward diffusion and $u_{xxxx}$ represents a dissipation term.

###  1.2 Methods, motivation and contributions

Note that the perturbed $(1+1)$-dimensional DLWE (3) is singular perturbation
system due to the introduction of $\varepsilon u_{xxxx}$. Considering the
geometric singular perturbations in the evolution equation is a effective
method to describe the real situation for the nonlinear waves. One of
important themes of shallow water waves is to study the existence of traveling
wave solutions. In this paper, we will apply the geometric singular
perturbation theory to study the traveling wave solutions of the singular
perturbed DLWE (3).

In fact, the theory of the geometric singular perturbation have been well
developed by (Fenichel [18], Jones [23], Szmolyan [30], Guckenheimer and
Holmes [19], Wiggins [33], Zhang [48] and Li et al. [24, 25, 26]). And it is a
powerful tool to solve the applications arising in nonlinear waves, biological
systems and other dynamic systems. Up till now, there are an extensive
literature applying GSP theory to study such systems, including delayed CH
equation (Du et al. [16]), CH Kuramoto-Sivashinsky equation (Du et al. [15]),
generalized CH equation (Qiu et al. [29]), perturbed BBM equation (Chen et al.
[3]), generalized BBM equation (Sun et al. [31], Zhu et al. [45]), the delayed
DP equation (Cheng and Li [7]), KdV equation (Derks and Gils [11], Ogama
[28]), generalized KdV equation (Chen et al. [4], Yan et al. [42], Zhang et
al. [43]), biological model (Chen and Zhang [9], Wang and Zhang [38]),
reaction-diffusion equation (Yang and Ni [41], Wu and Ni [39]), Belousov-
Zhabotinskii system (Du and Qiao [13]), MEMS model (et al. Iuorio et al [22]),
piecewise-smooth dynamical systems (Buzzi et al. [1]), perturbed Gardner
equation (Wen [34]), delayed Schrödinger equation (Xu et al. [40]), perturbed
mKdV and mK$(3,1)$ equation (Zhang et al. [43, 44]), generalized Keller-Segel
system (Du et al. [14]), and so on.

To apply the geometric singular perturbation approach to track invariant
manifolds of corresponding ordinary differential equations (ODEs), usually, it
is connected with the zeroes of the Melnikov functions associated to the
perturbed ODEs. However, it is not easy to determine the zeroes of the
Melnikov functions because of the complexity of the expression. To analyze the
Melnikov functions, usually, some auxiliary tools are needed. For examples,
one of the effective method is to detect the monotonicity of the ratio of
Abelian integral (for short, MRAI) associated to the Melnikov function. For
instance, Derks and Gils [11], Ogawa [28] computed the MRAI of perturbed KdV
equation. Chen et al. [3, 4, 2] detected the MRAI of perturbed BBM equation,
perturbed defocusing mKdV equation and perturbed generalized KdV equation by
the similar method. Du et al. [12] considered MRAI of a generalized Nizhnik-
Novikov-Veselov equation with diffusion term. Then they all proved the
existence of traveling wave solutions. Different from them, Sun et al. [32]
employed the Chebyshev criterion to detect the MRAI of a shallow water fluid
to analyze the coexistence of the solitary and periodic wave solutions. Sun
and Yu [31] also illustrated the existence and uniqueness of periodic waves of
a generalized BBM equation based on the same technique.

Different from applying the method of MRAI ([3, 4, 2, 12, 11, 28, 32, 31]) and
the method of Fredholm orthogonality (Du and Qiao [13]), in this paper, we
directly compute the explicit expression of Melnikov function to determine the
zeroes of the Melnikov function. The main purpose of this paper is to
investigate the existence of solitary and kink (anti-kink) wave solutions for
$(1+1)$-dimensional DLWE under a small singular perturbation. Firstly, we
introduce a new parameter ${\kappa}$ (horizontal ordinate of the saddle point
associated to the homoclinic orbit) to obtain the solitary wave solutions and
kink (anti-kink) wave solutions of the unperturbed $(1+1)$-dimensional DLWE.
Then, the geometric singular perturbation approach is used to construct
homoclinic or heteroclinic orbits by tracking invariant manifolds of
corresponding ODEs. Finally, by Melnikov method, we analyze conditions such as
the wave speed and parameter for the existence of solitary and kink (anti-
kink) wave solutions. The contributions and novelty of this paper is
summarized as follows:
(1) By GSP approach and the bifurcation analysis, we prove that the solitary
wave solution exists at a suitable wave speed $c$ for the bifurcation
parameter ${\kappa}\in(0,1-\frac{\sqrt{3}}{3})\cup(1+\frac{\sqrt{3}}{3},2)$,
while the kink and anti-kink wave solutions exist at a unique wave speed
$c^{*}=\sqrt{15}/3$ for $\kappa=0$ or $\kappa=2$. This indicates that the
solitary and kink (anti-kink) wave solutions can’t coexist at the same
bifurcation parameter value in $(1+1)$-dimensional DLWE.
(2) To prove the existence of solitary and kink (anti-kink) wave solutions, we
use the analytical Melnikov method and bifurcation theory. Not only the
analytical expression of Melnikov integral is directly obtained by hand (not
by mathematical software) for a perturbed PDE, but also the analytical
expression of the limit wave speed $c$ is directly obtained by hand (see
(56)-(40) in Section 5). Indeed, the calculations and expressions are verified
by the mathematical softwares (e.g., Maple).
(3) To obtained the analytical expressions, we use the following two
mathematical skills.
(i) We apply a factorization technique together with the handbook of integral
[20] to the calculate $I(\kappa)$ and $J(\kappa)$, which lead us to obtain
analytical Melnikov integral in section 5. We perform very detailed
computations and skillful mathematical analysis to obtain the expressions of
$I(\kappa)$ and $J(\kappa)$ (see (48)–(55) in Section 5).
(ii) $I{({\kappa})}$ and $J{({\kappa})}$ are integrated over a closed curve.
Usually, it is difficult to be calculated. We transfer them to an definite
integral by translating the time variable $\zeta$ into the state variable $z$
on homoclinic orbit. Then the up and lower limits of the definite integral
depend on the coordinates of the intersections between the homoclinic orbit
and $z$ axis. For example,

$I({\kappa})=\displaystyle\oint_{\Gamma({\kappa})}y^{2}d\zeta=\displaystyle\oint_{\Gamma({\kappa})}ydz=2\displaystyle\int_{z_{+}}^{{\kappa}}f(z)dz,$

where $f(z)$ is the expression related to the homoclinic orbit (see (48) in
Section 5).

(4) We introduce a new bifurcation parameter ${\kappa}$ (horizontal ordinate
of the saddle point associated to the homoclinic orbit) instead of the
integral constant $G$. This benefits to the factorization of the integrand
because there is a factor $(z-{\kappa})^{2}$. For example, when we compute
$I({\kappa})$, it will lead us to

$I({\kappa})=\cdots=\displaystyle\int_{z_{+}}^{{\kappa}}\sqrt{\frac{3}{4}(z-{\kappa})^{2}(z-z_{+})(z-z_{-})}dz=\cdots$

(see (48) in section 5).

(5) Numerical simulations are utilized to verify the theoretical results.

###  1.3 Outline of the paper

The outline of the paper is as follows. Some preliminaries including geometric
singular perturbation theory are introduced in section 2. Reduction of the
model by geometric singular perturbation theory and dynamical system method is
presented in Section 3. By introducing a bifurcation parameter ${\kappa}$, the
solitary and kink (anti-kink) wave solutions are determined in section 4. In
section 5, combining the GSP approach and Melnikov method, the existence of
solitary and kink (anti-kink) wave solutions for a perturbed
$(1+1)$-dimensional DLWE is investigated. In section 6, numerical simulations
are carried out to show the effectiveness of previous theoretical results.

##  2 Preliminaries

In this section, we introduce some known results on the theory of geometric
singular perturbation (see. e.g., [23]). Consider the system

$\begin{cases}x_{1}^{\prime}=f(x_{1},x_{2},\epsilon),\\\
x_{2}^{\prime}=\epsilon g(x_{1},x_{2},\epsilon),\end{cases}$ (4)

where ${}^{\prime}=\frac{d}{dt}$, $x_{1}\in\mathbb{R}^{n}$,
$x_{2}\in\mathbb{R}^{l}$ and $\epsilon$ is a positive real parameter,
$U\subseteq\mathbb{R}^{n+l}$ is open subset, and $I$ is an open subset of
$\mathbb{R}$, containing $0$. $f$ and $g$ are $C^{\infty}$ on a set $U\times
I$. Moreover, the $x_{1}$ (resp., $x_{2}$) variables are called fast (resp.,
slow) variables. Letting $\tau=\epsilon t$ which gives the following
equivalent system

$\begin{cases}\epsilon\dot{x_{1}}=f(x_{1},x_{2},\epsilon),\\\
\dot{x_{2}}=g(x_{1},x_{2},\epsilon),\end{cases}$ (5)

where $\cdot=\frac{d}{d\tau}$. We refer to $t$ (resp., $\tau$ ) as the fast
time scale or fast time (resp., slow time scale or slow time). Each of the
scalings is naturally associated with a limit as $\epsilon$ tend to zero.
These limits are respectively given by

$\begin{cases}x_{1}^{\prime}=f(x_{1},x_{2},0),\\\
x_{2}^{\prime}=0,\end{cases}$ (6)

and

$\begin{cases}0=f(x_{1},x_{2},0),\\\
x_{2}^{\prime}=g(x_{1},x_{2},0).\end{cases}$ (7)

System (6) is called the layer problem and system (7) is reduced system.

###### Definition 1

(see [23, 16, 7]) A manifold $M_{0}$ on which $f(x_{1},x_{2},0)=0$ is called a
critical manifold or slow manifold. A critical manifold $M_{0}$ is said to be
normally hyperbolic if the linearization of system (4) at each point in
$M_{0}$ has exactly $l$ eigenvalues on the imaginary axis $Re(\lambda)=0$.

###### Definition 2

(see [23, 16, 7]) A set $M$ is locally invariant under the flow of system (4)
if it has neighborhood $V$ so that no trajectory can leave $M$ without also
leaving $V$. In other words, it is locally invariant if for all $x_{1}\in M$,
$x_{1}\cdot[0,t]\subseteq V$ implies that $x_{1}\in M$,
$x_{1}\cdot[0,t]\subseteq M$, similarly with $[0,t]$ replaced by $[t,0]$ when
$t<0$, where $x_{1}\cdot[0,t]$ denotes the application of a flow after time
$t$ to the initial condition $x_{1}$.

###### Lemma 1

(see [18, 16, 7]) Let $M_{0}$ be a compact, normally hyperbolic critical
manifold given as a graph $\\{(x_{1},x_{2}):x_{1}=h^{0}(x_{2})\\}$. Then for
sufficiently small positive $\epsilon$ and any $0<r<+\infty$,

$\bullet$ there exists a manifold $M_{\epsilon}$, which is locally invariant
under the flow of system (4) and $C^{r}$ in $x_{1}$, $x_{2}$, $\epsilon$.
Moreover, $M_{\epsilon}$ is given as graph:

$M_{\epsilon}=\\{(x_{1},x_{2}):x_{1}=h^{\epsilon}(x_{2})\\}$

for some $C^{r}$ function $h^{\epsilon}(x_{2})$;

$\bullet$ $M_{\epsilon}$ possesses locally invariant stable and unstable
manifold $W^{s}(M_{\epsilon})$ and $W^{u}(M_{\epsilon})$ lying within
$O(\epsilon)$ and being $C^{r}$ diffeomorphic to the stable and unstable
manifold $W^{s}(M_{0})$ and $W^{u}(M_{0})$ of the critical manifold $M_{0}$;

$\bullet$ $W^{s}(M_{\epsilon})$ is partitioned by moving invariant
submanifolds $\digamma^{s}(p_{\epsilon})$, which are $O(\epsilon)$ close and
diffeomorphic to $\digamma^{s}(p_{0})$, with base point $p_{\epsilon}$
belonging to $M_{\epsilon}$. Moreover, they are $C^{r}$ with respect to $p$
and $\epsilon$. Moving invariance means the submanifold
$\digamma^{s}(p_{\epsilon})$ is mapped under the time $t$ flow to another
submanifold $\digamma^{s}(p_{\epsilon}\cdot t)$ whose base point is the time
$t$ evolution image of the taken base point $p_{\epsilon}$;

$\bullet$ the dynamics on $M_{\epsilon}$ is a regular perturbation of that
generated by system (7).

##  3 Reduction of the model by geometric singular perturbation theory and
dynamical system method

In this section, to consider the exact solutions of the unperturbed
$(1+1)$-dimensional DLWE (3), we firstly reduce the model by GSP theory and
dynamical system method. From the view of physical meanings, we only consider
the case of $c>0$ in this paper.

By introducing the following transformations:

$u(x,t)=\phi(\xi),\ \ \ v(x,t)=\psi(\xi),\ \ \ \xi=x-ct,$ (8)

we obtain

$\begin{array}[]{l}\frac{\partial u(x,t)}{\partial t}=-c\phi_{\xi},\ \
\frac{\partial v(x,t)}{\partial t}=-c\psi_{\xi},\\\
\frac{\partial^{2}u(x,t)}{\partial x^{2}}=\phi_{\xi\xi},\ \
\frac{\partial^{3}u(x,t)}{\partial x^{3}}=\phi_{\xi\xi\xi},\ \
\frac{\partial^{4}u(x,t)}{\partial x^{4}}=\phi_{\xi\xi\xi\xi},\end{array}$ (9)

where $\phi_{\xi}$ and $\psi_{\xi}$ are the first order derivative with
respect to $\xi$. And $\phi_{\xi\xi}$, $\phi_{\xi\xi\xi}$ and
$\phi_{\xi\xi\xi\xi}$ means the second, third and fourth order derivative with
respect to $\xi$, respectively.

Substituting (8) and (9) into (3), then system (3) is given by

$\begin{cases}c\phi_{\xi}=\phi\phi_{\xi}+\psi_{\xi},\\\
c\psi_{\xi}=\psi\phi_{\xi}+\phi\psi_{\xi}+\frac{1}{3}\phi_{\xi\xi\xi}-\varepsilon(\phi_{\xi\xi}+\phi_{\xi\xi\xi\xi}).\end{cases}$
(10)

Integrating both sides of the first equation of system (10) once with respect
to $\xi$ and letting the integration constant be zero yields

$\psi(\xi)=c\phi(\xi)-\frac{1}{2}\phi^{2}(\xi).$ (11)

We substitute (11) into the second equation of (10), then the coupled system
(10) becomes the following ODE:

$\frac{1}{3}\phi_{\xi\xi\xi}-\frac{3}{2}\phi^{2}\phi_{\xi}+3c\phi\phi_{\xi}-c^{2}\phi_{\xi}-\varepsilon(\phi_{\xi\xi}+\phi_{\xi\xi\xi\xi})=0.$
(12)

Integrating both sides on (12) once with respect to $\xi$ and rescaling
$\epsilon=3\varepsilon$, it follows that

$\phi_{\xi\xi}-\frac{3}{2}\phi^{3}+\frac{9}{2}c\phi^{2}-3c^{2}\phi-\epsilon(\phi_{\xi}+\phi_{\xi\xi\xi})=g,$
(13)

where $g$ is an integration constant ($g\in\mathbb{R}$).

Introducing new variables $\zeta=c\xi$, $z=\frac{\phi}{c}$ and
$G=\frac{g}{c^{3}}$, Eq. (13) is equivalent to

$-3z+\frac{9}{2}z^{2}-\frac{3}{2}z^{3}+\frac{d^{2}z}{d\zeta^{2}}-\epsilon(\frac{1}{c}\frac{dz}{d\zeta}+c\frac{d^{3}z}{d\zeta^{3}})=G.$
(14)

Obviously, Eq. (14) reduces to a three-dimensional system:

$\begin{cases}\frac{dz}{d\zeta}=y,\\\ \frac{dy}{d\zeta}=w,\\\ \epsilon
c\frac{dw}{d\zeta}=-3z+\frac{9}{2}z^{2}-\frac{3}{2}z^{3}+w-\frac{\epsilon}{c}y-G.\end{cases}$
(15)

where $\epsilon$ is a sufficiently small parameter such that $0<\epsilon\ll
1$. Therefore, the traveling wave solutions of Eq. (3) can be obtained by
studying the corresponding orbits of system (15).

Obviously, system (15) is a singularly perturbed system described as “slow
system”. Rescaling $\zeta=\epsilon\eta$, we have the following equivalent
“fast system”:

$\begin{cases}\frac{dz}{d\eta}=\epsilon y,\\\ \frac{dy}{d\eta}=\epsilon w,\\\
c\frac{dw}{d\eta}=-3z+\frac{9}{2}z^{2}-\frac{3}{2}z^{3}+w-\frac{\epsilon}{c}y-G.\end{cases}$
(16)

Let $\epsilon=0$ in system (15) and (16). Then, the corresponding reduced
system is given by

$\begin{cases}\frac{dz}{d\zeta}=y,\\\ \frac{dy}{d\zeta}=w,\\\
0=-3z+\frac{9}{2}z^{2}-\frac{3}{2}z^{3}+w-G,\end{cases}$ (17)

and the layer system is

$\begin{cases}\frac{dz}{d\eta}=0,\\\ \frac{dy}{d\eta}=0,\\\
c\frac{dw}{d\eta}=-3z+\frac{9}{2}z^{2}-\frac{3}{2}z^{3}+w-G.\end{cases}$ (18)

which admits a two-dimensional critical manifold as follows

$M_{0}=\\{(z,y,w)\in\mathbb{R}^{3}|w=3z-\frac{9}{2}z^{2}+\frac{3}{2}z^{3}+G\\}.$
(19)

Suppose $A$ to be the linearized matrix of system (18), then it is given by

$A=\left(\begin{array}[]{ccc}0&0&0\\\ 0&0&0\\\
\frac{-3+9z-\frac{9}{2}z^{2}}{c}&0&\frac{1}{c}\\\ \end{array}\right),$ (20)

that the eigenvalues of A are $0$, $0$ and $\frac{1}{c}$. Thus, $M_{0}$ is a
normally hyperbolic invariant manifold (see Definition 1 and Definition 2).
There exists a two-dimensional slow submanifold $M_{\epsilon}$ which is
$C^{1}$ $O(\epsilon)$ close to $M_{0}$ (see Lemma 1). The invariant
submanifold $M_{\epsilon}$ is represented by

$M_{\epsilon}=\\{(z,y,w)\in\mathbb{R}^{3}|w=3z-\frac{9}{2}z^{2}+\frac{3}{2}z^{3}+G+\epsilon[cy(\frac{9}{2}z^{2}-9z+3+\frac{1}{c^{2}})]+O(\epsilon^{2})\\}.$
(21)

The dynamical behaviour of slow system (15) or fast system (16) is governed by

$\begin{cases}\frac{dz}{d\zeta}=y,\\\
\frac{dy}{d\zeta}=3z-\frac{9}{2}z^{2}+\frac{3}{2}z^{3}+G+\epsilon[cy(\frac{9}{2}z^{2}-9z+3+\frac{1}{c^{2}})]+O(\epsilon^{2}).\end{cases}$
(22)

Obviously, the unperturbed system $(\ref{ODE5})\mid_{\epsilon=0}$ admits
homoclinic, heteroclinic and periodic orbits with the parameter $G$ taking
different values. By the bifurcation theory of planar dynamical systems (see
[27]), we have the following proposition of unperturbed system
$(\ref{ODE5})\mid_{\epsilon=0}$:

###### Proposition 1

(i) If $\mid G\mid>\frac{\sqrt{3}}{3}$, there exists a single saddle point
(see Fig.1 (a) or (g)).
(ii) If $\mid G\mid=\frac{\sqrt{3}}{3}$, there exist a saddle point and a cusp
point (see Fig.1 (b) or (f)).
(iii) If $G=0$, there exist two saddle points $(0,0)$ and $(2,0)$, and a
center point $(1,0)$. It has two heteroclinic orbits surrounding the center
point $(1,0)$ to two saddle points $(0,0)$ and $(2,0)$ (see Fig.1 (d)).
(iv) If $0<\mid G\mid<\frac{\sqrt{3}}{3}$, there exist two saddle points and a
center point. It has a homoclinic orbit surrounding the center point to one
saddle point (see Fig.1 (c) or (e)).

(a) $G<-\frac{\sqrt{3}}{3}$ | (b) $G=-\frac{\sqrt{3}}{3}$ | (c) $-\frac{\sqrt{3}}{3}<G<0$ | (d) $G=0$
---|---|---|---

(e) $0<G<\frac{\sqrt{3}}{3}$ | (f) $G=\frac{\sqrt{3}}{3}$ | (g) $G>\frac{\sqrt{3}}{3}$
---|---|---

Fig. 1 The bifurcation and phase portraits of unperturbed system
(22)$\mid_{\epsilon=0}$.

Assume that $z(\zeta)$ is a solution of Eq. (14) satisfying
$\lim\limits_{\zeta\longrightarrow\infty}z(\zeta)={\kappa}$, thus
$u(x,t)=cz(\zeta)=cz[c(x-ct)]$ is a solitary wave solution or a kink wave
solution of Eq. (14).

For $-\frac{\sqrt{3}}{3}<G\leq 0$ and $0\leq G<\frac{\sqrt{3}}{3}$, the
unperturbed system (22)$\mid_{\epsilon=0}$ admits a saddle point defined by
$S_{0}({\kappa},0)$. We know that
$G=-(3{\kappa}-\frac{9}{2}{{\kappa}}^{2}+\frac{3}{2}{{\kappa}}^{3})$, then we
obtain $1+\frac{\sqrt{3}}{3}<{\kappa}\leq 2$ for $G\in(-\frac{\sqrt{3}}{3},0]$
and $0\leq{\kappa}<1-\frac{\sqrt{3}}{3}$ for $G\in[0,\frac{\sqrt{3}}{3})$,
respectively. We can regard ${\kappa}$ as a bifurcation parameter to find the
homoclinic orbits of system (22). Therefore, we rewrite system (22), that is:

$\begin{cases}\frac{dz}{d\zeta}=y,\\\
\frac{dy}{d\zeta}=3z-\frac{9}{2}z^{2}+\frac{3}{2}z^{3}-(3{\kappa}-\frac{9}{2}{{\kappa}}^{2}+\frac{3}{2}{{\kappa}}^{3})+\epsilon[cy(\frac{9}{2}z^{2}-9z+3+\frac{1}{c^{2}})]+O(\epsilon^{2}).\end{cases}$
(23)

##  4 Traveling wave solutions of Eq. $(\ref{wuzhang1})\mid_{\varepsilon=0}$

In this section, we study the solitary and kink (anti-kink) wave solutions of
Eq. $(\ref{wuzhang1})\mid_{\varepsilon=0}$. By dynamical method, we consider
the traveling wave solution of system $(\ref{ODE6})\mid_{\epsilon=0}$ with
${\kappa}\in[0,1-\frac{\sqrt{3}}{3})\cup(1+\frac{\sqrt{3}}{3},2]$. The
unperturbed system $(\ref{ODE6})$ is of the form:

$\begin{cases}\frac{dz}{d\zeta}=y,\\\
\frac{dy}{d\zeta}=3z-\frac{9}{2}z^{2}+\frac{3}{2}z^{3}-(3{\kappa}-\frac{9}{2}{{\kappa}}^{2}+\frac{3}{2}{{\kappa}}^{3}).\end{cases}$
(24)

It is easy to obtain the first integral of system (24), it is given by

$H(z,y)=\frac{1}{2}y^{2}-\frac{3}{8}z^{4}+\frac{3}{2}z^{3}-\frac{3}{2}z^{2}+(\frac{3}{2}{{\kappa}}^{3}-\frac{9}{2}{{\kappa}}^{2}+3{\kappa})z.$
(25)

The homoclinic orbit $\Gamma({\kappa})$ to the saddle $({\kappa},0)$ is
defined by

$\begin{array}[]{rcl}H(z,y)&=&\frac{1}{2}y^{2}-\frac{3}{8}z^{4}+\frac{3}{2}z^{3}-\frac{3}{2}z^{2}+(\frac{3}{2}{{\kappa}}^{3}-\frac{9}{2}{{\kappa}}^{2}+3{\kappa})z\\\
&=&\frac{9}{8}{{\kappa}}^{4}-3{{\kappa}}^{3}+\frac{3}{2}{{\kappa}}^{2}.\end{array}$
(26)

And the heteroclinic orbits $\Upsilon({\kappa})_{\pm}$ (where $\pm$ represents
upper and lower branch curves) to the two saddle points $(0,0)$ and $(2,0)$,
namely, ${\kappa}=0$ or ${\kappa}=2$, the heteroclinic curve is determined by

$\begin{array}[]{rcl}H(z,y)=\frac{1}{2}y^{2}-\frac{3}{8}z^{4}+\frac{3}{2}z^{3}-\frac{3}{2}z^{2}=0.\end{array}$
(27)

###  4.1 Solitary wave solutions of Eq. $(\ref{wuzhang1})\mid_{\varepsilon=0}$

By (26) and Proposition 1 (_iv_), we can obtain the expression of $y$ as
follows:

$\begin{array}[]{rcl}y&=&\pm\sqrt{\frac{3}{4}z^{4}-3z^{3}+3z^{2}-(3{{\kappa}}^{3}-9{{\kappa}}^{2}+6{\kappa})z+\frac{9}{4}{{\kappa}}^{4}-6{{\kappa}}^{3}+3{{\kappa}}^{2}}\\\
&=&\pm\sqrt{\frac{3}{4}(z-{\kappa})^{2}(z-z_{+})(z-z_{-})},\end{array}$ (28)

where $z_{\pm}=-{\kappa}+2\pm\sqrt{-2{{\kappa}}^{2}+4{\kappa}}$. It can be
seen that $z_{-}<z_{+}<{\kappa}$ for $1+\frac{\sqrt{3}}{3}<{\kappa}<2$ (see
Fig. 1(c)) and ${\kappa}<z_{-}<z_{+}$ for $0<{\kappa}<1-\frac{\sqrt{3}}{3}$
(see Fig. 1(e)).

Since $\frac{dz}{d\zeta}=y$, it implies

$\frac{dz}{d\zeta}=\pm\sqrt{\frac{3}{4}(z-{\kappa})^{2}(z-z_{+})(z-z_{-})},$
(29)

which yields

$\zeta=\pm\int_{z_{+}}^{z}\frac{1}{\sqrt{\frac{3}{4}(s-{\kappa})^{2}(s-z_{+})(s-z_{-})}}ds,\
\ (1+\frac{\sqrt{3}}{3}<{\kappa}<2),$ (30)

and

$\zeta=\pm\int_{z_{-}}^{z}\frac{1}{\sqrt{\frac{3}{4}(s-{\kappa})^{2}(z_{+}-s)(z_{-}-s)}}ds,\
\ (0<{\kappa}<1-\frac{\sqrt{3}}{3}).$ (31)

Thus, for $1+\frac{\sqrt{3}}{3}<{\kappa}<2$, the parametric representation of
the dark solitary wave of system $(\ref{ODE6})|_{\epsilon=0}$ can be obtained
by (see Fig. 2(a))

$z(\zeta)=-\frac{2(6{{\kappa}}^{2}-12{\kappa}+4)}{2\textmd{cosh}(\frac{1}{2}\sqrt{18{{\kappa}}^{2}-36{\kappa}+12}\zeta)\sqrt{-2{{\kappa}}^{2}+4{\kappa}}+4{\kappa}-4}.$
(32)

Then, it corresponds to a dark solitary wave solution (see Fig. 2(b)) of Eq.
$(\ref{wuzhang1})\mid_{\varepsilon=0}$ is given by:

$u(x,t)=-\frac{2(6{{\kappa}}^{2}-12{\kappa}+4)c}{2\textmd{cosh}[\frac{1}{2}\sqrt{18{{\kappa}}^{2}-36{\kappa}+12}c(x-ct)]\sqrt{-2{{\kappa}}^{2}+4{\kappa}}+4{\kappa}-4}.$ (33) (a) Dark solitary wave | (b) Exact solution
---|---

Fig. 2 Dark solitary wave of system $(\ref{ODE6})|_{\epsilon=0}$ and exact
solution of Eq. $(\ref{wuzhang1})\mid_{\varepsilon=0}$ for $c=0.5$ and
${\kappa}=1.7$ .

For $0<{\kappa}<1-\frac{\sqrt{3}}{3}$, the expression of the bright solitary
wave of system $(\ref{ODE6})|_{\epsilon=0}$ can be obtained by (see Fig 3.(a))

$z(\zeta)=-\frac{2(6{{\kappa}}^{2}-12{\kappa}+4)}{4{\kappa}-4-2\textmd{cosh}(\frac{1}{2}\sqrt{18{{\kappa}}^{2}-36{\kappa}+12}\zeta)\sqrt{-2{{\kappa}}^{2}+4{\kappa}}}.$
(34)

The bright solitary wave solution of Eq.
$(\ref{wuzhang1})\mid_{\varepsilon=0}$ (see Fig 3.(b)) is of the form

$u(x,t)=-\frac{2(6{{\kappa}}^{2}-12{\kappa}+4)c}{4{\kappa}-4-2\textmd{cosh}\big{[}\frac{1}{2}\sqrt{18{{\kappa}}^{2}-36{\kappa}+12}c(x-ct)\big{]}\sqrt{-2{{\kappa}}^{2}+4{\kappa}}}.$
(35)

(a) Bright solitary wave | (b) Exact solution
---|---

Fig. 3 Bright solitary wave of system $(\ref{ODE6})|_{\epsilon=0}$ and exact
solution of Eq. $(\ref{wuzhang1})\mid_{\varepsilon=0}$ for $c=0.5$ and
${\kappa}=0.3$.

In summary, we give the theorem as follows.

###### Theorem 1

For any wave speed $c>0$ and
${\kappa}\in(0,1-\frac{\sqrt{3}}{3})\cup(1+\frac{\sqrt{3}}{3},2)$, Eq.
$(\ref{wuzhang1})\mid_{\varepsilon=0}$ has solitary wave solutions given by
(33) or (35).

###  4.2 Kink and anti-kink wave solutions of Eq.
$(\ref{wuzhang1})\mid_{\varepsilon=0}$

According (27) and Proposition 1 (_iii_), we can obtain the expression of $y$
as follows:

$\begin{array}[]{rcl}y=\pm\sqrt{\frac{3}{4}z^{4}-3z^{3}+3z^{2}}=\pm\sqrt{\frac{3}{4}(z-0)^{2}(2-z)^{2}}=\pm\frac{\sqrt{3}}{2}z(z-2).\end{array}$
(36)

Similarly, due to $\frac{dz}{d\zeta}=y$, we have

$\zeta=\pm\int_{1}^{z}\frac{1}{\frac{\sqrt{3}}{2}s(s-2)}ds.$ (37)

For $\kappa=0$ or $\kappa=0$, the expression of the kink or anti-kink wave of
system $(\ref{ODE6})|_{\epsilon=0}$ can be obtained by (see Fig. 4(a) and (c))

$z(\zeta)=-\frac{2c}{1+e^{\pm\sqrt{3}c\zeta}}.$ (38)

Therefore, we obtain the following kink and anti-kink wave solutions (see Fig.
4(b) and (b)):

$u(x,t)=-\frac{2c}{1+e^{\pm\sqrt{3}c(x-ct)}}.$ (39)

(a) Kink wave | (b) Kink wave solution
---|---

(c) Anti-kink wave | (d) Anti-kink wave solution
---|---

Fig. 4 Kink and anti-kink wave of system $(\ref{ODE6})|_{\epsilon=0}$. Kink
and anti-kink wave solutions of Eq. $(\ref{wuzhang1})\mid_{\varepsilon=0}$ for
$c=0.5$.

Based on above analysis, we have the following theorems:

###### Theorem 2

For any wave speed $c>0$ and $\kappa=0$ or $\kappa=2$, Eq.
$(\ref{wuzhang1})\mid_{\varepsilon=0}$ has kink or anti-kink solutions given
by (39).

##  5  Bifurcations of the travelling wave solutions for the perturbed Eq. (3)

In this section, before stating our main results on the existence of the
solitary and kink wave solutions for the perturbed Eq. (3), denote

$\begin{array}[]{rcl}c_{1}(\kappa)\Big{|}_{\kappa\in(0,1-\frac{\sqrt{3}}{3})}&=&\frac{\sqrt{15}}{3}\big{[}\big{(}6\textmd{ln}2{{\kappa}}^{3}+6{{\kappa}}^{3}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)-18\textmd{ln}2{{\kappa}}^{2}\\\
&&-18{{\kappa}}^{2}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)+12\textmd{ln}2{\kappa}\\\
&&+12{\kappa}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)+2\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}~{}\big{]}\\\
&&-\sqrt{3}{\kappa}\big{[}2{{\kappa}}^{2}\textmd{ln}(-\sqrt{2{\kappa}-{{\kappa}}^{2}})+3\textmd{ln}2{{\kappa}}^{2}-6{\kappa}\textmd{ln}(-\sqrt{2{\kappa}-{{\kappa}}^{2}})-9\textmd{ln}2{\kappa}\\\
&&+4\textmd{ln}(-\sqrt{2{\kappa}-{{\kappa}}^{2}})+6\textmd{ln}2\big{)}\big{/}\big{(}27{{\kappa}}^{4}\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}\\\
&&+30{{\kappa}}^{3}\textmd{ln}(-\sqrt{2{\kappa}-{{\kappa}}^{2}})+60{\kappa}\textmd{ln}(-\sqrt{2{\kappa}-{{\kappa}}^{2}})+15\textmd{ln}2{{\kappa}}^{3}\\\
&&-30{{\kappa}}^{3}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)-90{{\kappa}}^{2}\textmd{ln}(-\sqrt{2{\kappa}-{{\kappa}}^{2}})\\\
&&+90{{\kappa}}^{2}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)+30\sqrt{2}{\kappa}\sqrt{(3{{\kappa}}^{2}-6{\kappa}+2)^{3}}\\\
&&+294{{\kappa}}^{2}\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}-45\textmd{ln}2{{\kappa}}^{2}-198{{\kappa}}^{3}\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}\\\
&&-60{\kappa}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)+10\sqrt{2}\sqrt{(3{{\kappa}}^{2}-6{\kappa}+2)^{3}}\\\
&&-72{\kappa}\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+30\textmd{ln}2{\kappa}-18\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}\
\big{)}\big{]}^{\frac{1}{2}}.\end{array}$ (40)

and

$\begin{array}[]{rcl}c_{2}(\kappa)\Big{|}_{\kappa\in(1+\frac{\sqrt{3}}{3},2)}&=&\frac{\sqrt{15}}{3}\big{[}\big{(}6\textmd{ln}2{{\kappa}}^{3}+6{{\kappa}}^{3}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)-18\textmd{ln}2{{\kappa}}^{2}\\\
&&-18{{\kappa}}^{2}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)+12\textmd{ln}2{\kappa}\\\
&&+12{\kappa}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)+2\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}~{}\big{]}\\\
&&-\sqrt{3}{\kappa}\big{[}{{\kappa}}^{2}\textmd{ln}(2{\kappa}-{{\kappa}}^{2})+3\textmd{ln}2{{\kappa}}^{2}-3{\kappa}\textmd{ln}(2{\kappa}-{{\kappa}}^{2})-9\textmd{ln}2{\kappa}\\\
&&+2\textmd{ln}(2{\kappa}-{{\kappa}}^{2})+6\textmd{ln}2\big{)}\big{/}\big{(}27{{\kappa}}^{4}\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}\\\
&&+15{{\kappa}}^{3}\textmd{ln}(2{\kappa}-{{\kappa}}^{2})-198{{\kappa}}^{3}\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}\\\
&&-30{{\kappa}}^{3}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)+15\textmd{ln}2{{\kappa}}^{3}-45{{\kappa}}^{2}\textmd{ln}(2{\kappa}-{{\kappa}}^{2})\\\
&&+90{{\kappa}}^{2}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)+30\sqrt{2}{\kappa}\sqrt{(3{{\kappa}}^{2}-6{\kappa}+2)^{3}}\\\
&&+294{{\kappa}}^{2}\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}-45\textmd{ln}2{{\kappa}}^{2}+30{\kappa}\textmd{ln}(2{\kappa}-{{\kappa}}^{2})\\\
&&-60{\kappa}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)+10\sqrt{2}\sqrt{(3{{\kappa}}^{2}-6{\kappa}+2)^{3}}\\\
&&-72{\kappa}\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+30\textmd{ln}2{\kappa}-18\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}~{}\big{)}\big{]}^{\frac{1}{2}},\end{array}$
(41)

Now we are in a position to state our main results on the bifurcations of the
travelling waves for the perturbed Eq. (3).

###### Theorem 3

(i) If ${\kappa}\in(0,1-\frac{\sqrt{3}}{3})$, then for sufficiently small
$\epsilon$ ($0<\epsilon\ll 1$), there exists a wave speed
$\tilde{c}_{1}(\kappa,\epsilon)=c_{1}({\kappa})+O(\epsilon)$ such that system
(23) possesses a homoclinic orbit. Thus, Eq. (3) has a solitary wave solution
$u=u(x,t,{\kappa},\epsilon)$ with a wave speed $c=\tilde{c}_{1}$.

(ii) If ${\kappa}\in(1+\frac{\sqrt{3}}{3},2)$, then for sufficiently small
$\epsilon$ ($0<\epsilon\ll 1$), there exists a wave speed
$\tilde{c}_{2}(\kappa,\epsilon)=c_{2}({\kappa})+O(\epsilon)$ such that system
(23) possesses a homoclinic orbit. Thus, Eq. (3) has a solitary wave solution
$u=u(x,t,{\kappa},\epsilon)$ with a wave speed $c=\tilde{c}_{2}$.

(iii) If $\kappa=0$ or $\kappa=2$, then for sufficiently small $\epsilon$
($0<\epsilon\ll 1$), there exists a unique wave speed
$\tilde{c}_{3}(\epsilon)=\sqrt{15}/3+O(\epsilon)$ such that system (23) has a
heteroclinic orbit. Therefore, Eq. (3) has a pair of kink and anti-kink wave
solutions $u=u(x,t,\kappa,\epsilon)$ with a wave speed $c=\tilde{c}_{3}$.

Remark 1. Different kinds of travelling wave solutions exist in the different
bifurcation parameter regions. For DWLE, solitary wave and kink wave can not
coexist at a same wave speed.

Remark 2. The exact limit wave speed is given for each bifurcation parameter
regions.

In order to prove Theorem 3, we need to introduce some lemmas.

Poincaré map $d_{i}$ ($i=1,2$) is a powerful tool to study the bifurcation of
homoclinic or heteroclinic orbits (e.g., see [21, 33, 45, 34]).
Case a. It is defined by (see Fig. 5(a))

$d_{1}:A(h)\longrightarrow P(h,c,\varepsilon),$

and

$\begin{array}[]{rcl}d_{1}(h,c,\varepsilon)&=&\displaystyle\int_{A(h)}^{P(h,c,\varepsilon)}dH\\\
&=&\displaystyle\int_{-\infty}^{+\infty}(\frac{\partial H}{\partial
x}\frac{dx}{dt}+\frac{\partial H}{\partial
y}\frac{dy}{dt})\mid_{L_{(h,c,\varepsilon)}}dt,\end{array}$

where $A(h)$ is the initial point and $P(h,c,\varepsilon)$ is the mapping
point. $L(h,c,\varepsilon)$ and $L(h,c)$ are the perturbed and unperturbed
orbits, respectively. We have $\lim\limits_{\varepsilon\rightarrow
0}L(h,c,\varepsilon)=L(h,c)$. And $d_{1}(h,c,\varepsilon)$ can be Taylor
expanded in $\varepsilon$, one has

$d_{1}(h,c,\varepsilon)=\varepsilon M_{\textmd{hom}}(h,c)+O(\varepsilon^{2}).$

Case b. A heteroclinic orbit $L(h,c)$ connecting two saddle points $S_{1}$ and
$S_{2}$. Then, $L^{+}_{s}(h,c,\varepsilon)$ and $L^{-}_{u}(h,c,\varepsilon)$
represent the stable and unstable manifolds of perturbed heteroclinic orbit
$L(h,c,\varepsilon)$. Taking a point $P(h,c)\in L(h,c)$, we let $L^{*}$ be a
segment normal of $L(h,c)$ at point $P(h,c)$. For $0<\varepsilon\ll 1$, we
assume that $L^{+}_{s}(h,c,\varepsilon)$ and $L^{-}_{u}(h,c,\varepsilon)$
intersect the normal line $L^{*}$ transversally at points
$P^{+}_{s}(h,c,\varepsilon)$ and $P^{-}_{u}(h,c,\varepsilon)$ (see Fig. 5(b)).

Let

$\begin{array}[]{rcl}d_{2}(h,c,\varepsilon)&=&-\overrightarrow{n}\cdot\overrightarrow{P^{+}_{s}P^{-}_{u}},\end{array}$

where
$\overrightarrow{n}=\frac{(H_{z}(P(h,c)),H_{y}(P(h,c)))}{\mid(H_{y}(P(h,c)),-H_{z}(P(h,c)))\mid}$.
And $d_{2}(h,c,\varepsilon)$ can be Taylor expanded in $\varepsilon$, we have

$d_{2}(h,c,\varepsilon)=\varepsilon\cdot A\cdot
M_{\textmd{het}}(h,c)+O(\varepsilon^{2}),$

where $A$ is a constant.

Usually, $d_{i}(h,c,\varepsilon)$ ($i=1,2$) is used to measure the distance
between perturbed stable and unstable manifolds. Then, we say that
$M_{\textmd{hom}}(h,c)$ or $M_{\textmd{het}}(h,c)$ is the so called “first-
order” Melnikov integral.

(a) The type of periodic or homoclinic orbits | (b) The type of heteroclinic orbits
---|---

Fig. 5 Poincaré map

In this paper, “first-order” Melnikov integral is employed to detect the
persistence of the homoclinic and heteroclinic orbits under small
perturbation, respectively. Consequently, the existence of solitary and kink
wave solutions are proved under small perturbation.

Firstly, the homoclinic Melnikov function of system (23) is defined by

$\begin{array}[]{rcl}M_{\textmd{hom}}(c,{\kappa})=\displaystyle\oint_{\Gamma({\kappa})}(\frac{9}{2}z^{2}-9z+3+\frac{1}{c^{2}})y^{2}d\zeta=\frac{1}{c^{2}}I({\kappa})+J({\kappa}),\end{array}$
(42)

where $I({\kappa})=\displaystyle\oint_{\Gamma({\kappa})}y^{2}d\zeta$ and
$J({\kappa})=\displaystyle\oint_{\Gamma({\kappa})}(\frac{9}{2}z^{2}-9z+3)y^{2}d\zeta$.

###### Lemma 2

For any ${\kappa}\in(0,1-\frac{\sqrt{3}}{3})\cup(1+\frac{\sqrt{3}}{3},2)$,
there exists a positive root $c=c({\kappa})$ of
$M_{\textmd{hom}}(c,{\kappa})=0$. Moreover, $\frac{\partial
M_{\textmd{hom}}}{\partial c}\mid_{c=c({\kappa})}\neq 0$.

_Proof._ Obviously, we know that
$I({\kappa})=\displaystyle\oint_{\Gamma({\kappa})}y^{2}d\zeta>0$. Because
$J({\kappa})$ integrates along the homoclinic orbit $\Gamma({\kappa})$, then

$\begin{array}[]{rcl}y^{\prime\prime}=(\frac{9}{2}z^{2}-9z+3)z^{\prime}=(\frac{9}{2}z^{2}-9z+3)y,\end{array}$
(43)

which yields

$\begin{array}[]{rcl}dy^{\prime}=(\frac{9}{2}z^{2}-9z+3)z^{\prime}d\zeta=(\frac{9}{2}z^{2}-9z+3)yd\zeta.\end{array}$
(44)

By integration of parts, it is not difficult to have

$\begin{array}[]{rcl}J({\kappa})&=&\displaystyle\oint_{\Gamma({\kappa})}(\frac{9}{2}z^{2}-9z+3)y^{2}d\zeta=\displaystyle\oint_{\Gamma({\kappa})}ydy^{\prime}\\\
&=&-\displaystyle\int_{\mathbb{R}}(y^{\prime})^{2}d\zeta<0.\end{array}$ (45)

Hence, $-\frac{I({\kappa})}{J({\kappa})}>0$, there exists a single positive
root $c=c({\kappa})$ of $M_{\textmd{hom}}(c,{\kappa})=0$ that

$\begin{array}[]{rcl}c({\kappa})=\sqrt{-\frac{I({\kappa})}{J({\kappa})}}.\end{array}$
(46)

On the other hand, it follows from (42) that

$\begin{array}[]{rcl}\frac{\partial M_{\textmd{hom}}(c,{\kappa})}{\partial
c}\mid_{c=c({\kappa})}=-2\frac{1}{c^{3}}I({\kappa})=-2J({\kappa})\sqrt{-\frac{J({\kappa})}{I({\kappa})}}>0.\end{array}$
(47)

The proof is completed.

In fact, the analytical expressions of $I({\kappa})$ and $J({\kappa})$ can be
computed by dividing them into two cases.

Case (I): $1+\frac{\sqrt{3}}{3}<{\kappa}<2$.

Firstly, since $I{({\kappa})}$ and $J{({\kappa})}$ are integrated over a
closed curve, and the time variable $\zeta$ can be represented by the state
variable $z$ on homoclinic orbit. Thus, by (28), we have

$\begin{array}[]{rcl}I({\kappa})&=&\displaystyle\oint_{\Gamma({\kappa})}y^{2}d\zeta=\displaystyle\oint_{\Gamma({\kappa})}ydz\\\
&=&2\displaystyle\int_{z_{+}}^{{\kappa}}\sqrt{\frac{3}{4}z^{4}-3z^{3}+3z^{2}-(3{{\kappa}}^{3}-9{{\kappa}}^{2}+6{\kappa})z+\frac{9}{4}{{\kappa}}^{4}-6{{\kappa}}^{3}+3{{\kappa}}^{2}}dz\\\
&=&2\displaystyle\int_{z_{+}}^{{\kappa}}\sqrt{\frac{3}{4}(z-{\kappa})^{2}(z-z_{+})(z-z_{-})}dz\\\
&=&\sqrt{3}\big{[}\displaystyle\int_{z_{+}}^{{\kappa}}{\kappa}\sqrt{(z-z_{+})(z-z_{-})}dz-\displaystyle\int_{z_{+}}^{{\kappa}}z\sqrt{(z-z_{+})(z-z_{-})}dz\big{]}\\\
&=&\sqrt{3}\big{[}{\kappa}I_{1}({\kappa})-I_{2}({\kappa})\big{]},\end{array}$
(48)

where
$I_{1}({\kappa})=\displaystyle\int_{z_{+}}^{{\kappa}}\sqrt{(z-z_{+})(z-z_{-})}dz$
and
$I_{2}({\kappa})=\displaystyle\int_{z_{+}}^{{\kappa}}z\sqrt{(z-z_{+})(z-z_{-})}dz$.

Note that the formulas 2.261, 2.262 (1) and 2.262 (2) in handbook [20] read:

$\begin{array}[]{rcl}\displaystyle\int\frac{dx}{\sqrt{R}}&=&\frac{1}{\sqrt{c_{1}}}\textmd{ln}(2\sqrt{c_{1}R}+2c_{1}x+b),\
\ (c_{1}>0,\ 2c_{1}x+b>\sqrt{-\Delta},\ \Delta<0),\\\
\displaystyle\int\sqrt{R}dx&=&\frac{(2c_{1}x+b)\sqrt{R}}{4c_{1}}+\frac{\Delta}{8c_{1}}\displaystyle\int\frac{dx}{\sqrt{R}},\\\
\displaystyle\int
x\sqrt{R}dx&=&\frac{\sqrt{R^{3}}}{3c_{1}}-\frac{(2c_{1}x+b)b\sqrt{R}}{8c_{1}^{2}}-\frac{b\Delta}{16c_{1}^{2}}\displaystyle\int\frac{dx}{\sqrt{R}},\end{array}$
(49)

where $R=a+bx+c_{1}x^{2}$ and $\Delta=4ac_{1}-b^{2}$. Combining (48) and (49),
it follows that

$\begin{array}[]{rcl}I({\kappa})&=&\frac{\sqrt{3}}{3}\big{[}6\textmd{ln}2{{\kappa}}^{3}+6{{\kappa}}^{3}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)-18\textmd{ln}2{{\kappa}}^{2}\\\
&&-18{{\kappa}}^{2}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)+12\textmd{ln}2{\kappa}\\\
&&+12{\kappa}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)+2\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}~{}\big{]}\\\
&&-\sqrt{3}{\kappa}\big{[}{{\kappa}}^{2}\textmd{ln}(2{\kappa}-{{\kappa}}^{2})+3\textmd{ln}2{{\kappa}}^{2}-3{\kappa}\textmd{ln}(2{\kappa}-{{\kappa}}^{2})-9\textmd{ln}2{\kappa}\\\
&&+2\textmd{ln}(2{\kappa}-{{\kappa}}^{2})+6\textmd{ln}2\big{]}.\end{array}$
(50)

Secondly,

$\begin{array}[]{rcl}J({\kappa})&=&\displaystyle\oint_{\Gamma({\kappa})}(\frac{9}{2}z^{2}-9z+3)y^{2}d\zeta=\displaystyle\oint_{\Gamma({\kappa})}(\frac{9}{2}z^{2}-9z+3)ydz\\\
&=&2\displaystyle\int_{z_{+}}^{{\kappa}}(\frac{9}{2}z^{2}-9z+3)\scriptstyle\sqrt{\frac{3}{4}z^{4}-3z^{3}+3z^{2}-(3{{\kappa}}^{3}-9{{\kappa}}^{2}+6{\kappa})z+\frac{9}{4}{{\kappa}}^{4}-6{{\kappa}}^{3}+3{{\kappa}}^{2}}\displaystyle
dz\\\
&=&\sqrt{3}\displaystyle\int_{z_{+}}^{{\kappa}}(\frac{9}{2}z^{2}-9z+3)\sqrt{({\kappa}-z)^{2}(z-z_{+})(z-z_{-})}dz\\\
&=&\sqrt{3}\displaystyle\int_{z_{+}}^{{\kappa}}\big{[}-\frac{9}{2}z^{3}+(\frac{9}{2}{\kappa}+9)z^{2}-(9{\kappa}+3)z+3{\kappa}\big{]}\sqrt{(z-z_{+})(z-z_{-})}dz\\\
&=&\sqrt{3}\big{[}-\frac{9}{2}J_{1}({\kappa})+(\frac{9}{2}{\kappa}+9)J_{2}({\kappa})-(9{\kappa}+3)I_{2}({\kappa})+3{\kappa}I_{1}({\kappa})\big{]},\end{array}$
(51)

where
$J_{1}({\kappa})=\displaystyle\int_{z_{+}}^{{\kappa}}z^{3}\sqrt{(z-z_{+})(z-z_{-})}dz$,
$J_{2}({\kappa})=\displaystyle\int_{z_{+}}^{{\kappa}}z^{2}\sqrt{(z-z_{+})(z-z_{-})}dz$.

Similarly, by the formulas 2.262 (3) and 2.262 (4) in [20], they are expressed
by:

$\begin{array}[]{rcl}\displaystyle\int
x^{2}\sqrt{R}dx&=&(\frac{x}{4c_{1}}-\frac{5b}{24c_{1}^{2}})\sqrt{R^{3}}+(\frac{5b^{2}}{16c_{1}^{2}}-\frac{a}{4c_{1}})\frac{(2c_{1}x+b)\sqrt{R}}{4c_{1}}+(\frac{5b^{2}}{16c_{1}^{2}}-\frac{a}{4c_{1}})\frac{\Delta}{8c_{1}}\displaystyle\int\frac{dx}{\sqrt{R}},\\\
\displaystyle\int
x^{3}\sqrt{R}dx&=&(\frac{x^{2}}{5c_{1}}-\frac{7bx}{40c_{1}^{2}}+\frac{7b^{2}}{48c_{1}^{3}}-\frac{2a}{15c_{1}^{2}})\sqrt{R^{3}}-(\frac{7b^{3}}{32c_{1}^{3}}-\frac{3ab}{8c_{1}^{2}})\frac{(2c_{1}x+b)\sqrt{R}}{4c_{1}}\\\
&&-(\frac{7b^{3}}{32c_{1}^{3}}-\frac{3ab}{8c_{1}^{2}})\frac{\Delta}{8c_{1}}\displaystyle\int\frac{dx}{\sqrt{R}}.\end{array}$
(52)

Using (49), (51) and (52), we have

$\begin{array}[]{rcl}J({\kappa})&=&-\frac{\sqrt{3}}{5}\big{[}27{{\kappa}}^{4}\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+15{{\kappa}}^{3}\textmd{ln}(2{\kappa}-{{\kappa}}^{2})-198{{\kappa}}^{3}\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}\\\
&&-30{{\kappa}}^{3}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)+15\textmd{ln}2{{\kappa}}^{3}-45{{\kappa}}^{2}\textmd{ln}(2{\kappa}-{{\kappa}}^{2})\\\
&&+90{{\kappa}}^{2}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)+30\sqrt{2}{\kappa}\sqrt{(3{{\kappa}}^{2}-6{\kappa}+2)^{3}}\\\
&&+294{{\kappa}}^{2}\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}-45\textmd{ln}2{{\kappa}}^{2}+30{\kappa}\textmd{ln}(2{\kappa}-{{\kappa}}^{2})\\\
&&-60{\kappa}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)+10\sqrt{2}\sqrt{(3{{\kappa}}^{2}-6{\kappa}+2)^{3}}\\\
&&-72{\kappa}\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+30\textmd{ln}2{\kappa}-18\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}~{}\big{]}.\end{array}$
(53)

Case (II): $0<{\kappa}<1-\frac{\sqrt{3}}{3}$.

Same calculation method as case (I), we obtain

$\begin{array}[]{rcl}I({\kappa})&=&\displaystyle\oint_{\Gamma({\kappa})}y^{2}d\zeta=\displaystyle\oint_{\Gamma({\kappa})}ydz=2\displaystyle\int_{z_{-}}^{{\kappa}}\cdots
dz\\\
&=&\frac{\sqrt{3}}{3}\big{[}6\textmd{ln}2{{\kappa}}^{3}+6{{\kappa}}^{3}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)-18\textmd{ln}2{{\kappa}}^{2}\\\
&&-18{{\kappa}}^{2}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)+12\textmd{ln}2{\kappa}\\\
&&+12{\kappa}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)+2\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}~{}\big{]}\\\
&&-\sqrt{3}{\kappa}\big{[}2{{\kappa}}^{2}\textmd{ln}(-\sqrt{2{\kappa}-{{\kappa}}^{2}})+3\textmd{ln}2{{\kappa}}^{2}-6{\kappa}\textmd{ln}(-\sqrt{2{\kappa}-{{\kappa}}^{2}})-9\textmd{ln}2{\kappa}\\\
&&+4\textmd{ln}(-\sqrt{2{\kappa}-{{\kappa}}^{2}})+6\textmd{ln}2\big{]},\end{array}$
(54)

and

$\begin{array}[]{rcl}J({\kappa})&=&\displaystyle\oint_{\Gamma({\kappa})}(\frac{9}{2}z^{2}-9z+3)y^{2}d\zeta=\displaystyle\oint_{\Gamma({\kappa})}(\frac{9}{2}z^{2}-9z+3)ydz=2\displaystyle\int_{z_{-}}^{{\kappa}}\cdots
dz\\\
&=&-\frac{\sqrt{3}}{5}\big{[}27{{\kappa}}^{4}\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+30{{\kappa}}^{3}\textmd{ln}(-\sqrt{2{\kappa}-{{\kappa}}^{2}})+60{\kappa}\textmd{ln}(-\sqrt{2{\kappa}-{{\kappa}}^{2}})\\\
&&-30{{\kappa}}^{3}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)+15\textmd{ln}2{{\kappa}}^{3}-90{{\kappa}}^{2}\textmd{ln}(-\sqrt{2{\kappa}-{{\kappa}}^{2}})\\\
&&+90{{\kappa}}^{2}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)+30\sqrt{2}{\kappa}\sqrt{(3{{\kappa}}^{2}-6{\kappa}+2)^{3}}\\\
&&+294{{\kappa}}^{2}\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}-45\textmd{ln}2{{\kappa}}^{2}-198{{\kappa}}^{3}\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}\\\
&&-60{\kappa}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)+10\sqrt{2}\sqrt{(3{{\kappa}}^{2}-6{\kappa}+2)^{3}}\\\
&&-72{\kappa}\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+30\textmd{ln}2{\kappa}-18\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}\
\big{]}.\end{array}$ (55)

We plot the algebraic curve of $I({\kappa})$ with respect to ${\kappa}$ and
$J({\kappa})$ with respect to ${\kappa}$, respectively. (see Fig. 6 and Fig.
7)

(a) (${\kappa}$, $I(\kappa)$) | (b) (${\kappa}$, $J(\kappa)$)
---|---

Fig. 6 The algebraic curves of $I({\kappa})$ and $J({\kappa})$ with respect to
${\kappa}$ for $1+\frac{\sqrt{3}}{3}<{\kappa}<2$.

(a) (${\kappa}$, $I(\kappa)$) | (b) (${\kappa}$, $J(\kappa)$)
---|---

Fig. 7 The algebraic curves of $I({\kappa})$ and $J({\kappa})$ with respect to
${\kappa}$ for $0<{\kappa}<1-\frac{\sqrt{3}}{3}$.

We can see that $I({\kappa})>0$ and $J({\kappa})<0$ for
$1+\frac{\sqrt{3}}{3}<{\kappa}<2$ in Fig. 6, $I({\kappa})>0$ and
$J({\kappa})<0$ for $0<{\kappa}<1-\frac{\sqrt{3}}{3}$ in Fig. 7. It is shown
that a good agreement with the proof of Lemma 2.

Hence, we substitute (50) and (53) into (42), and the expression of the
homoclinic Melnikov function of system (23) for
$1+\frac{\sqrt{3}}{3}<{\kappa}<2$ can be rewritten as

$\begin{array}[]{rcl}&&M_{\textmd{hom}}(c,{\kappa})\Big{|}_{\kappa\in(1+\frac{\sqrt{3}}{3},2)}\\\
&=&\frac{1}{c^{2}}\frac{\sqrt{3}}{3}\big{[}6\textmd{ln}2{{\kappa}}^{3}+6{{\kappa}}^{3}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)-18\textmd{ln}2{{\kappa}}^{2}\\\
&&-18{{\kappa}}^{2}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)+12\textmd{ln}2{\kappa}\\\
&&+12{\kappa}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)+2\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}~{}\big{]}\\\
&&-\frac{1}{c^{2}}\sqrt{3}{\kappa}\big{[}{{\kappa}}^{2}\textmd{ln}(2{\kappa}-{{\kappa}}^{2})+3\textmd{ln}2{{\kappa}}^{2}-3{\kappa}\textmd{ln}(2{\kappa}-{{\kappa}}^{2})-9\textmd{ln}2{\kappa}\\\
&&+2\textmd{ln}(2{\kappa}-{{\kappa}}^{2})+6\textmd{ln}2\big{]}-\frac{\sqrt{3}}{5}\big{[}27{{\kappa}}^{4}\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}\\\
&&+15{{\kappa}}^{3}\textmd{ln}(2{\kappa}-{{\kappa}}^{2})-198{{\kappa}}^{3}\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}\\\
&&-30{{\kappa}}^{3}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)+15\textmd{ln}2{{\kappa}}^{3}-45{{\kappa}}^{2}\textmd{ln}(2{\kappa}-{{\kappa}}^{2})\\\
&&+90{{\kappa}}^{2}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)+30\sqrt{2}{\kappa}\sqrt{(3{{\kappa}}^{2}-6{\kappa}+2)^{3}}\\\
&&+294{{\kappa}}^{2}\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}-45\textmd{ln}2{{\kappa}}^{2}+30{\kappa}\textmd{ln}(2{\kappa}-{{\kappa}}^{2})\\\
&&-60{\kappa}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)+10\sqrt{2}\sqrt{(3{{\kappa}}^{2}-6{\kappa}+2)^{3}}\\\
&&-72{\kappa}\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+30\textmd{ln}2{\kappa}-18\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}~{}\big{]}.\end{array}$
(56)

It is easy to verify that for $\kappa\in(1+\frac{\sqrt{3}}{3},2)$, there is
$c=c_{1}(\kappa)$ such that
$M_{\textmd{hom}}(c,{\kappa})\Big{|}_{\kappa\in(1+\frac{\sqrt{3}}{3},2)}=0$.

Similarly, the expression of the homoclinic Melnikov function of system (23)
for $0<{\kappa}<1-\frac{\sqrt{3}}{3}$ is

$\begin{array}[]{rcl}&&M_{\textmd{hom}}(c,{\kappa})\Big{|}_{\kappa\in(0,1-\frac{\sqrt{3}}{3})}\\\
&=&\frac{1}{c^{2}}\frac{\sqrt{3}}{3}\big{[}6\textmd{ln}2{{\kappa}}^{3}+6{{\kappa}}^{3}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)-18\textmd{ln}2{{\kappa}}^{2}\\\
&&-18{{\kappa}}^{2}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)+12\textmd{ln}2{\kappa}\\\
&&+12{\kappa}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)+2\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}~{}\big{]}\\\
&&-\frac{1}{c^{2}}\sqrt{3}{\kappa}\big{[}2{{\kappa}}^{2}\textmd{ln}(-\sqrt{2{\kappa}-{{\kappa}}^{2}})+3\textmd{ln}2{{\kappa}}^{2}-6{\kappa}\textmd{ln}(-\sqrt{2{\kappa}-{{\kappa}}^{2}})-9\textmd{ln}2{\kappa}\\\
&&+4\textmd{ln}(-\sqrt{2{\kappa}-{{\kappa}}^{2}})+6\textmd{ln}2\big{]}-\frac{\sqrt{3}}{5}\big{[}27{{\kappa}}^{4}\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}\\\
&&+30{{\kappa}}^{3}\textmd{ln}(-\sqrt{2{\kappa}-{{\kappa}}^{2}})+60{\kappa}\textmd{ln}(-\sqrt{2{\kappa}-{{\kappa}}^{2}})+15\textmd{ln}2{{\kappa}}^{3}\\\
&&-30{{\kappa}}^{3}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)-90{{\kappa}}^{2}\textmd{ln}(-\sqrt{2{\kappa}-{{\kappa}}^{2}})\\\
&&+90{{\kappa}}^{2}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)+30\sqrt{2}{\kappa}\sqrt{(3{{\kappa}}^{2}-6{\kappa}+2)^{3}}\\\
&&+294{{\kappa}}^{2}\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}-45\textmd{ln}2{{\kappa}}^{2}-198{{\kappa}}^{3}\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}\\\
&&-60{\kappa}\textmd{ln}(\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+2{\kappa}-2)+10\sqrt{2}\sqrt{(3{{\kappa}}^{2}-6{\kappa}+2)^{3}}\\\
&&-72{\kappa}\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}+30\textmd{ln}2{\kappa}-18\sqrt{6{{\kappa}}^{2}-12{\kappa}+4}\
\big{]}.\end{array}$ (57)

It is easy to verify that for $\kappa\in(0,1-\frac{\sqrt{3}}{3})$, there is
$c=c_{2}(\kappa)$ such that
$M_{\textmd{hom}}(c,{\kappa})\Big{|}_{\kappa\in(0,1-\frac{\sqrt{3}}{3})}=0$.

Secondly, the Melnikov method is also employed to detect the existence of kink
(anti-kink) wave solutions under small perturbation. The heteroclinic Melnikov
function of system (23) is defined as follows:

$\begin{array}[]{rcl}M_{\textmd{het}}(c,\kappa)&=&\displaystyle\int_{\Upsilon(\kappa)}(\frac{9}{2}z^{2}-9z+3+\frac{1}{c^{2}})y^{2}d\zeta=\displaystyle\int_{-\infty}^{+\infty}(\frac{9}{2}z^{2}-9z+3+\frac{1}{c^{2}})y^{2}d\zeta\\\
&=&\displaystyle\int_{0}^{2}(\frac{9}{2}z^{2}-9z+3+\frac{1}{c^{2}})ydz\\\
&=&\frac{\sqrt{3}}{2}\displaystyle\int_{0}^{2}(\frac{9}{2}z^{2}-9z+3+\frac{1}{c^{2}})z(2-z)dz=\displaystyle-\frac{2\sqrt{3}(3c^{2}-5)}{15c^{2}}.\end{array}$
(58)

###### Lemma 3

For $\kappa=0$ or $\kappa=2$, $M_{het}(c,\kappa)=0$ has a positive root
$c^{\ast}=\frac{\sqrt{15}}{3}$. Moreover, $\frac{\partial
M_{\textmd{het}}}{\partial c}\mid_{c=c^{\ast}}\neq 0$.

_Proof._ Clearly, it directly follows form (58) that
$M_{\textmd{het}}(c,\kappa)=0$ has a positive root
$c^{\ast}=\frac{\sqrt{15}}{3}$. Meanwhile,

$\frac{\partial M_{\textmd{het}}}{\partial
c}=-\frac{4\sqrt{3}}{5c}+\frac{4\sqrt{3}(3c^{2}-5)}{15c^{3}},$

and

$\frac{\partial M_{\textmd{het}}}{\partial
c}\mid_{c=c^{\ast}}=-\frac{12\sqrt{5}}{25}\neq 0.$

Proof of the Theorem 3: By Lemma 2 and Lemma 3, Theorem 3 follows immediately.

##  6 Numerical analysis

Numerical simulations are employed to confirm the theoretical results derived
in previous sections. Here, maple software $18.0$ is used.

Firstly, we simulate the existence of solitary wave solution of perturbed
$(1+1)$-dimensional DLWE (3). Let ${\kappa}=1.7$ and $\epsilon=0.01$. By (50)
and (53), we obtain $I(1.7)\approx 0.0648882124$ and
$J(1.7)\approx-0.02790800120$ such that
$c(1.7)=\sqrt{-\frac{I(1.7)}{J(1.7)}}\approx 1.524819859$. Then, according to
Theorem 3, we take $c=c(1.7)$ and let the initial value be
$(z(0),y(0))=(z_{+},0)=(1.309950494,0)$ (red point) which the homoclinic orbit
would pass through. The phase portraits $(z,y)$ and time history curves
$(\zeta,z)$ of system (23) are plotted in Fig. 8.

(a) $c=c(1.7)$ | (b) $(\zeta,z)$
---|---

Fig. 8 (a) Phase portraits with $c=c(1.7)$, (b) Time history curves of
$(\zeta,z)$ for system (23) for $\epsilon=0.01$, ${\kappa}=1.7$ and initial
value $(z(0),y(0))=(z_{+},0)=(1.309950494,0)$.

One can see that the homoclinic orbit and dark solitary wave solution of
system (23) still persist after a small singular perturbation.

Secondly, set ${\kappa}=0.3$ and $\epsilon=0.01$. By (54) and (55), we have
$I(0.3)\approx 0.0648882229$ and $J(0.3)\approx-0.02790784012$ such that
$c(0.3)=\sqrt{-\frac{I(0.3)}{J(0.3)}}\approx 1.524824377$. Taking $c=c(0.3)$
and initial value to be $(z(0),y(0))=(z_{-},0)=(0.690049506,0)$ (red point)
which the homoclinic orbit would pass through. The phase portraits $(z,y)$ and
time history curves $(\zeta,z)$ of system (23) are plotted in Fig. 9.

(a) $c=c(0.3)$ | (b) $(\zeta,z)$
---|---

Fig. 9 (a) Phase portraits with $c=c(0.3)$, (b) Time history curves of
$(\zeta,z)$ for system (23) for $\epsilon=0.01$, ${\kappa}=0.3$ and initial
value $(z(0),y(0))=(z_{-},0)=(0.690049506,0)$.

We see that the homoclinic orbit and bright solitary wave solution of system
(23) still exist after a small singular perturbation which verifies the
Theorem 3 $(i)$.

Next, let us to see the existence of kink and anti-kink wave solution of
system (23). Giving $\kappa=0$ or $\kappa=2$, $\epsilon=0.01$ and initial
value $(z(0),y(0))=(1,\pm\frac{\sqrt{3}}{2})$ (red point). We draw the phase
portraits of system (23) after taking $c=\frac{\sqrt{15}}{3}$ in Fig. 10. Form
the Fig. 10 (a), it is shown that the heteroclinic orbits still persist.
Correspondingly, a pair of kink and anti-kink wave solutions exists (see Fig.
10 (b) and (c)), which is consistent with Theorem 3 $(ii)$.

(a) $c=\frac{\sqrt{15}}{3}$ | (b) Kink wave | (c) Anti-kink wave
---|---|---

Fig. 10 Phase portraits with $c=\frac{\sqrt{15}}{3}$ for $\epsilon=0.01$,
initial value $(z(0),y(0))=(1,\pm\frac{\sqrt{3}}{2})$.

##  7 Fund Acknowledgement

This work was jointly supported by the Zhejiang Provincial Natural Science
Foundation of China (No.LZ23A010001), National Natural Science Foundation of
China under Grant (No. 11671176, 11931016), Natural Science Foundation of
Fujian Province under Grant (No. 2021J011148), Fujian Province Young Middle-
Aged teachers education scientific research project (No. JAT210454) and
Teacher and Student Scientific Team Fund of Wuyi University (Grant No.
2020-SSTD-003).

##  8  Conflict of Interest

The authors declare that they have no conflict of interest.

##  9 Data Availability Statement

No data was used for the research in this article.

##  10  Author contributions

All the authors have same contributions to the paper.

## References

  * [1] C. A. Buzzi, P. R. da Silva, M. Teixeira, Slow-fast systems on algebraic varieties bordering piecewise-smooth dynamical systems, Bull. Sci. Math. 136 (2012) 444–462.
  * [2] A. Y. Chen, L. N. Guo, W. T. Huang, Existence of kink waves and periodic waves for a perturbed defocusing mKdV equation, Qual. Theory Dyn. Syst. 17 (2018) 495–517.
  * [3] A. Y. Chen, L. N. Guo, X. J. Deng, Existence of solitary waves and periodic waves for a perturbed generalized BBM equation. J. Differential Equations 261 (2016) 5324–5349.
  * [4] A. Y. Chen, C. Zhang, W. T. Huang, Monotonicity of limit wave speed of traveling wave solutions for a perturbed generalized KdV equation, Appl. Math. Lett. 121 (2021) 107381.
  * [5] C. Chen, X. Tang, S. Lou, Solutions of a $(2+1)$-dimensional dispersive long wave equation, Phys. Rev. E 66 (2002) 036605.
  * [6] C. Chen, X. Tang, S. Lou, Soliton excitations and periodic waves without dispersion relation in shallow water system, Chaos, Solitons and Fractals, 16 (2003) 27-35.
  * [7] F. F. Cheng, J. Li, Geometric singular perturbation analysis of Degasperis-Procesi equation with distributed delay, Discrete Contin. Dyn. Syst. 41 (2021) 967–985.
  * [8] J. Carr, J. K. Hale, Abelian integrals and bifurcation theory, J. Differential Equations 59 (1985) 413–436.
  * [9] X. F. Chen, X. Zhang. Dynamics of the predator-prey model with the Sigmoid functional response, Stud. Appl. Math. 147 (1) (2021) 300–318.
  * [10] Y. Chen, Q. Wang, A new general algebraic method with symbolic computation to construct new travelling wave solution for the $(1+1)$-dimensional dispersive long wave equation, Appl. Math. Comput. 168 (2005) 1189–1204.
  * [11] G. Derks, S. Gils, On the uniqueness of traveling waves in perturbed Korteweg-de Vries equations, Jpn. J. Ind. Appl. Math. 10 (1993) 413–430.
  * [12] Z. J. Du, X. J. Lin, S. S. Yu. Solitary wave and periodic wave for a generalized Nizhnik-Novikov-Veselov equation with diffusion term (in Chinese), Sci. Sin. Math. 50 (2020) 1–22.
  * [13] Z. J. Du, Q. Qiao, The dynamics of traveling waves for a nonlinear Belousov-Zhabotinskii system, J. Differential Equations 269 (2020) 7214–7230.
  * [14] Z. J. Du, J. Liu, Y. L. Ren, Traveling pulse solutions of a generalized Keller-Segel system with small cell diffusion via a geometric approach, J. Differential Equations 270 (2020) 1019–1042.
  * [15] Z. J. Du, J. Li, Geometric singular perturbation analysis to Camassa-Holm Kuramoto-Sivashinsky equation, J. Differential Equations 306 (2022) 418–438.
  * [16] Z. J. Du, J. Li, X. W. Li, The existence of solitary wave solutions of delayed Camassa-Holm equation via a geometric approach, J. Funct. Anal. 275 (2018) 988–1007.
  * [17] E. G. Fan, Extended tanh-function method and its applications to nonlinear equations, Phys. Lett. A 277 (2000) 212–218.
  * [18] N. Fenichel, Geometric singular perturbation theory for ordinary differential equations, J. Differential Equations 31 (1979) 53–98.
  * [19] J. Guckenheimer, P. Holmes, Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields, Sprigner-Verlag, New York, 1997, pp. 184–193.
  * [20] I. S. Gradshteyn, I. M. Ryzhik, Table of integrals, series, and products. Amsterdam: Elsevier/Academic Press, Amsterdam, Translation edited and with a preface by Daniel Zwillinger and Victor Moll. Eighth edition, 2015, pp. 94–95.
  * [21] M. A. Han, Bifurcation theory of limit cycles, Beijing: Science press, 2013, pp. 252–289.
  * [22] A. Iuorio, N. Popović, P. Szmolyan, Singular perturbation analysis of a regularized MEMS Model, SIAM J. Appl. Dyn. Syst. 18 (2019) 661–708.
  * [23] C.K.R. Jones, Geometric Singular Perturbation Theory Dynamical Systems, in: Lecture Notes Math., vol. 1609, Springer, Berlin, 1995, pp. 44–120.
  * [24] J. Li, K. Lu, P.W. Bates, Normally hyperbolic invariant manifolds for random dynamical systems, Trans. Amer. Math. Soc. 365 (2013) 5933–5966.
  * [25] J. Li, K. Lu, P. W. Bates, Invariant foliations for random dynamical systems, Discrete Contin. Dyn. Syst. 34 (2014) 3639-3666.
  * [26] J. Li, K. Lu, P.W. Bates, Geometric singular perturbation theory with real noise, J. Differential Equations 259 (2015) 5137-5167.
  * [27] J. B. Li, Singular Nonlinear Travelling Wave Equations: Bifurcations and Exact Solutions, Beijing: Science, 2013.
  * [28] T. Ogama, Travelling wave solutions to a perturbed Korteweg-de Vries equation, Hiroshima Math. J. 24 (1994) 401–422.
  * [29] H. M. Qiu, L. N. Zhong, J. H. Shen, Traveling waves in a generalized Camassa-Holm equation involving dual-power law nonlinearities, Commun. Nonlinear Sci. Numer. Simul. 106 (2022) 106106.
  * [30] P. Szmolyan, Transversal heteroclinic and homoclinic orbits in singular perturbation problems, J. Differential Equations 92 (1991) 252–281.
  * [31] X. B. Sun, P. Yu, Periodic traveling waves in a generalized BBM equation with weak backward diffusion and dissipation terms, Discrete Contin. Dyn. Syst. Ser. B 24 (2019) 965-987.
  * [32] X. B. Sun, W. T. Huang, J. N. Cai, Coexistence of the solitary and periodic waves in convecting shallow water fluid, Nonlinear Anal. Real World Appl. 53 (2020) 103067.
  * [33] S. Wiggins, Introduction to Applied Nonlinear Systems and Chaos, Springer, New York, 1990.
  * [34] Z. S. Wen, On existence of kink and antikink wave solutions of singularly perturbed Gardner equation, Math. Methods Appl. Sci. 43 (7) (2020) 4422–4427.
  * [35] J. Weiss, M. Tabor, Carnevale G. The Painlevé property for partial differential equations, J. Math. Phys. 24 (1983) 522–526.
  * [36] Q. Wang, Y. Chen, H. Q. Zhang, A new Jacobi elliptic function rational expansion method and its application to $(1+1)$-dimensional dispersive long wave equation, Chaos Soliton Fract. 23 (2005) 477–483.
  * [37] T. Wu, J. Zhang. On modeling nonlinear long wave, in: L.P. Cook, V. Roytbhurd, M. Tulin (Eds.), Mathematics is for Solving Problems, SIAM, 1996, pp. 233.
  * [38] C. Wang, X. Zhang. Canards, heteroclinic and homoclinic orbits for a slow-fast predator-prey model of generalized Holling type III, J. Differential Equations 267 (2019) 3397–3441.
  * [39] X. Wu, M. K. Ni, Existence and stability of periodic contrast structure in reaction-advection-diffusion equation with discontinuous reactive and convective terms, Commun. Nonlinear Sci. Numer. Simul. 91 (2020) 105457.
  * [40] C. H. Xu, Y. H. Wu, L. X. Tian, B. L. Guo, On kink and anti-kink wave solutions of Schrodinger equation with distributed delay, J. Appl. Anal. Comput. 8 (2018) 1385–1395.
  * [41] Q. Yang, M. K. Ni, Asymptotics of the solution to a stationary piecewise-smooth reaction-diffusion equation with a multiple root of the degenerate equation, Sci. China Math. 65 (2022) 291–308.
  * [42] W. F. Yan, Z. R. Liu, Y. Liang, Existence of solitary waves and periodic waves to a perturbed generalized KdV equation, Math. Model. Anal. 19 (2014) 537–555.
  * [43] L. J. Zhang, J. D. Wang, E. Shchepakina, V. Sobolev, New type of solitary wave solution with coexisting crest and trough for a perturbed wave equation, Nonlinear Dyn. 106 (202) 3479–3493.
  * [44] L. J. Zhang, M. A. Han, M. J. Zhang, C. M. Khalique, A new type of solitary wave solution of the mKdV equation under singular perturbations, Int. J. Bifurcat. Chaos 30 (2020) 1–14.
  * [45] K. Zhu, Y. H. Wu, Z, P. Yu, J. H. Shen, New solitary wave solutions in a perturbed generalized BBM equation. Nonlinear Dyn. 97 (2019) 2413–2423.
  * [46] X. Zheng, Y. Chen, H. Q. Zhang, Generalized extended tanh-function method and its application to $(1+1)$-dimensional dispersive long wave equation, Phys. Lett. A 311 (2003) 145–157.
  * [47] X. Zeng, D. S. Wang, A generalized extended rational expansion method and its application to (1+1)-dimensional dispersive long wave equation, Appl. Math. Comput. 212 (2) (2009) 296-304.
  * [48] X. Zhang, Homoclinic, heteroclinic and periodic orbits of singularly perturbed systems, Sci. China Math. 62 (2019) 1687–1704.


# A Deep Reinforcement Learning Environment for
Particle Robot Navigation and Object Manipulation

Jeremy Shen1∗, Erdong Xiao2∗, Yuchen Liu2, Chen Feng2†
https://ai4ce.github.io/DeepParticleRobot ∗ indicate equal contributions.1
Stuyvesant High School, New York, NY 10280, USA<EMAIL_ADDRESS>New York
University, Brooklyn, NY 11201, USA {ex276, yl5680<EMAIL_ADDRESS>Corresponding author<EMAIL_ADDRESS>
###### Abstract

Particle robots are novel biologically-inspired robotic systems where
locomotion can be achieved collectively and robustly, but not independently.
While its control is currently limited to a hand-crafted policy for basic
locomotion tasks, such a multi-robot system could be potentially controlled
via Deep Reinforcement Learning (DRL) for different tasks more efficiently.
However, the particle robot system presents a new set of challenges for DRL
differing from existing swarm robotics systems: the low degrees of freedom of
each robot and the increased necessity of coordination between robots. We
present a 2D particle robot simulator using the OpenAI Gym interface and
Pymunk as the physics engine, and introduce new tasks and challenges to
research the underexplored applications of DRL in the particle robot system.
Moreover, we use Stable-baselines3 to provide a set of benchmarks for the
tasks. Current baseline DRL algorithms show signs of achieving the tasks but
are yet unable to reach the performance of the hand-crafted policy. Further
development of DRL algorithms is necessary in order to accomplish the proposed
tasks.

##  I Introduction

Innovations in hardware and computing have sparked increased interest and
research in swarm robotics. Swarm robot systems provide a more adaptable,
robust, and scalable alternative to current state-of-the-art robots. While
traditional swarm robots, like drones, maintain a basic level of capabilities
in each individual, a new class of swarm robots with simpler hardware has been
proposed: the particle robot. Existing work in particle robots can be
classified into two categories: biological-inspired particle robots, like
viruses and sea urchins, that can operate independently [1], and statistical-
mechanics-inspired particle robots that are incapable of self-locomotion [2].
Both types of particle robots pose a collaborative control task, but in this
paper, we will focus on the statistical mechanics-inspired particle robot.
This type of robot is unable to complete any task on its own, so in order to
utilize the particle robot, we must control an extensive multi-robot system.

We present a Particle Robot simulator to study reinforcement learning methods
on this type of robots. The particle robot [2] we study in this simulator is a
disk-shaped robot only capable of changing its radius (increasing/decreasing
its size). It switches between two states: fully expanded (largest size) or
fully contracted (smallest size), by moving paddles attached to its central
body (see Figure 1). These paddles are not in contact with the ground meaning
individual particles are incapable of self-propulsion. However, as a part of a
particle robot system, the collective particle robots can leverage the
internal forces from the paddle movement between particle robots and the
external force from the friction from the ground acting on each particle
robot’s central body. It has been proven that if the particle robot system’s
expansion-contraction cycle follows a longitudinal wave, the particle robots’
motion will be in the opposite direction to the wave (see Figure 2). For
further detail on the dynamics of the particle robot system, please check
section III-A.

Figure 1: Upper left: Contracted single particle robot. Lower left: Expanded
single particle robot. Right: 56-particle robot system with mixed states
Figure 2: Longitudinal wave policy on a 25-particle robot system. The
direction of the expansion-contraction wave is downwards, and the motion of
the particle robots is upwards. Each image is an individual frame.

Control methods for the particle robot can be mathematically derived, like the
aforementioned “wave-policy”, but it would be a difficult and tedious task to
create a new control method for each new robot design and task. The appeal of
swarm robots is their adaptability, so it makes sense for particle robot
control to be easily repurposed and reapplied. This motivates the application
of deep reinforcement learning as a control method for particle robots. Since
the introduction of Deep Q-Networks (DQN) [3], new Deep Reinforcement Learning
(DRL) algorithms have become widespread, including Advantage Actor-Critic
(A2C) [4] and Proximal Policy Optimization (PPO) [5]. The combination of
traditional Q-Learning or Policy Gradient methods and the powerful modeling
ability of deep neural networks has improved the performance of DRL
algorithms, allowing them to tackle increasingly difficult problems. DRL has
become well suited for complex robot control problems, thus it is reasonable
to study if DRL methods could automatically find optimal control policies for
the particle robot system. Our particle robot simulator poses a novel swarm
robot environment with unique constraints on individual robot capabilities to
test and develop DRL methods.

As the first step of studying DRL for particle robots, we choose to use a
centralized control method to facilitate collaborative particle robots. We
consider this a “super-agent” setup, where there exists an outside controller
that possesses full observability of the world and the ability to send
simultaneous instructions to every particle robot. Although the particle robot
system seems naturally more suited for decentralized RL, and the super-agent
setup seems to make the problem easier due to the full observability and
coordination, we find our setup already non-trivial and challenging for
existing DRL baselines. Note that our environment supports the decentralized
setup, but in this work, we focus on thoroughly studying the super-agent
control before advancing to a more challenging problem.

Our simulator contains tasks varying in difficulty from point-to-point
navigation to object manipulation, and allows for varying the number and
states of the particle robots. By using the OpenAI gym interface, it is
possible to modify observation and action spaces, tasks and challenges, and
reward methods, and it is easy to apply reinforcement learning techniques on
the particle robots.

We also benchmark four baseline algorithms: (1) a hand-crafted control policy,
(2) DQN, (3) A2C, and (4) PPO, over four different tasks: (1) simple
navigation, (2) obstacle navigation, (3) navigation with unresponsive
particles, and (4) object manipulation.

In summary, our contributions are the following:

  * •

introduced a new deep reinforcement learning environment for particle robot
research.

  * •

implemented multiple, challenging tasks for the development of particle
robots.

  * •

benchmarked three existing deep reinforcement learning algorithms on our
tasks.

##  II Related Work

Particle-robot-like Swarm Robot Control. There are existing low-capability
robot systems similar to particle robots which are not controlled using deep
reinforcement learning. Research is focused on exploiting the randomness of
collective low-level actions to generate movement in larger robot systems. The
Smarticle [6] system implements statistical models and unsupervised learning
to create asymmetries in smarticle collisions by activating and deactivating
certain smarticles. Smarticle control was advanced with a more physics-
inspired model for self organization [7]. Deblais et al. [8] considers a
bounded collection of rod-like robots which can form mobile and deformable
superstructures. Boudet et al. [9] proposes a Langevin’s equation-inspired
dynamics model for this robot scaffold. These systems are collections of
robots bounded by a soft or hard container, creating parallels to fluid
behaviors. The particle robots in our environment are not container bounded
but use the magnetic attraction for a similar adhesive purpose as in [2].

Wang et al. [10] propose the Virtual Particle Exchange, a localization
algorithm that relies on simple sensors and calculations and can anonymize and
limit robot communications. Mayya et al. [11] considers equilibrium
thermodynamics in brushbot systems to take advantage of the differences in
robot densities. Yang et al. [12] creates vortex-like paramagnetic
nanoparticle swarms to track and control microbots. Like particle robots,
individual robots in these systems possess little to no self-locomotion
abilities. Unlike previous works that handcraft control policies from physics
rules, with our particle robot environment, we aim to use deep reinforcement
learning to discover explicit control policies. The robots used in our
simulator are inspired by particle robots [2], a disk-shaped, biological-cell-
inspired robot only capable of changing its size and without self-locomotion
abilities. Researchers hand-derived a sinusoidal-function control policy (an
expansion/contraction cycle determined from distance to goal position) for
large swarms of particle robots to complete various navigation tasks.

Deep Reinforcement Learning-based Control Method. Deep reinforcement learning
has become popular for solving multi-agent system control problems. Collision
avoidance between robots is a critical problem for multi-robot systems; Fan et
al. [13] and Long et al. [14] have developed reinforcement learning based
control methods for multi-robot navigation tasks by combining a collision
penalty term and the time cost term in their reward function, enabling the
robots to learn a time-optimal and collision-free trajectory. DRL is capable
of coordinating systems of robots that possess significantly more locomotive
abilities than our particle robots. Drone system control is a more difficult
multi-agent control problem which presents a larger state and action space
since drones can move freely in 3D space, but deep reinforcement learning is
still viable for cooperative multi-drones tasks. Real-time trajectory planning
of multi-drone systems [15] trains a large amount of drones (with limited
communication) using a PPO algorithm to find a collision-free trajectory. This
demonstrates the PPO algorithm’s capability of solving cooperative navigation
problems, which is similar to our proposed Particle Robot tasks. Wang et al.
[16] introduces an application of DRL in pattern formation: a multi-robot
system learns to configure in given patterns such as lines, circles, and more
complex shapes; these tasks could show the potential of using DRL to solve
similar cooperation tasks for particle robot systems. Moreover, Ding et al.
[17] shows that robots can learn a cooperative strategy in the multi-arms
manipulation task. Zhu et al. [18] discusses using DRL to solve the robot
flocking control problem (a migrating-birds-like robots system). Sun et al.
[19] give a systemic review of DRL in the robot motion planning field. DRL has
proven to be a successful method of discovering robot control policies, and it
is considered to have the potential to solve the particle robots’ control
problem.

Related Reinforcement Learning Environments. It is necessary for our
environment to preserve realism and maximize adaptability so that the control
methods developed in the simulator could also work with real-world particle
robots. Thus, our Deep Particle Robot simulator (DPR) features a continuous
world environment, centralized and decentralized communication structures,
navigation and manipulation tasks, and a physics engine (PhysE). Our
environment can be easily modified to add new robot designs, challenges, and
goals, and the Pymunk physics engines simulates realistic forces and
collisions. Compared to our environment, there are various existing
environments for multi-agent reinforcement learning, but none of them are
directly suitable for simulating particle robots (see Table I). Multi-Agent-
Reinforcement-Learning-Environment [20] offers a framework for multi-agent
problems like playing soccer, navigating mazes, and manipulating boxes. MAgent
[21] is one of the first environments to support millions of agents and offer
numerous variants of communication structures and tasks. Multi-agent-Particle-
Envs [22] has a physics engine collision handler and both a continuous and
discretized action space. However, these environments are simplified to a
discretized state space and non-continuous world, limiting simulator realism.
Multi-Agent-Emergence-Environments [23] is solely a hide-and-seek game and
every agent within a team shares its observations. Playground [24], an
environment for the multi-agent game Pommerman, emphasizes coordination and
strategy but uses a grid world and a binary action space. There exists
specialized simulators for video games like Neural-MMOs [25], SMAC [26] for
Starcraft II, and Google Research Football [27] for soccer, and while these
simulators are used for multi-agent reinforcement learning research, they are
unable to be applied in particle robots.

TABLE I: Related Multi-agent Reinforcement Learning Environments

| Statea | Commb | PhysEc | Navd | Manipe | #Agentsf | PRg
---|---|---|---|---|---|---|---
[20] | D | C | $\times$ | ✓ | ✓ |  $<$10 | $\times$
[21] | D | D | $\times$ | ✓ | $\times$ |  $>$100 | $\times$
[22] | D | D | 2D | ✓ | $\times$ | 10-100 | $\times$
[23] | C | C | 3D | ✓ | ✓ |  $<$10 | $\times$
[24] | D | D | $\times$ | $\times$ | $\times$ |  $<$10 | $\times$
[25] | C | C+D | $\times$ | $\times$ | ✓ |  $<$100 | $\times$
[26] | C | D | $\times$ | ✓ | $\times$ |  $<$100 | $\times$
[27] | C | C+D | $\times$ | $\times$ | ✓ | 10-100 | $\times$
Ours | C | C+D | 2D | ✓ | ✓ |  $<$100 | ✓

  * a

Continuous state (C) or discretized state/grid world (D)

  * b

Centralized (C) or decentralized (D) communication structure

  * c

2D, 3D, or No Physics Engine

  * d

Navigation Task

  * e

Manipulation Task

  * f

Number of Agents Supported

  * g

Low-level locomotion-inhibited robotics environment

##  III Environment Setup

Our Particle Robots environment uses Pymunk, a 2D Physics library, to simulate
the environment, and Pygame for visualization. The environment also uses the
OpenAI Gym interface, allowing for easy reinforcement learning integration.

###  III-A Robot Design

The basic Particle Robot consists of a circular central body with sixteen
paddles attached to it. Each paddle has a range of 90 degrees, from being
tangent to being perpendicular to the central body, and is powered with a
simple motor. A second set of paddles is attached to the end of the first set
of paddles which rotate to remain parallel to the surface of the central body
(see Figure 3).

The paddles of each particle robot are lined with magnets in order to exert an
attractive force on its neighboring particle robots. The magnetic force is
necessary to keep the particle robot system together.

The original paper of Particle Robot [2] has shown that the force from the
paddles, friction force, and magnetic force all must lie within a specific
range but has not provided definite values:

$F_{f,i}<F_{paddles/magnetic}<\sum_{j=1}^{n}F_{f,j},\quad\forall i,$ (1)

where $F_{f,i}$ is the magnitude of the friction force on the $i$-th particle
robot, $F_{paddles}$ is the force the paddles exert, and $F_{magnetic}$ is the
magnetic force.

Figure 3: Anatomy of Particle Robot

###  III-B Observations

Observations are comprised of any combination of information derived from a
top-down camera. The top-down camera is centered over the goal, and the
position, velocity (linear and rotational), distance to goal, and state of
each particle robot and any objects can be captured. It is also possible to
directly use each Pygame frame as the observation.

###  III-C Actions

Individual robot. The basic particle robot has two states, fully contracted
and fully expanded. Actions cannot be overridden, so when the particle robot’s
paddles are in motion, all proposed actions will be voided. This discretized
action space is used for our experiments, but it is possible to use a
continuous action space to set an exact angle for the paddle’s state, allowing
the particle robot to partially expand.

Super-agent. Our super-agent setup treats the individual particle robots
within its system as another degree of freedom. Each particle robot is
assigned an index, and the action space can be a list of binary states (multi-
binary action space) or a decimal integer (discrete action space) to be
converted into a binary value. Each particle robot will take its corresponding
action based on its index.

###  III-D Reward Design Motivation

New reward functions can easily be added to the Gym interface. The reward
functions are the same for any type of agent: the center of mass of the
Particle Robot system, individual particle robots, or manipulated objects.

We present a piecewise reward function that is used in our experiments where
$\mathbf{G}$ is the coordinate of the goal, $\mathbf{A}_{t}$ is the coordinate
of the Agent’s position at time $t$, and $\theta$ is the angle between
$\mathbf{G}-\mathbf{A}_{t-1}$ and $\mathbf{A}_{t}-\mathbf{A}_{t-1}$:

$r_{t}=\begin{cases}||\mathbf{A}_{t}-\mathbf{A}_{t-1}||\frac{d\cos(\theta)}{||\mathbf{G}-\mathbf{A}_{t}||},\quad&||\mathbf{G}-\mathbf{A}_{t}||\leq
d;\\\
||\mathbf{A}_{t}-\mathbf{A}_{t-1}||\cos(\theta),\quad&||\mathbf{G}-\mathbf{A}_{t}||>d.\\\
\end{cases}$ (2)

Note that $||\mathbf{A}_{t}-\mathbf{A}_{t-1}||\cos(\theta)$ is the length of
the agent’s displacement vector (previous position to current position)
projected onto the ideal vector (previous position to goal position). This is
the distance the agent travels in the correct direction. When the agent is
within a threshold distance $d$, the magnitude of the reward is scaled
inversely proportional to the distance. This helps the model converge as the
agent finalizes its approach to the goal.

Figure 4: Left Simple navigation task. Middle Obstacle navigation task: gray
blocks are the gates. Right Object manipulation task: orange ball is the
target object. The goal for all three tasks is the green square in the lower
left

###  III-E Navigation Task

####  III-E1 Problem formulation

For an N-particle robot system. The state space at time step t, where
$0<t<T=2,500$:

$s_{t}=(\bm{x}_{t},\bm{y}_{t},\bm{\dot{x}}_{t},\bm{\dot{y}}_{t}),$ (3)

where $\bm{x}_{t}$ and $\bm{y}_{t}$ are the x-coordinates vector and
y-coordinates vector respectively of the particle robots.

The action space is N-dimensional binary, where each particle robot can only
contract (0) or expand (1) and the action of the entire system is the joint
action of all individual particle robots.

The reward $r_{t}$ takes on the same form as equation (2), where the agent is
the center of mass of the entire particle robot system. The score function is:

$J(\theta)=\sum_{k=0}^{T-1}r_{k+1}\gamma^{T-k-1},$ (4)

and our goal is to find a policy function that maximizes the score:

$\pi(s;\theta)=\mathop{\arg\max}_{\theta}J(\theta).$ (5)

####  III-E2 Navigation Variations

Simple navigation. This is the plainest task designed to demonstrate the
Particle Robots ability to perform basic locomotion. The Particle Robots are
initialized in a grid at a set start position, and the task is complete when
the center of mass of the Particle Robot system reaches the goal position (see
Figure 4).

Obstacle navigation. The obstacle navigation task test the Particle Robots’
ability to respond to external forces. A gate is set up between the start and
goal locations such that the opening is greater than the width of one Particle
Robot and less than the width of the entire Particle Robot system (see Figure
4). The shape of the particle robot system will be deformed when going through
the gate.

Navigation with unresponsive particle robots. The navigation task can be made
more challenging by ”killing” a fixed number of Particle Robots. Any action
sent to an unresponsive Particle Robot is ignored, and the unresponsive
Particle Robot will remain in a fully contracted position. Unresponsive
Particle Robots still have magnetic attractions, meaning these particle robots
will travel with the system. This demonstrates the particle robot system’s
ability to continue operating despite some faulty robots.

###  III-F Object Manipulation Task:

####  III-F1 Problem formulation

For an N-particle robot system with a target object. At time step t, where
$0<t<T=2,500$:

$s_{t}=(\bm{x}_{t},x_{t}^{o},\bm{y}_{t},y_{t}^{o},\bm{\dot{x}}_{t},\dot{x}_{t}^{o},\bm{\dot{y}}_{t},\dot{y}_{t}^{o}).$
(6)

Note that this state space is nearly identical to the state space introduced
in equation (3), except that the position and velocity of the target object
are also included, denoted with the superscript “o”. The action space of the
object manipulation task is identical to that of the navigation tasks.

The reward $r_{t}$ takes on the same form as equation (2), where the agent is
the manipulated object.

The score function for the object manipulation task is identical to that of
the navigation tasks (reference equation (4)), and our goal is to find a
policy function of the same form as equation (5).

There is significant overlap between the problem formulations of the
navigation tasks (section III-E1) and the object manipulation task (section
III-F1). The key difference is that the reward of the navigation tasks is
dependent on the particle robots, while the reward of the manipulation task is
solely dependent on the manipulated object.

####  III-F2 Manipulation Task

The object manipulation task tests the Particle Robots ability to apply an
external force on a non-Particle Robot object (see Figure 4). The Particle
Robots are required to move an object from a start position to a goal
position. The shape, dimensions, and size of this object can be changed to
simulate different real world objects (ball, brick, door, etc.). To simplify
the problem, the object is initialized adjacent to and in between the particle
robot system and the goal location. The difficulty of this task is that the
Particle Robots must learn and maintain the contact point between the system
and the object.

###  III-G Single-agent vs. Multi-agent Environment Setup

The proposed tasks remain identical between the single and multi-agent setups,
but the control methods for the Particle Robots change.

The single-agent environment regards the $N$ agents system as a super-agent
with $N$ degrees of freedom, therefore the action space is the union of the
action space for a single agent in the multi-agent setup. This assumption
could simplify the problem as a single agent reinforcement learning, and there
are already many existing baseline algorithms to help solve the problem.
However, this assumption also naturally limits the potential of the multi
robots system, for example, the number of particle robots must remain
constant.

There are many possible configurations for the multi-agent environment. For
our experiments, we use a network architecture where the observation and
reward is shared among all the particle robots. It is possible to configure
the environment with other structures such as a fully decentralized,
competitive setting.

##  IV Experiments

###  IV-A Baseline Methods

We benchmarked the following four baseline methods for particle robot control.

Hand-crafted algorithm: We implemented the wave policy discussed in section
III-A and shown in equation (1) and figure 2. The wave policy function returns
a sequence of actions that simulates one full expansion/contraction cycle.

DQN, A2C, PPO: We used stable-baselines3[28] to implement a discrete setting
for DQN, A2C, and PPO.

###  IV-B Evaluation Metrics

Each algorithm is given a fixed number of steps to run (2500), and we evaluate
a model’s performance given the following results.

Total distance traveled: Compared with the net displacement, the total
distance the particle robot system travels gives us information about how
efficient the system’s motion is.

Net Displacement: The total displacement of the particle robot system
(distance between initial and final coordinates) gives us insight into the
system’s locomotive abilities. This does not take into consideration the
direction of motion.

Projected Displacement: The displacement vector is projected onto the ideal
trajectory (vector pointing from start to goal). This is the distance the
system moves in the correct direction.

$Projection=||\bm{A}_{f}-\bm{A}_{0}||\cos(\theta)$ (7)

Where, in equation (7), $\bm{A}_{t}$ is the agent’s position at step
$t\in[0,f]$ and $\theta$ is the angle between the agent’s displacement vector
and the ideal trajectory.

###  IV-C Significant Hyperparameter Motivations

Learning rate. The initial learning rate must be large to take advantage of
the early exploration steps but the learning rate must eventually decrease in
order to converge to an optimal policy which resulted in us using a linear
schedule learning rate (from 0.0075 to 0).

Network architecture. After experimenting with different architectures, we
discovered that a Q-network/Policy network architecture of [64, 128, 128, 128,
64] (4 hidden layers) resulted in the quickest convergence to an optimal
policy.

Training length. We chose 10,000,000 to be our total timesteps during
training. This ensures that by the end of the training, the value loss and
policy gradient loss (if it exists) converges to 0.

###  IV-D Results

---
---

---
---

Figure 5: Sample trajectories for (starting upper-left, clockwise) Simple
Navigation, Obstacle Navigation, Unresponsive Particles, and Object
Manipulation tasks. Each figure simulates the trajectory of the four baseline
models all starting from the same point.

---
---

---
---

---
---

---
---

---
---

---
---

Figure 6: Benchmark results for all baseline models (handcrafted baseline,
DQN, A2C, PPO): median performance plotted in bar graphs and maximum and
minimum performance plotted in error bars.

Navigation tasks. From Figure 6, the DRL algorithms are outperformed
drastically by the handcrafted algorithm. The Projected Displacement shows
that in eight out of nine navigation tasks, the particle robots move less than
one unit in the correct direction. However in Row 2, we see that the net
displacement of the particle robots in the simple navigation task for DQN,
A2C, and PPO is a significant amount, 10.8%, 9.1% , and 6.6% of the
handcrafted algorithm respectively. Figure 5 qualitatively corroborates these
findings; the particle robot systems move significantly but in the incorrect
direction during trial runs. From Row 1, we also notice that the total
distance traveled by the handcrafted algorithm is nearly identical to the DRL
algorithms. This suggests that existing DRL algorithms are beginning to solve
a self-propulsion problem but mainly struggle with navigation.

The handcrafted algorithm’s performance significantly decreases as the tasks
become more challenging. This is evident in both columns 2 and 3, where the
handcrafted algorithm is 75% as effective in obstacle navigation and 56% as
effective with unresponsive particle robots. However DRL algorithms do not
have this clear drop-off with harder tasks. Using the net displacement metric,
DQN and PPO performed better in the obstacle navigation than the simple
navigation task. A2C and PPO performed better with unresponsive particles.
PPO’s net displacement performance relative to the handcrafted algorithm
increased from 9.1% in the simple navigation task to 25.4% in the obstacle
navigation task and 26.5% with unresponsive particles. In potential real-world
particle robot systems, DRL algorithms posses a greater ability to adapt to
unforeseen challenges than the handcrafted algorithm.

Manipulation tasks. From the last column of Figure 6, we notice that DQN and
PPO are 9% as effective as the handcrafted algorithm, and A2C is 34% as
effective. A2C was the best performing DRL model in the manipulation task. For
each policy, the evaluation metrics: total distance, net displacement, and
projected displacement, are all similar to each other. The particle robots are
only able to push the target object, meaning any oscillatory motion by the
particle robots is not reflected in the target object’s motion. In contrast to
the navigation tasks, the particle robots are able to move and navigate the
object ($\frac{\text{net displacement}}{\text{projected displacement}}>0.97$
for all DRL models, meaning the target object is being pushed mostly in the
correct direction), but DRL still lags behind handcrafted control.

##  V Conclusions

Current deep reinforcement learning algorithms are far from matching the
performance of a handcrafted control policy. These DRL algorithms struggle
with navigation but possess limited locomotive abilities. The most encouraging
result from the DRL algorithms was the performance in more challenging
settings was similar to the simple navigation task, indicating a high level of
adaptability.

A potential reason for the unsatisfying performance of the particle robots is
that the probability of the super-agent guessing a good action during training
is very low. Therefore, the training is heavily dependent on the random search
of the super-agent’s action space, which is inefficient. In the future, we
will study imitation learning on the handcrafted wave policy and examine
supervised learning approach.

We hope to stimulate research into deep reinforcement learning for swarm
robots by proposing a suite of navigation and construction tasks and
benchmarking the performance of three state-of-the-art DRL algorithms and a
handcrafted policy. We plan on introducing new particle robot designs with
similar functions to study how mixed particle robot systems operate. We also
plan to expand our environment into 3d and improve robot realism. Although the
super-agent assumption is still a challenging task, we plan on benchmarking
the environment in a decentralized multi-agent setting in order to fully
extract the advantages of a low-level swarm robot system.

Acknowledgment. This research is supported by the NSF CPS program under
CMMI-1932187.

## References

  * Mateos [2020] L. A. Mateos, “Particle robots a new species of hybrid bio-inspired robotics,” _arXiv:2003.08289_ , 2020.
  * Li et al. [2019] S. Li, R. Batra, D. Brown, H.-D. Chang, N. Ranganathan, C. Hoberman, and D. R. . H. Lipson, “Particle robotics based on statistical mechanics of loosely coupled components,” _Nature 567, 361–365_ , 2019.
  * Mnih et al. [2015] V. Mnih, K. Kavukcuoglu, D. Silver, A. A. Rusu, J. Veness, M. G. Bellemare, A. Graves, M. Riedmiller, A. K. Fidjeland, G. Ostrovski, S. Petersen, C. Beattie, A. Sadik, I. Antonoglou, H. King, D. Kumaran, D. Wierstra, and S. L. . D. Hassabis, “Human-level control through deep reinforcement learning,” _Nature 518, 529–533_ , 2015.
  * Mnih et al. [2016] V. Mnih, A. P. Badia1, M. Mirza, A. Graves, T. Harley, T. P. Lillicrap, D. Silver, and K. Kavukcuoglu, “Asynchronous methods for deep reinforcement learning,” _arXiv:1602.01783v2_ , 2016.
  * Schulman et al. [2017] J. Schulman, F. Wolski, P. Dhariwal, A. Radford, and O. Klimov, “Proximal policy optimization algorithms,” _arXiv:1707.06347v2_ , 2017.
  * Savoie et al. [2019] W. Savoie, T. A. Berrueta, Z. Jackson, A. Pervan, R. Warkentin, S. Li, T. D. Murphey, and K. W. . D. I. Goldman, “A robot made of robots: Emergent transport and control of a smarticle ensemble,” _Science Robotics, Volume 4, Issue 34_ , 2019.
  * Chvykov et al. [2021] P. Chvykov, T. A. Berrueta, A. Vardhan, W. Savoie, A. Samland, T. D. Murphey, K. Wiesenfeld, D. I. Goldman, J. L., and England, “Low rattling: A predictive principle for self-organization in active collectives,” _Science, Vol 371, Issue 6524 • pp. 90-95_ , 2021.
  * Deblais et al. [2018] A. Deblais, T. Barois, T. Guérin, P.-H. Delville, R. Vaudaine, J. S. Lintuvuori, J.-F. Boudet, J.-C. Baret, and H. Kellay, “Boundaries control collective dynamics of inertial self-propelled robots,” _Physical Review Letters, 120 (18), pp.188002 (1-5)_ , 2018.
  * Boudet et al. [2021] J. F. Boudet, J. Lintuvuori, C. Lacouture, T. Barois, A. Deblais, K. Xie, S. Cassagnere, B. Tregon, D. B. Bruckner, J. C. Baret, and H. Kellay, “From collections of independent, mindless robots to flexible, mobile and directional superstructures,” _SCIENCE ROBOTICS Vol 6, Issue 56_ , 2021.
  * Wang et al. [2020] J. Wang, T. Wang, W. Wang, X. Dong, and Y. Wang, “Distributed localization without direct communication inspired by statistical mechanics,” _arXiv:2006.02658v1_ , 2020.
  * Mayya et al. [2019] S. Mayya, G. Notomista, D. Shell, S. Hutchinson, and M. Egerstedt, “Non-uniform robot densities in vibration driven swarms using phase separation theory,” in _2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2019, pp. 4106–4112.
  * Yang et al. [2020] L. Yang, J. Yu, and L. Zhang, “Statistics-based automated control for a swarm of paramagnetic nanoparticles in 2-d space,” _IEEE Transactions on Robotics_ , vol. 36, no. 1, pp. 254–270, 2020.
  * Fan et al. [2018] T. Fan, P. Long, W. Liu, and J. Pan, “Fully distributed multi-robot collision avoidance via deep reinforcement learning for safe and efficient navigation in complex scenarios,” _arXiv:1808.03841v1 [cs.RO] 11 Aug 2018_ , 2018.
  * Long et al. [2018] P. Long, T. Fan, X. Liao, W. Liu, H. Zhang, and J. Pan, “Towards optimally decentralized multi-robot collision avoidance via deep reinforcement learning,” _International Conference on Robotics and Automation (ICRA), pp. 6252-6259_ , 2018.
  * Arya Kumar [2019] A. Z. Arya Kumar, “Using reinforcement learning for real-time trajectory planning of aerial multi agent systems,” _TJHSST Computer Systems Research Lab, Senior Research_ , 2019.
  * Wang et al. [2019] J. Wang, J. Cao, M. Stojmenovic, M. Zhao, J. Chen, and S. Jiang, “Pattern-rl: Multi-robot cooperative pattern formation via deep reinforcement learning,” _International Conference on Machine Learning and Applications (ICMLA)_ , 2019\.
  * Ding et al. [2020] G. Ding, J. J. Koh, K. Merckaert, B. Vanderborght, M. M. Nicotra, C. Heckman, A. Roncone, and L. Chen, “Distributed reinforcement learning for cooperative multi-robot object manipulation,” _arXiv:2003.09540_ , 2020.
  * ZHU et al. [2020] P. ZHU, W. DAI, W. YAO, J. MA, Z. ZENG, and H. LU, “Multi-robot flocking control based on deep reinforcement learning,” _IEEE Access, vol. 8, pp. 150397-150406, 2020_ , 2020.
  * SUN et al. [2021] H. SUN, W. ZHANG, R. YU, and Y. ZHANG, “Motion planning for mobile robots—focusing on deep reinforcement learning: A systematic review,” _IEEE Access, vol. 9, pp. 69061-69081, 2021_ , 2021.
  * Jiang [2018] S. Jiang, “Multi-agent-reinforcement-learning-environment,” https://github.com/Bigpig4396/Multi-Agent-Reinforcement-Learning-Environment, 2018\.
  * Zheng et al. [2017] L. Zheng, J. Yang, H. Cai, W. Zhang, J. Wang, and Y. Yu, “Magent: A many-agent reinforcement learning platform for artificial collective intelligence,” _arXiv:1712.00600_ , 2017.
  * Lowe et al. [2020] R. Lowe, Y. Wu, A. Tamar, J. Harb, P. Abbeel, and I. Mordatch, “Multi-agent actor-critic for mixed cooperative-competitive environments,” _arXiv:1706.02275v4_ , 2020.
  * Baker et al. [2020] B. Baker, I. Kanitscheider, T. Markov, Y. Wu, G. Powell, B. McGrew, and I. Mordatch, “Emergent tool use from multi-agent autocurricula,” _arXiv:1909.07528_ , 2020.
  * pla [2019] “Playground,” https://github.com/MultiAgentLearning/playground, 2019.
  * Suarez [2019] J. Suarez, “Neural mmos,” https://github.com/openai/neural-mmo, 2019.
  * Samvelyan [2019] M. Samvelyan, “Starcraft ii,” https://github.com/oxwhirl/smac, 2019.
  * foo [2019] “Google research football,” https://github.com/google-research/football, 2019.
  * Raffin et al. [2019] A. Raffin, A. Hill, M. Ernestus, A. Gleave, A. Kanervisto, and N. Dormann, “Stable baselines3,” https://github.com/DLR-RM/stable-baselines3, 2019\.