Two-dimensional fractional algebraic internal solitary waves model and its solution

Abstract

With the development of theory and advancement of scientific research, fractional calculus has begun to be considered as a method for the study of physical systems. In this work, based on the basic system of equations for internal solitary waves, we have derived two-dimensional Benjamin–Ono equation. Then, the integer-order two-dimensional Benjamin–Ono equation is transformed into the time-fractional Benjamin–Ono equation. To study the properties of the algebraic internal solitary waves, we discuss the conservation laws of the new model. Also, using the Hirota bilinear method, the derived new model is solved. Finally, we explore the characteristics of motion of the algebraic internal solitary waves with the help of the multi-soliton solutions.

Introduction

Nonlinear partial differential equations have been used to describe a wide range of physics phenomena as a model for the evolution and interaction of nonlinear waves.^1–8 Particularly, internal solitary waves are a commonly occuring feature along the belt of the continental edge with suitable stratification, current, and topography. The research on the internal waves is a research hot spot and various applications have been found in ocean engineering, global environment, and underwater acoustic exploration.^9–15 Algebraic solitary wave theory is an important part of the study of nonlinear partial differential equation. This type of equation admits exact solutions representing “algebraic” solitary waves.^16–18 In recent years, the application of fractional differential equations has attracted increasing attention in many different areas.^19–26 Researchers have discovered that the derivatives and integrals of fractional order models are suitable for describing various physical phenomena.

Conservation laws are a mathematical formulation which indicates that the total amount of a certain physical quantity remains the same during the evolution of a physical system. Conservation laws provide much information as regards systems simulated by differential equations. ²⁷ Many fractional conservation laws were obtained in the past two decades.^26–32

Many solution methods have been found and used to solve the fractional order equation, for instance, the homotopy perturbation method, ³³ Hirota bilinear method,^34,35 computational approach, ³⁶ and others.^37–44 Various phenomena can be explained via the application of the solutions given by the above methods.^45–47 Approximation and numerical methods such as variational iteration and adomian decomposition have been used because most fractional differential equations do not have exact analytic solutions. In recent years, Odibat and Momani have realized the solution method of fractional order nonlinear ordinary differential equation. They use the homotopy perturbation method to derive analytic approximate solution of the nonlinear partial differential equation of a fractional order system.

Since the concept of soliton was first put forward, the solution of integral partial differential development equation has been one of the most important research topics.^48–51 For fractional partial differential equations, how to obtain the soliton solutions of a nonlinear wave equation is also of great research value. Recently, the Hirota bilinear method has been used to solve fractional models. This method is one direct and efficient way for obtaining the multi-solitary wave solutions. In this article, we obtain soliton solutions for the new model using the Hirota bilinear method.

This article is organized as follows. Section “Derivation of the 2D BO equation” is devoted to the derivation of the two-dimensional (2D) Benjamin–Ono (BO) equation which can describe the internal wave propagation in a deep stratified fluid. In section “Formulation of the time-fractional BO equation,” the resultant BO equation is converted into the time-fractional BO equation. In section “Conservation laws of the time-fractional BO equation,” conservation laws of the time-fractional BO equation are discussed. In section “Multi-soliton solutions for the time-fractional BO equation,” we seek multiple-soliton solutions of the equation using the simplified Hirota bilinear method. And then, based on the analytic solutions of the BO equation, the propagation of the two-solitary waves are discussed and illustrated. Section “Conclusion” contains the results and discussion of this work.

Derivation of the 2D BO equation

We shall begin by considering an inviscid, incompressible, density-stratified fluid with boundary conditions of infinite depth. Initially, we shall suppose that the flow is three dimensional (3D) and can be described by the spatial coordinates $(x,y,z)$ where x and y are horizontal and z is vertical. Also, we suppose that the fluid has a free surface at $z=0$ and the density $\rho (z)$ varies with depth in the upper layer but remains uniform in the lower layer. Here, in the basic state, the fluid has density ${\rho }_0(z)$ and a corresponding pressure $p_0(z)$ such that $d{p}_0/\mathrm{dz}=-g{\rho }_0$ describes the basic hydrostatic equilibrium. Then, in standard notation, the equations of motion relative to this basic state are

$$\{\begin{array}{c}\frac{\partial \rho }{\partial t}+u\frac{\partial \rho }{\partial x}+v\frac{\partial \rho }{\partial y}+w\frac{\partial \rho }{\partial z}=0\hfill \\ \multicolumn{1}{c}{\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0}\\ \multicolumn{1}{c}{\rho (\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z})=-\frac{\partial p}{\partial x}}\\ \multicolumn{1}{c}{\rho (\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z})=-\frac{\partial p}{\partial y}}\\ \multicolumn{1}{c}{\rho (\frac{\partial w}{\partial t}+u\frac{\partial w}{\partial x}+v\frac{\partial w}{\partial y}+w\frac{\partial w}{\partial z})=-\frac{\partial p}{\partial z}-g\rho }\end{array}$$

(1)

where t is the time variable, g is the gravitational acceleration, u, v, and w are the three fluid velocity components in the x-, y-, and z-directions, respectively, $\rho$ is the density, and p is the fluid pressure.

The boundary conditions are

