A Melnikov analysis and chaotic dynamics of low-frequency dissipative ion-acoustic solitons in multicomponent plasmas

Abstract

This study investigates the nonlinear dynamics of low-frequency ion-acoustic waves in a four-component unmagnetized plasma consisting of warm adiabatic ions, κ-distributed superthermal electrons and positrons, and isothermal electrons. Using a reductive perturbation approach, we derive a damped, forced modified Korteweg-de Vries-Burgers equation to describe the evolution of compressive ion-acoustic solitary waves (IASWs). The properties of these IASWs are significantly influenced by key plasma parameters, including the temperature ratio between superthermal electrons and positrons, the density ratio of superthermal to isothermal electrons, the temperature ratio of warm ions to superthermal electrons, and the superthermal index. Further, we apply Melnikov’s analytical method to identify criteria for the onset of chaos when an external periodic force is introduced. Supported by numerical simulations including time series, phase portraits, and Poincaré sections, our analysis reveals that the system can exhibit regular periodic oscillations, quasi-periodic regimes, and chaotic dynamics depending on the strength and frequency of the external perturbation. These results offer new insights into the transition to chaos in dissipative plasmas and may help interpret nonlinear wave phenomena in environments where collisions and external forcing are present.

I. Introduction

Recently, several satellite missions and measurements^1,2 have revealed that electrons and positrons coexist in an electron-positron-ion (e-p-i) plasma in both space (like the early universe,^3–5 Saturn’s magnetosphere,^6,7 solar wind,^8–11 neutron star polar regions, ¹² and active galactic nuclei,^13,14 etc.) and Laboratory environments (like semiconductor plasmas, ¹³ laser irradiation, ³ magnetic confinement systems, ¹³ and hot cathode discharge ³ ). To comprehend the fundamental characteristics of e-p-i plasmas, numerous authors have investigated ion-acoustic waves (IAWs),^15–17 electron-acoustic waves (EAWs), and positron-acoustic waves (PAWs), ⁹ along with their associated nonlinear structures, including solitons,^18–20 rogue waves, ⁸ shocks, and double layers. Specifically, Mouhammadoul et al. ⁸ studied PAWs in e-p-i magnetoplasmas, while Paul et al.^21,22 examined the stability of these waves.

Observations from Cassini CAPS2 ²³ and Voyager PLS2 ²⁴ in Saturn’s magnetosphere revealed the existence of two-temperature electrons (TTE) (cold and hot), a finding later confirmed by multiple satellite missions such as FAST Auroral Snapshot. ²⁵ These electrons follow a superthermal κ-distribution rather than the conventional Maxwellian distribution, a fact that numerous authors have incorporated in their investigations of nonlinear electrostatic wave propagation in space plasmas.^26–28 The superthermality of plasma species is characterized by the superthermal parameter in the κ-distribution, where small κ-values indicate a substantial deviation from thermal equilibrium, while large κ-values approach a Maxwellian distribution. Alinejad et al. ⁶ studied ion-acoustic (IA) solitary waves in a plasma model with cold ions and TTE, confirming the existence of both compressive and rarefactive solitary structures. Similar investigations were conducted by other researchers,^29–31 who analyzed IA solitons theoretically and numerically in the presence of TTE. In a plasma system containing inertial cold ions and inertialess TTE, Panwar et al. ³ established IA cnoidal waves and discovered that the superthermality of a cold electron increases the amplitude of the cnoidal wave (CW). Sing et al. studied ion thermal effects on IAWs in magnetized superthermal plasma via the Sagdeev potential method, reporting good agreement with Viking satellite observations in the auroral region. ³² Subsequent investigations have demonstrated that superthermal effects significantly increase the modulational instability of electrostatic wave packets in e-p-i plasmas.^33,34

Recently, the study of the chaotic behavior of (un)magnetized plasmas has attracted considerable interest.^7,35,36 The dynamical features of dust IAWs (DIAWs) in dusty plasma were investigated by Saha et al. ³⁷ Employing the bifurcation theory of planar dynamical systems, they derived the soliton and periodic solutions of these waves in the absence of a periodic external perturbation. However, when a periodic external perturbation is introduced into the plasma system, DIAWs have also been shown to exhibit quasi-periodic and chaotic behavior. Saha et al. ³⁸ looked at the presence of homoclinic and periodic orbits as well as the nonlinear behavior of IAWs. They also showed that numerous quasi-periodic and chaotic features can be produced by adding a periodic disturbance to this model. The coexistence of periodic, quasi-periodic, and chaotic behaviors for the IAWs has been obtained for suitable parameters. ³⁹ Pradhan et al. ⁴⁰ constructed the modified Korteweg-de Vries (mKdV) equation for EAWs in quantum plasma and used the concept of phase plane analysis to get a variety of solutions, including periodic and super-periodic waves. Moreover, they illustrated the existence of multiple equilibrium points, which are used to examine different solutions. Additionally, they investigate how the stationary mKdV equation affects the emergence of quasi-periodic, super-periodic, and chaotic regimes.

