Computational Analysis of Electrostatic Potential Waves in Rocket Engine Plasma Using Advanced Zakharov–Kuznetsov Equation Solvers

A fundamental aspect of contemporary aerospace research is the characterization of plasma dynamics in the exhaust plumes of high-performance rocket engines [1], [2], [3]. Engine stability and electromagnetic interference are directly impacted by complicated electrostatic phenomena that arise from the interaction between nonlinearity and dispersion as ionized gases are released via the nozzle [4], [5]. The electrostatic potential, which acts as a normalized representation of the voltage inside the plasma field, is essential to these dynamics [6], [7], [8]. To model these multi-dimensional fluctuations, this study utilizes the modified Zakharov–Kuznetsov equation [9], a governing nonlinear partial differential equation, formulated as follows:

where, $\alpha$ is the electrostatic potential (normalized voltage in the plasma); $\tau$ denotes the temporal coordinate, and $\mathscr{E}$, $\vartheta$, and $\omega$ denote the three spatial coordinates. The non-linear fluid steepening is controlled by the coefficients $A$ and $B$, which are mainly influenced by changes in background ion density and electron temperature $T_e$. The electronic Debye length ($\gamma_D$) and ion Larmor radius $\rho_i$, which are inversely proportional to the engine's magnetic confinement field strength $B_0$, determine the dispersion coefficients $C$ and $D$, which stop unconstrained steepening [10], [11], [12]. To resolve this completely, this study presents a comprehensive summary table (Table 1) detailing these variables under typical high-density plasma propulsion environments. The parameter ranges are adopted from representative plasma conditions reported in previous studies.

Despite the mathematical robustness of the modified Zakharov–Kuznetsov equation, obtaining exact analytical solutions remains a significant computational challenge due to the higher-order dispersive terms [13], [14], [15], [16]. This study introduces an advanced computational solver based on a specific wave transformation:

where, $\beta_1, \beta_2, \beta_3$ the direction-cosine components of the wave in the three spatial directions represent the direction of the wave in the rocket nozzle; $\rho$ denotes the wave frequency, unrelated to $\rho_{\mathrm{i}}$ where $\rho_{\mathrm{i}}$ is ion Larmor radius (physical parameter); and $\gamma$ denotes the transformed traveling-wave coordinate. A nonlinear ordinary differential equation is created by transforming the nonlinear partial differential equation that is particular to the coordinates of a rocket nozzle. Using the homogeneous balancing principle [17], [18], [19] and an ansatz-based approach [20], this study obtains a complete set of computational solutions that characterize the evolution of the electrostatic potential of the plasma.

Applying the chain rule of partial differentiation to project Eq. (1) onto the one-dimensional coordinates field, the individual derivative transformations can be obtained as follows:

$ \alpha_\tau=-\rho \alpha^{\prime}, \alpha_{\mathscr{E}}=\beta_1 \alpha^{\prime}, \alpha_{\mathscr{E} \mathscr{E} \mathscr{E}}=\beta_1^3 \alpha^{\prime \prime \prime}, \alpha_{\mathscr{E} \vartheta \vartheta}=\beta_1 \beta_2^2 \alpha^{\prime \prime \prime}, \alpha_{\mathscr{E} \omega \omega}=\beta_1 \beta_3^2 \alpha^{\prime \prime \prime} $

where, the prime represents the total derivative with respect to $\gamma\left(\alpha^{\prime}=d \alpha / {d} \gamma\right)$. Substituting these expressions back into Eq. (1) yields:

$ -\rho \alpha^{\prime}+A \beta_1 \alpha \alpha^{\prime}+B \beta_1 \alpha^2 \alpha^{\prime}+C \beta_1^3 \alpha^{\prime \prime \prime}+D \beta_1\left(\beta_2^2+\beta_3^2\right) \alpha^{\prime \prime \prime}=0 $

Factoring out the shared spatial derivative $\beta_1$ from the dispersive terms allows us to group the combined dispersion coefficient as follows:

$ -\rho \alpha^{\prime}+A \beta_1 \alpha \alpha^{\prime}+B \beta_1 \alpha^2 \alpha^{\prime}+\beta_1\left[C \beta_1^2+D\left(\beta_2^2+\beta_3^2\right)\right] \alpha^{\prime \prime \prime}=0 $

