Exact solitary wave solutions and linear stability of the β -time-fractional Gardner equation in shallow water dynamics

Abstract

This study investigates the nonlinear dynamics of soliton structures governed by the fractional Gardner equation incorporating the $\beta$ -time derivative. The model is of significant importance because it generalizes classical nonlinear evolution equations to capture memory effects that arise in dispersive media, plasma physics, and shallow-water waves. Using the improved modified extended tanh function method, a variety of exact analytical solutions are derived, including bright, dark, and singular solitons, as well as periodic wave profiles. The influence of the fractional-order parameter $\beta$ on the amplitude and stability of the resulting waveforms is analyzed in detail. Linear stability analysis is also performed to identify the stable regions of soliton propagation. Numerical simulations and three-dimensional plots confirm the validity of the obtained analytical results and illustrate the impact of $\beta$ on the soliton structure. The findings reveal that decreasing $\beta$ enhances the amplitude and steepness of the soliton, demonstrating the strong memory-dependent behavior of the system. Compared with previous works on the classical and time-fractional Gardner equations, the present results extend analytical solution classes by employing the $\beta$ -fractional derivative, providing a deeper physical interpretation of fractional-order wave interactions. This novel approach contributes to a better understanding of nonlinear wave phenomena in fractional dispersive systems and bridges the gap between analytical and numerical frameworks.

Introduction

The investigation of nonlinear wave equations has played a pivotal role in advancing our understanding of wave propagation phenomena across a broad spectrum of disciplines, including shallow water dynamics, plasma physics, nonlinear optics, and geophysical flows.^1–5 Nonlinear partial differential equations (PDEs) such as the Korteweg-de Vries (KdV) equation, the modified KdV (mKdV) equation, and their generalizations have been successfully applied to describe solitons, cnoidal waves, and other nonlinear structures.^6,7 Among these models, the Gardner equation, which combines features of both the KdV and mKdV equations, stands out as a versatile framework for analyzing solitary waves influenced by both quadratic and cubic nonlinearities. This equation, sometimes referred to as the combined KdV–mKdV equation, has been widely used to capture a wider range of nonlinear interactions compared to its classical counterparts.^8–10 Historically, the Gardner equation emerged as a hybrid of the KdV and mKdV equations to capture the interaction of quadratic and cubic nonlinearities in shallow water systems.^11–13 Its extensions have been widely studied in plasma physics, optics, and fluid dynamics. Recent developments in fractional calculus have led to fractional Gardner models,^14,15 which incorporate temporal effects. In recent years, fractional modeling has been widely applied to nonlinear wave dynamics in both plasma physics and optical media. For instance, the computational analysis of hyper-geometric soliton waves in plasma physics was explored using the auxiliary equation method, ¹⁶ while assorted optical soliton structures of a nonlinear fractional model with the $\beta$ -derivative were investigated in optical fibers. ¹⁷ These studies demonstrate the versatility of fractional operators in describing memory-dependent phenomena. Motivated by these developments, the present work employs the $\beta$ -fractional derivative to obtain new exact analytical solutions of the Gardner equation, thereby extending the applicability of fractional calculus to nonlinear dispersive systems.

In the context of shallow water wave dynamics, the Gardner equation provides an effective description of unidirectional weakly nonlinear long waves where higher-order nonlinear effects cannot be neglected. Such models are essential for understanding the evolution of solitons, bores, and internal waves in coastal and oceanic systems. However, many physical systems inherently exhibit memory effects that classical integer-order PDEs cannot fully capture. To address this limitation, researchers have turned to fractional calculus, which generalizes differentiation and integration to noninteger orders.^18,19 Fractional derivatives allow for more accurate representations of processes with hereditary properties, anomalous diffusion, and long-range temporal correlations.^20,21

One particularly useful generalization is the $\beta$ -time-fractional derivative, which introduces an additional degree of freedom into the governing equations. ²² This fractional parameter $\beta$ not only enriches the mathematical structure of the equation but also has a direct impact on the physical characteristics of wave propagation, such as amplitude modulation, soliton width, and wave stability. ²³ Incorporating fractional operators into the Gardner equation results in a $\beta$ -time-fractional Gardner equation, which has emerged as a promising model for studying shallow water waves under more realistic physical conditions. This fractional framework provides deeper insight into nonlinear dispersive systems and their sensitivity to fractional-order effects.

Despite its theoretical appeal, solving nonlinear fractional PDEs remains a formidable challenge. Classical analytical techniques often fall short, motivating the development of advanced symbolic and computational methods. Among the various approaches, methods based on hyperbolic function expansions—such as the tanh method and its extensions—have proven particularly effective in deriving closed-form solutions. In this context, the improved modified extended tanh function (IMETF) method offers a systematic and generalized procedure for constructing exact traveling wave solutions. In contrast to its classical formulations, the IMETF method integrates enhanced balancing procedures and generalized transformation frameworks, allowing the derivation of a broader spectrum of analytical solutions. These include solitary wave structures (such as dark solitons), periodic solutions represented through Jacobi elliptic functions, quasi-periodic patterns expressed via Weierstrass elliptic functions, and singular waveforms.