$$\{\begin{array}{c}\frac{\partial w}{\partial z}-\frac{\mathrm{gw}}{c^2}=0,(z=0)\hfill \\ \multicolumn{1}{c}{w\to 0,(z\to \infty )}\end{array}$$

(2)

and the density is defined as

$$\rho =\{\begin{array}{c}\rho (z),(-2h_0\le z\le 0)\hfill \\ \multicolumn{1}{c}{{\rho }_{\infty }(\mathrm{constant}),(z\le -2h_0)}\end{array}$$

(3)

Now we consider the flow in the upper layer $-2h_0\le z\le 0$ and introduce the following independent stretched variables

$$\tau ={ϵ}^{2}t,\xi =ϵk(x-\mathrm{ct}),\eta ={ϵ}^{\frac{3}{2}}y,\zeta =z$$

(4)

where $ϵ$ is a small parameter characterizing the strength of nonlinearity, k is the wavenumber factor, and c is the wave velocity. Thus, we can obtain

$$\begin{array}{c}\multicolumn{1}{c}{\frac{\partial }{\partial t}=-ϵ\mathrm{ck}\frac{\partial }{\partial \xi }+{ϵ}^2\frac{\partial }{\partial \tau },\frac{\partial }{\partial x}=ϵk\frac{\partial }{\partial \xi },\frac{\partial }{\partial y}={ϵ}^{\frac{3}{2}}\frac{\partial }{\partial \eta },}\\ \multicolumn{1}{c}{\frac{\partial }{\partial z}=\frac{\partial }{\partial \zeta }}\end{array}$$

(5)

The variables are expanded in the following form

$$\{\begin{array}{c}c=c_0+c_1ϵ+c_2{ϵ}^2+\cdots \hfill \\ \multicolumn{1}{c}{k=1+ϵk_1+{ϵ}^2{k}_2+\cdots }\\ \multicolumn{1}{c}{u(\xi ,\eta ,\zeta ,\tau )=ϵu_1(\xi ,\eta ,\zeta ,\tau )+{ϵ}^2{u}_2(\xi ,\eta ,\zeta ,\tau )+\cdots }\\ \multicolumn{1}{c}{v(\xi ,\eta ,\zeta ,\tau )={ϵ}^{\frac{3}{2}}{v}_1(\xi ,\eta ,\zeta ,\tau )+{ϵ}^{\frac{5}{2}}{v}_2(\xi ,\eta ,\zeta ,\tau )+\cdots }\\ \multicolumn{1}{c}{w(\xi ,\eta ,\zeta ,\tau )={ϵ}^2{w}_1(\xi ,\eta ,\zeta ,\tau )+{ϵ}^3{w}_2(\xi ,\eta ,\zeta ,\tau )+\cdots }\\ \multicolumn{1}{c}{p(\xi ,\eta ,\zeta ,\tau )=p_0(\zeta )+ϵp_1(\xi ,\eta ,\zeta ,\tau )+{ϵ}^2{p}_2(\xi ,\eta ,\zeta ,\tau )+\cdots }\\ \multicolumn{1}{c}{\rho (\xi ,\eta ,\zeta ,\tau )={\rho }_0(\zeta )+ϵ{\rho }_1(\xi ,\eta ,\zeta ,\tau )+{ϵ}^2{\rho }_2(\xi ,\eta ,\zeta ,\tau )+\cdots }\end{array}$$

(6)

where $c_0$ is the linear long wave speed.

Substituting equations (5) and (6) into equation (1), we can obtain the approximate equations for $\epsilon$ in the following form

$${ϵ}^0:\frac{d{p}_0}{d\zeta }=-{\rho }_{0}g$$

(7)

$${ϵ}^1:\{\begin{array}{c}{c}_0\frac{\partial {\rho }_1}{\partial \xi }=w_1\frac{d{\rho }_0}{d\zeta }\hfill \\ \multicolumn{1}{c}{\frac{\partial u_1}{\partial \xi }+\frac{\partial w_1}{\partial \zeta }=0}\\ \multicolumn{1}{c}{-{\rho }_0{c}_0\frac{\partial u_1}{\partial \xi }+\frac{\partial p_1}{\partial \xi }=0}\\ \multicolumn{1}{c}{-{\rho }_0{c}_0\frac{\partial v_1}{\partial \xi }+\frac{\partial p_1}{\partial \eta }=0}\\ \multicolumn{1}{c}{\frac{\partial p_1}{\partial \zeta }=-{\rho }_{1}g}\end{array}$$

(8)

$${ϵ}^2:\{\begin{array}{c}-c_0\frac{\partial {\rho }_2}{\partial \xi }+w_2\frac{d{\rho }_0}{d\zeta }=Q_1\hfill \\ \multicolumn{1}{c}{\frac{\partial u_2}{\partial \xi }+\frac{\partial w_2}{\partial \zeta }=Q_2}\\ \multicolumn{1}{c}{-{\rho }_0{c}_0\frac{\partial u_2}{\partial \xi }+\frac{\partial p_2}{\partial \xi }=Q_3}\\ \multicolumn{1}{c}{-{\rho }_0{c}_0\frac{\partial v_2}{\partial \xi }+\frac{\partial p_2}{\partial \eta }=Q_4}\\ \multicolumn{1}{c}{\frac{\partial p_2}{\partial \zeta }=-{\rho }_{2}g}\end{array}$$

(9)

where

$$\{\begin{array}{c}{Q}_1=(c_1+c_0{k}_1-u_1)\frac{\partial {\rho }_1}{\partial \xi }-w_1\frac{\partial {\rho }_1}{\partial \zeta }-\frac{\partial {\rho }_1}{\partial \tau }\hfill \\ \multicolumn{1}{c}{Q_2=-k_1\frac{\partial u_1}{\partial \xi }-\frac{\partial v_1}{\partial \eta }}\\ \multicolumn{1}{c}{Q_3={\rho }_1{c}_0\frac{\partial u_1}{\partial \xi }+{\rho }_0(c_1+c_0{k}_1-u_1)\frac{\partial u_1}{\partial \xi }-{\rho }_0{w}_1\frac{\partial u_1}{\partial \zeta }-k_1\frac{\partial p_1}{\partial \xi }-{\rho }_0\frac{\partial u_1}{\partial \tau }}\\ \multicolumn{1}{c}{Q_4={\rho }_1{c}_0\frac{\partial v_1}{\partial \xi }+{\rho }_0(c_1+c_0{k}_1-v_1)\frac{\partial v_1}{\partial \xi }-{\rho }_0{w}_1\frac{\partial v_1}{\partial \zeta }-{\rho }_0\frac{\partial v_1}{\partial \tau }}\end{array}$$

(10)

Eliminating ${\rho }_1,p_1,u_1,\mathrm{and}{v}_1$ from equation (8), we can obtain the governing equation and the boundary condition of $w_1(\xi ,\eta ,\zeta ,\tau )$

$$\{\begin{array}{c}\frac{\partial }{\partial \zeta }\left({\rho }_0\frac{\partial w_1}{\partial \zeta }\right)-\frac{g{w}_1}{{(c_0)}^2}\frac{d{\rho }_0}{d\zeta }=0\hfill \\ \multicolumn{1}{c}{\frac{\partial w_1}{\partial \zeta }-\frac{g{w}_1}{{(c_0)}^2}=0,\zeta =0}\end{array}$$

(11)

To solve the above equations, we suppose that the solution to equation (11) is written as follows

$$w_1(\xi ,\eta ,\zeta ,\tau )=-\phi (\zeta )\frac{\partial f(\xi ,\eta ,\zeta ,\tau )}{\partial \xi }$$

(12)

where f denotes the wave amplitude which acquires

$$f\to 0(\xi \to \infty ,\eta \to \infty )$$

(13)

Substituting equations (12) and (13) into equation (11), we can obtain the governing equation and the boundary condition of $\phi (\zeta )$

$$\{\begin{array}{c}\frac{d}{d\zeta }({\rho }_0\frac{d\phi }{d\zeta })-\frac{g\phi }{{(c_0)}^2}\frac{d{\rho }_0}{d\zeta }=0\hfill \\ \multicolumn{1}{c}{\frac{d\phi }{d\zeta }-\frac{g\phi }{{(c_0)}^2}=0,\zeta =0}\end{array}$$

(14)

Also, by equation (12), we can rewrite ${\rho }_1,p_1,u_1,\mathrm{and}{v}_1$ as

$$\{\begin{array}{c}{u}_1=f\frac{d\phi }{d\zeta }\hfill \\ \multicolumn{1}{c}{v_1=\frac{d\phi }{d\zeta }\underset{0}{\overset{\xi }{\int }}\frac{\partial f}{\partial \eta }d\xi }\\ \multicolumn{1}{c}{{\rho }_1=-\frac{f\phi }{c_0}\frac{d{\rho }_0}{d\zeta }}\\ \multicolumn{1}{c}{p_1=c_0{\rho }_{0}f\frac{d\phi }{d\zeta }}\end{array}$$

(15)

Eliminating ${\rho }_2,p_2,u_2,\mathrm{and}{v}_2$ from equation (9), we can obtain the governing equation and the boundary condition of $w_2(\xi ,\eta ,\zeta ,\tau )$

$$\{\begin{array}{c}\frac{\partial }{\partial \zeta }({\rho }_0\frac{\partial w_2}{\partial \zeta })-\frac{g{w}_2}{{(c_0)}^2}\frac{d{\rho }_0}{d\zeta }=J(f,{\rho }_0,\phi )\hfill \\ \multicolumn{1}{c}{\frac{\partial w_2}{\partial \zeta }-\frac{g{w}_2}{{(c_0)}^2}=-\frac{2g{c}_1{w}_1}{c_0},\zeta =0}\end{array}$$

(16)

where

$$\begin{array}{c}J(f,{\rho }_0,\phi )=\frac{1}{c_0}\frac{\partial Q_3}{\partial \zeta }-\frac{g{Q}_1}{{(c_0)}^2}+\frac{\partial ({\rho }_0{Q}_2)}{\partial \zeta }\hfill \end{array}$$

(17)

By equation (15), we can rewrite $J(f,{\rho }_0,\phi )$ as

$$\begin{array}{c}J(f,{\rho }_0,\phi )=b_1{f}_{\tau }+b_2{f}_{\xi }+b_{3}f{f}_{\xi }+b_4\underset{0}{\overset{\xi }{\int }}\frac{{\partial }^{2}f}{\partial {\eta }^2}d\xi \hfill \end{array}$$

(18)

where

$$\{\begin{array}{c}{b}_1=-\frac{2({\rho }_0\phi \prime )\prime }{c_0}\hfill \\ \multicolumn{1}{c}{b_2=\frac{2c_1({\rho }_0\phi \prime )\prime }{c_0}}\\ \multicolumn{1}{c}{b_3=-\frac{3{\rho \prime }_0{(\phi \prime )}^2+2{\rho }_0\phi \prime \phi ″-2{\rho \prime }_0\phi \phi ″-2{\rho }_0\phi \phi ″\prime }{c_0}}\\ \multicolumn{1}{c}{b_4=-({\rho }_0\phi \prime )\prime }\end{array}$$

(19)

By the solvable condition of $f,\phi$ in equations (14) and (16)

$$\begin{array}{c}[{\rho }_0(\phi {w\prime }_2-\phi \prime w_2)]{|}_0^{\zeta =-2h_0}=\underset{0}{\overset{-2h_0}{\int }}J(f,{\rho }_0,\phi )\phi (\zeta )d\zeta \hfill \end{array}$$

(20)

we obtain

$$\begin{array}{l}[-c_0{\rho }_0(\phi {w^{\prime }}_2-{\phi }^{\prime }{w}_2)]{|}_{\zeta }=-2h_0=a_1{f}_{\tau }+a_2{f}_{\xi }+a_{3}f{f}_{\xi }\\ +a_4\underset{0}{\overset{\xi }{\int }}\frac{{\partial }^{2}f}{\partial {\eta }^2}d\xi \end{array}$$

(21)

where

$$\{\begin{array}{c}{a}_1=\underset{0}{\overset{-2h_0}{\int }}2({\rho }_0\phi \prime )\prime \phi d\zeta \hfill \\ \multicolumn{1}{c}{a_2=\underset{0}{\overset{-2h_0}{\int }}-2c_1({\rho }_0\phi \prime )\prime \phi d\zeta +c_1\frac{2g{\rho }_0(0){\phi }^2(0)}{{(c_0)}^2}}\\ \multicolumn{1}{c}{a_3=\underset{0}{\overset{-2h_0}{\int }}[3{\rho \prime }_0{(\phi \prime )}^2+2{\rho }_0\phi \prime \phi ″-2{\rho \prime }_0\phi \phi ″-2{\rho }_0\phi \phi ″\prime ]\phi d\zeta }\\ \multicolumn{1}{c}{a_4=\underset{0}{\overset{-2h_0}{\int }}{c}_0({\rho }_0\phi \prime )\prime \phi d\zeta }\end{array}$$

(22)

Next, we consider the flow in the lower layer $z\le -2h_0$ and introduce the following independent stretched variables

$$\tau ={ϵ}^{2}t,\xi =ϵk(x-\mathrm{ct}),\eta ={ϵ}^{2}y,\zeta =ϵk(z+2h_0)$$

(23)

and thus we can obtain

$$\frac{\partial }{\partial t}=-ϵ\mathrm{ck}\frac{\partial }{\partial \xi }+{ϵ}^2\frac{\partial }{\partial \tau },\frac{\partial }{\partial x}=ϵk\frac{\partial }{\partial \xi },\frac{\partial }{\partial y}={ϵ}^2\frac{\partial }{\partial \eta },\frac{\partial }{\partial z}=ϵk\frac{\partial }{\partial \zeta }$$

(24)

The variables are expanded in the following form

$$\{\begin{array}{c}c=c_0+c_1ϵ+c_2{ϵ}^2+\cdots \hfill \\ \multicolumn{1}{c}{k=1+ϵk_1+{ϵ}^2{k}_2+\cdots }\\ \multicolumn{1}{c}{u(\xi ,\eta ,\zeta ,\tau )={ϵ}^2{U}_1(\xi ,\eta ,\zeta ,\tau )+{ϵ}^3{U}_2(\xi ,\eta ,\zeta ,\tau )+\cdots }\\ \multicolumn{1}{c}{v(\xi ,\eta ,\zeta ,\tau )={ϵ}^2{V}_1(\xi ,\eta ,\zeta ,\tau )+{ϵ}^3{V}_2(\xi ,\eta ,\zeta ,\tau )+\cdots }\\ \multicolumn{1}{c}{w(\xi ,\eta ,\zeta ,\tau )={ϵ}^2{W}_1(\xi ,\eta ,\zeta ,\tau )+{ϵ}^3{W}_2(\xi ,\eta ,\zeta ,\tau )+\cdots }\\ \multicolumn{1}{c}{p(\xi ,\eta ,\zeta ,\tau )=p_0(\zeta )+{ϵ}^2{P}_1(\xi ,\eta ,\zeta ,\tau )+{ϵ}^2{P}_2(\xi ,\eta ,\zeta ,\tau )+\cdots }\end{array}$$

(25)

Substituting equations (24) and (25) into equation (1), we can obtain the following equation

$$\frac{{\partial }^2{W}_1}{\partial {\xi }^2}+\frac{{\partial }^2{W}_1}{\partial {\zeta }^2}=0$$

(26)

and the boundary condition of $W_1$

$$\{\begin{array}{c}{W}_1(\xi ,\eta ,\zeta ,\tau ){|}_{\zeta \to 0}\to -\frac{\partial f}{\partial \xi }\phi (z){|}_{z\to -2h_0}\hfill \\ \multicolumn{1}{c}{W_1(\xi ,\eta ,\zeta ,\tau ){|}_{\zeta \to \infty }\to 0}\end{array}$$

(27)

The solution to equation (26) with boundary conditions is the principal value of Cauchy integral

$$W_1(\xi ,\eta ,\zeta ,\tau )=\frac{1}{\pi }P.V\underset{-\infty }{\overset{+\infty }{\int }}-\frac{\partial f(\xi \prime ,\eta ,\tau )}{\partial \xi \prime }\frac{\zeta }{{\zeta }^2+{(\xi -\xi \prime )}^2}d\xi \prime$$

(28)

Taking the derivative with respect to $\zeta$ on both sides with the above equation and letting $\zeta \to 0$ , we obtain

$$\frac{\partial W_1(\xi ,\eta ,0,\tau )}{\partial \zeta }=\frac{1}{\pi }\frac{{\partial }^2}{\partial {\xi }^2}\underset{-\infty }{\overset{+\infty }{\int }}-\frac{f(\xi \prime ,\eta ,\tau )}{\xi -\xi \prime }d\xi \prime$$

(29)

Matching the value of the lower layer and the upper layer at $z=-2h_0$ , we can obtain

$$\{\begin{array}{c}{W}_1(\xi ,\eta ,\zeta ,\tau ){|}_{\zeta \to 0}=-\frac{\partial f}{\partial \xi }\phi (z){|}_{z\to -2h_0}\hfill \\ \multicolumn{1}{c}{\phi \prime (z){|}_{z\to -2h_0}=0}\\ \multicolumn{1}{c}{\frac{\partial w_2}{\partial z}{|}_{z\to -2h_0}=\frac{\partial W_1}{\partial \zeta }{|}_{\zeta \to 0}=\frac{{\partial }^2}{\partial {\xi }^2}\aleph (f){|}_{z\to -2h_0}}\\ \multicolumn{1}{c}{\phi ″(z){|}_{z\to -2h_0}=0}\end{array}$$

(30)

where

$$\aleph (f)=\frac{1}{\pi }\underset{-\infty }{\overset{+\infty }{\int }}\frac{f(\xi \prime ,\eta ,\tau )}{\xi -\xi \prime }d\xi \prime$$

(31)

is the famous Hilbert transformation. Then setting $\phi (-2h_0)=1$ and substituting equation (30) into equation (21), we can obtain

$$\begin{array}{c}{a}_1{f}_{\tau }+a_2{f}_{\xi }+a_{3}f{f}_{\xi }+a_4\underset{0}{\overset{\xi }{\int }}\frac{{\partial }^{2}f}{\partial {\eta }^2}d\xi +c_0{\rho }_{\infty }\frac{{\partial }^2}{\partial {\xi }^2}\aleph (f)=0\hfill \end{array}$$

(32)

Now setting $\phi (0)=0$ and $A_1=a_2/a_1,A_2=a_3/a_1,$ $A_3=c_0{\rho }_{\infty }/a_1,A_4=a_4/a_1$ , we obtain

$$\begin{array}{c}{f}_{\tau }+A_1{f}_{\xi }+A_{2}f{f}_{\xi }+A_3\frac{{\partial }^2}{\partial {\xi }^2}\aleph (f)+A_4\underset{0}{\overset{\xi }{\int }}\frac{{\partial }^{2}f}{\partial {\eta }^2}d\xi =0\hfill \end{array}$$

(33)

which is the 2D BO equation.

Formulation of the time-fractional BO equation

In section “Derivation of the 2D BO equation,” we have obtained the 2D BO equation. In this section, the BO equation will be converted to the time-fractional BO equation.

According to the BO equation (33), assuming $f(\xi ,\eta ,\tau )=g_{\xi }(\xi ,\eta ,\tau )$ where $g(\xi ,\eta ,\tau )$ is a potential function, the potential equation of equation (33) can be written in the following form

$$g_{\xi \tau }+A_1{g}_{\xi \xi }+A_2{g}_{\xi }{g}_{\xi \xi }+A_3[\aleph (f){]}_{\xi \xi }+A_4\underset{0}{\overset{\xi }{\int }}\frac{{\partial }^{3}g}{\partial {\eta }^2\partial \xi }d\xi =0$$

(34)

Then, the functional of the potential equation (34) can be described as

$$\begin{array}{c}J(g)=\underset{R}{\int }d\xi \underset{R}{\int }d\eta \underset{T}{\int }d\tau g(\xi ,\eta ,\tau )[B_1{g}_{\xi \tau }+B_2{A}_1{g}_{\xi \xi }\\ +B_3{A}_2{g}_{\xi }{g}_{\xi \xi }+B_4{A}_3{[\aleph (f)]}_{\xi \xi }+B_5{A}_4\underset{0}{\overset{\xi }{\int }}\frac{{\partial }^{3}g}{\partial {\eta }^2\partial \xi }d\xi ]\end{array}$$

(35)

where $B_i$ , $i=1,2,3,4,5$ , are the Lagrangian multipliers which can be obtained later. Using integration by parts for equation (35) and taking $g_{\xi }{|}_R=g_{\eta }{|}_R=g_{\tau }{|}_T=0$ , we obtain

$$\begin{array}{c}\multicolumn{1}{c}{J(g)=\underset{R}{\int }d\xi \underset{R}{\int }d\eta \underset{T}{\int }d\tau [-B_1{g}_{\xi }{g}_{\tau }-B_2{A}_1{g}_{\xi }^2-\frac{1}{2}{B}_3{A}_2{g}_{\xi }^3}\\ \multicolumn{1}{c}{+B_4{A}_{3}g{[\aleph (f)]}_{\xi \xi }-B_5{A}_4{g}_{\eta }\underset{0}{\overset{\xi }{\int }}\frac{{\partial }^{2}g}{\partial \eta \partial \xi }d\xi ]}\end{array}$$

(36)

Using the variation of the above function, we obtain

$$\begin{array}{c}\multicolumn{1}{c}{2B_1{g}_{\xi \tau }+2B_2{A}_1{g}_{\xi \xi }+3B_3{A}_2{g}_{\xi }{g}_{\xi \xi }+B_4{A}_3[\aleph (f){]}_{\xi \xi }}\\ \multicolumn{1}{c}{+B_5{A}_4\underset{0}{\overset{\xi }{\int }}\frac{{\partial }^{3}g}{\partial {\eta }^2\partial \xi }d\xi =0}\end{array}$$

(37)

Comparing equation (34) with equation (37), we obtain the following Lagrangian multipliers

$$B_1=\frac{1}{2},B_2=\frac{1}{2},B_3=\frac{1}{3},B_4=1,B_5=1$$

(38)

Therefore, the Lagrangian form of the BO equation is given by

$$\begin{array}{c}\multicolumn{1}{c}{L(g,g_{\xi },g_{\tau },g_{\eta })=-\frac{1}{2}{g}_{\xi }{g}_{\tau }-\frac{1}{2}{A}_1{g}_{\xi }^2-\frac{1}{6}{A}_2{g}_{\xi }^3}\\ \multicolumn{1}{c}{+A_{3}g[\aleph (f){]}_{\xi \xi }-A_4{g}_{\eta }\underset{0}{\overset{\xi }{\int }}\frac{{\partial }^{2}g}{\partial \eta \partial \xi }d\xi }\end{array}$$

(39)

Similarly, the Lagrangian form of the time-fractional BO equation is given by

$$\begin{array}{c}F(g,g_{\xi },D_{\tau }^{\alpha }g,g_{\eta })=-\frac{1}{2}{g}_{\xi }{D}_{\tau }^{\alpha }g-\frac{1}{2}{A}_1{g}_{\xi }^2-\frac{1}{6}{A}_2{g}_{\xi }^3+A_{3}g{[\aleph (f)]}_{\xi \xi }\hfill \\ \multicolumn{1}{c}{-A_4{g}_{\eta }\underset{0}{\overset{\xi }{\int }}\frac{{\partial }^{2}g}{\partial \eta \partial \xi }d\xi }\end{array}$$

(40)

where the fractional derivative $D_{\tau }^{\alpha }$ indicates the modified Riemann–Liouville (mRL) fractional derivative

$$D_z^{\gamma }h(z)=\frac{1}{\Gamma (1-\gamma )}\frac{\partial }{\partial z}\underset{a}{\overset{z}{\int }}\frac{[h(\zeta )-h(a)]}{{(z-\zeta )}^{\gamma }}d\zeta ,0\le \gamma <1$$

(41)

Thus, the functional of the time-fractional BO equation can be obtained as

$$\begin{array}{c}{J}_F(g)=\underset{R}{\int }d\xi \underset{R}{\int }d\eta \underset{R}{\int }d\tau F(g,g_{\xi },D_{\tau }^{\alpha }g,g_{\eta })\hfill \end{array}$$

(42)

According to Agrawal’s method, the variation of the functional (equation (42)) can be written as

$$\begin{array}{l}\delta J_F(g)=\underset{R}{\int }d\xi \underset{R}{\int }d\eta \underset{T}{\int }d\tau \\ \left[\frac{\partial F}{\partial D_{\tau }^{\alpha }g}\delta (D_{\tau }^{\alpha }g)+\frac{\partial F}{\partial g}\delta g+\frac{\partial F}{\partial g_{\xi }}\delta g_{\xi }+\frac{\partial F}{\partial g_{\eta }}\delta g_{\eta }\right]\end{array}$$

(43)

The fractional integration by parts satisfies the following rules

$$\underset{a}{\overset{b}{\int }}{(d\tau )}^{j}h(\tau )=j\underset{a}{\overset{b}{\int }}d\tau {(b-\tau )}^{j}h(\tau )$$

(44)

and

$$\begin{array}{c}\multicolumn{1}{c}{\underset{a}{\overset{b}{\int }}{(d\tau )}^{j}h(z)D_z^{j}g(z)=\Gamma (1+j)}\\ \multicolumn{1}{c}{[g(z)h(z){|}_a^b-\underset{a}{\overset{b}{\int }}{(\mathrm{dz})}^{j}g(z)D_z^{j}h(z)]}\end{array}$$

(45)

Using the fractional integration by parts for equation (43), we can obtain

$$\begin{array}{l}\delta J_F(g)=\underset{R}{\int }d\xi \underset{R}{\int }d\eta \underset{T}{\int }d\tau \\ \left[-D_{\tau }^{\alpha }\left(\frac{\partial F}{\partial D_{\tau }^{\alpha }}\right)+\frac{\partial F}{\partial g}-\frac{\partial }{\partial \xi }\left(\frac{\partial F}{\partial g_{\xi }}\right)-\frac{\partial }{\partial \eta }\left(\frac{\partial F}{\partial g_{\eta }}\right)\right]\delta g\end{array}$$

(46)

Optimizing the variation (equation (46)) and setting $\delta J_F(g)=0$ , we can obtain

$$\begin{array}{c}-D_{\tau }^{\alpha }\left(\frac{\partial F}{\partial D_{\tau }^{\alpha }}\right)+\frac{\partial F}{\partial g}-\frac{\partial }{\partial \xi }\left(\frac{\partial F}{\partial g_{\xi }}\right)-\frac{\partial }{\partial \eta }\left(\frac{\partial F}{\partial g_{\eta }}\right)=0\hfill \end{array}$$

(47)

Substituting equation (40) into equation (47), we have

$$D_{\tau }^{\alpha }{g}_{\xi }+A_1{g}_{\xi \xi }+A_2{g}_{\xi }{g}_{\xi \xi }+A_3[\aleph (f){]}_{\xi \xi }+A_4\underset{0}{\overset{\xi }{\int }}\frac{{\partial }^{3}g}{\partial {\eta }^2\partial \xi }d\xi =0$$

(48)

Substituting $g_{\xi }(\xi ,\eta ,\tau )=f(\xi ,\eta ,\tau )$ into equation (48), we can obtain

$$\begin{array}{c}{D}_{\tau }^{\alpha }f+A_1{f}_{\xi }+A_{2}f{f}_{\xi }+A_3\frac{{\partial }^2}{\partial {\xi }^2}\aleph (f)+A_4\underset{0}{\overset{\xi }{\int }}\frac{{\partial }^{2}f}{\partial {\eta }^2}d\xi =0\hfill \end{array}$$

(49)

which is the time-fractional BO equation.

Conservation laws of the time-fractional BO equation

In this section, we will discuss conservation laws of the time-fractional BO equation to learn about the properties of the new model (49). First, we assume

$$\begin{array}{c}f,f_{\xi },f_{\xi \xi }\to 0,|\xi |\to \infty \hfill \\ \multicolumn{1}{c}{f,f_{\eta }\to 0,|\eta |\to \infty }\end{array}$$

(50)

Integrating equation (49) with respect to $\xi$ from $-\infty$ to $+\infty$ yields

$$E_1=\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}{D}_{\tau }^{\alpha -1}\mathrm{fd}\xi d\eta ,\frac{\partial E_1}{\partial \tau }=0$$

(51)

where $E_1$ means the mass of the algebraic internal solitary waves and equation (51) indicates that the mass of the algebraic internal solitary waves is conserved.

Next, we multiply equation (49) by $f(\xi ,\eta ,\zeta )$ and integrate it with respect to $\xi$ from $-\infty$ to $+\infty$ . Using the property of the mRL fractional derivative

$$D_z^{\alpha }[f(x)g(x)]=f(x)D_z^{\alpha }g(x)+g(x)D_z^{\alpha }f(x)$$

(52)

and the properties of the Hilbert operator

$$\{\begin{array}{c}\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}f(X,Y)\aleph [f(X,Y)]dXdY=0\\ \underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}f(X,Y)\aleph [g(X,Y)]dXdY\\ =\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}g(X,Y)\aleph [f(X,Y)]dXdY\end{array}$$

(53)

where $f(X,Y)\to 0,g(X,Y)\to 0,(|X|\to \infty ,|Y|\to \infty )$ , we have

$$E_2=\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}\frac{1}{2}{D}_{\tau }^{\alpha -1}{f}^{2}d\xi d\eta ,\frac{\partial E_2}{\partial \tau }=0$$

(54)

where $E_2$ means the momentum of the internal gravity solitary waves and equation (54) indicates that the momentum of the algebraic internal solitary waves is conserved.

Now, multiplying equation (49) by $f_2-(A_4/A_2)\aleph (f_{\xi })$ , we have

$$\begin{array}{l}\left[f_2-\frac{A_4}{A_2}\aleph (f_{\xi })\right]\\ \left[D_{\tau }^{\alpha }f+A_1{f}_{\xi }+A_{2}f{f}_{\xi }+A_3\frac{{\partial }^2}{\partial {\xi }^2}\aleph (f)+A_4\underset{0}{\overset{\xi }{\int }}\frac{{\partial }^{2}f}{\partial {\eta }^2}d\xi \right]=0\end{array}$$

(55)

Deriving equation (49) about $\xi$ and multiplying $f_{\xi }+(A_4/A_2)\aleph (f)$ , we obtain

$$\begin{array}{l}\left[f_{\xi }+\frac{A_4}{A_2}\aleph (f)\right]\frac{\partial }{\partial \xi }\\ \left[D_{\tau }^{\alpha }f+A_1{f}_{\xi }+A_{2}f{f}_{\xi }+A_3\frac{{\partial }^2}{\partial {\xi }^2}\aleph (f)+A_4\underset{0}{\overset{\xi }{\int }}\frac{{\partial }^{2}f}{\partial {\eta }^2}d\xi \right]=0\end{array}$$

(56)

Deriving equation (49) about $\eta$ and multiplying $f_{\eta }$ , we obtain

$$\begin{array}{c}{f}_{\eta }\frac{\partial }{\partial \eta }[D_{\tau }^{\alpha }f+A_1{f}_{\xi }+A_{2}f{f}_{\xi }+A_3\frac{{\partial }^2}{\partial {\xi }^2}\aleph (f)+A_4\underset{0}{\overset{\xi }{\int }}\frac{{\partial }^{2}f}{\partial {\eta }^2}d\xi ]=0\hfill \end{array}$$

(57)

According to the above equations, we obain

$$\begin{array}{c}{E}_3=\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}[\frac{1}{3}{D}_{\tau }^{\alpha -1}(f^3)+\frac{1}{2}{D}_{\tau }^{\alpha -1}(f_{\xi }^2)\hfill \\ +\frac{1}{2}{D}_{\tau }^{\alpha -1}(f_{\eta }^2)+\frac{A_4}{A_2}{D}_{\tau }^{\alpha -1}[f_{\xi }\aleph (f)]]d\xi d\eta ,\frac{\partial E_3}{\partial \tau }=0\hfill \end{array}$$

(58)

where $E_3$ means the energy of the algebraic internal solitary waves and equation (58) indicates that the energy of the algebraic internal solitary waves is conserved.

Multi-soliton solutions for the time-fractional BO equation

In this section, using the simplified Hirota bilinear method, we seek multiple-soliton solutions of equation (49) and the propagation of the two-solitary waves based on the obtained analytic results.

First, we introduce the following fractional transforms

$$T=\frac{p{\tau }^{\alpha }}{\Gamma (1+\alpha )},D_{\tau }^{\alpha }f=p{f}_{\tau }$$

(59)

where p is a constant. Using the above transformations, we can convert equation (49) into

$$\begin{array}{c}{f}_T(\xi ,\eta ,T)+A_1{f}_{\xi }(\xi ,\eta ,T)+A_{2}f(\xi ,\eta ,T)f_{\xi }(\xi ,\eta ,T)\hfill \\ +A_3\frac{{\partial }^2}{\partial {\xi }^2}\aleph [f(\xi ,\eta ,T)]+A_4\underset{0}{\overset{\xi }{\int }}\frac{{\partial }^{2}f(\xi ,\eta ,T)}{\partial {\eta }^2}d\xi =0\hfill \end{array}$$

(60)

We assume that the solution of equation (60) has the form

$$f(\xi ,\eta ,T)=\mathrm{iR}{\left[\ln \frac{s^{*}(\xi ,\eta ,T)}{s(\xi ,\eta ,T)}\right]}_{\xi }$$

(61)

where

$$\{\begin{array}{c}s(\xi ,\eta ,T)\propto \underset{n=1}{\overset{N}{\Pi }}[\xi -{\xi }_n(T)]\hfill \\ \multicolumn{1}{c}{\mathrm{Im}{\xi }_n>0,n=1,2,\dots ,N}\end{array}$$

(62)

where $s^{*}$ is the complex conjugate of s and ${\xi }_n(n=1,2,\dots ,N)$ are all complex functions for T.

Substituting equation (61) into equation (60), we can obtain

$$R=\frac{2A_3}{A_2}$$

(63)

Hence, equation (61) can be written as

$$f(\xi ,\eta ,T)=\frac{i2A_3}{A_2}{\left[\ln \frac{s^{*}(\xi ,\eta ,T)}{s(\xi ,\eta ,T)}\right]}_{\xi }$$

(64)

Then using formula

$$P\frac{1}{\eta -\xi }=\frac{1}{2}\underset{\epsilon \to +0}{lim}(\frac{1}{\eta -\xi +i\epsilon }+\frac{1}{\eta -\xi -i\epsilon })$$

(65)

we have

$$\begin{array}{c}\aleph [f(\xi ,\eta ,T)]=-\frac{2A_3}{A_2}(\sum _{n=1}^N\frac{1}{\xi -{\xi }_n(T)}+\sum _{n=1}^N\frac{1}{{\xi }^{*}-{\xi }_n^{*}(T)})\hfill \\ \multicolumn{1}{c}{=-\frac{2A_3}{A_2}(\frac{S_{\xi }(\xi ,\eta ,T)}{S(\xi ,\eta ,T)}+\frac{S_{\xi }^{*}(\xi ,\eta ,T)}{S^{*}(\xi ,\eta ,T)})}\end{array}$$

(66)

Substituting equations (62), (64), and (66) into equation (60), the internal wave BO equation can be transformed to the bilinear form

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{Im}(S_T^{*}S)+A_1{I}_m(S_{\xi }^{*}S)=A_2{S}_{\xi }^{*}{S}_{\xi }-A_3\mathrm{Re}(S_{\xi \xi }^{*}S)}\\ \multicolumn{1}{c}{-A_4\mathrm{Re}\underset{0}{\overset{\xi }{\int }}{S}_{\eta \eta }^{*}\mathrm{Sd}\xi }\end{array}$$

(67)

Based on the bilinear equation (67) and using the bilinear method, the N-solitary wave solution is given in the matrix form $S_N=\det (H_{N\times N})$ where H is an $N\times N$ matrix, the entries of which are

$$H_{\mathrm{nm}}=\{\begin{array}{c}i{\theta }_n+1,n=m\hfill \\ \multicolumn{1}{c}{\frac{2{(e_n{e}_m)}^{\frac{1}{2}}}{e_n-e_m},n\ne m}\end{array}$$

(68)

with

$${\theta }_n(\xi ,\eta ,T)=e_N(\xi -\eta -e_{n}T-{\omega }_n)$$

(69)

where $e_n,{\omega }_n(n=1,2,\dots ,N)$ are arbitrary constants and $e_n\ne e_m$ when $n\ne m$ .

When $N=1$ , we have

$$S_1=i{\theta }_1+1=i{e}_1(\xi -\eta -e_{1}T-{\omega }_1)$$

(70)

Substituting equation (70) into equation (64), we obtain the one-solitary wave solution

$$f_1=\frac{4A_3{e}_1}{A_2({\theta }_1^2+1)}$$

(71)

When $N=2$ , we have

$$S_2=-{\theta }_1{\theta }_2+i({\theta }_1+{\theta }_2)+e_{12}$$

(72)

where

$$e_{\mathrm{nm}}={\left(\frac{e_n+e_m}{e_n-e_m}\right)}^2$$

(73)

Substituting equation (72) into equation (64), we obtain the two-solitary wave solution

$$f_2(\xi ,\eta ,T)=\frac{4A_3}{A_2}\frac{e_1{e}_2{({\theta }_1+{\theta }_2)}^2-2(e_1+e_2)(e_{12}-{\theta }_1{\theta }_2)}{{(e_{12}-{\theta }_1{\theta }_2)}^2+{({\theta }_1+{\theta }_2)}^2}$$

(74)

Now, we rewrite equation (74) as

$$f_2(\xi ,\eta ,\tau )=\frac{4A_3}{A_2}\frac{e_1{e}_2{({\theta }_1+{\theta }_2)}^2-2(e_1+e_2)(e_{12}-{\theta }_1{\theta }_2)}{{(e_{12}-{\theta }_1{\theta }_2)}^2+{({\theta }_1+{\theta }_2)}^2}$$

(75)

where

$${\theta }_n(\xi ,\eta ,\tau )=e_n(\xi -\eta -\frac{e_n{\tau }^{\alpha }}{\Gamma (1+\alpha )}-{\omega }_n)$$

(76)

Based on the above conclusion, we can give the images of the two-soliton solution with different fractional orders $\alpha$ .

As we can see from Figure 1, when $0<\alpha <1/2$ , with the increase of the independent variable, there will be multiple increasing wave peaks, and the propagation speed is slow. When $\alpha =1/2$ , there will be two obvious wave peaks and the propagation speed of soliton waves will be faster. When $1/2<\alpha <1$ , there is only a single wave peak and the propagation speed is the fastest. In conclusion, the change of $\alpha$ has a great impact on the waves. With the continuous increase of the fractional order $\alpha$ , the propagation form of soliton waves changes from multiple peaks to a single peak, and the propagation speed increases continuously.

Figure 1.

Plots for the two-soliton solution with different fractional orders $\alpha$ .

In addition to single- and double-soliton wave solutions, there may be three-soliton wave solutions to the equation,^52–54 but for the fractional order equation in this article, we only study the first two solutions, and whether and how the three-soliton solutions exist will be the focus of our future research.

Conclusion

In this work, using multi-scale analysis and the perturbation method, we have obtained the 2D BO equation based on the basic system of equations. Furthermore, based on the integer-order model and using the semi-inverse method and the fractional variational principle, we obtain the 2D fractional BO equation which opens the door to the study of internal solitary waves. Also, we study the conservation laws of the new model using the Riemann–Liouville fractional derivative. Then, using the Hirota bilinear method, we obtain the soliton solutions of the new fractional model. Based on the multi-soliton solutions, we study the characteristics of motion of internal solitary waves.