To simplify the expression, this study defines the net multi-dimensional engine dispersion factor as $\beta=\beta_1\left[C \beta_1^2+D\left(\beta_2^2+\beta_3^2\right)\right]$. This condenses the equation to:

$ -\rho \alpha^{\prime}+A \mathrm{Q}_1 \alpha \alpha^{\prime}+B \beta_1 \alpha^2 \alpha^{\prime}+\beta \alpha^{\prime \prime \prime}=0 $

This study then integrates this equation once with respect to $\gamma$, setting the constant of integration to zero under localized boundary conditions $\left(\alpha, \alpha^{\prime}, \alpha^{\prime \prime}, \rightarrow 0\right.$ as $\left.|\gamma| \rightarrow \infty\right)$:

$ -\rho \alpha+\frac{1}{2} A \beta_1 \alpha^2+\frac{1}{3} B \beta_1 \alpha^3+\beta \alpha^{\prime \prime}=0 $

By identifying the localized wave frequency term relative to nozzle velocity as $\wp=-\rho$, now called wave velocity relative to exhaust speed, this study successfully arrives at the final operational computational ordinary differential equation:

where, $\wp$ denotes the wave velocity relative to the exhaust speed, and Plain $\beta$, defined as $\beta=\beta_1 [{C} \beta_1{ }^2+{D}(\beta_2{ }^2+\beta_3{ }^2)]$ the “combined dispersive/engine dispersion factor” represents the combined dispersive effects of the engine’s magnetic field.

The primary objective of this investigation is to categorize the resulting wave structures into three physically distinct regimes as follows:

For soliton regimes ($t >$ 0): These exhibit non-dispersive, highly localized profiles. For example, the bright soliton ($C_1$) has a structural width that is precisely scaled as $1 /(\sqrt{\mathrm{t}})$ and an explicit, stable peak amplitude that is proportional to $\sqrt{t}$. This suggests that the energy packet is compressed by stronger magnetic fields (which raise $t$), raising peak concentration while reducing spatial footprints.

For periodic instabilities ($t <$ 0): In contrast to solitons, these do not show localized degradation. Over the whole spatial range from $-\infty$ to $+\infty$, there is a constant variation in the amplitude. As the nozzle cross-section gets smaller, the corrugated structures’ spatial frequency scales linearly with $\sqrt{(-t)}$, causing high-frequency oscillations and dense wave trains.

For rational and transition states ($t$ = 0): These serve as the precise limits in mathematics. They have an unlimited effective spatial breadth and no periodic footprint, exhibiting algebraic decay proportional to $1 /(\gamma)$, which represents a smooth, transient spatial dissipation into background noise. Through this computational analysis, this study provides a theoretical foundation for predicting electrostatic potential distributions, offering critical insights for the optimization of next-generation plasma-based propulsion systems.

This study considers the following function:

where, $C(\mathscr{M}, h)$ is an unknown function, and $\mathscr{H}$ is a polynomial in $C(\mathscr{M}, h)$.

Step 1: The following transformation is used:

where, $\beta$ and $F$ are constants to be determined in $\gamma=\beta\left(\mathscr{M}-{Fh}+\gamma_0\right)$ a wave-number placeholder in the abstract method, never tied back to the $\beta_1 / \beta_2 / \beta_3$ or the dispersion factor above and $\gamma_0$ is an arbitrary constant. From Eqs. (4) and (5), the following equation can be derived:

Step 2: The ansatz method is considered as the form [2]:

where, $N=\left(\frac{\mathfrak{I}^{\prime}}{\mathfrak{I}}+\frac{\mathscr{Q}}{2}\right),\left|B_{-S}\right|+\left|B_S\right| \neq 0$ and $\mathfrak{I}=\mathfrak{I}(\gamma)$ satisfies the equation.

where, $B_k( \pm 1, \pm 2, \ldots \ldots \ldots, \pm S), \mathcal{L}$ and $\sigma$ are coefficient constants. Implementing the homogeneous balance principle in Eq. (6), the positive integer $S$ can be determined. The following equation can be derived from Eq. (8):