External perturbations have an impact on plasma environments such as planetary magnetospheres or turbulent solar wind.^41–43 It is evident that a flexible, fast waveform generator may provide such an external force. Sen et al. ³⁵ derived the forced Korteweg-de Vries (fKdV) equation using a sophisticated plasma model. They showed that one of the primary causes of the density issue was the external periodic force term. A solitary wave under the effect of a periodic force was observed by Ali et al. ³⁶ The fKdV equation was recently developed by Ghorui et al., ⁷ who obtained analytical solutions in the form of solitary waves in quantum plasma with the presence of external disturbance. They also examined how the different physical plasma parameters affect the propagation of nonlinear waves. On the other hand, extensive theoretical investigations have been conducted to investigate dissipative soliton dynamics in collisional plasmas. Alinejad et al. examined the evolution of dissipative DIAWs in a dusty plasma and found that the associated dissipative rogue wave profile is extremely sensitive to changes in the dissipation coefficient. ⁴⁴ El-Tantawy analyzed the nonlinear dynamics of soliton collisions in electronegative plasma and reported that variations in plasma parameters markedly influence the phase shifts and trajectories of the colliding solitons. ⁴⁵ Albawi et al. studied dissipative IAWs in anisotropic rotating magnetoplasmas, showing that higher collision frequencies produce shorter and broader dissipative solitons.^46,47 Despite extensive experimental and theoretical research, many fundamental issues regarding dissipative IA solitary waves (IASWs) at a crucial parameter combination involving the collisional frequency and an external periodic force remain unresolved.

Numerous numerical and analytical methods are now used to investigate the characteristics of nonlinear complex systems.^48–52 For example, the multiple-scale approach is recommended for studying nonlinear resonances, and the averaging method is typically used to determine the amplitude response of self-excited systems, such as the van der Pol oscillator. Strong numerical tools like phase portraits, Poincaré sections, and the Lyapunov exponent in conjunction with the bifurcation diagram are typically employed to investigate the chaotic behavior of nonlinear systems.5 ⁵³ Despite these advances, there are not many analytical markers available to identify or forecast chaos in nonlinear systems. One of the effective methods used recently to identify chaos is the Melnikov criterion.^54–59 The Melnikov approach has been used by several authors^60–63 to determine the critical levels at which solitonic trajectories become sensitive to disturbances. These studies show that Melnikov analysis provides a bridge between theory and experimental diagnostics by directly relating the structure of differential equations to the observed chaotic dynamics. Recent research has confirmed the applicability of Melnikov’s method for diagnosing chaos in a wide range of dynamical systems, from vibration impact oscillators subject to noise ⁶⁴ to fractional systems with nonlinear amortization. ⁶⁵ These studies demonstrate that the Melnikov criterion remains a universal tool for characterizing transitions to chaotic regimes, encompassing a broad class of plasma systems. Our method is embedded in this dynamics, specifically applying to nonlinear structures in superthermal and dissipative plasmas. The requirement is to examine the unperturbed Hamiltonian system with a regularly periodic invariant set where stable and unstable manifolds may intersect. The Melnikov theory has also been employed in numerous nonlinear systems, including biology, epidemiology, and engineering.^66,67

To date, no research has been published that demonstrates how Melnikov’s theory of chaotic behavior and ion-neutral collisions affects IASWs within the context of the damped forced modified Korteweg-de Vries-Burgers (DFMKdVB) equation. We are moving forward with this work as a result. Thus, our aim to study the impact of the collision frequency, the ion kinematic viscosity, and the external periodic force on the properties of IASWs. On the other hand, we are going to perform the Melnikov criteria to identify the region for which the regular and irregular state appear with respect to the strength and the frequency of the external periodic perturbation. Using the indicators of chaos, the periodic, quasi-periodic, and chaotic patterns of IAWs within the context of the DFMKdVB equation are shown.

This work is structured as follows: The model equations are given in Section II. In Section III, we derive the DFMKdVB equation with its corresponding soliton solution. In Section IV, we investigate the dynamical system associated with the DFMKdVB equation, the bifurcation theory, and the Melnikov analysis. Section V analyzes the quasi-periodic and chaotic motions of the perturbed system. Section VI summarizes the key findings of our investigation.

II. Basic fluid equations

Warm adiabatic ions (charge Z₊ e, mass m₊), κ-distributed superthermal electrons (charge −e, mass $m_{e_1}$ ), isothermal electrons (charge −e, mass $m_{e_2}$ ), and κ-distributed superthermal positrons (charge +e, mass m _p ) make up the four constituents of an unmagnetized, collisionally dissipative plasma model that we examine. The system maintains quasi-neutrality in equilibrium, expressed by the relation ϑ₁ + ϑ₂ − ϑ₃ = 1, where the dimensionless density ratios are defined as ϑ₁ = n_e10/Z₊ n₊₀, ϑ₂ = n_e20/Z₊ n₊₀, and ϑ₃ = n_p0/Z₊ n₊₀. Here, n₊₀, n_e10, n_e20, and n_p0 stand for the unperturbed number densities of ions, superthermal electrons, isothermal electrons, and positrons, respectively. We modify the fundamental normalized fluid equations by adding ion viscosity η₊, ion-neutral collisional damping μ _i , and a spatially and temporally source term S (x, t) in Poisson’s equation to account for realistic dissipative events and external field contributions. The source term may originate from various sources, such as a flexible high-speed waveform generator, ⁶⁸ or from physical processes including photoemission and secondary electron emission induced by debris objects (e.g., broken solar panels, metals fragments, or insulation materials) traversing Earth’s ionosphere. ⁶⁹ These extra terms are driven by astrophysical and laboratory settings where external charge disturbances and frictional dissipation are not negligible ⁷⁰ :