The analytical determination of such solutions is not only of mathematical interest but also of physical significance. Exact soliton and elliptic solutions provide benchmark cases for validating numerical simulations and approximations, as well as deepening our understanding of the interplay between nonlinearity and dispersion in fractional-order systems.^24,25 Furthermore, the classification of wave solutions into solitary, periodic, and singular types reflects the diverse behaviors encountered in real-world wave propagation scenarios, particularly in shallow water environments where nonlinear and dispersive effects coexist.^26,27 In recent years, a variety of analytical and semi-analytical techniques have been developed to derive exact or approximate solutions for nonlinear and fractional-order differential equations. Among these are the enhanced direct algebraic method, the enhanced Kudryashov method, and the new projective Riccati equation method, which have shown strong capability in handling nonlinear dispersive systems. For instance, recent work ²⁸ demonstrates the application of hybrid analytical approaches to fractional models, further highlighting the diversity of available solution frameworks. In comparison, the improved modified extended tanh method (IMETM) employed in this study offers a compact and systematic algebraic structure that allows for the derivation of multiple solution families under a unified scheme.

The aim of the present study is to systematically analyze the dimensionless $\beta$ -time-fractional derivative form of the unperturbed Gardner equation using the IMETF. By constructing a rich family of exact solutions, this work highlights how fractional calculus and modern symbolic methods can be combined to explore nonlinear wave dynamics more comprehensively. The solutions obtained here are expected to shed light on the role of fractional parameters in shaping the properties of nonlinear waves and to contribute to a growing body of literature on fractional nonlinear wave equations.

Governing equation and model formulation

The unperturbed form of the Gardner-type shallow water wave equation, which accounts for both quadratic and cubic nonlinearities along with higher-order dispersive effects, can be expressed in its classical integer-order form as ²⁹ :

$$u_t+(a_{1}u+a_5{u}^2)u_x+a_2{u}_{xxx}+a_3{u}_{xxt}+a_4{u}_{xtt}=0$$

(1)

where $u(x,t)$ denotes the wave profile, and $a_1,a_2,a_3,a_4,a_5$ are real constants characterizing the effects of nonlinearity and dispersion. The nonlinear terms $(a_{1}u+a_5{u}^2)u_x$ represent quadratic–cubic nonlinear interactions, while the higher-order derivatives $u_{xxx},u_{xxt},u_{xtt}$ capture dispersive and mixed spatio-temporal effects.

To investigate the influence of memory effects and anomalous temporal dynamics, we generalize the above model by replacing the classical first-order time derivative $q_t$ with the $\beta$ -time-fractional derivative $D_t^{\beta }u(x,t)$ , yielding the fractional model:

$$D_t^{\beta }u+(a_{1}u+a_5{u}^2)u_x+a_2{u}_{xxx}+a_3{u}_{xxt}+a_4{u}_{xtt}=0,0<\beta \le 1$$

(2)

which will serve as the governing equation for the present analysis.

The fractional derivative $D_t^{\beta }f(t)$ is defined as^30–33:

$$D^{\beta }g(\tau )=\underset{h\to 0}{lim}\frac{g(\tau +h{(\tau +\frac{1}{\mathrm{\Gamma }(\beta )})}^{1-\beta })-g(\tau )}{h},\tau >0,0<\beta \le 1$$

(3)

The fractional derivative $D_t^{\beta }$ used in this study satisfies several important properties, which are summarized below:

$$\begin{array}{rll}{D}_t^{\beta }(af(t)+bg(t))& =a{D}_t^{\beta }f(t)+b{D}_t^{\beta }g(t),\mathrm{\forall }a,b\in \mathbb{R},& \text{( linearity) }\\ D_t^{\beta }(f(t)g(t))& =f(t)D_t^{\beta }g(t)+g(t)D_t^{\beta }f(t),& \text{( product rule) }\\ D_t^{\beta }\left(\frac{f(t)}{g(t)}\right)& =\frac{g(t)D_t^{\beta }f(t)-f(t)D_t^{\beta }g(t)}{{[g(t)]}^2},g(t)\ne 0,& \text{( quotient rule) }\\ D_t^{\beta }f(t)& ={(t+\frac{1}{\mathrm{\Gamma }(\beta )})}^{1-\beta }\frac{\mathrm{d}f(t)}{\mathrm{d}t},& \text{( fundamental definition) }\\ D_t^{\beta }(f(g(t)))& ={(g^{\prime }(t))}^{\beta }{D}_t^{\beta }f(g(t)),& \text{( fractional chain rule) }\end{array}$$

These properties show that the $\beta$ -time-fractional derivative is a local fractional operator that preserves the fundamental algebraic structure of classical calculus, including linearity and product-type rules. Unlike nonlocal fractional derivatives, it does not involve convolution kernels or historical memory effects. Instead, it modifies the classical first-order derivative through a time-dependent scaling factor ${(t+(1/\mathrm{\Gamma }(\beta )))}^{1-\beta }$ , which enables a smooth transition between integer-order and fractional-order dynamics while maintaining analytical tractability.

Thus, the parameter $\beta$ acts as a control parameter for fractional temporal dynamics: it modifies the dispersive and nonlinear behaviors of the Gardner equation and enriches the physical description of shallow water wave propagation. The novelty of this work lies in extending the Gardner framework to a $\beta$ -time-fractional form and systematically deriving exact analytical solutions using an IMETM adapted for fractional operators. Furthermore, a new stability analysis is performed to elucidate the influence of the fractional order $\beta$ on the temporal and structural stability of nonlinear shallow-water waves.