where, $t=\left(\frac{\mathscr{D}^2-4 \partial}{4}\right)$ and $t$ is calculated by $\mathscr{L}$ and $\partial$. Therefore, $N$ satisfies Eq. (9), which admits five types of solutions.

If $t>$ 0, then we get:

$ \begin{aligned} & N=\sqrt{t} \tanh (\sqrt{t} \gamma) ; \\ & N=\sqrt{t} \operatorname{coth}(\sqrt{t} \gamma) ; \end{aligned} $

If $t$ = 0, then we get:

$ N=\frac{1}{\gamma} ; $

If $t<$ 0, then we get:

$ \begin{gathered} N=-\sqrt{-t} \tanh (\sqrt{-t} \gamma) ; \\ N=\sqrt{-t} \operatorname{coth}(\sqrt{-t} \gamma) . \end{gathered} $

Step 3: The left-hand side of Eq. (6) is transformed into a polynomial in $N$ by applying Eqs. (7), (6), and (9) and grouping all terms of the same order of $N$ together. Equating each polynomial's coefficient to zero, a set of algebraic equations can be obtained, which can be solved to find the values of $B_k, k= \pm 1, \pm 2, \ldots \ldots, \pm S, \mathcal{L}, \partial$. Finally, this study obtains the general solutions of Eq. (8) from $B_k, \mathcal{L}, \partial$.

Setting $S$ = 1, the following equation can be derived from Eq. (7):

After Eq. (10) is inserted into Eq. (3), the coefficient of $N$ is obtained, and the resulting system is solved. The three sets below can be found.

Using $B_1=\sqrt{\frac{-6 \beta}{B \beta_1}}, B_{-1}=t \sqrt{\frac{-6 \beta}{B \beta_1}}$, and $B_0=-\frac{A}{2 B}$

When $t>$ 0, the following soliton solutions (hyperbolic tanh and hyperbolic coth) are as follows:

$ C_1(\gamma)=-\frac{A}{2 B}+\sqrt{\frac{-6 \beta}{B \beta_1}}(\sqrt{t} \tanh (\sqrt{t} \gamma))+t \sqrt{\frac{-6 \beta}{B \beta_1}}(\sqrt{t} \tanh (\sqrt{t} \gamma))^{-1} $

$ C_2(\gamma)=-\frac{A}{2 B}+\sqrt{\frac{-6 \beta}{B \beta_1}}(\sqrt{t} \operatorname{coth}(\sqrt{t} \gamma))+t \sqrt{\frac{-6 \beta}{B \beta_1}}(\sqrt{t} \operatorname{coth}(\sqrt{t} \gamma))^{-1} $

When $t$ = 0, the rational solution is as follows:

$ C_3(\gamma)=-\frac{A}{2 B}+\sqrt{\frac{-6 \beta}{B \beta_1}}\left(\frac{1}{\gamma}\right) $

When $t<$ 0, and the periodic solution (tan and cot) are as follows:

$ C_4(\gamma)=-\frac{A}{2 B}+\sqrt{\frac{-6 \beta}{B \beta_1}}(\sqrt{-t} \tan (\sqrt{-t} \gamma))+t \sqrt{\frac{-6 \beta}{B \beta_1}}(\sqrt{-t} \tan (\sqrt{-t} \gamma))^{-1} $

$ C_5(\gamma)=-\frac{A}{2 B}+\sqrt{\frac{-6 \beta}{B \beta_1}}(-\sqrt{-t} \cot (\sqrt{-t} \gamma))+t \sqrt{\frac{-6 \beta}{B \beta_1}}(-\sqrt{-t} \cot (\sqrt{-t} \gamma))^{-1} $

Using $B_1=-\sqrt{\frac{-6 \beta}{B \beta_1}}, B_{-1}=-t \sqrt{\frac{-6 \beta}{B \beta_1}}$, and $B_0=-\frac{A}{2 B}$

When $t>$ 0, the following soliton solutions (hyperbolic tanh and hyperbolic coth) are as follows:

$ C_6(\gamma)=-\frac{A}{2 B}-\sqrt{\frac{-6 \beta}{B \beta_1}}(\sqrt{t} \tanh (\sqrt{t} \gamma))-t \sqrt{\frac{-6 \beta}{B \beta_1}}(\sqrt{t} \tanh (\sqrt{t} \gamma))^{-1} $

$ C_7(\gamma)=-\frac{A}{2 B}-\sqrt{\frac{-6 \beta}{B \beta_1}}(\sqrt{t} \operatorname{coth}(\sqrt{t} \gamma))-t \sqrt{\frac{-6 \beta}{B \beta_1}}(\sqrt{t} \operatorname{coth}(\sqrt{t} \gamma))^{-1} $

When $t$ = 0, the rational solution is as follows:

$ C_8(\gamma)=-\frac{A}{2 B}-\sqrt{\frac{-6 \beta}{B \beta_1}}\left(\frac{1}{\gamma}\right) $

When $t<$ 0, the periodic solutions (tan and cot) are as follows:

$ C_9(\gamma)=-\frac{A}{2 B}-\sqrt{\frac{-6 \beta}{B \beta_1}}(\sqrt{-t} \tan (\sqrt{-t} \gamma))-t \sqrt{\frac{-6 \beta}{B \beta_1}}(\sqrt{-t} \tan (\sqrt{-t} \gamma))^{-1} $

$ C_{10}(\gamma)=-\frac{A}{2 B}-\sqrt{\frac{-6 \beta}{B \beta_1}}(-\sqrt{-t} \cot (\sqrt{-t} \gamma))-t \sqrt{\frac{-6 \beta}{B \beta_1}}(-\sqrt{-t} \cot (\sqrt{-t} \gamma))^{-1} $

The soliton tanh and coth (opposite sign ${B}_{-1}$) are as follows:

$ C_{11}(\gamma)=-\frac{A}{2 B}+\sqrt{\frac{-6 \beta}{B \beta_1}}(\sqrt{t} \tanh (\sqrt{t} \gamma))-t \sqrt{\frac{-6 \beta}{B \beta_1}}(\sqrt{t} \tanh (\sqrt{t} \gamma))^{-1} $

$ C_{12}(\gamma)=-\frac{A}{2 B}+\sqrt{\frac{-6 \beta}{B \beta_1}}(\sqrt{t} \operatorname{coth}(\sqrt{t} \gamma))-t \sqrt{\frac{-6 \beta}{B \beta_1}}(\sqrt{t} \operatorname{coth}(\sqrt{t} \gamma))^{-1} $

The periodic tan and cot (opposite sign ${B}_{-1}$) are as follow:

$ \mathrm{C}_{13}(\gamma)=-\frac{A}{2 B}+\sqrt{\frac{-6 \beta}{B \beta_1}}(\sqrt{-t} \tan (\sqrt{-t} \gamma))-t \sqrt{\frac{-6 \beta}{B \beta_1}}(\sqrt{-t} \tan (\sqrt{-t} \gamma))^{-1} $

$ C_{14}(\gamma)=-\frac{A}{2 B}-\sqrt{\frac{-6 \beta}{B \beta_1}}(\sqrt{-t} \cot (\sqrt{-t} \gamma))+t \sqrt{\frac{-6 \beta}{B \beta_1}}(\sqrt{-t} \cot (\sqrt{-t} \gamma))^{-1} $

The rational inverse solution is as follows:

$ C_{15}(\gamma)=-\frac{A}{2 B}+t \sqrt{\frac{-6 \beta}{B \beta_1}} \gamma $

A crucial understanding of the plasma dynamics of a rocket engine is provided by discussing the 15 three-dimensional graphical representations of the electrostatic potential. These graphs show how dispersion, magnetic fields, and nonlinearity interact to impact wave propagation.

The analyses of the solutions are as follows:

• Analysis of Soliton Solutions (Figure 1, Figure 2, Figure 3, Figure 4, Figure 5, Figure 6)