$\begin{array}{ll}\hfill {\partial }_t{n}_{+}+{\partial }_x\left(n_{+}{u}_{+}\right)& =0,\hfill \end{array}$

(1a)

$\begin{array}{ll}\hfill {\partial }_t{u}_{+}+u_{+}{\partial }_x{u}_{+}+\varpi n_{+}{\partial }_x{n}_{+}& =-{\mu }_i{u}_{+}+{\eta }_{+}{\partial }_x^2{u}_{+}-{\partial }_x\varphi ,\hfill \end{array}$

(1b)

$\begin{array}{ll}\hfill {\partial }_x^2\varphi & ={\vartheta }_1{n}_{e_1}+{\vartheta }_2{n}_{e_2}-{\vartheta }_3{n}_p-n_{+}+S.\hfill \end{array}$

(1c)

Here, n₊ stands for the number density of warm adiabatic ions, normalized by the unperturbed ions densities n₊₀, whereas u₊₀ represents their corresponding speed, normalized by the ion sound speed $C_{+}=\sqrt{Z_{+}{k}_B{T}_{e_1}/m_{+}}$ . The electrostatic wave potential ϕ is normalized by the quantity $k_B{T}_{e_1}/e$ . The ratio of warm ion temperature to superthermal electron temperature is defined by $\varpi =3T_{+}/Z_{+}{T}_{e_1}$ . The expression of ion pressure for one dimension is $P_{+}=\mathrm{\Xi }{(N_{+}/n_{+0})}^2$ , where Ξ = k _B N₊ T₊. The Debye length (inverse plasma frequency) scales the space (time) variable ${\lambda }_{D+}=\sqrt{k_B{T}_{e_1}/4\pi e^2{Z}_{+}{n}_{+0}}$ $\left({\omega }_{p+}=\sqrt{4\pi e^2{Z}_{+}^2{n}_{+0}/m_{+}}\right)$ .

The number density equations for superthermal electrons $(n_{e_1})$ , isothermal electrons $(n_{e_2})$ and superthermal positrons (n _p ) are as follows^71,72:

$\begin{array}{ll}\hfill n_{e_1}& ={\left(1-\frac{2\varphi }{2\kappa -3}\right)}^{-\kappa +1/2},\hfill \end{array}$

(2a)

$\begin{array}{ll}\hfill n_{e_2}& ={\exp }^{\varrho \varphi },\hfill \end{array}$

(2b)

$\begin{array}{ll}\hfill n_p& ={\left(1+\frac{2\rho \varphi }{2\kappa -3}\right)}^{-\kappa +1/2},\hfill \end{array}$

(2c)

where the temperature ratio between superthermal electrons and positron is given by $\rho =T_{e_1}/T_p$ . The κ-distribution requires κ > 3/2 to ensure physical meaningfulness, where κ → 3/2 indicate a significant non-Maxwellian tail, while κ → ∞ recovers the Maxwellian limit. Applying Taylor series expansion to Eqs. (2a)–(2c) and substituting these expressions into Eq. (1c), yields the following modified Poisson equation:

$\begin{array}{c}{\partial }_x^2\varphi ={\vartheta }_1+{\vartheta }_2-{\vartheta }_3+S-n_{+}+{\alpha }_1\varphi +{\alpha }_2{\varphi }^2+{\alpha }_3{\varphi }^3,\hfill \end{array}$

(3)

where the coefficients α₁, α₂, and α₃ incorporate all relevant plasma parameters through the following expressions:

$\begin{array}{ll}\hfill {\alpha }_1& =\frac{{\vartheta }_2\varrho (2\kappa -3)+(2\kappa -1)({\vartheta }_1+{\vartheta }_3\rho )}{2\kappa -3},\hfill \end{array}$

(4a)

$\begin{array}{ll}\hfill {\alpha }_2& =\frac{{\vartheta }_2{\varrho }^2{(2\kappa -3)}^2+(4\kappa -1)({\vartheta }_1-{\vartheta }_3{\rho }^2)}{2{(2\kappa -3)}^2},\hfill \end{array}$

(4b)

$\begin{array}{ll}\hfill {\alpha }_3& =\frac{{\vartheta }_2{\varrho }^3{(2\kappa -3)}^3+(4\kappa -1)(2\kappa +3)({\vartheta }_1+{\vartheta }_3{\rho }^3)}{6{(2\kappa -3)}^3}.\hfill \end{array}$

(4c)

III. Derivation of DFMKdVB equation and its solution