The organization of this study is as follows: Section “The employed technique” outlines the proposed methodology in brief. Section “Applying to studied model” applies the method to derive exact solutions for the examined model. Section “Analysis and visualization of the derived solution profiles” presents graphical representations of selected solutions to highlight the characteristics of the propagating wave. Finally, the last section provides the conclusion of the work.

The employed technique

In this section, we present the IMETM, an analytical approach widely applied to derive exact solutions of nonlinear evolution equations, including fractional-order and higher-order wave models. IMETM refines the classical tanh technique by incorporating auxiliary parameters and generalized transformations, which enhance its ability to manage strong nonlinearities that appear in various physical systems. Due to its flexibility and robustness, IMETM has become a preferred method in soliton theory, capable of generating diverse solutions such as dark solitons, singular solutions, and Jacobi elliptic solutions under different nonlinear frameworks.^34,35

Consider a general nonlinear partial differential equation (NLPDE):

$$V(u,D_t^{\beta }u,u_x,u_{xx},\dots )=0$$

(4)

To apply the IMETM to equation (4), the following steps are carried out:

• Step (1):

Transformation. We reduce equation (4) into an ordinary differential equation (ODE) through the wave transformation:

$$u(x,t)=W(z)e^{i\varphi },z=hx-\frac{k{(\frac{1}{\mathrm{\Gamma }(\beta )}+t)}^{\beta }}{\beta },\varphi =x-\frac{\omega {(\frac{1}{\mathrm{\Gamma }(\beta )}+t)}^{\beta }}{\beta }$$

(5)

where $k$ denotes the wave speed. Substitution of (5) transforms equation (4) into:

$$T(W,W^{\prime },W^{″},W^{{}^{‴}},W^{(4)},\dots )=0$$

(6)

Step (2):

Ansatz for the solution. The solution of the reduced ODE is assumed in the form:

$$M(z)=\sum _{j=0}^N{a}_j{\mathrm{{\rm Y}}}^j(z)+\sum _{j=-1}^{-N}{b}_{-j}{\mathrm{{\rm Y}}}^j(z)$$

(7)

where $\mathrm{{\rm Y}}(z)$ satisfies the auxiliary equation

$${\mathrm{{\rm Y}}}^{\prime }(z)=\sqrt{d_0+d_1\mathrm{{\rm Y}}(z)+d_2{\mathrm{{\rm Y}}}^2(z)+d_3{\mathrm{{\rm Y}}}^3(z)+d_4{\mathrm{{\rm Y}}}^4(z)}$$

(8)

Step (3):

Balancing principle. The parameter $N$ in (7) is determined by applying the homogeneous balance rule between the nonlinear and derivative terms.

Step (4):

Substitution. By substituting equations (7) and (8) into equation (6) and equating the coefficients of identical powers of $\mathrm{{\rm Y}}(z)$ , a corresponding system of nonlinear algebraic equations is obtained.

Step (5):

Solving the system. The resulting algebraic equations are solved using symbolic computation software such as Mathematica V.13.1 on a laptop core i5, 8 GB of RAM, and a cache memory of 6 MB, yielding the coefficients $a_j,b_j$ and the parameters $h,k,\omega$ .

Step (6):

Parameter specification. By assigning appropriate values to $d_0,d_1,d_2,d_3,d_4$ , different types of solutions can be constructed. For example:

• Case 1:

For $d_0=d_1=d_3=0$ ,

$$\mathrm{{\rm Y}}(z)=\sqrt{-\frac{d_2}{d_4}}\text{sech}\left(\sqrt{d_2}z\right),d_2>0,d_4<0$$

$$\mathrm{{\rm Y}}(z)=\sqrt{-\frac{d_2}{d_4}}\text{sec}\left(\sqrt{-d_2}z\right),d_2<0,d_4>0$$

Case 2:

For $d_1=d_3=0$ ,

$$\mathrm{{\rm Y}}(z)=\sqrt{-\frac{d_2{m}^2}{d_4(1+m^2)}}\text{sn}\left(\sqrt{\frac{-d_2}{1+m^2}}z\right),d_2<0,d_4>0,d_0=\frac{d_2^2{m}^2}{d_4(1+m^2{)}^2}$$

Case 3:

For $d_2=d_4=0$ ,

$$\mathrm{{\rm Y}}(z)=\mathrm{\wp }(\frac{\sqrt{d_3}}{2}z;-\frac{4d_1}{d_3},-\frac{4d_0}{d_3})$$

Case 4:

For $d_3=d_4=0$ ,

$$\mathrm{{\rm Y}}(z)=\exp \left(\sqrt{d_2}z\right)-\frac{d_1}{2d_2},d_0=\frac{d_1^2}{4d_2},d_2>0$$

Case 5:

For $d_1=d_0=0$ ,

$$\mathrm{{\rm Y}}(z)=\frac{1}{2}\sqrt{\frac{d_2}{d_4}}(\tanh \left(\frac{\sqrt{d_2}z}{2}\right)+1),d_3=2\sqrt{d_2{d}_4},d_2>0$$

Step (7):

Final construction. The constants $a_j,b_{-j}$ obtained in Step (5) are substituted back into equation (7), combined with the explicit forms of $\mathrm{{\rm Y}}(z)$ , to generate closed-form solutions of equation (4).