These graphs, derived from the hyperbolic tangent and cotangent auxiliary solutions, depict solitary waves where $t >$ 0. These “ridges” represent stable energy packets that do not split as they travel through the exhaust plume. Figure 1 and Figure 3 (tanh) represent “bright” or “dark” solitons, showing localized increases or decreases in potential. Figure 2 and Figure 4 (coth) exhibit unique peaks, indicating zones of high-intensity electrostatic concentration. Figure 5 and Figure 6 (mixed) show the interaction of $N$ and $N^{-1}$ terms, resulting in multi-peaked or complex solitary structures.

• Analysis of Periodic Solutions (Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12)

These surfaces show continuous fluctuations or “wave trains” in the plasma potential when $t <$ 0. The oscillatory form of the potential is based on trigonometric functions. The recurring corrugated surface in Figure 7 and Figure 9 (tan) indicates frequent periodic instability. Figure 8 and Figure 10 (cot) show sharp periodic spikes that correspond to fast voltage fluctuations in the rocket nozzle. Figure 11 and Figure 12 (mixed) show modulated oscillations in which wave interference causes the amplitude of the oscillations to change throughout the plasma field.

• Analysis of Rational and Growth Solutions (Figure 13, Figure 14, Figure 15)

These take place by direct linear transformation or at the crucial border where $t$ = 0. These solutions show transition states in which waves do not form stable solitons or oscillate. Figure 13 and Figure 14 (rational) show a single, localized “lump” that decays algebraically, signifying a transient potential disturbance that fades into the surrounding plasma. Figure 15 (linear growth) shows a non-oscillatory potential ramp and is frequently linked to a steady-state shift in the background ionization of the plume.

Computational research leads to 15 different topological solutions to the modified Zakharov–Kuznetsov problem. Depending on the type of electrostatic potential $\alpha$, these solutions are divided into three physical regimes. The proposed approach is innovative because it uses a generalized auxiliary ansatz equation based on a comprehensive discriminant $t=\left(\frac{\mathcal{L}^2-4\partial}{4}\right)$ in conjunction with a unified multi-dimensional wave transformation. This blended composition offers a number of special benefits.

• Simultaneous triple-regime capture: The proposed framework simultaneously extracts 15 different analytical solutions covering three completely different physical topologies (hyperbolic solitons, periodic trigonometric wave trains, and rational algebraic lumps) from a single algebraic system, rather than performing separate calculations for different wave states.

• Asymmetric and mixed structural modes: Unlike standard setups, the solver in this study introduces “Set C” mixed coefficient parameters, yielding complex hybrid models (such as Solutions $C_{11}$ and $C_{12}$) that characterize multi-peaked wave structures arising from intense local plasma-wall boundaries.

• Direct mapping to plasma confinement: The transformation in this study quickly applies the analytical results to plasma diagnostic interpretations by directly connecting the mathematical coordinates to the geometry of a converging-diverging rocket nozzle.

Solitary waves are supported by the plasma in the regime where the discriminant $t$ is positive. These physically correspond to localized electrostatic energy packets that move through the exhaust plume without changing form.

• Bright solitons: These appear as high-potential peaks in Solutions 1, 3, and 5. These correspond to areas of localized electron compression in a rocket nozzle, which may be a sign of significant heat strain on the nozzle walls.

• Dark solitons: These are “holes” or localized potential dips (Solutions 2, 4). When the plasma density briefly falls below the equilibrium condition, they are indicative of rarefactive ion waves.

• Singular solitons (Solutions 6, 7): The “cusp-like” peaks seen in this study are suggestive of wave-breaking occurrences, in which the dispersion is dominated by nonlinearity, frequently resulting in localized turbulence.

The solutions change into periodic trigonometric functions when $t$ is negative. Continuous wave trains are represented by these solutions.

• Harmonic oscillations (Solutions 8–13): These solutions describe the “hum” or high-frequency instabilities that are frequently found in plasma plumes or Hall effect rockets.

• Corrugated potential surfaces: These solutions display regular, repeated ridges in their three-dimensional visualizations. These are periodic variations in the thrust vector in the context of propulsion, which, if not reduced by magnetic shielding, can cause vibration and structural fatigue.

The boundary between solitary and periodic behavior is represented by the rational answers (Answers 14 and 15). These are usually “lump” solutions, which are transient, localized disruptions. These could simulate “sparking” or brief arcs in the rocket engine brought on by abrupt changes in the ion Larmor radius or variations in the surrounding magnetic field.