Let’s use the reductive perturbation method through the stretched coordinates to investigate the IASWs propagating in the unmagnetized plasma. In order to explain the hierarchy of elements supporting the nonlinear dispersive dynamics, we introduce specific coordinates that combine short variations along the propagation direction and even shorter modulations in the transverse directions as: ξ = ϵ (l₁ x + l₂ y + l₃ z − v _p t), τ = ϵ³ t. The small parameter ϵ ≪ 1 represents the weak nonlinearity and small amplitude of the wave. The coefficients l₁, l₂, l₃ account for oblique propagation in an unmagnetized plasma. Since the system is isotropic, the direction vector (l₁, l₂, l₃) satisfies $l_1^2+l_2^2+l_3^2=1$ , ensuring that the derivation is general and not restricted to one-dimensional propagation along a coordinate axis. This choice guarantees the existence of coherent stable structures of nonlinear evolution equations, as established in relativistically enhanced environments and in magnetized plasmas, where the compensation of nonlinearity, dispersion, and transverse diffraction precisely requires this separation of electrons.^73,74 Additionally, this scaling makes it possible to analyze transversal stability without altering the longitudinal equilibrium’s order, making comparisons with reconstructions based on physical learning and recent numerical diagnostics in multi-species magnetized environments easier.^75,76 The reductive perturbation approach extends the independent variables in the following ways ⁷⁷ :

$\begin{array}{ll}\hfill n_{+}& =1+ϵn_{+}^{(1)}+{ϵ}^2{n}_{+}^{(2)}+{ϵ}^3{n}_{+}^{(3)}+\dots ,\hfill \end{array}$

(5a)

$\begin{array}{ll}\hfill u_{+}& =ϵu_{+}^{(1)}+{ϵ}^2{u}_{+}^{(2)}+{ϵ}^3{u}_{+}^{(3)}+\dots ,\hfill \end{array}$

(5b)

$\begin{array}{ll}\hfill \varphi & =ϵ{\varphi }_1+{ϵ}^2{\varphi }_2+{ϵ}^3{\varphi }_3+\dots ,\hfill \end{array}$

(5c)

$\begin{array}{ll}\hfill S& ={ϵ}^3{S}_3,\hfill \end{array}$

(5d)

$\begin{array}{ll}\hfill {\eta }_{+}& =ϵ\eta .\hfill \end{array}$

(5e)

The solutions of first order ϵ, read:

$\begin{array}{c}{n}_{+}^{(1)}={\varphi }_1{\alpha }_1,u_{+}^{(1)}=v_p{\mathrm{\Phi }}_1{\alpha }_1,\hfill \end{array}$

(6)

and the phase velocity expression is as follow: $v_p^2=\varpi {\alpha }_1+1/{\alpha }_1$ .

Additionally, following a few algebraic steps, the second-order perturbed values, $n_{+}^{(2)}$ , $u_{+}^{(2)}$ , have been expressed in the following manner for the next-order in ϵ:

$\begin{array}{c}{n}_{+}^{(2)}={{\varphi }_1}^2{\alpha }_2+{\varphi }_2{\alpha }_1,u_{+}^{(2)}=\left(-{{\alpha }_1}^2{v}_p+{\alpha }_2{v}_p\right){{\varphi }_1}^2+v_p{\alpha }_1{\varphi }_2.\hfill \end{array}$

(7)

Moreover, collecting the next order series of ϵ, we obtain:

$\begin{array}{c}{n}_{+}^{(3)}=2{\varphi }_1{\varphi }_2{\alpha }_2+{\varphi }_3{\alpha }_1-\frac{{\partial }^2}{\partial {\xi }^2}{\varphi }_1-S_3.\hfill \end{array}$

(8)

The DFMKdVB equation is obtained by substituting Eqs. (6)–(8) in Poisson’s equation and after algebraic manipulation, we have:

$\frac{\partial }{\partial \tau }\varphi +A{\varphi }^2\frac{\partial }{\partial \xi }\varphi +B\frac{{\partial }^3}{\partial {\xi }^3}\varphi +C\frac{{\partial }^2}{\partial {\xi }^2}\varphi +D\varphi =B\frac{\partial }{\partial \xi }{S}_3,$

(9)

where the following is the expression for this equation’s coefficients:

$\begin{array}{c}A=\frac{-6{\alpha }_1^3{v}_p^2+9{\alpha }_1{\alpha }_2{v}_p^2+3\varpi {\alpha }_1{\alpha }_2}{2v_p{\alpha }_1},B=\frac{v_p^2-\varpi }{2v_p{\alpha }_1},C=-\frac{\eta }{2},D=\frac{{\mu }_i}{2}.\hfill \end{array}$

(10)

Recent studies have demonstrated that external periodic forces can significantly modify the behavior of nonlinear waves in plasma systems. Various forcing terms have been proposed to model different physical scenarios, including trigonometric, hyperbolic, and Gaussian functions, which may represent effects from space debris or other perturbations. ⁷⁸ Motivated by these findings, we consider a spatially linear and temporally periodic forcing term of the form S₃ = f₀ ξ/B cos(ω τ), where f₀ represents the forcing amplitude, ω is the modulation frequency, and τ denotes the normalized time variable. Accordingly, we obtain:

$\frac{\partial }{\partial \tau }\varphi +A{\varphi }^2\frac{\partial }{\partial \xi }\varphi +B\frac{{\partial }^3}{\partial {\xi }^3}\varphi +C\frac{{\partial }^2}{\partial {\xi }^2}\varphi +D\varphi =f_{0}cos(\omega \tau ),$

(11)

In general, when the Burgers term C is added to the Korteweg-de Vries-Burgers (KdVB) equation, the system of equations becomes dissipative and the total energy is not conservative. As a result, it becomes more challenging to construct the exact solution of the KdVB problem. Numerous computational techniques have been created to address different types of nonlinear differential equations. We refer to the equation above as the DFMKdVB equation. The equation (11), identified as the mKdV equation, admits a solution in a frame moving with phase speed M₀. By neglecting the impact of viscosity, external force, and collision, the corresponding soliton solution reads^78,79:

$\varphi ={\varphi }_{0}sech\left(\frac{\xi -M_0\tau }{W_0}\right),$

(12)

where, ${\varphi }_0=\sqrt{6M_0/A}$ , and $W_0=\sqrt{B/M_0}$ , stand for the amplitude and the width of the wave respectively. Here, M₀ represents the constant phase speed of the IASW in the absence of damping and external forcing. In our present study, we assumed for our numerical simulation that $T_{e_1}=T_p=(10\sim 1000)T_{e_2}$ , $T_{+}=0.1T_{e_2}$ , Z = 1 − 20.^70–72

The profile of the IASWs solution to the mKdV equation is displayed in Figure 1 for the relevant plasma parameters. Panel 1(a) illustrates the impact of the ratio between the temperature of superthermal electrons and that of superthermal positrons. It is evident that the amplitude of the IASWs grows by increasing the value of ρ. This is because hotter electrons increase the wave’s compressive strength by making a greater contribution to the pressure balance and to nonlinear steepening. According to the physical interpretation, the Sagdeev pseudo-potential framework favors higher-amplitude solitons because electron pressure dominates, creating deeper potential wells. Both relativistic and non-relativistic plasma models have demonstrated this behavior. ⁸⁰ In panel 1(b), the behavior of the IASWs with respect to the ratio of superthermal electron density to warm ion density is observed. It can be observed that the amplitude of the IASWs increases, while their width decreases, as this density ratio increases. It can be inferred that the plasma system’s nonlinearity increases with the ratio of warm-ion density to superthermal-electron density. Stronger electrostatic potential wells support higher-amplitude IASWs. However, the wave profile tends to become sharper due to the increased nonlinearity, which results in lower soliton widths. This makes the plasma more sensitive to disturbances, enabling energy to focus into more intense, tightly isolated formations. The influence of the temperature ratio between warm ions and superthermal electrons is shown in panel 1(c). From this panel, one can see that the amplitude of the IASWs decreases as ϖ increases. This reveals that the ions are becoming increasingly energetic than the electrons as the temperature ratio between warm ions and superthermal electrons increases. In IASWs, ions provide the inertia, while electrons usually support the pressure. The impact of the superthermal parameter on IASWs is shown in panel 1(d). This panel shows that the wave’s amplitude decreases while its width increases slowly as the superthermal parameter is increased. This indicates that the dispersive effects are enhanced and the plasma medium’s nonlinearity is reduced when both superthermal electrons and positrons are present. Because the increasing population of energetic particles smoothes out potential gradients and lowers the electrostatic trapping required for sharp solitary structures, this dual superthermality results in IASWs that are broader and have lower amplitudes.

Figure 1.

The profile of IASWs (ϕ) obtained from solution provided in Eq. (12), for (a) ϑ₁ = 0.6, ϑ₂ = 0.3, κ = 2, ϖ = 0.05, for (b) ϑ₂ = 0.3, κ = 2, ρ = 1.1, ϖ = 0.05, for (c) ϑ₁ = 0.6, ϑ₂ = 0.3, κ = 2, ρ = 1.1, for (d) ϑ₁ = 0.6, ϑ₂ = 0.3, ρ = 1.1, ϖ = 0.05. Along with M₀ = 0.1.

The soliton solution which depends on the viscosity, collision and forced term is obtained in the form^78,79:

${\varphi }_1={\varphi }_0(\tau )sech\left(\frac{\xi -M(\tau )\tau }{W(\tau )}\right),$

(13)

where ${\varphi }_0(\tau )=\sqrt{6M(\tau )/A}$ , and $W(\tau )=\sqrt{B/M(\tau )}$ , stand for the amplitude and the width of the wave respectively. The expression of the time-dependent velocity is [details in Appendix A]:

$M(\tau )=\frac{3BD{M}_0{e}^{4D\tau }}{M_{0}C(e^{4D\tau }-1)+4D}+\frac{3A{f}_0\left(\frac{4}{3}Ccos(\omega \tau )+\omega sin(\omega \tau )\right)}{4\pi (2D^2+{\omega }^2)}.$

(14)

We now present our numerical analysis of how key physical parameters influence IASW solutions of the DMFKdVB equation (11). Figure 2(a) demonstrates the impact of collision frequency μ _i on wave propagation, holding all other parameters constant as in Figure 1. The results reveal an apparent damping effect: as μ _i increases, we observe a systematic reduction in solitary wave amplitude while maintaining the overall wave structure. This damping phenomenon arises from multiple collisions between plasma particles, which enhance energy dissipation in the plasma model. Similarly, considering the same parameters as in Figure 1, it is seen in Figure 2(b) that the amplitude of the IASW decreases with the enhancement of the viscosity coefficient η. In Figure 2(c), it is shown that as the force f₀ grows, the IASWs’ amplitude increases and their width diminishes. Figure 2(d) shows the fluctuation of the solitary wave solution for IASWs associated with the DMFKdVB equation (11) for various frequencies ω of the external periodic force. As the frequency of the external periodic force rises, both the IASWs’ amplitude and width decrease.

Figure 2.

