Analytical study of (3 + 1)-dimensional fractional ultralow-frequency dust acoustic waves in a dual-temperature plasma

Abstract

Ultralow-frequency dust acoustic waves in plasma are an important branch of study in plasma physics. Due to the important application of dust acoustic waves in the space environment, industrial processing, and the laboratory, they have aroused great interest in physics and other disciplines. In this paper, dust acoustic waves in dust plasma containing double-temperature ions are studied. According to the equations of dust plasma motion considering the higher order perturbation factors, the (3 + 1)-dimensional Kadomtsev–Petviashvili equation is derived by using multiscale analysis and reduced perturbation methods. The result is a three-dimensional model describing the propagation of waves in space. Compared with the low-dimensional model, this model is more closely aligned with the actual physical situation. Further, to better comprehend the propagation features of dust acoustic waves in plasmas, with the help of the semi-inverse method and the fractional variational principle, the (3 + 1)-dimensional time-space-fractional Kadomtsev–Petviashvili equation is deduced. At the same time, based on the Lie symmetry analysis method, the conservation laws of the time-fractional Kadomtsev–Petviashvili equation are discussed, and the conserved quantities of dust acoustic waves in dual-temperature plasma are obtained. Finally, by virtue of the fractional subequation method, three types of exact analytical solutions of fractional Kadomtsev–Petviashvili equations are given. Drawing on the chirp concept in optical soliton communication, the chirp effect of ultralow-frequency dust acoustic waves is discussed, and the effects of dust temperature, ion temperature, and fractional order on dust acoustic waves in a dual-temperature plasma are studied using the analytical solutions.

Introduction

In 1929, Langmuir and Tonks introduced the word “plasma” for the first time to indicate the fourth state of matter. Plasma is an ionized gas-like substance comprising atoms and atoms that are deprived by partial electrons and positive and negative electrons generated by ionization. Humans gained an understanding of plasma much later than they did for the other three states of matter. Since the 1950s, plasma physics has matured with the help of controlled thermal fusion research and space technology. In the late 1970s, plasma physics became a new independent branch of physics recognized by the physics community. In recent years, plasma physics has become an important area for human beings to understand and control changes in the global environment and to maintain global communications. For this reason, plasma physics has received increasingly more attention.

To date, various characteristics of ordinary plasma have been extensively studied. However, with the rapid development of science and technology, complex plasmas have emerged. A special aspect of dust plasma compared to ordinary plasma is the presence of dust particles. It is precisely because of this feature that dust plasma exhibits new characteristics that are quite different from ordinary plasma and presents many new physical phenomena. Dust plasma is widely used in the field of cosmology, industrial processing, laboratories, etc. Due to its importance, many scholars have been involved in dust plasma physics since the mid-1990s. In the following decade, dust plasma research developed rapidly, forming a new branch in the plasma discipline.

Collective fluctuations and instability are among the most studied topics in dust plasma physics. In addition to the fluctuation pattern in ordinary plasma, dust plasma also has ultralow-frequency dust acoustic waves oscillating on a time scale of the dust particle movement. This property is because the presence of dust particles can significantly alter the nature and behavior of the plasma into which they are immersed, affecting various plasma fluctuation patterns and producing new, very-low-frequency fluctuation patterns. The dust sound wave was first predicted by Rao et al. ¹ in 1990. Barkan first confirmed it experimentally in 1995. Since then, interest in the study of dust sound waves and instability has intensified.

In an earlier study, Ghosh et al. ² used a (1 + 1)-dimensional model to study the effect of nonadiabatic dust charge changes in nonlinear dust ion acoustic waves in collision-free dust plasmas. Kadomtsev and Petviashvili ³ attempted to describe solitons in (2 + 1) dimensional systems by applying the Kadomtsev–Petviashvili (KP) equation. Subsequently, Duan ⁴ used the reduced perturbation method to derive the (2 + 1)-dimensional KP equation of the nonmagnetized hot dust plasma, which indicates that the nonlinear dust acoustic waves in the hot dust plasma are stable. Gill et al. ⁵ considered the variable dust charge and two-temperature ions to derive the (2 + 1)-dimensional KP equation for the dust plasma in the same way. Pakzad ⁶ studied the propagation of nonlinear waves in a warm dust plasma with variable dust charge, two-temperature ions, and nonthermal electrons based on the KP equation. Saini et al. ⁷ analyzed and numerically discussed the characteristics of dust acoustic solitary waves in dust plasma systems composed of dust fluids and superheated electrons and ions. However, in reality, waves propagate in three-dimensional space. One- and two-dimensional models are limited in solving some practical problems involving dust acoustic waves, which require us to consider a higher dimensional model.

On the other hand, with the development of nonlinear science,^8,9 fractional calculus theory has received increasingly more attention. In particular, fractional differential equations abstracted from practical problems have become research hot spots.^10,11 However, in the study of some practical problems, most studies have established integer-order models^12–15; fractional-order models are few in number. Due to the presence of dispersion and/or dissipation forces, the physical process is nonconservative in the actual plasma system. The classical processing of forces is achieved with an integer-order difference equation, which means that these conservative descriptions are not convenient for dealing with nonconservative physical processes. Therefore, it is necessary to establish a fractional-order model to study dust acoustic waves.

The conservation law^16–18 plays a very important role in the study of nonlinear physical phenomena. The symmetry analysis of fractional differential equations^19,20 has aroused great interest among researchers. The Lie symmetry analysis method was first introduced by the Norwegian mathematician Lie ²¹ at the end of the 19th century. Noether’s theorem ²² famously establishes a relationship between the conservation laws of symmetric and differential equations. Previously, researchers studied the conservation laws of integer-order equations. Only a few scholars have discussed the conservation laws of fractional-order equations.

The solution of nonlinear partial differential equations^23–26 has become an indispensable method for mathematics and physics researchers to study nonlinear problems. As fractional differential equations are gradually used to solve problems in optical, thermal, and mechanical systems, the solution of fractional differential equations^27,28 has attracted increasingly more attention. The methods for solving fractional differential equations^29–31 have also gradually matured, such as the extended Kudryaskov method, ³² the $\exp (-\phi (\xi ))$ method, ³³ and the monotone iterative method. ³⁴ In this paper, the fractional subequation method ³⁵ is adopted to derive the exact solutions of the (3 + 1)-dimensional time-space-fractional Kadomtsev–Petviashvili (TSF-KP) equation. The acquisition of the exact analytical solution lays the theoretical foundation for further research on the properties of dust acoustic waves.

The rest of this article is organized as follows. In the next section, based on the dust particle motion equation, the (3 + 1)-dimensional integer-order KP equation is derived using multiscale analysis and perturbation expansion methods.^36,37 By virtue of the semi-inverse method and the fractional variation principle,^38,39 the integer-order KP equation is transformed into a TSF-KP equation. Then, based on the Lie symmetric analysis method, the conservation laws of the time-fractional KP equation are discussed, and the conserved quantities of the dust acoustic waves during propagation are obtained. Next, the three types of exact analytical solutions of the fractional KP equations are given by the fractional subequation method. Finally, with the help of the exact solutions of the (3 + 1)-dimensional TSF-KP equation, the chirp effect of dust acoustic waves is studied. The effects of dust temperature, ion temperature, ion number density, and fractional order on the characteristics of dust acoustic waves are discussed.

Derivation of the (3 + 1)-dimensional KP equation

To study ultralow-frequency dust sound waves in dual-temperature dust plasma, we assume that the dust plasma is nonmagnetized, that the sound wave propagates along the x direction, and that the high-order transverse perturbations in the y and z directions are weak. The dual-temperature dust plasma studied in this paper consists of three parts: negatively charged dust particles, low-temperature ions, and high-temperature ions. Dust plasma is electrically neutral under the equilibrium conditions of $n_{il0}+n_{ih0}=Z_d{n}_{d0}$ , where Z_d is the number of charges of the undisturbed dust particles, and $n_{il0},n_{ih0}$ , and $n_{d0}$ represent the number density of the low-temperature ions, high-temperature ions, and dust particles, respectively. The dynamic behavior of the dust particles is determined by the continuity equation and the momentum equation, and the dimensionless governing equation is as follows

$$\{\begin{array}{l}\frac{\partial n_d}{\partial t}+\frac{\partial n_d{u}_d}{\partial x}+\frac{\partial n_d{v}_d}{\partial y}+\frac{\partial n_d{w}_d}{\partial z}=0,\hfill \\ \frac{\partial u_d}{\partial t}+u_d\frac{\partial u_d}{\partial x}+v_d\frac{\partial u_d}{\partial y}+w_d\frac{\partial u_d}{\partial z}+\frac{\sigma }{n_d}\frac{\partial p_d}{\partial x}=\frac{\partial \varphi }{\partial x},\hfill \\ \frac{\partial v_d}{\partial t}+u_d\frac{\partial v_d}{\partial x}+v_d\frac{\partial v_d}{\partial y}+w_d\frac{\partial v_d}{\partial z}+\frac{\sigma }{n_d}\frac{\partial p_d}{\partial y}=\frac{\partial \varphi }{\partial y},\hfill \\ \frac{\partial w_d}{\partial t}+u_d\frac{\partial w_d}{\partial x}+v_{{\hspace{0.17em}}_d}\frac{\partial w_d}{\partial y}+w_d\frac{\partial w_d}{\partial z}+\frac{\sigma }{n_d}\frac{\partial p_d}{\partial z}=\frac{\partial \varphi }{\partial z},\hfill \\ \frac{\partial p_d}{\partial t}+u\frac{\partial p_d}{\partial x}+v\frac{\partial p_d}{\partial y}+w\frac{\partial p_d}{\partial z}+\gamma p_d\left(\frac{\partial u}{\partial t}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}\right)=0,\hfill \\ \frac{{\partial }^2\varphi }{\partial x^2}+\frac{{\partial }^2\varphi }{\partial y^2}+\frac{{\partial }^2\varphi }{\partial z^2}=n_d-n_{il}-n_{ih}\hfill \end{array}$$

(1)

where the particle number density n_d is normalized by $n_{d0}$ , the pressure p_d is normalized by $Z_d{n}_{d0}{T}_d$ , T_d is the dust temperature, the spatial coordinates (x, y, z) are normalized by the equivalent Debye length $\lambda ={\left[\frac{T_{eff}}{4\pi e^2\left(n_{il0}+n_{ih0}\right)}\right]}^{\frac{1}{2}}$ , the time t is normalized by the reciprocal of the equivalent dust plasma frequency ${\omega }_{pd}^{-1}={\left(\frac{m_d}{4\pi Z_d^2{n}_{d0}{e}^2}\right)}^{\frac{1}{2}}$ , the speeds $\left(u_d,v_d,w_d\right)$ are normalized by the equivalent speed of sound $C_d={\left(\frac{Z_d{T}_{eff}}{m_d}\right)}^{\frac{1}{2}}$ , and the potential $\varphi$ is normalized by $\frac{T_{eff}}{e}$ . $T_{eff}=\frac{T_l{T}_h}{\nu T_l+\mu T_h}$ is the equivalent temperature, ν and μ are the number density of the nondimensionalized high- and low-temperature ions, respectively, T_l and T_h are the temperature of the high- and low-temperature ions, and the ratio between the dust temperature and equivalent temperature is $\sigma =\frac{T_d}{T_{eff}}$ in the case of adiabatic γ = 3.

The number densities of the low- and high-temperature ions satisfy the Boltzmann distribution

$$n_{il}=\mu \exp \left(\frac{-\varphi }{\mu +\nu \beta }\right),n_{ih}=\nu \exp \left(\frac{-\beta \varphi }{\mu +\nu \beta }\right)$$

(2)

where $\beta =\frac{T_l}{T_h}$ and $\mu +\nu =1$ .

To obtain the KP equation, multiscale analysis and perturbation expansion are used. The independent variables are stretched as

$$\xi ={\mathrm{ϵ}}^{\frac{1}{2}}\left(x-v_0\right),\eta =\mathrm{ϵ}y,\zeta =\mathrm{ϵ}z,\tau ={\mathrm{ϵ}}^{\frac{3}{2}}t$$

(3)