It is worth highlighting the advantages of IMETM compared with other analytical approaches. Hirota’s bilinear method is highly effective in producing multi-soliton solutions for integrable systems, but is less suitable for fractional models. Similarly, the inverse scattering transform provides a rigorous framework for soliton construction in integrable PDEs, but its computational complexity limits its applicability to fractional or nonintegrable equations. In contrast, IMETM offers notable strengths: it is versatile in generating a wide spectrum of solutions (including solitons, Jacobi elliptic, Weierstrass elliptic, and exponential solutions), it reduces the problem to solvable algebraic systems, and it integrates fractional derivatives naturally through the wave transformation. These features establish IMETM as a robust and efficient tool for exploring the nonlinear dynamics of fractional wave equations, including the $\beta$ -time-fractional Gardner model studied in this work.

Applying to studied model

Assuming the following, our target is to obtain solutions for equation (2):

$$u(x,t)=Q(z),\text{where}z=hx-\frac{k{(\frac{1}{\mathrm{\Gamma }(\beta )}+t)}^{\beta }}{\beta }$$

(9)

This transformation maps the original space–time-dependent solution $u(x,t)$ into a profile function $Q(z)$ that evolves in a traveling wave coordinate $z$ , with a time-dependent oscillatory phase $\varphi$ . The parameters involved are interpreted as follows:

• $h$ : A spatial scaling parameter controlling the stretching or compression of the soliton profile in space.

• $k$ : A velocity-related parameter associated with the temporal evolution of the wave in the coordinate frame. It governs the propagation speed in the $z$ direction after the fractional transformation.

• $\beta$ : The order of the beta time-fractional derivative, with $0<\beta \le 1$ . It encapsulates the memory effect and temporal nonlocality of the medium.

• $Q(z)$ : The wave profile in the transformed frame.

This transformation is essential for reducing the fractional partial differential equation into a more tractable form amenable to exact solution techniques such as the IMETM.

By substituting equation (9) into equation (2), the fractional NLPDE is transformed into an ODE with full derivatives, resulting in the following expression:

$$a_2{h}^3{Q}^{(3)}(z)-a_3{h}^{2}k{Q}^{(3)}(z)+a_{4}h{k}^2{Q}^{(3)}(z)+h(a_{5}Q{(z)}^2+a_{1}Q(z))Q^{\prime }(z)-k{Q}^{\prime }(z)=0$$

(10)

To apply the proposed method, it is first necessary to determine the integer value of $N$ . Using the balancing principle in equation (10), and equating the terms $Q^2(z)Q^{\prime }(z)$ with $Q^{(3)}(z)$ , we obtain $N=1$ .

The general solution of the derived ODE can be expressed as:

$$Q(z)=s_0+s_1\mathrm{{\rm Y}}(z)+\frac{s_2}{\mathrm{{\rm Y}}(z)}$$

(11)

where $s_0,s_1,s_2$ are constants determining background offset and soliton amplitude scaling. By substituting equation (11) together with equation (8) into equation (10), and then equating the coefficients of $\lambda (z)$ to zero, a system of nonlinear algebraic equations is obtained as follows:

$$\begin{array}{rl}& 3d_{3}h{s}_1(a_2{h}^2+k(a_{4}k-a_{3}h))+a_{1}h{s}_1^2+2a_{5}h{s}_0{s}_1^2=0\\ & 6d_{4}h{s}_1(a_2{h}^2+k(a_{4}k-a_{3}h))+a_{5}h{s}_1^3=0\\ & -6a_2{d}_0{h}^3{s}_2+6a_3{d}_0{h}^{2}k{s}_2-6a_4{d}_{0}h{k}^2{s}_2-a_{5}h{s}_2^3=0\\ & d_{2}h{s}_1(a_2{h}^2+k(a_{4}k-a_{3}h))+a_{5}h{s}_2{s}_1^2+h{s}_0{s}_1(a_5{s}_0+a_1)-k{s}_1=0\\ & -3a_2{d}_1{h}^3{s}_2+3a_3{d}_1{h}^{2}k{s}_2-3a_4{d}_{1}h{k}^2{s}_2-a_{1}h{s}_2^2-2a_{5}h{s}_0{s}_2^2=0\\ & -a_2{d}_2{h}^3{s}_2+a_3{d}_2{h}^{2}k{s}_2-a_4{d}_{2}h{k}^2{s}_2-a_{5}h{s}_1{s}_2^2-a_{5}h{s}_0^2{s}_2-a_{1}h{s}_0{s}_2+k{s}_2=0\end{array}$$

These equations are subsequently solved using Mathematica V.13.1 on a laptop core i5, 8 GB of RAM, and a cache memory of 6 MB, and the corresponding analytical results are summarized as follows:

• Case(1):

  $d_0=d_1=d_3=0$

  $$s_1=-\frac{\sqrt{6}\sqrt{a_5{d}_4(-h)s_0^2-d_{4}k}}{\sqrt{a_5}\sqrt{d_2}\sqrt{h}},s_2=0,a_1=-2a_5{s}_0,d_2=\frac{a_{5}h{s}_0^2+k}{h(a_2{h}^2-a_{3}hk+a_4{k}^2)}$$
  Then a bright soliton solution is obtained:
  $$\begin{array}{r}u(x,t)=s_0+\sqrt{\frac{6(a_{5}h{s}_0^2+k)}{a_{5}h}}\text{sech}\left(\frac{\sqrt{\frac{a_{5}h{s}_0^2+k}{a_2{h}^3+hk(a_{4}k-a_{3}h)}}(\beta hx-k{(\frac{1}{\mathrm{\Gamma }(\beta )}+t)}^{\beta })}{\beta }\right)\end{array}$$

  (12)
  and a singular periodic solution is obtained :
  $$u(x,t)=s_0-\sqrt{6}\sqrt{\frac{a_{5}h{s}_0^2+k}{a_{5}h}}\text{sec}\left(\frac{\sqrt{\frac{-(a_{5}h{s}_0^2+k)}{a_2{h}^3+hk(a_{4}k-a_{3}h)}}(\beta hx-k{(\frac{1}{\mathrm{\Gamma }(\beta )}+t)}^{\beta })}{\beta }\right)$$

• Case (2):

  $d_1=d_3=0,d_0=\frac{d_2^2{m}^2}{d_4(1+m^2{)}^2},0<m\le 1$

  $$s_1=\frac{\sqrt{6}\sqrt{a_5{d}_4(-h)s_0^2-d_{4}k}}{\sqrt{a_5}\sqrt{d_2}\sqrt{h}},s_2=0,a_1=-2a_5{s}_0,d_2=\frac{a_{5}h{s}_0^2+k}{h(a_2{h}^2-a_{3}hk+a_4{k}^2)}$$
  then a Jacobi periodic solution is obtained:
  $$\begin{array}{rl}& u(x,t)=\\ & \frac{\sqrt{\frac{6m^2(a_{5}h{s}_0^2+k)}{m^2+1}}\text{sn}\left((hx-k{(t+\frac{1}{\mathrm{\Gamma }(\beta )})}^{\beta })\sqrt{-\frac{h{a}_5{s}_0^2+k}{h(m^2+1)(a_2{h}^2-k{a}_{3}h+k^2{a}_4)}}|m\right)}{\sqrt{h{a}_5}}\\ & +s_0\end{array}$$

  (13)

• Case 3:

  $d_2=d_1=d_4=0$

  $$\begin{array}{rl}{s}_1& =0,d_0=\frac{a_5{d}_{1}h{s}_0{s}_2}{2(a_{5}h{s}_0^2+k)},s_0=\frac{-\sqrt{h}\sqrt{a_1^{2}h+4a_{5}k}-a_{1}h}{2a_{5}h}\\ s_2& =\frac{3(a_3{d}_1{h}^{3/2}(-k)\sqrt{a_1^{2}h+4a_{5}k}+a_2{d}_1{h}^{5/2}\sqrt{a_1^{2}h+4a_{5}k}+a_4{d}_1\sqrt{h}{k}^2\sqrt{a_1^{2}h+4a_{5}k})}{a_1^{2}h+4a_{5}k}\end{array}$$
  then a Weierstrass elliptic solution is obtained:
  $$\begin{array}{rl}u(x,t)=& \frac{3d_1\sqrt{h}(a_2{h}^2+k(a_{4}k-a_{3}h))}{\sqrt{a_1^{2}h+4a_{5}k}\mathrm{\wp }(\frac{1}{2}(hx-k(t+\frac{1}{\mathrm{\Gamma }(\beta )}{)}^{\beta })\sqrt{d_3};}\\ & -\frac{4d_1}{d_3},\frac{6h(a_2{h}^2+k(k{a}_4-h{a}_3))a_5(\sqrt{h{a}_1^2+4k{a}_5}+\sqrt{h}{a}_1)d_1^2}{\sqrt{h{a}_1^2+4k{a}_5}(h{a}_1^2+\sqrt{h}\sqrt{h{a}_1^2+4k{a}_5}{a}_1+4k{a}_5)d_3})\\ & -\frac{\frac{\sqrt{a_1^{2}h+4a_{5}k}}{\sqrt{h}}+a_1}{2a_5}\end{array}$$

  (14)

• Case 4:

  $d_3=d_4=0$

  $$\begin{array}{rl}{s}_0& =\frac{1}{6}(\frac{\sqrt{3}\sqrt{a_5{d}_{2}h{s}_2^2-12d_{0}k}}{\sqrt{a_5}\sqrt{d_0}\sqrt{h}}+\frac{3\sqrt{d_0{d}_2}{s}_2}{d_0}),s_1=0,d_0=-\frac{a_5{s}_2^2}{6(a_2{h}^2-a_{3}hk+a_4{k}^2)}\\ d_2& =\frac{a_1^2(-h)-4a_{5}k}{2a_{5}h(a_2{h}^2-a_{3}hk+a_4{k}^2)}\end{array}$$
  then an exponential solution is obtained:
  $$u(x,t)=\frac{1}{\exp \left(\frac{\sqrt{-\frac{a_1^{2}h+4a_{5}k}{a_{5}h(a_2{h}^2+k(a_{4}k-a_{3}h))}}(hx-k{(\frac{1}{\mathrm{\Gamma }(\beta )}+t)}^{\beta })}{\sqrt{2}}\right)+\frac{a_5}{\sqrt{3}\sqrt{\frac{a_1^{2}h+4a_{5}k}{h}}}}-\frac{\sqrt{3}\sqrt{\frac{a_1^{2}h+4a_{5}k}{h}}}{2a_5}+\frac{a_1}{2a_5}$$

  (15)