The profile of IASWs (ϕ₁) obtained from solution provided in Eq. (13), for (a) f₀ = 1.5, ω = 0.5, η = 0.1, for (b) f₀ = 1.5, ω = 0.5, μ _i = 0.1, for (c) ω = 0.5, η = 0.1, μ _i = 0.1, for (d) f₀ = 1.5, η = 0.1, μ _i = 0.1. Along with ϑ₁ = 0.6, ϑ₂ = 0.3, κ = 2, ρ = 1.1, ϖ = 0.05, and τ = 10.

The fluctuation of the width and amplitude of IASWs corresponding to the DMFKdVB equation (11) due to the impact of the external periodic force is shown in Figure 3. It is observed from Figure 3(a) that, increasing the value of f₀ leads to diminish the width of IASWs, whereas in Figure 3(b) the width of IASWs is an increasing function of the frequency of the periodic force. Moreover, Figure 3(c) displays the dependence of the amplitude of IASWs solution of the DMFKdVB (11) on the strength f₀ of the periodic force. It is noticed that when the strength f₀ of the periodic force increased, the amplitude of IASWs increased. Figure 3(d) shows the variation of the amplitude of IASWs with respect to the frequency ω. It is seen that the amplitude is a decreasing function of the frequency of the external periodic force. In general, as f₀ increases, the width decreases over time, while the amplitude increases. This result indicates that a stronger external force leads to greater energy injection into the system. In contrast, higher forcing frequency ω reduces the soliton amplitude and increases its width, indicating that rapid oscillations of the external drive tend to suppress nonlinearity and enhance dispersive effects, producing smoother, weaker soliton structures. Moreover, for small times, the width (amplitude) decreases (increases) with time, but beyond a certain time, the opposite occurs. Overall, the figure highlights that external forcing acts as a control parameter for IASWs: strong, low-frequency forcing leads to sharp, localized solitons, while high-frequency forcing produces weaker, smoother structures.

Figure 3.

Evolution of the width and amplitude of IASWs against τ for various values of f₀ and ω. For (a) and (c) ω = 0.7, for (b) and (d) f₀ = 0.5. Along with ϑ₁ = 0.6, ϑ₂ = 0.3, κ = 2, ρ = 1.1, ϖ = 0.05, η = 0.1, and μ _i = 0.1.

IV. Dynamical system of DMFKdV-Burger equation

We convert the DMFKdVB Eq. (9) to the system of traveling waves by using these transformations χ = (ξ − u₀ τ), ϕ(ξ, τ), S₃(ξ, τ) = S₃(χ), and assuming that D = 0 thus, the DMFKdVB Eq. (9) becomes:

$-u_0\frac{\mathrm{d}}{\mathrm{d}\chi }\varphi +A{\varphi }^2\frac{\mathrm{d}}{\mathrm{d}\chi }\varphi +B\frac{{\mathrm{d}}^3}{\mathrm{d}{\chi }^3}\varphi +C\frac{{\mathrm{d}}^2}{\mathrm{d}{\chi }^2}\varphi -B\frac{d}{d\chi }{S}_3=0.$

(15)

Then, after integrating Eq. (15) once as a function of χ and choosing the source term in the form S₃ (χ) = f₀/B cos (ωχ), the following dynamical system can be easily found:

$\left\{\begin{array}{cc}\hfill \frac{\mathrm{d}\mathrm{\varphi }}{\mathrm{d}\chi }& =Y,\hfill \\ \hfill \frac{\mathrm{d}\mathrm{Y}}{\mathrm{d}\chi }& =\frac{u_0\varphi }{B}-\frac{A{\varphi }^3}{3B}-\frac{CY}{B}+f_0\cos (\omega \chi ).\hfill \end{array}\right.$

(16)

IV.1. The Hamiltonian system

The Hamiltonian system according to Eq. (16) by neglecting the Burgers’s term and the external force, reads:

$\left\{\begin{array}{cc}\hfill \frac{\mathrm{d}\mathrm{\varphi }}{\mathrm{d}\chi }& =Y,\hfill \\ \hfill \frac{\mathrm{d}\mathrm{Y}}{\mathrm{d}\chi }& =\alpha \varphi -\beta {\varphi }^3.\hfill \end{array}\right.$

(17)

The Hamiltonian function with respect to system (17) is obtained as: H (ϕ, Y) = Y²/2 + dV(ϕ)/dϕ, where V(ϕ) = −α ϕ²/2 + 4 β ϕ⁴ stand for the potential energy of the system, and α = u₀/B, and β = A/3 B.

When the conditions for a double-well potential are satisfied, a homoclinic orbit emerges in the system’s phase space (see Figure 4), revealing three equilibrium points: two centers and a saddle point. The homoclinic orbit encloses the centers, supporting periodic waves, while the saddle point allows for soliton solutions. By setting H (ϕ, Y) = 0, the homoclinic trajectory can be determined, and solving for the resulting displacement and differentiating to find velocity.

$\left(\genfrac{}{}{0ex}{}{{\varphi }_h}{Y_h}\right)=\left(\genfrac{}{}{0ex}{}{\sqrt{\frac{2\alpha }{\beta }}sech(\sqrt{\alpha }\chi )}{-\sqrt{\frac{2\alpha }{\beta }}\sqrt{\alpha }sech(\sqrt{\alpha }\chi )tanh(\sqrt{\alpha }\chi )}\right),$

(18)

Figure 4.

Phase portrait in panel (a) and potential energy in panel (b) for the unperturbed dynamical system (17). Along with ϑ₁ = 0.6, ϑ₂ = 0.3, κ = 2, ϖ = 0.05, ρ = 1.1, and u₀ = 0.03.

IV.2. Melnikov analysis

We apply here the Melnikov method to analytically and numerically determine the threshold conditions for the emergence of periodic regimes and horseshoe chaos in the system under perturbation. Following conventional approaches, the perturbed Hamiltonian equation can be expressed as follows^54–56:

$\frac{\partial }{\partial t}\left(\genfrac{}{}{0ex}{}{\varphi }{Y}\right)=F\left(\varphi ,Y\right)+ϵG\left(\varphi ,Y,\chi -{\chi }_0\right),$

(19)

here, the expression of the vectors fields F (ϕ, Y) and G (ϕ, Y, χ − χ₀) take the form:

$\left\{\begin{array}{cc}\hfill F\left(\varphi ,Y\right)& =\left(\genfrac{}{}{0ex}{}{Y}{\frac{u_0\varphi }{B}-\frac{A{\varphi }^3}{3B}}\right)\hfill \\ \hfill \\ \hfill G\left(\varphi ,Y,\chi -{\chi }_0\right)& =\left(\genfrac{}{}{0ex}{}{0}{-\frac{CY}{B}+f_{0}cos(\omega (\chi -{\chi }_0))}\right).\hfill \end{array}\right.$

(20)

For ϵ = 0, the homoclinic orbit connects to the unstable (saddle) point, where the eigenvalues of the linearized system are real and of opposite sign. When ϵ ≠ 0, chaos may arise due to the transverse intersection of the stable and unstable manifolds, a condition detectable via the Melnikov theorem.

The core idea of this approach is to derive a function that quantifies the separation between the stable and unstable manifolds near one or two saddle points in the perturbed system. When the Melnikov function associated with this method vanishes for specific bifurcation parameters, the stable and unstable manifolds intersect transversally in the Poincaré section, away from the saddle point. As a result, they will create a kind of Smale horseshoe map that will be chaotic. According to Smale-Birkhoff theorem,^57–59 the presence of such an intersection causes a gradual shift toward chaotic dynamics. The Melnikov theorem states that the Melnikov function M (χ₀) provides the distance between the stable and unstable manifolds as follows:

$M({\chi }_0)={\int }_{-\infty }^{\infty }F\left(\varphi ,Y\right)\wedge G\left(\varphi ,Y,\chi -{\chi }_0\right)\mathrm{d}\chi .$

(21)

We assume that, M (χ₀) = σ I₁ + f₀ I₂, with σ = −C/B, $I_1={\int }_{-\infty }^{\infty }{Y_h}^2\mathrm{d}\chi$ , and $I_2={\int }_{-\infty }^{\infty }{Y}_{h}cos(\omega (\chi -{\chi }_0))\mathrm{d}\chi$ .

Following the assessment of these basic integrals (refer to Appendix B), the homoclinic Melnikov function is obtained. Let’s examine the saddle point’s invariant manifold intersections. It is well recognized that chaos can only exist when certain crossings occur. Since the distance between the perturbed stable and unstable manifolds in the Poincaré section can be determined by the Melnikov function theory, respecting the homoclinic loops under perturbations involves that, at χ₀, if M (χ₀) has a simple zero, then a homoclinic bifurcation occurs, indicating the potential for chaotic behavior. Then, the following represents the general required condition for which the invariant manifolds intersect:

$f_0=\mid \frac{4\sigma \sqrt{{\alpha }^3}}{3\pi \omega \sqrt{2\beta }sech\left(\frac{\pi \omega }{2\sqrt{\alpha }}\right)}\mid .$

(22)

According to Melnikov, the following are prerequisites for chaos:

$f_0\ge f_{0_{c}r}=\mid \frac{4\sigma \sqrt{{\alpha }^3}}{3\pi \omega \sqrt{2\beta }sech\left(\frac{\pi \omega }{2\sqrt{\alpha }}\right)}\mid ,$

(23)

where f _0cr stand for the threshold function.

When ϵ ≠ 0, non-periodic motion emerges as predicted by the Melnikov criterion, which requires the gap between stable and unstable manifolds to vanish. The threshold value of f _0cr for the emergence of a transverse intersection between the perturbed and unperturbed manifolds is determined by the criterion in Eq. (23). Such a circumstance is referred to as essential to chaos’s existence. Figure 5 plots the threshold condition in the space (f₀, ω). The system exhibits periodic motion for (f₀, ω) taken below the lower bound line, but the upper domain shows potentially chaotic motion.

Figure 5.

Critical amplitude f_0cr for the emergence of non periodic motion as a function of the frequency, ϑ₁ = 0.6, ϑ₂ = 0.3, κ = 2, ρ = 1.1, ϖ = 0.05, η = 0.1, and μ _i = 0.1.

V. Quasi-periodic and chaotic oscillations

Here, to confirm the analytical prediction according to the Melnikov criteria, we examine how the system responds to the external periodic perturbation f₀ in terms of its strength in Figures 6–8. So, we are going to select the values of f₀ and ω in Figure 5 that satisfy the threshold condition. The perturbed dynamical system given by Eq. (16) is integrated numerically using a fourth-order Runge-Kutta method with a fixed time step of Δχ = 0.01. Initial conditions is chosen near the unperturbed homoclinic orbit: (ϕ (0), Y (0)) = (ϕ _h (0) + 10⁻⁴, Y _h (0)), where (ϕ _h , Y _h ) is given by Eq. (18). Figure 6 illustrates the profile of the IAWs, for f₀ = 0.01 and ω = 0.5. For this value of f₀ and ω, the time series plot (Figure 6(a)) and phase space (Figure 6(b)) display a periodic oscillation. The Poincaré section (Figure 6(c)) shows one point, which confirms the existence of periodic motion. On the other hand, for f₀ = 0.3 and ω = 0.75, the time series plot (Figure 7(a)) display oscillations with incommensurable frequencies. A torus’s surface is heavily filled with phase space oscillations (Figure 7(b)). Furthermore, a collection of highly populated points makes up the Poincaré section of Figure 7(c), which further illustrates the motion’s quasi-periodic nature.

Figure 6.

Periodic regime of IAWs for (f₀, ω) = (0.01, 0.5), and ϑ₁ = 0.6, ϑ₂ = 0.3, κ = 2, ρ = 1.1, ϖ = 0.05, η = 0.01, and μ _i = 0.1.

Figure 7.

Quasi-periodic regime of IAWs for (f₀, ω) = (0.3, 0.75), and ϑ₁ = 0.6, ϑ₂ = 0.3, κ = 2, ρ = 1.1, ϖ = 0.05, η = 0.01, and μ _i = 0.1.

Figure 8.

Chaotic regime of IAWs for (f₀, ω) = (0.3, 0.5), ϑ₁ = 0.6, ϑ₂ = 0.3, κ = 2, ρ = 1.1, ϖ = 0.05, η = 0.01, and μ _i = 0.1.

The findings shown in Figure 8 are obtained for parameters ω = 0.5 and f₀ = 0.3. The time series plot of Figure 8(a) demonstrates erratic oscillations and nearly random behavior. The points are widely distributed and display random oscillations with no discernible structure, as seen in the corresponding Poincaré section in Figure 8(c). These features characterize a chaotic pattern. Thus, our analytical predictions agree with the numerical results. With the same settings as in Figure 8(a)–(c). A similar outcome was also attained in Refs. 81, where the authors conducted a thorough investigation of PAWs.

VI. Conclusion

This work has presented a comprehensive investigation of low-frequency, dissipative ion-acoustic solitary waves (IASWs) in a multi-component plasma composed of warm adiabatic ions, superthermal (κ-distributed) electrons and positrons, and isothermal electrons. By applying the reductive perturbation method, we derived a damped forced modified Korteweg-de Vries-Burgers (DFMKdVB) equation, which governs the nonlinear dynamics of compressive solitary waves in the presence of collisional dissipation, kinematic viscosity, and an external periodic perturbation. Our analytical and numerical results lead to the following key conclusions:

1. The DFMKdVB model appropriately describes low-frequency dissipative ion-acoustic waves in the studied plasma. Its coefficients, which depend explicitly on fundamental plasma parameters (ϑ₁, ϑ₂, κ, ρ, ϖ), provide a self-consistent link between the microscopic plasma conditions and the observed macroscopic wave behavior.

2. The characteristics of IASWs are sensitively controlled by plasma parameters. The superthermal index κ, temperature ratios (ρ, ϖ), and density ratios significantly modify the wave’s amplitude and width, as demonstrated analytically through the soliton solution (Eq. (12)) and confirmed numerically in Figure 1.

3. Dissipation and external forcing are competing influences the shape of wave evolution. Ion-neutral collisions μ _i and kinematic viscosity η systematically dampen wave amplitude (Figure 2(a) and (b)), while an external periodic force can amplify it, depending on its strength f₀ and frequency ω (Figure 2(c) and (d)).

4. The transition from ordered to chaotic dynamics is predictable via Melnikov analysis. By applying Melnikov’s method to the perturbed Hamiltonian form of the DFMKdVB equation, we derived an explicit analytical threshold f_0cr(ω) (Eq. (23)) for the onset of homoclinic chaos. This prediction, visualized in Figure 5, was rigorously validated by numerical diagnostics:

     Forcing parameters (f₀, ω) below this threshold yielded periodic oscillations (Figure 6). Parameters near the threshold led to quasi-periodic motion (Figure 7). Exceeding the threshold reliably induced chaotic dynamics, characterized by a fractal Poincaré section and erratic time series (Figure 8).

In summary, this study elucidates the intricate balance between nonlinearity, dispersion, dissipation, and external forcing in a complex plasma environment. The synthesis of analytical perturbation theory, Melnikov criteria, and phase-space diagnostics provides a robust framework for analyzing chaotic regimes in dissipative, driven plasma systems. The strong agreement between the analytically predicted chaos threshold and the numerically observed dynamical transitions confirms the reliability of this approach. The results of this study offer valuable insights for interpreting nonlinear wave phenomena in both laboratory experiments and space plasma settings, such as planetary magnetospheres or turbulent solar wind, where superthermal distributions, collisional effects, and external perturbations are prevalent.