For engine optimization, the electrostatic potential's sensitivity to the governing coefficients ($A$, $B$, $C$, and $D$) is crucial:

• Dispersion modulation ($C$, $D$): The solitons broaden as the Debye length grows due to the stronger dispersive effects. This suggests that the electrostatic potential is more “spread out” with lower plasma densities, lowering the possibility of localized wall erosion.

• Nonlinear effects ($A$, $B$): Narrower, higher-amplitude “bright” solitons are produced when the nonlinearity is increased by higher electron temperatures ($T_e$). This implies that high-intensity electrostatic spikes are more likely to occur in hotter plumes.

The physical insights drawn from the wave topologies in this study are summarized as follows:

• Plume confinement and erosion mitigation: High-intensity localized potential spike locations are captured by the bright solitary solitons (Solutions 6 and 7). These spikes show areas of high electron compression in an operational nozzle, which can accelerate plasma-wall erosion and produce localized heat hotspots. These mathematical profiles can be used by engineers to map and forecast wall damage.

• Magnetic field regulation: The proposed parametric model shows that the discriminant $t$, which is highly dependent on coefficient $D$ (the ion Larmor radius effect), determines the key structural transition boundary. The findings offer a direct theoretical method for magnetic nozzle management because the Larmor radius is inversely proportional to the external magnetic field strength ($B_0$). Operators can drive the system into stable, predictable solitary energy packets ($t >$ 0) and out of chaotic, high-frequency periodic instability regimes ($t <$ 0) by actively manipulating magnetic field coils.

• Aviation communication protection: High-frequency plasma oscillations are modeled by the periodic corrugated surfaces seen when $t <$ 0. During engine firings, these variations cause significant electromagnetic interference, which results in communication blackout windows. Engineers can create focused magnetic shielding techniques by using the proposed analytical framework to pinpoint the precise thrust and density thresholds that cause these instabilities.

The nonlinear dynamics of electrostatic potential waves within rocket engine plasma plumes, modeled by the modified Zakharov–Kuznetsov equation, were effectively investigated in this study using a reliable computer solver. This study obtained fifteen precise analytical solutions describing the intricate behavior of ionized gases in magnetized nozzle environments by combining an advanced wave transformation with an ansatz-based methodology. For solitary wave solutions ($t >$ 0), the study divided these solutions into three different physical regimes. Understanding localized heat flow on nozzle surfaces requires the demonstration of stable, localized energy packets (bright, dark, and solitary solitons) that can move through the exhaust without dissipating. For periodic oscillations ($t <$ 0), wave trains and high-frequency instabilities that indicate periodic potential fluctuations were found; these are crucial for identifying electromagnetic interference in communication systems. For rational solutions ($t$ = 0), the transition between stable and unstable plasma states is represented by captured transient disturbances. The parametric analysis showed that the Debye length and the ion Larmor radius had a significant impact on the morphology of these waves, particularly their amplitude and width. These results imply that engineers can successfully regulate the change from chaotic instabilities to stable energy transport by varying the external magnetic field of a plasma thruster. In the end, this computational framework offers a theoretical standard for designing and refining next-generation propulsion systems based on plasma. In order to improve the forecast accuracy in a variety of aerospace situations, future work will concentrate on expanding this model to incorporate collisional effects and non-thermal ion distributions.

The highlights of this study are as follows:

• Advanced modeling: Successfully modeled three-dimensional electrostatic potential waves in magnetized rocket plasma using the modified Zakharov–Kuznetsov equation.

• Novel computational approach: Implemented a unified wave transformation and ansatz-based solver to derive 15 distinct exact analytical solutions.

• Diverse wave topologies: Characterized stable bright/dark solitons, periodic wave trains, and rational “lump” solutions within the exhaust plume.

• Parametric sensitivity: Quantified the influence of ion Larmor radius and Debye length on wave stability and energy localization.

• Propulsion application: Provided a theoretical framework for mitigating plasma instabilities and optimizing magnetic nozzle efficiency in aerospace systems.