• Case 5:

  $d_0=d_1=0,d_4>0,d_3=2\sqrt{d_2{d}_4}$

  $$\begin{array}{rl}{s}_0& =\frac{1}{6}(\frac{\sqrt{3}\sqrt{a_5{d}_{2}h{s}_1^2-12d_{4}k}}{\sqrt{a_5}\sqrt{d_4}\sqrt{h}}+\frac{3\sqrt{d_2{d}_4}{s}_1}{d_4}),s_2=0,d_2=\frac{a_1^2(-h)-4a_{5}k}{2a_{5}h(a_2{h}^2-a_{3}hk+a_4{k}^2)}\\ d_4& =-\frac{a_5{s}_1^2}{6(a_2{h}^2-a_{3}hk+a_4{k}^2)}\end{array}$$
  then a dark solitary solution is obtained:
  $$\begin{array}{rl}& u(x,t)=\\ & \frac{1}{6}(3\sqrt{3}\sqrt{\frac{a_1^{2}h+4a_{5}k}{a_5^{2}h}}(\tanh \left(\frac{\sqrt{-\frac{a_1^{2}h+4a_{5}k}{a_{5}h(a_2{h}^2+k(a_{4}k-a_{3}h))}}(hx-k{(\frac{1}{\mathrm{\Gamma }(\beta )}+t)}^{\beta })}{2\sqrt{2}}\right)+1)\\ & -\frac{3\sqrt{3}\sqrt{\frac{a_1^{2}h+4a_{5}k}{h}}}{a_5}+\frac{3a_1}{a_5})\end{array}$$

  (16)

Analysis and visualization of the derived solution profiles

Figure 1 presents a graphical representation of the computed results, emphasizing the characteristic features and physical behaviors of the selected solution profiles. In particular, it depicts the bright soliton given by equation (12), with the parameter values $h=1.18$ , $s_0=1.53$ , $a_2=0.01$ , $a_3=1.58$ , $a_4=-1.61$ , $a_5=-1.1$ , $k=2.54$ , and $t=13$  s.

Figure 1.

Spatial evolution of the bright soliton solution corresponding to equation (12) for various fractional orders $\beta$ . The soliton maintains a localized structure and stable propagation, demonstrating the balance between dispersion and nonlinearity. (a) $\beta =0.55$ , (b) $\beta =0.7$ , (c) $\beta =1$ , and (d) two-dimensional (2D) profile of the bright soliton at various $\beta$ .

Figure 2 presents a periodic solution corresponding to equation (13), plotted using the parameters $h=1.88$ , $s_0=0.605$ , $a_2=0.01$ , $a_3=0.56$ , $a_4=0.41$ , $a_5=1.1$ , $m=0.762$ , $k=0.72$ , and $t=2.562$  s.

Figure 2.

Periodic wave solution of equation (13) illustrated for different fractional values $\beta$ . The modulation pattern reflects the influence of fractional dynamics on wave periodicity and amplitude. (a) $\beta =0.5$ , (b) $\beta =0.8$ , (c) $\beta =1$ , and (d) two-dimensional (2D) profile of the Jaccobi soliton at various $\beta$ .

Figure 3 shows a dark soliton solution associated with equation (16), under the values $h=-2.74$ , $a_1=-0.82$ , $a_2=-0.08$ , $a_3=-0.78$ , $a_4=-0.64$ , $a_5=-0.2$ , $k=-0.7$ , and $t=10$  s.

Figure 3.

Dark soliton profile obtained from equation (16) for selected $\beta$ values. The reduction in amplitude with decreasing $\beta$ demonstrates the dissipative effect introduced by fractional-order memory. (a) $\beta =0.5$ , (b) $\beta =0.8$ , (c) $\beta =1$ , and (d) two-dimensional (2D) profile of the dark soliton at various $\beta$ .

These wave structures exhibit remarkable stability, maintaining their form and propagation speed across extended spatial domains. Their persistence arises from a delicate interplay between dispersive spreading and nonlinear self-modulation. While dispersion acts to stretch and diffuse the wave profile, nonlinear effects counteract this tendency by inducing a self-compression mechanism. When these two influences are in equilibrium, the resulting soliton travels with minimal deformation or attenuation, preserving its localized structure throughout propagation.

Comparative discussion and novelty justification

To underline the novelty of the obtained results, it is instructive to compare them with previous findings on the Gardner equation and its fractional extensions. For the classical Gardner model,^11–13 analytical solutions were derived only for the integer-order case ( $\beta =1$ ), where soliton amplitudes and widths remain constant. Recent fractional formulations, such as those reported by Ghanbari ¹⁵ and Islam et al., ¹⁴ primarily relied on approximate or numerical approaches and did not yield closed-form analytical expressions. In contrast, the present study provides exact analytical $\beta$ -fractional soliton and periodic solutions using the IMETM, thus extending the available solution classes beyond prior fractional models. Furthermore, when compared with the recent numerical approach of Chen et al., ³⁶ who employed the exponential time-differencing Runge–Kutta (ETDRK) scheme for the fractional Gardner equation, the analytical expressions obtained here offer direct insight into how the fractional order $\beta$ modifies soliton amplitude, width, and stability. This highlights the originality of the present work in bridging analytical and numerical frameworks within fractional nonlinear dynamics.

Linear stability analysis