where ϵ represents a small parameter of nonlinear strength, v₀ is the speed of the nonlinear dust sound waves, which is dimensionless by C_d. According to equation (3) we have

$$\begin{array}{l}\frac{\partial }{\partial t}={\mathrm{ϵ}}^{\frac{3}{2}}\frac{\partial }{\partial \tau }-v_0{\epsilon }^{\frac{1}{2}},\frac{\partial }{\partial x}={\mathrm{ϵ}}^{\frac{1}{2}}\frac{\partial }{\partial \xi },\frac{\partial }{\partial y}=\mathrm{ϵ}\frac{\partial }{\partial \eta },\frac{\partial }{\partial z}=\mathrm{ϵ}\frac{\partial }{\partial \zeta }\hfill \end{array}$$

(4)

The dependent variables $n_d,u_d,w_d$ , and $\varphi$ are expanded as follows

$$\{\begin{array}{l}{n}_d=1+\mathrm{ϵ}{n}_{d1}+{\mathrm{ϵ}}^2{n}_{d2}+{\mathrm{ϵ}}^3{n}_{d3}+\cdots ,\hfill \\ u_d=u_0+\mathrm{ϵ}{u}_{d1}+{\mathrm{ϵ}}^2{u}_{d2}+{\mathrm{ϵ}}^3{u}_{d3}+\cdots ,\hfill \\ v_d={\mathrm{ϵ}}^{\frac{3}{2}}{v}_{d1}+{\mathrm{ϵ}}^{\frac{5}{2}}{v}_{d2}+{\mathrm{ϵ}}^{\frac{7}{2}}{v}_{d3}+\cdots ,\hfill \\ w_d={\mathrm{ϵ}}^{\frac{3}{2}}{w}_{d1}+{\mathrm{ϵ}}^{\frac{5}{2}}{w}_{d2}+{\mathrm{ϵ}}^{\frac{7}{2}}{w}_{d3}+\cdots ,\hfill \\ p_d=1+\mathrm{ϵ}{p}_{d1}+{\mathrm{ϵ}}^2{p}_{d2}+{\mathrm{ϵ}}^3{p}_{d3}+\cdots ,\hfill \\ \varphi =\mathrm{ϵ}{\varphi }_1+{\mathrm{ϵ}}^2{\varphi }_2+{\mathrm{ϵ}}^3{\varphi }_3+\cdots \hfill \end{array}$$

(5)

Substituting equations (4) and (5) into equation (1), we obtain the following equation

$$\{\begin{array}{l}{\mathrm{ϵ}}^{\frac{3}{2}}\frac{\partial n_d}{\partial \tau }-{\mathrm{ϵ}}^{\frac{1}{2}}{v}_0\frac{\partial n_d}{\partial \xi }+{\mathrm{ϵ}}^{\frac{1}{2}}\frac{\partial n_d{u}_d}{\partial \xi }+\mathrm{ϵ}\frac{\partial n_d{v}_d}{\partial \eta }+\mathrm{ϵ}\frac{\partial n_d{w}_d}{\partial \zeta }=0,\hfill \\ {\mathrm{ϵ}}^{\frac{3}{2}}\frac{\partial u_d}{\partial \tau }-{\mathrm{ϵ}}^{\frac{1}{2}}{v}_0\frac{\partial u_d}{\partial \xi }+{\mathrm{ϵ}}^{\frac{1}{2}}{u}_d\frac{\partial u_d}{\partial \xi }+\mathrm{ϵ}{v}_d\frac{\partial u_d}{\partial \eta }+\mathrm{ϵ}{w}_d\frac{\partial u_d}{\partial \zeta }+{\mathrm{ϵ}}^{\frac{1}{2}}\frac{\sigma }{n_d}\frac{\partial p_d}{\partial \xi }={\mathrm{ϵ}}^{\frac{1}{2}}\frac{\partial \varphi }{\partial \xi },\hfill \\ {\mathrm{ϵ}}^{\frac{3}{2}}\frac{\partial v_d}{\partial \tau }-{\mathrm{ϵ}}^{\frac{1}{2}}{v}_0\frac{\partial v_d}{\partial \xi }+{\mathrm{ϵ}}^{\frac{1}{2}}{u}_d\frac{\partial v_d}{\partial \xi }+\mathrm{ϵ}{v}_d\frac{\partial v_d}{\partial \eta }+\mathrm{ϵ}{w}_d\frac{\partial v_d}{\partial \zeta }+\mathrm{ϵ}\frac{\sigma }{n_d}\frac{\partial p_d}{\partial \eta }=\mathrm{ϵ}\frac{\partial \varphi }{\partial \eta },\hfill \\ {\mathrm{ϵ}}^{\frac{3}{2}}\frac{\partial w_d}{\partial \tau }-{\mathrm{ϵ}}^{\frac{1}{2}}{v}_0\frac{\partial w_d}{\partial \xi }+{\mathrm{ϵ}}^{\frac{1}{2}}{u}_d\frac{\partial w_d}{\partial \xi }+\mathrm{ϵ}{v}_d\frac{\partial w_d}{\partial \eta }+\mathrm{ϵ}{w}_d\frac{\partial w_d}{\partial \zeta }+\mathrm{ϵ}\frac{\sigma }{n_d}\frac{\partial p_d}{\partial \zeta }=\mathrm{ϵ}\frac{\partial \varphi }{\partial \zeta },\hfill \\ {\mathrm{ϵ}}^{\frac{3}{2}}\frac{\partial p_d}{\partial \tau }-{\mathrm{ϵ}}^{\frac{1}{2}}{v}_0\frac{\partial p_d}{\partial \xi }+{\mathrm{ϵ}}^{\frac{1}{2}}{u}_d\frac{\partial p_d}{\partial \xi }+\mathrm{ϵ}{v}_d\frac{\partial p_d}{\partial \eta }+\mathrm{ϵ}{w}_d\frac{\partial p_d}{\partial \zeta }\hfill \\ +\gamma p_d\left({\mathrm{ϵ}}^{\frac{1}{2}}\frac{\partial u_d}{\partial \xi }+\mathrm{ϵ}\frac{\partial v_d}{\partial \eta }+\mathrm{ϵ}\frac{\partial w_d}{\partial \zeta }\right)=0,\hfill \\ \mathrm{ϵ}\frac{{\partial }^2\varphi }{\partial x^2}+{\mathrm{ϵ}}^2\frac{{\partial }^2\varphi }{\partial y^2}+{\mathrm{ϵ}}^2\frac{{\partial }^2\varphi }{\partial z^2}=n_d+C_1\varphi +C_2{\varphi }^2+C_3{\varphi }^3\hfill \end{array}$$

(6)

According to the different power expansion of ϵ, the following equation is obtained

$$\mathrm{ϵ}:\{\begin{array}{l}-v_0\frac{\partial n_1}{\partial \xi }+\frac{\partial u_1}{\partial \xi }=0,\hfill \\ -v_0\frac{\partial u_1}{\partial \xi }+\sigma \frac{\partial p_1}{\partial \xi }=\frac{\partial {\varphi }_1}{\partial \xi },\hfill \\ -v_0\frac{\partial p_1}{\partial \xi }+\gamma \frac{\partial u_1}{\partial \xi }=0,\hfill \\ 0=n_1+{\varphi }_1\hfill \end{array}$$

(7)

$${\mathrm{ϵ}}^2:\{\begin{array}{l}\frac{\partial n_1}{\partial \tau }-v_0\frac{\partial n_2}{\partial \xi }+\frac{\partial n_1{u}_1}{\partial \xi }+\frac{\partial u_2}{\partial \xi }+\frac{\partial v_1}{\partial \eta }+\frac{\partial w_1}{\partial \zeta }=0,\hfill \\ \frac{\partial u_1}{\partial \tau }-v_0\frac{\partial u_2}{\partial \xi }-v_0{n}_1\frac{\partial u_1}{\partial \xi }+u_1\frac{\partial u_1}{\partial \xi }+\sigma \frac{\partial p_2}{\partial \xi }=\frac{\partial {\varphi }_2}{\partial \xi }+n_1\frac{\partial {\varphi }_1}{\partial \xi },\hfill \\ \frac{\partial p_1}{\partial \tau }-v_0\frac{\partial p_2}{\partial \xi }+u_1\frac{\partial p_1}{\partial \xi }+\gamma p_1\frac{\partial u_1}{\partial \xi }+\gamma \frac{\partial u_2}{\partial \xi }+\gamma \frac{\partial v_1}{\partial \eta }+\gamma \frac{\partial w_1}{\partial \zeta }=0,\hfill \\ \frac{{\partial }^2{\varphi }_1}{\partial x^2}=n_2+{\varphi }_2-\frac{\mu +\nu {\beta }^2}{2{\left(\mu +\nu \beta \right)}^2}{\varphi }_1^2\hfill \end{array}$$

(8)

The following equation is obtained under the low-order approximation of ε

$$\{\begin{array}{l}{n}_1=-{\varphi }_1,u_1=-v_0{\varphi }_1,p_1=-\gamma {\varphi }_1,v_0^2=1+\sigma \gamma ,\hfill \\ \frac{\partial v_1}{\partial \xi }=-v_0\frac{\partial {\varphi }_1}{\partial \eta },\frac{\partial w_1}{\partial \xi }=-v_0\frac{\partial {\varphi }_1}{\partial \zeta }\hfill \end{array}$$

(9)

Under the high-order approximation of ε, we obtain the (3 + 1) dimensional KP equation

$$\frac{\partial }{\partial \xi }\left(\frac{\partial {\varphi }_1}{\partial \tau }+a_1{\varphi }_1\frac{\partial {\varphi }_1}{\partial \xi }+a_2\frac{{\partial }^3{\varphi }_1}{{\partial }^3\xi }\right)+a_3\left(\frac{{\partial }^2{\varphi }_1}{{\partial }^2\eta }+\frac{{\partial }^2{\varphi }_1}{{\partial }^2\zeta }\right)=0$$

(10)

where

$$\{\begin{array}{l}{a}_1=-\frac{1}{2v_0}\left[1+\frac{3}{2}{v}_0^2+\frac{\gamma }{2}+\frac{{\gamma }^2\sigma }{2}-\frac{\mu +\nu {\beta }^2}{{\left(\mu +\nu \beta \right)}^2}\right],\hfill \\ a_2=\frac{1}{2v_0},a_3=\frac{v_0}{2}\hfill \end{array}$$

(11)

Letting ${\varphi }_1\left(\xi ,\eta ,\zeta ,\tau \right)=A\left(\xi ,\eta ,\zeta ,\tau \right)$ , equation (10) can be rewritten as follows

$$\frac{\partial }{\partial \xi }\left(\frac{\partial A}{\partial \tau }+a_{1}A\frac{\partial A}{\partial \xi }+a_2\frac{{\partial }^{3}A}{{\partial }^3\xi }\right)+a_3\left(\frac{{\partial }^{2}A}{{\partial }^2\eta }+\frac{{\partial }^{2}A}{{\partial }^2\zeta }\right)=0$$

(12)

Remark: a₁ is the nonlinear coefficient, and a₂ and a₃ are dispersion coefficients. Assuming there is no high-order lateral disturbance, equation (12) is the KdV equation. Compared with the lower dimensional model, equation (12) is more in line with the actual situation.

Derivation of the (3 + 1)-dimensional TSF-KP equation

Based on the above section, we obtain the (3 + 1)-dimensional KP equation. However, with the development of scientific research, the integer-order model is insufficient for studying the problems in the real world. Fractal calculus ⁴⁰ and fractional calculus^41-42 have become hot spots in mathematics and engineering. In 1998, He studied the analytical method of fractional percolation model of porous media, ⁴³ and we can see the significance of studying fractional percolation model. To deeply study the ultralow-frequency dust acoustic waves in dust plasma, in this section, we use the semi-inverse method and the fractional variation principle to derive the (3 + 1)-dimensional TSF-KP equation.

The fractional derivatives are defined as follows:

Definition 1. ³⁹ The Riemann–Liouville fractional derivative operator of a function $A(\xi ,\eta ,\zeta ,\tau )$ is defined as

$$D_t^{n}A=\{\begin{array}{ll}\frac{{\partial }_{m}A}{\partial t^m}\hfill & n=m\in N,\hfill \\ \frac{1}{\Gamma \left(m-n\right)}\frac{{\partial }_m}{\partial t^m}{\displaystyle {\int }_0^t\frac{A(\tau ,x)}{{\left(t-\tau \right)}^{n+1-m}}}\text{d}\tau ,\hfill & m-1<n<m,m\in N\hfill \end{array}$$

(13)

Definition 2. ³⁹ The modified Riemann–Liouville fractional derivative operator of a function $A(\xi ,\eta ,\zeta ,\tau )$ is defined as

$$D_t^{n}f\left(t\right)=\frac{1}{\Gamma \left(1-n\right)}\frac{d}{dt}\left[{\displaystyle {\int }_a^t\frac{f\left(\tau \right)-f\left(a\right)}{{\left(t-\tau \right)}^n}}\text{d}\tau \right],0\le n<1$$

(14)

In addition to the above definitions, He’s fractional derivative and Caputo’s definition are also frequently used and their definition are as follows:

Definition 3. ⁴⁴ He’s fractal derivative

$$\frac{Du}{D{x}^{\alpha }}=\Gamma (1+\alpha )\underset{\Delta x=x_1-x_2\to L}{\lim }\frac{u(x_1)-u(x_2)}{{(x_1-x_2)}^{\alpha }}$$

(15)

Definition 4. ⁴⁴ Caputo’s definition

$$D_x^{\alpha }f(x)=\frac{1}{\Gamma (n-\alpha )}{\displaystyle {\int }_0^x{(x-t)}^{n-\alpha -1}}\frac{{\text{d}}^{n}f(t)}{\text{d}{t}^n}\text{d}t$$

(16)

Equation (12) can be written as

$$A_{\tau }+a_{1}A{A}_{\xi }+a_{2}A{A}_{\xi \xi \xi }+a_3{D}^{-1}\left(A_{\eta \eta }+A_{\zeta \zeta }\right)=0$$

(17)

Letting $A(\xi ,\eta ,\zeta ,\tau )=B_{\xi }(\xi ,\eta ,\zeta ,\tau ),B(\xi ,\eta ,\zeta ,\tau )$ is a potential function, and the potential equation of the (3 + 1)-dimensional TSF-KP equation is as follows

$$B_{\xi \tau }+a_1{B}_{\xi }{B}_{\xi \xi }+a_2{B}_{\xi \xi \xi \xi }+a_3{B}_{\eta \eta }+a_3{B}_{\zeta \zeta }=0$$

(18)

The functional of equation (18) is as follows

$$J\left(B\right)={\displaystyle \underset{R}{\iiint }\text{d}}\xi \text{d}\eta \text{d}\zeta {\displaystyle {\int }_T\text{d}}\tau \left[B\left(b_1{B}_{\xi \tau }+b_2{a}_1{B}_{\xi }{B}_{\xi \xi }+b_3{a}_2{B}_{\xi \xi \xi \xi }+b_4{a}_3{B}_{\eta \eta }+b_5{a}_3{B}_{\zeta \zeta }\right)\right]$$

(19)

where $b_i(i=1,2,3,4,5)$ are Lagrangian multipliers.

Using the partial integration method for equation (19), and letting $B_{\xi }{|}_R=B_{\eta }{|}_R=B_{\zeta }{|}_R=B_{\tau }{|}_T=B_{\xi \xi }{|}_R=0$ , we obtain

$$J(B)={\displaystyle \underset{R}{\iiint }\text{d}}\xi \text{d}\eta \text{d}\zeta {\displaystyle {\int }_T\text{d}}\tau \left[-b_1{B}_{\xi }{B}_{\tau }-\frac{1}{2}{b}_2{a}_1{B}_{\xi }^3+b_3{a}_2{(B_{\xi \xi })}^2+b_4{a}_3{(B_{\eta })}^2+b_5{a}_3{(B_{\zeta })}^2)\right]$$

(20)

Using the variational method for equation (20) the following equation is obtained

$$\begin{array}{l}F\left(\xi ,\eta ,,\tau ,B,B_{\tau },B_{\xi },B_{\xi \xi },B_{\xi \eta },B_{\xi \zeta }\right)\hfill \\ =\frac{\partial F}{\partial B}-\frac{\partial }{\partial \tau }\left(\frac{\partial F}{\partial B_{\tau }}\right)-\frac{\partial }{\partial \xi }\left(\frac{\partial F}{\partial B_{\xi }}\right)-\frac{\partial }{\partial \eta }\left(\frac{\partial F}{\partial B_{\eta }}\right)-\frac{\partial }{\partial \zeta }\left(\frac{\partial F}{\partial B_{\zeta }}\right)+\frac{{\partial }^2}{\partial {\xi }^2}\left(\frac{\partial F}{\partial B_{\xi \xi }}\right)\hfill \\ =2b_1{B}_{\xi \tau }+3b_2{a}_1{B}_{\xi }{B}_{\xi \xi }+2b_3{a}_2{B}_{\xi \xi \xi \xi }+2b_4{a}_3{B}_{\eta \eta }+2b_5{a}_3{B}_{\zeta \zeta }=0\hfill \end{array}$$

(21)

Equation (21) is equivalent to equation (18). By comparing the coefficients, we obtain the Lagrangian multipliers $b_i(i=1,2,3,4,5)$ as follows

$$b_1=\frac{1}{2}, b_2=\frac{1}{3}, b_3=\frac{1}{2}, b_4=\frac{1}{2}, b_5=\frac{1}{2}$$

(22)

Therefore, the Lagrangian form of the (3 + 1)-dimensional integer-order KP equation is as follows

$$L\left(B_{\tau },B_{\xi },B_{\eta },B_{\zeta },B_{\xi \xi }\right)=-\frac{1}{2}{B}_{\tau }{B}_{\xi }-\frac{1}{6}{a}_1{\left(B_{\xi }\right)}^3+\frac{1}{2}{a}_2{\left(B_{\xi \xi }\right)}^2-\frac{1}{2}{a}_3{\left(B_{\eta }\right)}^2-\frac{1}{2}{a}_3{\left(B_{\zeta }\right)}^2$$

(23)

Similarly, the Lagrangian form of the (3 + 1)-dimensional TSF-KP equation is as follows

$$\begin{array}{c}F\left(D_{\tau }^{\omega }B,D_{\xi }^{\alpha }B,D_{\eta }^{\beta }B,D_{\zeta }^{\gamma }B,D_{\xi }^{\alpha \alpha }B\right)=-\frac{1}{2}{D}_{\tau }^{\omega }B{D}_{\xi }^{\alpha }B-\frac{1}{6}{a}_1{(D_{\xi }^{\alpha }B)}^3+\frac{1}{2}{a}_2{(D_{\xi }^{\alpha \alpha }B)}^2\\ -\frac{1}{2}{a}_3{(D_{\eta }^{\beta }B)}^2-\frac{1}{2}{a}_3{(D_{\zeta }^{\gamma }B)}^2\end{array}$$

(24)

where $D_{\xi }^{\alpha \alpha }B=D_{\xi }^{\alpha }(D_{\xi }^{\alpha }B)$ . Thus, the function of the (3 + 1)-dimensional TSF-KP equation is as follows

$$\begin{array}{l}{J}_F(B)={\displaystyle \underset{R}{\iiint }{(\text{d}\xi )}^{\alpha }}{(\text{d}\eta )}^{\beta }{(\text{d}\zeta )}^{\gamma }{\displaystyle {\int }_T{(\text{d}\tau )}^{\omega }}F\left(D_{\tau }^{\omega }B,D_{\xi }^{\alpha }B,D_{\eta }^{\beta }B,D_{\zeta }^{\gamma }B,D_{\xi }^{\alpha \alpha }B\right)\hfill \end{array}$$

(25)

According to Agrawal’s method, ⁴⁴ equation (25) changes to the following form

$$\begin{array}{l}\delta J_F\left(B\right)={\displaystyle {\int }_R{\left(\text{d}\xi \right)}^{\alpha }}{\displaystyle {\int }_R{\left(\text{d}\eta \right)}^{\alpha }}{\displaystyle {\int }_R{\left(\text{d}\zeta \right)}^{\gamma }}{\displaystyle {\int }_T{\left(\text{d}\tau \right)}^{\omega }}[\left(\frac{\partial F}{\partial D_{\tau }^{\omega }B}\right)\delta D_{\tau }^{\omega }B+\left(\frac{\partial F}{\partial D_{\xi }^{\alpha }B}\right)\delta D_{\xi }^{\alpha }B\\ +\left(\frac{\partial F}{\partial D_{\xi }^{\alpha \alpha }B}\right)\delta D_{\xi }^{\alpha \alpha }+\left(\frac{\partial F}{\partial D_{\eta }^{\beta }B}\right)\delta D_{\eta }^{\beta }B+\left(\frac{\partial F}{\partial D_{\zeta }^{\gamma }B}\right)\delta D_{\zeta }^{\gamma }B]\end{array}$$

(26)

Using fractional integration by parts

$${\displaystyle {\int }_a^b{\left(\text{d}\tau \right)}^i}f\left(x\right)D_x^{i}g\left(x\right)=\Gamma \left(1+i\right)\left[g\left(x\right)f\left(x\right){|}_a^b-{\displaystyle {\int }_a^b{\left(\text{d}x\right)}^i}g\left(x\right)D_x^{i}f\left(i\right)\right], f\left(x\right),g\left(x\right)\in \left[a,b\right]$$

(27)

Then, from equation (26), we obtain the following equation

$$\begin{array}{l}\delta J_F\left(B\right)={\displaystyle {\int }_R{\left(\text{d}\xi \right)}^{\alpha }}{\displaystyle {\int }_R{\left(\text{d}\eta \right)}^{\alpha }}{\displaystyle {\int }_R{\left(\text{d}\zeta \right)}^{\gamma }}{\displaystyle {\int }_T{\left(\text{d}\tau \right)}^{\omega }}[-D_{\tau }^{\omega }\left(\frac{\partial F}{\partial D_{\tau }^{\omega }B}\right)-D_{\xi }^{\alpha }\left(\frac{\partial F}{\partial D_{\xi }^{\alpha }B}\right)\\ -D_{\eta }^{\beta }\left(\frac{\partial F}{\partial D_{\eta }^{\beta }B}\right)-D_{\zeta }^{\gamma }\left(\frac{\partial F}{\partial D_{\zeta }^{\gamma }B}\right)+D_{\xi }^{\alpha \alpha }\left(\frac{\partial F}{\partial D_{\xi }^{\alpha \alpha }B}\right)]\end{array}$$

(28)

Letting $\delta J_F(B)=0$ , the Euler–Lagrange equation of the (3 + 1)-dimensional TSF-KP equation can be written as

$$-D_{\tau }^{\omega }\left(\frac{\partial F}{\partial D_{\tau }^{\omega }B}\right)-D_{\xi }^{\alpha }\left(\frac{\partial F}{\partial D_{\xi }^{\alpha }B}\right)-D_{\eta }^{\beta }\left(\frac{\partial F}{\partial D_{\eta }^{\beta }B}\right)-D_{\zeta }^{\gamma }\left(\frac{\partial F}{\partial D_{\zeta }^{\gamma }B}\right)+D_{\xi }^{\alpha \alpha }\left(\frac{\partial F}{\partial D_{\xi }^{\alpha \alpha }B}\right)=0$$

(29)

Substituting equation (24) into equation (29), we have

$$D_{\tau }^{\omega }{D}_{\xi }^{\alpha }B+a_1(D_{\xi }^{\alpha }B)D_{\xi }^{\alpha \alpha }B+a_2{D}_{\xi }^{\alpha \alpha \alpha \alpha }B+a_3{D}_{\eta }^{\beta \beta }B+a_3{D}_{\zeta }^{\gamma \gamma }B=0$$

(30)

Letting $D_{\xi }^{\alpha }B(\xi ,\eta ,\zeta ,\tau )=A(\xi ,\eta ,\zeta ,\tau )$ and substituting it into equation (30), we obtain the following equation

$$D_{\tau }^{\omega }A+a_{1}A{D}_{\xi }^{\alpha }A+a_2{D}_{\xi }^{\alpha \alpha \alpha }A+{\displaystyle \int a_3}{D}_{\eta }^{\beta \beta }A+a_3{D}_{\zeta }^{\gamma \gamma }A=0$$

(31)

Finding the fractional derivative of equation (31), we obtain

$$D_{\xi }^{\alpha }\left(D_{\tau }^{\omega }A+a_{1}A{D}_{\xi }^{\alpha }A+a_2{D}_{\xi }^{\alpha \alpha \alpha }A\right)+a_3{D}_{\eta }^{\beta \beta }A+a_3{D}_{\zeta }^{\gamma \gamma }A=0$$

(32)

Remark: Equation (31) is the (3 + 1)-dimensional TSF-KP equation. When $\omega =\alpha =\beta =\gamma =1$ , equation (30) is the integer-order (3 + 1)-dimensional KP equation. Therefore, the fractional model can better describe the propagation of dust sound waves in space.

Conservation laws of (3 + 1) dimensional time-fractional KP equation

To further study the properties of fractional-order models, we study the symmetry and conservation laws of the (3 + 1)-dimensional time-fractional KP equations.

Lie symmetry analysis

The (3 + 1)-dimensional time-fractional KP equation is as follows

$$D_{\tau }^{\omega }A+a_{1}A\frac{\partial A}{\partial \xi }+a_{2}A\frac{{\partial }^{3}A}{\partial {\xi }^3}+a_3{D}^{-1}\left(\frac{{\partial }^{2}A}{\partial {\eta }^2}+\frac{{\partial }^{2}A}{\partial {\zeta }^2}\right)=0$$

(33)

Equation (31) is transformed into the form of a fractional partial differential equation as shown below

$$\begin{array}{l}{D}_{\tau }^{\gamma }A=Q\left(\xi ,\eta ,\zeta ,\tau ,A,A_{\xi \xi \xi },D_{\tau }^{\omega }A,A_{\eta \eta },A_{\zeta \zeta },\cdots \right), \omega >0\hfill \end{array}$$

(34)

We consider that equation (4) is invariant under the following one-parameter Lie group of point transformations

$$\{\begin{array}{l}\overline{\xi }=\xi +\mathrm{ϵ}X(\xi ,\eta ,\zeta ,\tau ,A)+O({\mathrm{ϵ}}^2),\hfill \\ \overline{\eta }=\eta +\mathrm{ϵ}Y(\xi ,\eta ,\zeta ,\tau ,A)+O({\mathrm{ϵ}}^2),\hfill \\ \overline{\zeta }=\zeta +\mathrm{ϵ}Z(\xi ,\eta ,\zeta ,\tau ,A)+O({\mathrm{ϵ}}^2),\hfill \\ \overline{\tau }=\tau +\mathrm{ϵ}T(\xi ,\eta ,\zeta ,\tau ,A)+O({\mathrm{ϵ}}^2),\hfill \\ \overline{A}=A+\mathrm{ϵ}\psi (\xi ,\eta ,\zeta ,\tau ,A)+O({\mathrm{ϵ}}^2),\hfill \\ D_{\tau }^{\omega }\overline{A}\to D_{\tau }^{\omega }A+\mathrm{ϵ}{\psi }_{\omega }^{\tau }+O({\mathrm{ϵ}}^2),\hfill \\ \frac{\partial \overline{A}}{\partial \xi }\to \frac{\partial A}{\partial \xi }+\mathrm{ϵ}{\psi }_{\xi }+O({\mathrm{ϵ}}^2),\hfill \\ \frac{{\partial }^3\overline{A}}{\partial {\xi }^3}\to \frac{{\partial }^{3}A}{\partial {\xi }^3}+\mathrm{ϵ}{\psi }_{\xi \xi \xi }+O({\mathrm{ϵ}}^2),\hfill \\ \frac{{\partial }^2\overline{A}}{\partial {\eta }^2}\to \frac{{\partial }^{2}A}{\partial {\eta }^2}+\mathrm{ϵ}{\psi }_{\eta \eta }+O({\mathrm{ϵ}}^2),\hfill \\ \frac{{\partial }^2\overline{A}}{\partial {\zeta }^2}\to \frac{{\partial }^{2}A}{\partial {\zeta }^2}+\mathrm{ϵ}{\psi }_{\zeta \zeta }+O({\mathrm{ϵ}}^2)\hfill \end{array}$$

(35)

where $\mathrm{ϵ}\ll 1$ is the group parameter, ξ, ζ, ψ, τ, η are the infinitesimals of the transformations. The expressions ${\eta }_{\omega }^T,{\eta }_{\alpha }^X,{\eta }_{\alpha }^{XXX},{\eta }_{\alpha ,\gamma }^{XYY},{\eta }_{\alpha ,\gamma }^{XZZ}$ are as follows

$$\{\begin{array}{l}{\psi }_{\omega }^{\tau }=D_{\tau }^{\omega }(\psi )+X{D}_{\tau }^{\omega }(A_{\xi })-D_{\tau }^{\omega }(X{A}_{\xi })+Y{D}_{\tau }^{\omega }(A_{\eta })-D_{\tau }^{\omega }(Y{A}_{\eta })+Z{D}_{\tau }^{\omega }(A_{\zeta })\hfill \\ -D_{\tau }^{\omega }(Z{A}_{\zeta })+D_{\tau }^{\omega }(D_{\tau }(T)A)-D_{\tau }^{\gamma +1}(\tau A)+T{D}_{\tau }^{\gamma +1}(A),\\ {\psi }_{\xi }=D_{\xi }(\psi )-A_{\xi }{D}_{\xi }(X)-A_{\eta }{D}_{\xi }(Y)-A_{\zeta }{D}_{\xi }(Z)-A_{\tau }{D}_{\xi }(T),\hfill \\ {\psi }_{\xi \xi \xi }=D_{\xi }({\psi }^{\xi \xi })-A_{\xi \xi \xi }{D}_{\xi }(X)-A_{\xi \xi \eta }{D}_{\xi }(Y)-A_{\xi \xi \zeta }{D}_{\xi }(Z)-A_{\xi \xi \tau }{D}_{\xi }(T),\hfill \\ {\psi }_{\eta \eta }=D_{\eta }({\psi }^{\eta })-A_{\xi \eta }{D}_{\eta }(X)-A_{\eta \eta }{D}_{\eta }(Y)-A_{\eta \zeta }{D}_{\eta }(Z)-A_{\eta \tau }{D}_{\eta }(T),\hfill \\ {\psi }_{\zeta \zeta }=D_{\zeta }({\psi }^{\zeta })-A_{\xi \zeta }{D}_{\zeta }(X)-A_{\eta \zeta }{D}_{\zeta }(Y)-A_{\zeta \zeta }{D}_{\zeta }(Z)-A_{\zeta \tau }{D}_{\zeta }(T)\hfill \end{array}$$

(36)

where $D_{\tau }^{\omega }$ is the total fractional derivative operator. $D_{\tau },D_{\xi },D_{\eta }$ , and $D_{\zeta }$ are the total differentiation of τ, ξ, η, and ζ, respectively, and their definitions are as follows

$$\{\begin{array}{l}{D}_{\tau }=\frac{\partial }{\partial \tau }+A_{\tau }\frac{\partial }{\partial \tau }+A_{\tau \tau }\frac{\partial }{\partial A_{\tau }}+A_{\xi \tau }\frac{\partial }{\partial A_{\xi }}+A_{\eta \tau }\frac{\partial }{\partial A_{\eta }}+A_{\zeta \tau }\frac{\partial }{\partial A_{\zeta }}+\cdots ,\hfill \\ D_{\xi }=\frac{\partial }{\partial \xi }+A_{\xi }\frac{\partial }{\partial A}+A_{\xi \xi }\frac{\partial }{\partial A_{\xi }}+A_{\tau \xi }\frac{\partial }{\partial A_{\tau }}+A_{\eta \xi }\frac{\partial }{\partial A_{\eta }}+A_{\zeta \xi }\frac{\partial }{\partial A_{\zeta }}+\cdots .\hfill \\ D_{\eta }=\frac{\partial }{\partial \eta }+A_{\eta }\frac{\partial }{\partial A}+A_{\eta \eta }\frac{\partial }{\partial A_{\eta }}+A_{\tau \eta }\frac{\partial }{\partial A_{\tau }}+A_{\xi \eta }\frac{\partial }{\partial A_{\xi }}+A_{\zeta \eta }\frac{\partial }{\partial A_{\zeta }}+\cdots .\hfill \\ D_{\zeta }=\frac{\partial }{\partial \zeta }+A_{\zeta }\frac{\partial }{\partial A}+A_{\zeta \zeta }\frac{\partial }{\partial A_{\zeta }}+A_{\tau \zeta }\frac{\partial }{\partial A_{\tau }}+A_{\eta \zeta }\frac{\partial }{\partial A_{\eta }}+A_{\xi \zeta }\frac{\partial }{\partial A_{\xi }}+\cdots \hfill \end{array}$$

(37)

Here, we give the definition of the generalized Leibnitz rule ⁴⁵ as follows

$$D_T^{\omega }\left(f\left(t\right)g\left(t\right)\right)={\displaystyle \sum _{n=0}^{\infty }\left(\begin{array}{c}\omega \\ n\end{array}\right)}{D}_t^{\omega -n}f(t)D_t^{n}g(t),\omega >0$$

(38)

where

$$\left(\begin{array}{c}\omega \\ n\end{array}\right)=\frac{{\left(-1\right)}^{n-1}\omega \Gamma \left(n-\omega \right)}{\Gamma \left(1-\omega \right)\Gamma \left(n+1\right)}$$

(39)

Applying the generalized Leibnitz rule, we have

$$\begin{array}{c}{\psi }_{\omega }^{\tau }=D_{\tau }^{\omega }(\psi )-\omega D_{\tau }^{\omega }(T)\frac{{\partial }^{\omega }A}{\partial {\tau }^{\omega }}-{\displaystyle \sum _{n=1}^{\infty }\left(\begin{array}{l}\omega \hfill \\ n\hfill \end{array}\right)}{D}_{\tau }^n(X)D_{\tau }^{\omega -n}{A}_{\xi }-{\displaystyle \sum _{n=1}^{\infty }\left(\begin{array}{l}\omega \hfill \\ n\hfill \end{array}\right)}{D}_{\tau }^n(Y)·D_{\tau }^{\omega -n}{A}_{\eta }\\ -{\displaystyle \sum _{n=1}^{\infty }\left(\begin{array}{l}\omega \hfill \\ n\hfill \end{array}\right)}{D}_{\tau }^n(Z)D_{\tau }^{\omega -n}{A}_{\zeta }-{\displaystyle \sum _{n=1}^{\infty }\left(\begin{array}{l}\omega \hfill \\ n+1\hfill \end{array}\right)}{D}_{\tau }^{n+1}(T)D_{\tau }^{\omega -n}A\end{array}$$

(40)

The chain rule for a compound function ⁴⁶ is defined as

$$\frac{d^{m}f\left(g\left(t\right)\right)}{d{t}^m}={\displaystyle \sum _{k=0}^m{\displaystyle \sum _{r=0}^k\left(\begin{array}{c}k\\ r\end{array}\right)}}\frac{1}{k!}\left[-g{\left(t\right)}^{k-r}\right]\frac{d^{k}f\left(g\left(t\right)\right)}{d{t}^k}$$

(41)

According to equations (38) and (41), when f(t) = 1, we can obtain

$$D_{\tau }^{\omega }\psi =\frac{{\partial }^{\omega }\psi }{\partial {\tau }^{\omega }}+{\psi }_A\frac{{\partial }^{\omega }A}{\partial {\tau }^{\omega }}-A\frac{{\partial }^{\omega }{\psi }_A}{\partial {\tau }^{\omega }}+{\displaystyle \sum _{n=1}^{\infty }\left(\begin{array}{c}\omega \\ n\end{array}\right)}\frac{{\partial }^n{\psi }_A}{\partial {\tau }^n}{D}_{\tau }^{\omega -n}A+Ra$$

(42)

where

$$Ra={\displaystyle \sum _{n=2}^{\infty }{\displaystyle \sum _{m=2}^n{\displaystyle \sum _{k=2}^m{\displaystyle \sum _{r=0}^{k-1}\left[\left(\begin{array}{l}\omega \hfill \\ n\hfill \end{array}\right)\left(\begin{array}{l}n\hfill \\ m\hfill \end{array}\right)\left(\begin{array}{l}k\hfill \\ r\hfill \end{array}\right)\frac{1}{k!}\frac{{\tau }^{n-\omega }}{\Gamma \left(n+1-\omega \right)}{\left(-A\right)}^r\frac{{\partial }^A}{\partial {\tau }^A}{\left(A\right)}^{k-r}\frac{{\partial }^{n-m+k}\psi }{\partial {\tau }^{n-m}\partial A^k}\right]}}}}$$

(43)

Therefore, equation (40) can be written as follows

$$\begin{array}{c}{\psi }_{\omega }^{\tau }=\frac{{\partial }^{\omega }\psi }{\partial {\tau }^{\omega }}+\left({\psi }_A-\omega D_{\tau }\left(T\right)\right)\frac{{\partial }^{\omega }A}{\partial {\tau }^{\omega }}-A\frac{{\partial }^{\omega }{\psi }_A}{\partial {\tau }^{\omega }}+{\displaystyle \sum _{n=1}^{\infty }\left[\left(\begin{array}{l}\omega \hfill \\ n\hfill \end{array}\right)\frac{{\partial }^{\omega }{\psi }_A}{\partial {\tau }^{\omega }}-\left(\begin{array}{l}\omega \hfill \\ n+1\hfill \end{array}\right)D_{\tau }^{n+1}(T)\right]D_{\tau }^{\omega -n}A}\\ -{\displaystyle \sum _{n=1}^{\infty }\left(\begin{array}{l}\omega \hfill \\ n\hfill \end{array}\right)}\left[D_{\tau }^n(X)D_{\tau }^{\omega -n}(A_{\xi })+D_{\tau }^n(Y)D_{\tau }^{\omega -n}(A_{\eta })+D_{\tau }^n(Z)D_{\tau }^{\omega -n}{A}_{\zeta }\right]+Ra\end{array}$$

(44)

Through the Lie symmetry theory, we obtain the infinitesimal generators M as follows

$$M=X\frac{\partial }{\partial \xi }+Y\frac{\partial }{\partial \eta }+Z\frac{\partial }{\partial \zeta }+T\frac{\partial }{\partial \tau }+\psi \frac{\partial }{\partial A}$$

(45)

The invariance criteria of system equation (34) are as follows

$$\{\begin{array}{c}P{r}^{(n)}M(\Delta ){|}_{\Delta =0}=0,n=1,2,3,\cdots ,\\ \Delta =D_{\tau }^{\omega }A+a_{1}A\frac{\partial A}{\partial \xi }+a_2\frac{{\partial }^{3}A}{\partial {\xi }^3}+a_3{D}^{-1}\left(\frac{{\partial }^{2}A}{\partial {\eta }^2}+\frac{{\partial }^{2}A}{\partial {\zeta }^2}\right)\end{array}$$

(46)

According to equations (44) and (46), we have

$$\begin{array}{c}P{r}^{(4)}M(\Delta )=T\frac{{\partial }^{\omega }}{\partial {\tau }^{\omega }}+X\frac{\partial }{\partial \xi }+Y\frac{\partial }{\partial \eta }+Z\frac{\partial }{\partial \zeta }+\psi \frac{\partial }{\partial A}+{\psi }_{\tau }^{\omega }\frac{\partial }{\partial D_{\tau }^{\gamma }A}\\ +{\psi }_{\xi }\frac{\partial }{\partial A_{\xi }}+{\psi }_{\xi \xi \xi }\frac{\partial }{\partial A_{\xi \xi \xi }}+{\psi }_{\eta \eta }\frac{\partial }{\partial A_{\eta \eta }}+{\psi }_{\zeta \zeta }\frac{\partial }{\partial A_{\zeta \zeta }}\end{array}$$

(47)

Using the second prolongation of equation (33), we can obtain the following invariance criteria

$${\psi }_{\omega }^{\tau }+a_1\psi \frac{\partial }{\partial \xi }+a_1{\psi }_{\xi }A+a_2{\psi }_{\xi \xi \xi }+a_3{D}^{-1}\left({\psi }_{\eta \eta }+{\psi }_{\zeta \zeta }\right)=0$$

(48)

By substituting equations (36) and (44) into equation (48) and setting each individual coefficient equal to 0, the following set of equations can be obtained

$$\{\begin{array}{l}\left(\begin{array}{l}\omega \hfill \\ n\hfill \end{array}\right)\frac{{\partial }^{\omega }{\psi }_A}{\partial {\tau }^{\omega }}-\left(\begin{array}{l}\omega \hfill \\ n+1\hfill \end{array}\right)D_{\tau }^{n+1}(\tau )=0,\hfill \\ X_A=X_{\tau }=0,\hfill \\ Y_A=Y_{\tau }=Y_{\xi }=0,\hfill \\ Z_A=Z_{\tau }=Z_{\xi }=0,\hfill \\ {\psi }_A-\omega T_{\tau }=0,\hfill \\ {\psi }_A-a_1{X}_{\xi }=0,\hfill \\ {\psi }_A-a_1{Y}_{\eta }=0,\hfill \\ {\psi }_A-a_1{Z}_{\zeta }=0\hfill \end{array}$$

(49)

By solving equation (49), we have following infinitesimals

$$\{\begin{array}{l}\psi =c_{1}A,X=\frac{c_1\xi }{a_1}+c_2,Y=\frac{c_1\eta }{a_3}+c_3,\hfill \\ Z=\frac{c_1\zeta }{a_3}+c_4,T=\frac{c_1\tau }{\omega }+c_5\hfill \end{array}$$

(50)

Therefore, the infinitesimal symmetric Lie algebra of equation (33) is spanned by the following infinitesimal generators

$$\{\begin{array}{l}{M}_1=\frac{\partial }{\partial \tau },\hfill \\ M_2=\frac{\partial }{\partial \xi },\hfill \\ M_3=\frac{\partial }{\partial \eta },\hfill \\ M_4=\frac{\partial }{\partial \zeta },\hfill \\ M_5=\frac{\xi }{a_1}\frac{\partial }{\partial \xi }+\frac{\eta }{a_3}\frac{\partial }{\partial \eta }+\frac{\zeta }{a_3}\frac{\partial }{\partial \zeta }+\frac{\tau }{\omega }\frac{\partial }{\partial \tau }-A\frac{\partial }{\partial A}\hfill \end{array}$$

(51)

Conservation laws

Based on the Lie symmetry analysis method, we study the conservation law of the (3 + 1)-dimensional time-fractional KP equation. The conservation law of equation (33) satisfies the following equation

$$D_{\tau }(C^{\tau })+D_{\xi }(C^{\xi })+D_{\eta }(C^{\eta })+D_{\zeta }(C^{\zeta })=0$$

(52)

where $C^{\tau },C^{\xi },C^{\eta }$ , and $C^{\zeta }$ are the conserved vectors.

A formal Lagrangian of equation (33) can be written in the following form

$$L=s\left(\xi ,\eta ,\zeta ,\tau \right)\left[(D_{\tau }^{\omega }A+a_{1}A\frac{\partial A}{\partial \xi }+a_2\frac{{\partial }^{3}A}{\partial {\xi }^3}+a_3{D}^{-1}\left(\frac{{\partial }^{2}A}{\partial {\eta }^2}+\frac{{\partial }^{2}A}{\partial {\zeta }^2}\right)\right]$$

(53)

where $s(\xi ,\eta ,\zeta ,\tau )$ is a new dependent variable. According to equation (53), an action integral is defined as follows

$${\displaystyle {\hspace{0.17em}}_{\int }^R{\displaystyle {\hspace{0.17em}}_{\int }^R{\displaystyle {\hspace{0.17em}}_{\int }^R{\displaystyle {\hspace{0.17em}}_{\int }^{T}L}}}}\left(\xi ,\eta ,\zeta ,\tau ,D_{\tau }^{\omega }A,A_{\xi },A_{\xi \xi \xi },A_{\eta \eta },A_{\zeta \zeta }\right)\text{d}\xi \text{d}\eta \text{d}\zeta \text{d}\tau$$

(54)

The Euler–Lagrangian operator is defined as

$$\begin{array}{ll}\frac{\delta }{\delta A}=\hfill & \frac{\partial }{\partial A}+{\left(D_{\tau }^{\omega }\right)}^{*}\frac{\partial }{\partial D_{\tau }^{\omega }A}+D_{\xi }\frac{\partial }{\partial A_{\xi }}-D_{\xi }^3\frac{\partial }{\partial A_{\xi \xi \xi }}-D_{\eta }^2\frac{\partial }{\partial A_{\eta \eta }}-D_{\zeta }^2\frac{\partial }{\partial A_{\zeta \zeta }}\hfill \end{array}$$

(55)

where ${(D_{\tau }^{\omega })}^{*}$ is the adjoint operator of $D_{\tau }^{\omega }$ , which is defined as follows

$${\left(D_{\tau }^{\omega }\right)}^{*}={\left(-1\right)}^n{I}_p^{n-\omega }\left(D_{\tau }^n\right){=}_{\tau }^C{D}_p^{\omega }$$

(56)

where $I_p^{n-\omega }$ is the right-sided fractional integral operator and ${\hspace{0.17em}}_{\tau }^C{D}_p^{\omega }$ is the right-sided Caputo fractional differential operator.

Therefore, we take the adjoint equation of equation (33) as the Euler–Lagrange equation, given by

$$F^{*}=\frac{\delta L}{\delta A}=0$$

(57)

and the above equation can be written as

$$F^{*}=a_1{A}_{\xi }+{\left(D_{\tau }^{\omega }\right)}^{*}s-a_1{D}_{\xi }\left(As\right)-a_2{D}_{\xi }^{3}s-a_3{D}^{-1}\left(D_{\eta }^{2}s+D_{\zeta }^{2}s\right)$$

(58)

The Lie characteristic function W is given by

$$W=\psi -T{A}_{\tau }-X{A}_{\xi }-Y{A}_{\eta }-Z{A}_{\zeta }$$

(59)

According to equation (51), we can obtain

$$\{\begin{array}{l}{W}_1=-A_{\xi },\hfill \\ W_2=-A_{\eta },\hfill \\ W_3=-A_{\zeta },\hfill \\ W_4=-A_{\tau },\hfill \\ W_5=-A-\frac{\xi }{a_1}{A}_{\xi }-\frac{\eta }{a_3}{A}_{\eta }-\frac{\zeta }{a_3}{A}_{\zeta }-\frac{\tau }{\omega }{A}_{\tau }\hfill \end{array}$$

(60)

According to the Riemann–Liouville fractional derivative, the component of conserved vectors is defined as

$$C^{\tau }=TL+{\displaystyle \sum _{k=0}^{n-1}{\left(-1\right)}_0^k}{D}_T^{\omega -1-k}\left(W_i\right)D_{\tau }^k\frac{\partial L}{\partial \left({\hspace{0.17em}}_0{D}_{\tau }^{\omega }A\right)}-{\left(-1\right)}^{n}J\left(W_i,D_{\tau }^n\frac{\partial L}{\partial \left({\hspace{0.17em}}_0{D}_{\tau }^{\omega }A\right)}\right)$$

(61)

where $J(·)$ is defined as

$$J\left(a,b\right)=\frac{1}{\Gamma \left(n-\omega \right)}{\displaystyle {\int }_0{\displaystyle {\int }_{\tau }\frac{f\left(T,\xi ,\eta ,\zeta \right)g\left(\mu ,\xi ,\eta ,\zeta \right)}{{\left(\mu -T\right)}^{\omega +1-n}}}}\text{\hspace{0.17em}d}\mu \text{d}T$$

(62)

$C^i(i=1,2,3)$ is defined as

$$\begin{array}{l}{C}^i={\rho }^{i}L+W_{\gamma }\left[\frac{\partial L}{\partial A_i}-D_i\left(\frac{\partial L}{\partial A_{ij}}\right)+D_i{D}_k\left(\frac{\partial L}{\partial A_{ijk}}\right)-\cdots \right]+D_j\left(W_{\gamma }\right)\left[\frac{\partial L}{\partial A_{ij}}-D_k\frac{\partial L}{\partial A_{ijk}}+\cdots \right]\\ +D_j{D}_k\left(W_{\gamma }\right)\left(\frac{\partial L}{\partial A_{ijk}}-\cdots \right)+\cdots \end{array}$$

(63)

where ${\rho }^1=X,{\rho }^2=Y,{\rho }^3=Z$ .

Using W₅ as an example to calculate the conserved vector of equation (31) by the above definition, we have

$$\{\begin{array}{l}{C}^{\tau }=TL{+}_0{D}_T^{\omega -1}\left(W_5\right)\frac{\partial L}{{\partial }_0{D}_{\tau }^{\gamma }A}+J\left(W_5,D_{\tau }\frac{\partial L}{\partial D_{\tau }^{\gamma }A}\right)\hfill \\ =s_0{D}_{\tau }^{\omega -1}\left(W_5\right)+J\left(W_5,s_{\tau }\right)\hfill \\ C^1=XL+W_5\left[\frac{\partial L}{\partial A_{\xi }}+D_{\xi }{D}_{\xi }\left(\frac{\partial L}{\partial A_{\xi \xi \xi }}\right)\right]+D_{\xi }\left(W_5\right)\left[-D_{\xi }\left(\frac{\partial L}{\partial A_{\xi \xi \xi }}\right)\right],\hfill \\ C^2=YL+W_5\left[-D_{\eta }\left(\frac{\partial L}{\partial A_{\eta \eta }}\right)\right]+D_{\eta }\left(W_5\right)\left[\frac{\partial L}{\partial A_{\eta \eta }}\right],\hfill \\ C^3=ZL+W_5\left[-D_{\zeta }\left(\frac{\partial L}{\partial A_{\zeta \zeta }}\right)\right]+D_{\zeta }\left(W_5\right)\left[\frac{\partial L}{\partial A_{\zeta \zeta }}\right]\hfill \end{array}$$

(64)

The above equations can be expressed as

$$\{\begin{array}{l}{C}^{\tau }=s_0{D}_{\tau }^{\omega -1}\left(-A-\frac{\xi }{a_1}{A}_{\xi }-\frac{\eta }{a_3}{A}_{\eta }-\frac{\zeta }{a_3}{A}_{\zeta }-\frac{\tau }{\omega }{A}_{\tau }\right)\hfill \\ +J\left[\left(-A-\frac{\xi }{a_1}{A}_{\xi }-\frac{\eta }{a_3}{A}_{\eta }-\frac{\zeta }{a_3}{A}_{\zeta }-\frac{\tau }{\omega }{A}_{\tau }\right),s_{\tau }\right],\hfill \\ C^1=XL+\left(-A-\frac{\xi }{a_1}{A}_{\xi }-\frac{\eta }{a_3}{A}_{\eta }-\frac{\zeta }{a_3}{A}_{\zeta }-\frac{\tau }{\omega }{A}_{\tau }\right)\left[\frac{\partial L}{\partial A_{\xi }}+D_{\xi }{D}_{\xi }\left(\frac{\partial L}{\partial A_{\xi \xi \xi }}\right)\right]\hfill \\ +D_{\xi }\left(-A-\frac{\xi }{a_1}{A}_{\xi }-\frac{\eta }{a_3}{A}_{\eta }-\frac{\zeta }{a_3}{A}_{\zeta }-\frac{\tau }{\omega }{A}_{\tau }\right)\left[-D_{\xi }\left(\frac{\partial L}{\partial A_{\xi \xi \xi }}\right)\right],\hfill \\ C^2=YL+\left(-A-\frac{\xi }{a_1}{A}_{\xi }-\frac{\eta }{a_3}{A}_{\eta }-\frac{\zeta }{a_3}{A}_{\zeta }-\frac{\tau }{\omega }{A}_{\tau }\right)\left[-D_{\eta }\left(\frac{\partial L}{\partial A_{\eta \eta }}\right)\right]\hfill \\ +D_{\eta }\left(-A-\frac{\xi }{a_1}{A}_{\xi }-\frac{\eta }{a_3}{A}_{\eta }-\frac{\zeta }{a_3}{A}_{\zeta }-\frac{\tau }{\omega }{A}_{\tau }\right)\left[\frac{\partial L}{\partial A_{\eta \eta }}\right],\hfill \\ C^3=ZL+\left(-A-\frac{\xi }{a_1}{A}_{\xi }-\frac{\eta }{a_3}{A}_{\eta }-\frac{\zeta }{a_3}{A}_{\zeta }-\frac{\tau }{\omega }{A}_{\tau }\right)\left[-D_{\zeta }\left(\frac{\partial L}{\partial A_{\zeta \zeta }}\right)\right]\hfill \\ +D_{\zeta }\left(-A-\frac{\xi }{a_1}{A}_{\xi }-\frac{\eta }{a_3}{A}_{\eta }-\frac{\zeta }{a_3}{A}_{\zeta }-\frac{\tau }{\omega }{A}_{\tau }\right)\left[\frac{\partial L}{\partial A_{\zeta \zeta }}\right]\hfill \end{array}$$

(65)

Remark: We have obtained a conservation law of W₅. Similarly, we can calculate the conservation laws of W₁, W₂, W₃, W₄, which are omitted here.

Solution of (3 + 1)-dimensional TSF-KP equation

In this section, we use the fractional subequation method to solve the (3 + 1)-dimensional TSF-KP equation and obtain three types of exact solutions. ⁴⁷

Introduction of method

Jumarie’s modified Riemann–Liouville derivative is defined as follows

$$D_{\tau }^{\omega }f(\tau )=\{\begin{array}{l}\frac{1}{\Gamma (1-\omega )}{\displaystyle {\int }_0^{\tau }{(\tau -T)}^{-\omega -1}}[f(T)-f(0)]\mathrm{d}T\omega <0,\hfill \\ \frac{1}{\Gamma (1-\omega )}{\displaystyle {\int }_0^{\tau }{(\tau -T)}^{-\omega }}[f(T)-f(0)]\mathrm{d}T0<\omega <1,\hfill \\ {[f^{\omega -n}(\tau )]}^{n}n\le \omega <n+1,n\ge 1\hfill \end{array}$$

(66)

Some properties of the modified Riemann–Liouville derivative ⁴⁸ are as follows

$$\begin{array}{l}{D}_{\tau }^{\omega }{\tau }^{\gamma }=\frac{\Gamma (\gamma +1)}{\Gamma (\gamma +1-\omega )}{\tau }^{\gamma -\omega }\gamma >0,\hfill \\ D_{\tau }^{\omega }(f(\tau )g(\tau ))=g(\tau )D_{\tau }^{\omega }f(\tau )+f(\tau )D_{\tau }^{\omega }g(\tau ),\hfill \\ D_{\tau }^{\omega }f[g(\tau )]=f\prime g[g(\tau )]D_{\tau }^{\omega }g(\tau )=D_{\tau }^{\omega }f[g(\tau )]{(g\prime (\tau ))}^{\omega }\hfill \end{array}$$

(67)

Step 1

Consider a fractional differential equation given by

$$P(A,A_{\xi },A_{\tau }\tau ,D_{\xi }^{\alpha },D_{\tau }^{\alpha })=0,0<\alpha \le 1$$

(68)

where $D_{\xi }^{\alpha }$ and $D_{\tau }^{\alpha }$ are Jumarie’s modified Riemann–Liouville derivatives, and $A=A(\xi ,\tau )$ is an unknown function.

Step 2

By using the following traveling wave transformation

$$A(\xi ,\tau )=A(\rho ),\rho =k\xi +c\tau$$

(69)

where k and c are constants, equation (68) is simplified to the following nonlinear fractional ordinary differential equation

$$P\left(A,kA\prime ,cA\prime ,k^{\alpha }{D}_{\rho }^{\alpha }A,c^{\alpha }{D}_{\rho }^{\alpha }A,\cdots \right)$$

(70)

Step 3

Supposing that equation (68) has the solution form

$$A(\rho )={\displaystyle \sum _{i=0}^n{a}_i}{\varphi }^i$$

(71)

where $a_i(i=1,2,3,\dots ,n)$ are constant and $\varphi =\varphi (\rho )$ satisfies the following Riccati equation

$$D_{\rho }^{\alpha }\varphi =\sigma +{\varphi }^2$$

(72)

where σ is a constant. The following solution of equation (72) is obtained by using the generalized Exp-function method via Mittag-Leffler functions

$$\varphi (\rho )=\{\begin{array}{ll}-\sqrt{-\sigma }{\tanh }_{\alpha }(\sqrt{-\sigma }\rho ),\hfill & \sigma <0,\hfill \\ -\sqrt{-\sigma }{\coth }_{\alpha }(\sqrt{-\sigma }\rho ),\hfill & \sigma <0,\hfill \\ \sqrt{\sigma }{\tan }_{\alpha }(\sqrt{\sigma }\rho ),\hfill & \sigma >0,\hfill \\ \sqrt{\sigma }{\tan }_{\alpha }(\sqrt{\sigma }\rho ),\hfill & \sigma >0,\hfill \\ -\frac{\Gamma (1+\alpha )}{{\rho }^{\alpha }\alpha +\kappa },\hfill & \kappa =const\hfill \end{array}$$

(73)

where the generalized hyperbolic and trigonometric functions are defined as

$$\varphi (\rho )=\{\begin{array}{l}{\sinh }_{\alpha }(x)=\frac{E_{\alpha }(x^{\alpha })-E_{\alpha }(-x^{\alpha })}{2},\hfill \\ {\cosh }_{\alpha }(x)=\frac{E_{\alpha }(x^{\alpha })+E_{\alpha }(-x^{\alpha })}{2},\hfill \\ {\tanh }_{\alpha }(x)=\frac{\sin {\mathrm{h}}_{\alpha }(x)}{\cos {\mathrm{h}}_{\alpha }(x)},{\coth }_{\alpha }(x)=\frac{{\cosh }_{\alpha }(x)}{{\sinh }_{\alpha }(x)},\hfill \\ {\sin }_{\alpha }(x)=\frac{E_{\alpha }(i{x}^{\alpha })-E_{\alpha }(-i{x}^{\alpha })}{2i},\hfill \\ {\cos }_{\alpha }(x)=\frac{E_{\alpha }(i{x}^{\alpha })+E_{\alpha }(-i{x}^{\alpha })}{2i},\hfill \\ {\tan }_{\alpha }(x)=\frac{{\sin }_{\alpha }(x)}{{\cos }_{\alpha }(x)},{\cot }_{\alpha }(x)=\frac{{\cos }_{\alpha }(x)}{{\sin }_{\alpha }(x)}\hfill \end{array}$$

(74)

where $E_{\alpha }(x^{\alpha })$ is the Mittag-Leffler function, which is defined as follows

$$E_{\alpha }(z)={\displaystyle \sum _{k=0}^{\infty }\frac{z^k}{\Gamma (1+k\alpha )}}$$

(75)

Step 4

Bringing equations (71) and (72) into equation (70) to obtain a polynomial for $\varphi (\rho )$ and setting all coefficients of ${\varphi }^i$ (i = 0, 1, 2, …) equal to zero, a set of overdetermined nonlinear algebraic equations for $b_i(i=0,1,2,\dots ,n)$ , k, and c are obtained.

Step 5

By solving the equations in step 4, the constants $b_i(i=0,1,2,\dots ,n)$ , k, and c are obtained. Substituting these constants and the solutions of equation (72) into equation (71), we can obtain the explicit solutions of equation (68).

Application of the method

The (3 + 1) dimensional TSF-KP equation is as follows

$$D_{\xi }^{\alpha }{D}_{\tau }^{\alpha }A+a_1{(D_{\xi }^{\alpha }A)}^2+a_{1}A{D}_{\xi }^{\alpha \alpha }A+a_2{D}_{\xi }^{\alpha \alpha \alpha \alpha }A+a_3{D}_{\eta }^{\alpha \alpha }A+a_3{D}_{\zeta }^{\alpha \alpha }A=0$$

(76)

First, through the following traveling wave transformation

$$A(\xi ,\eta ,\zeta ,\tau )=A(\rho ),\rho =k\xi +l\eta +m\zeta +c\tau$$

(77)

Equation (76) changes into a nonlinear fractional-order ordinary differential equation

$$k^{\alpha }{c}^{\alpha }{D}_{\rho }^{\alpha \alpha }A+a_1{\left(k^{\alpha }{D}_{\rho }^{\alpha }A\right)}^2+a_1{k}^{\alpha }A{D}_{\rho }^{\alpha \alpha }A+a_2{k}^{4\alpha }{D}_{\rho }^{\alpha \alpha \alpha \alpha }A+a_3{l}^{2\alpha }{D}_{\rho }^{\alpha \alpha }A+a_3{m}^{2\alpha }{D}_{\rho }^{\alpha \alpha }A=0$$

(78)

We suppose equation (76) has the following solution form

$$A(\rho )=b_0+b_1\varphi (\rho )$$

(79)

where $\varphi (\rho )$ satisfies the following equation

$$D_{\rho }^{\alpha }=\sigma +{\varphi }^2$$

(80)

Substituting equations (77), (79), and (80) into equation (78) and then setting all coefficients of ${\varphi }^i$ (i = 0, 1, 2, 3, 4) to zero, we obtain the following set of algebraic equations for $b_0,b_1,k,c$

$$\varphi (\rho )=\{\begin{array}{l}{\varphi }^0:b_1^2{a}_1{k}^{2\alpha }{\sigma }^2+6b_1{a}_2{k}^{4\alpha }{\sigma }^2=0,\hfill \\ {\varphi }^1:2b_1{k}^{\alpha }{c}^{\alpha }\sigma +2b_0{b}_1{a}_1{k}^{2\alpha }\sigma +4b_1{a}_2{k}^{4\alpha }{\sigma }^2+2b_1{a}_3{l}^{2\alpha }\sigma +2b_1{a}_3{m}^{2\alpha }\sigma =0,\hfill \\ {\varphi }^2:2b_1^2{a}_1{k}^{2\alpha }\sigma +2b_1^2{a}_1{k}^{2\alpha }\sigma +6b_1{a}_2{k}^{4\alpha }\sigma +18b_1{a}_2{k}^{4\alpha }\sigma =0,\hfill \\ {\varphi }^3:2b_1{k}^{\alpha }{c}^{\alpha }+2b_0{b}_1{a}_1{k}^{2\alpha }+4b_1{a}_2{k}^{4\alpha }\sigma +2b_1{a}_3{l}^{2\alpha }+2b_1{a}_3{m}^{2\alpha }=0,\hfill \\ {\varphi }^4:b_1^2{a}_1{k}^{2\alpha }+18b_1{a}_2{k}^{4\alpha }+2b_1^2{a}_1{k}^{2\alpha }=0\hfill \end{array}$$

(81)

Solving equation (81), we obtain

$$\varphi (\rho )=\{\begin{array}{l}{b}_0=b_0,\hfill \\ b_1=-\frac{6a_2}{a_1}{k}^{2\alpha },\hfill \\ \sigma =\frac{-a_3{l}^{2\alpha }-a_3{m}^{2\alpha }-k^{\alpha }{c}^{\alpha }-b_0{a}_1{k}^{2\alpha }}{2a_2{k}^{4\alpha }}\hfill \end{array}$$

(82)

The exact solutions of the three types of equation (76) are obtained as follows:

(I) When $\frac{-a_3{l}^{2\alpha }-a_3{m}^{2\alpha }-k^{\alpha }{c}^{\alpha }-b_0{a}_1{k}^{2\alpha }}{2a_2{k}^{4\alpha }}<0$ , two generalized hyperbolic function solutions are written as

$$\begin{array}{l}A=b_0+\frac{6a_2}{a_1}{k}^{2\alpha }\sqrt{\frac{a_3{l}^{2\alpha }+a_3{m}^{2\alpha }+k^{\alpha }{c}^{\alpha }+b_0{a}_1{k}^{2\alpha }}{2a_2{k}^{4\alpha }}}\hfill \\ {\tanh }_{\alpha }\left[\sqrt{\frac{a_3{l}^{2\alpha }+a_3{m}^{2\alpha }+k^{\alpha }{c}^{\alpha }+b_0{a}_1{k}^{2\alpha }}{2a_2{k}^{4\alpha }}}\left(k\xi +l\eta +m\zeta +c\tau \right)\right],\hfill \\ A=b_0+\frac{6a_2}{a_1}{k}^{2\alpha }\sqrt{\frac{a_3{l}^{2\alpha }+a_3{m}^{2\alpha }+k^{\alpha }{c}^{\alpha }+b_0{a}_1{k}^{2\alpha }}{2a_2{k}^{4\alpha }}}\hfill \\ {\coth }_{\alpha }\left[\sqrt{\frac{a_3{l}^{2\alpha }+a_3{m}^{2\alpha }+k^{\alpha }{c}^{\alpha }+b_0{a}_1{k}^{2\alpha }}{2a_2{k}^{4\alpha }}}\left(k\xi +l\eta +m\zeta +c\tau \right)\right]\hfill \end{array}$$

(83)

(II) When $\frac{-a_3{l}^{2\alpha }-a_3{m}^{2\alpha }-k^{\alpha }{c}^{\alpha }-b_0{a}_1{k}^{2\alpha }}{2a_2{k}^{4\alpha }}>0$ , two generalized trigonometric function solutions are written as

$$\begin{array}{l}A=b_0-\frac{6a_2}{a_1}{k}^{\alpha }\sqrt{\frac{a_3{l}^{2\alpha }+a_3{m}^{2\alpha }+k^{\alpha }{c}^{\alpha }+b_0{a}_1{k}^{2\alpha }}{2a_2{k}^{4\alpha }}}\hfill \\ {\tan }_{\alpha }\left[\sqrt{\frac{-a_3{l}^{2\alpha }-a_3{m}^{2\alpha }-k^{\alpha }{c}^{\alpha }-b_0{a}_1{k}^{2\alpha }}{2a_2{k}^{4\alpha }}}(k\xi +l\eta +m\zeta +c\tau )\right],\hfill \\ A=b_0+\frac{6a_2}{a_1}{k}^{2\alpha }\sqrt{\frac{a_3{l}^{2\alpha }+a_3{m}^{2\alpha }+k^{\alpha }{c}^{\alpha }+b_0{a}_1{k}^{2\alpha }}{2a_2{k}^{4\alpha }}}\hfill \\ {\cot }_{\alpha }\left[\sqrt{\frac{a_3{l}^{2\alpha }+a_3{m}^{2\alpha }+k^{\alpha }{c}^{\alpha }+b_0{a}_1{k}^{2\alpha }}{2a_2{k}^{4\alpha }}}(k\xi +l\eta +m\zeta +c\tau )\right]\hfill \end{array}$$

(84)

(III) When $\frac{-a_3{l}^{2\alpha }-a_3{m}^{2\alpha }-k^{\alpha }{c}^{\alpha }-b_0{a}_1{k}^{2\alpha }}{2a_2{k}^{4\alpha }}=0$ , the rational solution is written as

$$A=b_0+\frac{6a_2}{a_1}{k}^{2\alpha }\frac{\Gamma (1+\alpha )}{{(k\xi +l\eta +m\zeta +c\tau )}^{\alpha }+\mathrm{\varrho }}$$

(85)

In addition to the above methods, the Exp-function method ⁴⁹ is often used to solve the differential equation. Here we briefly explain this approach.

Consider the fractional differential equation (equation (68)). Using the fractional complex transform

$$\begin{array}{l}A(\xi ,\eta ,\zeta ,\tau )=U(\mathrm{\varrho }),\hfill \\ \mathrm{\varrho }=\frac{k{\xi }^{\alpha }}{\Gamma (1+\alpha )}+\frac{l{\eta }^{\beta }}{\Gamma (1+\beta )}+\frac{m{\zeta }^{\gamma }}{\Gamma (1+\gamma )}-\frac{c{\tau }^{\omega }}{\Gamma (1+\omega )}\hfill \end{array}$$

(86)

He and Li ⁵⁰ proposed the fractional complex transform; He expounded the physical basis in his recent review articles.

The fractional chain rule is defined as

$$\begin{array}{c}{D}_{\tau }^{\omega }A={\sigma }_{\tau }\frac{dU}{d\mathrm{\varrho }}{D}_{\tau }^{\omega }\mathrm{\varrho },D_{\xi }^{\omega }A={\sigma }_{\xi }\frac{dU}{d\mathrm{\varrho }}{D}_{\xi }^{\omega }\mathrm{\varrho },\\ D_{\eta }^{\omega }A={\sigma }_{\eta }\frac{dU}{d\mathrm{\varrho }}{D}_{\eta }^{\omega }\mathrm{\varrho },D_{\zeta }^{\omega }A={\sigma }_{\zeta }\frac{dU}{d\mathrm{\varrho }}{D}_{\zeta }^{\omega }\mathrm{\varrho }\end{array}$$

(87)

where ${\sigma }_{\tau },{\sigma }_{\xi },{\sigma }_{\eta },{\sigma }_{\zeta }$ are named as sigma indexes. We can choose ${\sigma }_{\tau }={\sigma }_{\xi }={\sigma }_{\eta }={\sigma }_{\zeta }=h$ , where h is a constant.

Substituting equations (84) and (87) into equation (68), we can rewrite equation (68) as follows

$$Q(U,U^{\prime },U^{\prime \prime },U^{\prime \prime \prime },\dots )=0$$

(88)

Suppose equation (86) has the following form of solution

$$U(\mathrm{\varrho })=\frac{{\displaystyle {\sum }_d^{n-s}{a}_n}\exp (n\mathrm{\varrho })}{{\displaystyle {\sum }_q^{m-p}{b}_m}\exp (m\mathrm{\varrho })}$$

(89)

where p, q, s, and d are positive integers, and a_n and b_m are unknown constants. In order to determine the value of p, q, s, and d, the lowest order linear term of equation (86) is balanced with the lowest order nonlinear term.

So both of these methods translate the fractional-order equation, but the difference is fractional subequation method turns the fractional-order equation into nonlinear fractional ODE, the Exp-function method turns the fractional-order equation into nonlinear ODE. The two methods assume different forms of the solution, so some new solutions can be found in the exponential method.

Remark: As α = 1, the solution obtained above can be used as the solution of the integer-order KP equation. The acquisition of the exact solution lays a theoretical foundation for further study of the dust sound waves in the dual-temperature dust plasma.

Chirp effect of dust acoustic waves

In optical soliton communication, the chirp effect refers to the phenomenon of the center wavelength shifting during pulse transmission. Drawing on that idea, in this section, we study the chirp effect of the dispersion and nonlinearity of the KP equation describing dust acoustic waves.

For the TSF-KP equation, taking the fractional order as α = 1, the equation can be rewritten as

$${(A_{\tau }+a_{1}A{A}_{\xi }+a_2{A}_{\xi \xi \xi })}_{\xi }+a_3(A_{\eta \eta }+A_{\zeta \zeta })=0$$

(90)

where a₁ is the nonlinear coefficient, and a₂ and a₃ are dispersion coefficients.

According to equation (83), the initial waveform is

$$A=\frac{6a_2}{a_1}{k}^2\sqrt{\frac{a_3{l}^2+a_3{m}^2+kc}{2a_2{k}^4}}\tanh \left[\sqrt{\frac{a_3{l}^2+a_3{m}^2+kc}{2a_2{k}^4}\left(k\xi +l\eta +m\zeta \right)}\right]$$

(91)

First, we consider the effect of the nonlinearity on the dust sound waves separately. The equation is written as

$${(A_{\tau }+a_{1}A{A}_{\xi })}_{\xi }=0$$

(92)

Investigating the condition of time τ from 0 to $\Delta \tau$ , where $\Delta \tau$ is an infinitesimal variable, and applying equation (91), the approximate solution of equation (92) is

$$\begin{array}{l}A\left(\Delta \tau ,\xi ,\eta ,\zeta \right)=\frac{6a_2}{a_1}{k}^2\sqrt{\frac{a_3{l}^2+a_3{m}^2+kc}{2a_2{k}^4}}\tanh \left[\sqrt{\frac{a_3{l}^2+a_3{m}^2+kc}{2a_2{k}^4}}\left(k\xi +l\eta +m\zeta \right)\right]\\ \times \exp \{-6\frac{a_3{l}^2+a_3{m}^2+kc}{k}\Delta \tau +6\frac{a_3{l}^2+a_3{m}^2+kc}{k}\\ \times {\tanh }^2\left[\sqrt{\frac{a_3{l}^2+a_3{m}^2+kc}{2a_2{k}^4}}\left(k\xi +l\eta +m\zeta \right)\right]\}\Delta \tau \end{array}$$

(93)

The phase of the wave is

$$\begin{array}{c}{\phi }_N=-6\frac{a_3{l}^2+a_3{m}^2+kc}{k}\Delta \tau +6\frac{a_3{l}^2+a_3{m}^2+kc}{k}{\tanh }^2\left[\sqrt{\frac{a_3{l}^2+a_3{m}^2+kc}{2a_2{k}^4}}\left(k\xi +l\eta +m\zeta \right)\right]\Delta \tau \end{array}$$

(94)

The chirp caused by the nonlinearity is

$$\begin{array}{l}\Delta {\vartheta }_N=-\nabla {\phi }_N=6\frac{a_3{l}^2+a_3{m}^2+kc}{k}\sqrt{\frac{a_3{l}^2+a_3{m}^2+kc}{2a_2{k}^4}}2\tanh [\sqrt{\frac{a_3{l}^2+a_3{m}^2+kc}{2a_2{k}^4}}(k\xi +l\eta +m\zeta )]\\ \times \{a1-{\tanh }^2[\sqrt{\frac{a_3{l}^2+a_3{m}^2+kc}{2a_2{k}^4}}(k\xi +l\eta +m\zeta )\left]\right\}\Delta \tau k\\ +6\frac{a_3{l}^2+a_3{m}^2+kc}{k}\sqrt{\frac{a_3{l}^2+a_3{m}^2+kc}{2a_2{k}^4}}2\tanh [\sqrt{\frac{a_3{l}^2+a_3{m}^2+kc}{2a_2{k}^4}}(k\xi +l\eta +m\zeta )]\\ \times \{a1-{\tanh }^2[\sqrt{\frac{a_3{l}^2+a_3{m}^2+kc}{2a_2{k}^4}}(k\xi +l\eta +m\zeta )\left]\right\}\Delta \tau l+6\frac{a_3{l}^2+a_3{m}^2+kc}{k}\\ \times \sqrt{\frac{a_3{l}^2+a_3{m}^2+kc}{2a_2{k}^4}}2\tanh [\sqrt{\frac{a_3{l}^2+a_3{m}^2+kc}{2a_2{k}^4}}(k\xi +l\eta +m\zeta )]\\ \times \{a1-{\tanh }^{}[\sqrt{\frac{a_3{l}^2+a_3{m}^2+kc}{2a_2{k}^4}}(k\xi +l\eta +m\zeta )\left]\right\}\Delta \tau m\end{array}$$

(95)

Second, we consider the effect of dispersion on dust sound waves separately. The equation is as follows

$${(A_{\tau }+a_2{A}_{\xi \xi \xi })}_{\xi }+a_3(A_{\eta \eta }+A_{\zeta \zeta })=0$$

(96)

Equation (96) can be expressed as

$$A_{\tau }=-a_2{A}_{\xi \xi \xi }-a_3{D}^{-1}(A_{\eta \eta }+A_{\zeta \zeta })$$

(97)

Using the same method, we can obtain the approximate solution of equation (97) as follows

$$\begin{array}{l}A(\Delta \tau ,\xi ,\eta ,\zeta )=\frac{6a_2}{a_1}{k}^2\sqrt{\frac{a_3{l}^2+a_3{m}^2+kc}{2a_2{k}^4}}\tanh \left[\sqrt{\frac{a_3{l}^2+a_3{m}^2+kc}{2a_2{k}^4}}\left(k\xi +l\eta +m\zeta \right)\right]\\ \times \exp \{4\frac{{(a_3{l}^2+a_3{m}^2+kc)}^{\frac{3}{2}}}{{(2a_2)}^{\frac{3}{2}}{k}^3}\tanh \left[\sqrt{\frac{a_3{l}^2+a_3{m}^2+kc}{2a_2{k}^4}}(k\xi +l\eta +m\zeta )\right]\Delta \tau \\ +3\frac{{(a_3{l}^2+a_3{m}^2+kc)}^{\frac{3}{2}}}{{(2a_2)}^{\frac{3}{2}}{k}^3}{\tanh }^3\left[\sqrt{\frac{a_3{l}^2+a_3{m}^2+kc}{2a_2{k}^4}}(k\xi +l\eta +m\zeta )\right]\Delta \tau \\ +\frac{a_1{a}_3{l}^2}{a_2^{\frac{3}{2}}{k}^3}\sqrt{\frac{a_3{l}^2+a_3{m}^2+kc}{2}}\tanh \left[\sqrt{\frac{a_3{l}^2+a_3{m}^2+kc}{2a_2{k}^4}}(k\xi +l\eta +m\zeta )\right]\Delta \tau \\ +\frac{a_1{a}_3{m}^2}{a_2^{\frac{3}{2}}{k}^3}\sqrt{\frac{a_3{l}^2+a_3{m}^2+kc}{2}}\tanh \left[\sqrt{\frac{a_3{l}^2+a_3{m}^2+kc}{2a_2{k}^4}}(k\xi +l\eta +m\zeta )\right]\Delta \tau \}\end{array}$$

(98)

The phase of the wave is

$$\begin{array}{l}{\phi }_D=4\frac{{(a_3{l}^2+a_3{m}^2+kc)}^{\frac{3}{2}}}{{(2a_2)}^{\frac{3}{2}}{k}^3}\tanh [\sqrt{\frac{a_3{l}^2+a_3{m}^2+kc}{2a_2{k}^4}}(k\xi +l\eta +m\zeta )]\Delta \tau \\ +3\frac{{(a_3{l}^2+a_3{m}^2+kc)}^{\frac{3}{2}}}{{(2a_2)}^{\frac{3}{2}}{k}^3}{\tanh }^3[\sqrt{\frac{a_3{l}^2+a_3{m}^2+kc}{2a_2{k}^4}}(k\xi +l\eta +m\zeta )]\Delta \tau \\ +\frac{a_1{a}_3{l}^2}{a_2^{\frac{3}{2}}{k}^3}\sqrt{\frac{a_3{l}^2+a_3{m}^2+kc}{2}}\tanh [\sqrt{\frac{a_3{l}^2+a_3{m}^2+kc}{2a_2{k}^4}}(k\xi +l\eta +m\zeta )]\Delta \tau \\ +\frac{a_1{a}_3{m}^2}{a_2^{\frac{3}{2}}{k}^3}\sqrt{\frac{a_3{l}^2+a_3{m}^2+kc}{2}}\tanh [\sqrt{\frac{a_3{l}^2+a_3{m}^2+kc}{2a_2{k}^4}}(k\xi +l\eta +m\zeta )]\Delta \tau \end{array}$$

(99)

The chirp caused by the dispersion is

$$\Delta {\vartheta }_D=-\nabla {\phi }_D$$

(100)

$\Delta {\vartheta }_D$ are explained in detail in the later section (see Appendix 1).

According to equations (99) and (100), the whole chirp is as follows

$$\Delta {\vartheta }_S=\Delta {\vartheta }_N+\Delta {\vartheta }_D$$

(101)

Figures 1 to 3 show the change in total chirp when the velocity v₀ of the dust acoustic wave takes different values. It can be seen that the nonlinearity gradually increases and disperses as the velocity v₀ of the dust acoustic wave increases.

Figure 1.

The total chirp of dust sound waves when $v_0=0.5$ .

Figure 2.

The total chirp of dust sound waves when $v_0=1$ .

Figure 3.

The total chirp of dust sound waves when $v_0=1.5$ .

Study on the nature of dust acoustic waves

Figures 4 and 5 show the variation of the solitary wave amplitude with number density when σ and β take different values. For a given σ, the solitary amplitude will decrease as β increases. For a given β, the solitary amplitude increases as σ increases. When σ and β are given, the amplitude of the solitary wave will increase at first and then decrease as the number density of the low-temperature ions increases.

Figure 4.

The amplitude of the solitary wave changes with μ when σ = 3 and $\beta =3,6,9$ .

Figure 5.

The amplitude of the solitary wave changes with μ when β = 3 and $\sigma =4,5,6$ .

Figure 6 shows the variation of the solitary wave amplitude with the number density of the low-temperature ions when α is taken at different values. As the fractional order is increased, the solitary wave amplitude will be reduced.

Figure 6.

The amplitude of the solitary wave changes with μ when $\alpha =0.8,0.9,1$ .

Discussion

According to the equation of motion for the dust particles, the (3 + 1)-dimensional KP equation is obtained by using multiscale analysis and a reduced perturbation method, and then, a new (3 + 1) equation is given by applying the semi-inverse method and the fractional variation principle. The dimensional TSF-KP equation is more suitable for describing the propagation of dust sound waves of plasma in space. Next, we study the conservation laws and exact solutions of the fractional KP equations and use the exact solutions to study the enthalpy effect of dust sound waves and the effect of dust sound velocity on the total enthalpy. Finally, the effects of dust temperature, ion temperature, ion number density, and fractional order on the acoustic characteristics of dust are discussed.

The concept of fractional derivatives is more suitable for modeling real-world problems than integer derivatives. However, it has certain features that lead to difficulties when applying to real-world problems. We suspect that this might be due to a different definition or the definition itself, Therefore, based on the Riemann–Liouville fractional derivative operator and Ji-Huan He’s fractal derivative, we propose possible forms of modification

1.

$$D_t^{n}f(t)=\frac{1}{\Gamma (1-n)}\frac{\text{d}}{\text{d}t}\left[{\displaystyle {\int }_a^t\frac{k_{1}f(\tau )-k_{2}f(a)}{{(k_{3}t-k_4\tau )}^n}}\text{d}\tau \right],0\le n<1$$

(102)

where $k_1,k_2,k_3$ , and k₄ are constants.

2.

$$D_t^{n}f(t)=\frac{1}{\Gamma (1-n)}\frac{\text{d}}{\text{d}t}\left[{\displaystyle {\int }_a^t\frac{f{(\tau )}^{k_1}-f{(a)}^{k_2}}{{(t^{k_3}-{\tau }^{k_4})}^n}}\text{d}\tau \right],0\le n<1$$

(103)

where $k_1,k_2,k_3$ , and k₄ are constants.

3.

$$\frac{Du}{D{x}^{\alpha }}=\Gamma (1+\alpha )\underset{\Delta x=x_1-x_2\to L}{\lim }\frac{k_{1}u(x_1)-k_{2}u(x_2)}{{(k_3{x}_1-k_4{x}_2)}^{\alpha }}$$

(104)

where $k_1,k_2,k_3$ , and k₄ are constants

4.

$$\frac{Du}{D{x}^{\alpha }}=\Gamma (1+\alpha )\underset{\Delta x=x_1-x_2\to L}{\lim }\frac{u{(x_1)}^{k_1}-u{(x_2)}^{k_2}}{{(x_1^{k_3}-x_2^{k_4})}^{\alpha }}$$

(105)

where $k_1,k_2,k_3$ , and k₄ are constants.

When $k_x=1(x=1,2,3,4)$ , the definition of the above amendments was the Riemann–Liouville splitting derivative operator and He’s fractal derivative. ⁵¹ These parameters have potential value for understanding the application of fractional order in dust plasma when $k_x=1(x=1,2,3,4)$ is taken as different values. For example, if $\frac{k_4}{k_3}$ can be used to express the wave number in the dust plasma direction. Can the definition of the above form better explain the problem of dust acoustic waves in dust plasma? In future research, we will continue to study such issues.