In this section, we investigate the linear stability of the $\beta$ -time-fractional Gardner equation derived earlier. The purpose of this analysis is to determine the conditions under which small perturbations to the steady-state or soliton solutions remain bounded or grow in time. This analysis follows standard approaches used in the stability study of fractional and nonlinear dispersive systems as discussed in Pakzad et al. ³⁷ and Li et al. ³⁸ These works provide the theoretical foundation of fractional differentiation in dynamic systems with memory effects. The linear stability analysis can be summarized as the following algorithmic sequence: (i) introduce a small perturbation to the steady-state solution; (ii) linearize the governing fractional PDE by neglecting higher-order perturbation terms; (iii) assume harmonic normal modes of the form $\omega \sim e^{i(kx-\omega t)}$ ; (iv) substitute the Fourier representations into the linearized equation to derive the dispersion relation; and (v) analyze the real and imaginary parts of the complex frequency $\omega$ to determine stability. This procedure is shown in the flowchart in Figure 4.

Figure 4.

Flowchart of the linear stability analysis procedure for the $\beta$ -time-fractional Gardner equation.

Linearized model

Starting from the governing model

$$D_t^{\beta }q+(a_{1}q+a_5{q}^2)q_x+a_2{q}_{xxx}+a_3{q}_{xxt}+a_4{q}_{xtt}=0,0<\beta \le 1$$

(17)

we introduce a small perturbation around the steady background $q_0$ :

$$u(x,t)=u_0+\epsilon \varphi (x,t),0<\epsilon \ll 1$$

(18)

Substituting into (17) and neglecting nonlinear terms of $\mathcal{O}({\epsilon }^2)$ yields the linearized fractional evolution equation:

$$D_t^{\beta }\varphi +(a_1+2a_5{q}_0){\varphi }_x+a_2{\varphi }_{xxx}+a_3{\varphi }_{xxt}+a_4{\varphi }_{xtt}=0$$

(19)

Normal mode analysis

We assume a harmonic perturbation of the form:

$$\varphi (x,t)=\exp [i(kx-\omega t)]$$

(20)

where $k$ is the wave number and $\omega$ is the complex frequency. This harmonic perturbation approach has been successfully used in several recent studies to analyze fractional and nonlinear dispersive equations.^39,40 Using the Fourier representations of derivatives:

$${\varphi }_x=ik\varphi ,{\varphi }_{xxx}=-i{k}^3\varphi ,{\varphi }_{xxt}=i{k}^2\omega \varphi ,{\varphi }_{xtt}=-ik{\omega }^2\varphi ,D_t^{\beta }\varphi =(-i\omega {)}^{\beta }\varphi$$

the substitution of these expressions into equation (19) leads to the general dispersion relation:

$$(-i\omega {)}^{\beta }+ik(a_1+2a_5{q}_0)-i{a}_2{k}^3+i{a}_3{k}^2\omega -i{a}_{4}k{\omega }^2=0$$

(21)

Classical limit $(\beta =1)$

For the classical time derivative, $D_t^1=\mathrm{d}/\mathrm{d}t$ , equation (21) reduces to

$$a_{4}k{\omega }^2-(1-a_3{k}^2)\omega +(a_2{k}^3-k(a_1+2a_5{q}_0))=0$$

(22)

Equation (22) is a quadratic relation in $\omega$ , whose explicit solution is given by

$$\omega =\frac{(1-a_3{k}^2)\pm \sqrt{(1-a_3{k}^2{)}^2-4a_{4}k(a_2{k}^3-k(a_1+2a_5{q}_0))}}{2a_{4}k}$$

(23)

The discriminant

$$\mathrm{\Delta }=(1-a_3{k}^2{)}^2-4a_{4}k(a_2{k}^3-k(a_1+2a_5{q}_0))$$

governs the system’s response:

• $\mathrm{\Delta }>0$ : $\omega$ is real $\Rightarrow$ neutrally stable wave propagation.

• $\mathrm{\Delta }<0$ : $\omega$ is complex $\Rightarrow$ exponential growth or decay (instability).

Hence, the sign of $\mathrm{Im}(\omega )$ determines the linear stability of the base state:

$$\mathrm{Im}(\omega )<0\Rightarrow \text{stable},\mathrm{Im}(\omega )>0\Rightarrow \text{unstable}$$

Fractional case $(0<\beta <1)$

When $0<\beta <1$ , the fractional derivative introduces a memory-dependent term:

$$(-i\omega {)}^{\beta }={\omega }^{\beta }{e}^{-i\pi \beta /2}$$

Substituting this into equation (21) yields the fractional dispersion relation:

$${\omega }^{\beta }=i{e}^{i\pi \beta /2}[a_{4}k{\omega }^2-a_3{k}^2\omega +(a_2{k}^3-k(a_1+2a_5{q}_0))]$$

(24)

The complex exponential $e^{-i\pi \beta /2}$ introduces a phase lag between temporal evolution and forcing, modifying both amplitude and frequency of perturbations. For $\beta <1$ , perturbations decay algebraically in time as $t^{-\beta }$ instead of exponentially, reflecting long-memory damping effects inherent to fractional dynamics.

The effective growth rate can be estimated as

$${\mathrm{\Gamma }}_{\text{eff}}\propto |\mathrm{Im}(\omega ){|}^{\beta }\cos \left(\frac{\pi \beta }{2}\right)$$

indicating that smaller $\beta$ values suppress the instability growth rate, enhancing temporal stability.

Summary of stability conditions

See Table 1 for details.

Table 1.

Summary of dispersion relations and corresponding stability control parameters.

Case	Dispersion relation	Dominant stability control
$\beta =1$ (classical)	$a_{4}k{\omega }^2-(1-a_3{k}^2)\omega +(a_2{k}^3-k(a_1+2a_5{q}_0))=0$	Sign of discriminant $\mathrm{\Delta }$
$0<\beta <1$ (fractional)	${\omega }^{\beta }=i{e}^{i\pi \beta /2}[a_{4}k{\omega }^2-a_3{k}^2\omega +(a_2{k}^3-k(a_1+2a_5{q}_0))]$	Fractional damping and phase shift

Graphical representation of stability characteristics

Figure 5 illustrates the dispersion and growth-rate behavior of the $\beta$ -time-fractional Gardner equation for several fractional orders $\beta =1.0,0.8,0.6,0.4$ using the parameter set $a_1=-0.82$ , $a_2=-0.08$ , $a_3=-0.78$ , $a_4=-0.64$ , $a_5=-0.2$ , and $q_0=1$ . The left panel shows the variation of the real part of the complex frequency, $\mathrm{Re}(\omega )$ , with respect to the wavenumber $k$ , which determines the phase velocity and dispersion characteristics of the propagating modes. The right panel presents the imaginary component, $\mathrm{Im}(\omega )$ , which quantifies the temporal growth or attenuation rate of perturbations.

Figure 5.

Dispersion and growth-rate profiles of the $\beta$ -time-fractional Gardner equation for several fractional orders $\beta$ . (Left) Real part of $\omega$ representing the dispersion relation. (Right) Imaginary part of $\omega$ illustrating the perturbation growth rate. A smaller fractional order $\beta$ reduces $\mathrm{Im}(\omega )$ , indicating enhanced temporal stability due to fractional damping and long-memory effects.

It is evident that the real part $\mathrm{Re}(\omega )$ remains negative for all $\beta$ values, indicating forward-propagating dispersive waves. The magnitude of $\mathrm{Re}(\omega )$ slightly decreases with smaller $\beta$ , revealing that fractional-order dynamics reduce the effective phase speed due to memory effects in time.

More importantly, the imaginary part $\mathrm{Im}(\omega )$ , which governs the linear growth rate, shows a clear stabilizing trend as $\beta$ decreases. For the classical case $\beta =1$ , perturbations exhibit the largest growth rate. When $\beta$ is reduced to 0.8, 0.6, and 0.4, the growth rate progressively declines, implying enhanced damping of small disturbances. This behavior confirms that the fractional time derivative introduces a memory effect that dissipates energy more efficiently, thereby stabilizing the nonlinear wave structures.

Conclusion

Recently, several studies have explored fractional generalizations of the Gardner equation to capture memory effects in nonlinear dispersive media. Numerical approaches have been particularly effective for such models. For instance, Chen et al. ³⁶ applied the ETDRK method to investigate nonlinear wave dynamics in the fractional Gardner equation, demonstrating high numerical stability and accuracy. Their results provide valuable computational support for analytical developments of fractional-order Gardner-type systems. While the ETDRK method in Chen et al. ³⁶ provides efficient numerical approximations, the present study focuses on developing closed-form analytical solutions using a fractional $\beta$ -derivative framework, allowing a deeper insight into the qualitative structure and stability of the obtained waves. In this study, we have successfully applied the IMETF method to obtain exact analytical solutions for the dimensionless $\beta$ -time-fractional Gardner equation, which models shallow water wave dynamics. A variety of solutions were derived, including bright and dark solitons, Jacobi elliptic periodic solutions, Weierstrass elliptic functions, exponential solutions, and singular solutions. These results demonstrate the effectiveness and versatility of the IMETF method in handling nonlinear fractional PDEs. The obtained fractional soliton and periodic solutions extend earlier integer-order results ¹¹ by incorporating the $\beta$ -time derivative, which introduces memory-dependent modulation. The stability findings further show that the fractional order enhances temporal damping compared to the classical Gardner equation, offering a new physical interpretation of fractional wave dynamics.

The inclusion of the fractional derivative of order $\beta$ ( $0<\beta \le 1$ ) significantly influences the behavior of the obtained wave solutions, affecting their amplitude, width, and periodicity. This underscores the importance of fractional calculus in capturing memory effects and nonlocal temporal dynamics in physical systems, offering a more realistic modeling approach compared to classical integer-order models.

Graphical representations of selected solutions further illustrate the propagation characteristics and stability of these waves under various parameter conditions. The balance between nonlinearity and dispersion in maintaining soliton structures is clearly evident in the results.

The validity of the obtained analytical solutions was confirmed by direct substitution into the original fractional Gardner equation, ensuring that each satisfies the governing model identically. Moreover, numerical simulations and three-dimensional surface plots were generated for different values of $\beta$ , showing excellent agreement with the analytical results. This cross-verification provides strong confidence in the correctness, stability, and physical consistency of the derived solutions.

This work not only contributes to the theoretical understanding of fractional nonlinear wave equations but also provides a framework for future studies on more complex systems. The proposed method can be extended to other nonlinear fractional models in fluid dynamics, plasma physics, and optical communications.