Nonlinear evolution, propagation, electron-trapping, and damping effects of ion-acoustic solitons using fully kinetic PIC simulations

Nonlinear waves in a plasma emerge from the collective interplay between charged particles and the self-consistent electromagnetic field within the plasma. Some examples of these waves include Alfvén waves, Langmuir waves, and ion acoustic waves.^1–3 The study of nonlinear wave phenomena in plasma constitutes a pivotal aspect of plasma physics, given its broad relevance in fields such as astrophysics and space physics. Ion acoustic waves exhibiting nonlinear behavior, particularly in the form of solitons, have garnered extensive research interest due to their ubiquity and significance across a broad spectrum of plasma applications.^4–28 The predominant approach in these studies has involved analytical and numerical studies grounded in fluid model descriptions of the plasma environment. The reductive perturbation method for weak nonlinearities⁶ and the fully nonlinear Sagdeev approach^5,29 stand out as the classic methodologies. Using the reductive perturbation method, Washimi and Taniuti demonstrated that, when considering small but finite amplitude perturbations in such plasmas, the evolution of ion acoustic waves can be described by the integrable Korteweg–de Vries (KdV) equation, which admits a distinct class of soliton solutions as exact traveling wave solutions.^6,30

These soliton structures arise from a precise balance between the opposing effects of nonlinear steepening and dispersive broadening of the wave. Specifically, in the KdV equation, the quadratic nonlinear term accounts for wave steepening, while the third-order dispersion term contributes to pulse spreading. Solitons possess several unique properties stemming from the mathematical characteristics and infinite number of conservation laws intrinsic to the KdV equation.¹⁹ Notable properties include the ability to propagate over long distances without degradation of their shape or amplitude, as well as the remarkable phenomenon of preserving individual identities upon colliding with other solitons. Ikezi et al.¹² first experimentally detected the formation and propagation of solitons in an argon plasma and since then solitons have been commonly observed in both laboratory and space plasma environments.^31–34 Recent investigations into the advancements in inertial confinement fusion,³⁵ the complementary heating method in fusion devices,³⁶ and techniques for accelerating charged particles³⁷ have demonstrated close ties to the dynamics of solitons. A broad spectrum of research endeavors has been undertaken concerning solitons, encompassing their formation,³⁸ propagation,^39,40 reflection,⁴¹ collisions,⁴² decay,^43,44 wave instabilities,^45,46 etc.

One emerging area of interest is the idea of using precursor solitons to detect sub-centimeter orbital debris. Current ground-based methods and on-orbit detectors cannot detect such small debris, which can still prove fatal to spacecraft. Sen et al. numerically solved the forced Korteweg–de Vries (fKdV) equation and demonstrated the existence of precursor solitons created by orbiting charged debris.^47,48 Tiwari and co-workers presented similar results from solving the full set of fluid equations.^49,50 Jaiswal et al. experimentally validated the existence of these waves.⁵¹ Truitt and Hartzell used a similar approach to Sen by numerically solving the fKdV equation for varying space regions. They showed that precursor waves with detectable amplitudes could propagate indefinitely without changing shape or size.⁵² 

However, these results were generated using the reductive perturbation method. It is essential to understand whether deviating from the fundamental assumptions underlying the fKdV equation, such as the small amplitude perturbation assumption or the idealized fluid model, can significantly influence the characteristics, formation, propagation, and stability of the nonlinear solution. To gain insights into this, a first-principles study initializing KdV solitons in a one-dimensional (1D) plasma and analyzing their properties seems appropriate. Therefore, a fully kinetic approach with minimal assumptions, such as a particle-in-cell (PIC) simulation, can be employed to investigate this phenomenon. This research motivation underpins the present study.

The study of solitons using various numerical methods has been previously conducted by a number of authors. Nopoush and Abbasi⁵³ employed hybrid particle-in-cell (PIC) simulations to explore the generation of ion acoustic solitary waves (IASWs). Dharodi et al.⁵⁴ demonstrated the formation of precursor solitons using fully kinetic PIC simulations and their dependence on the source shape. Kakad et al.²⁰ examined the formation of ion acoustic solitary wave chains using fluid simulations. Qi et al.⁵⁵ utilized hybrid PIC simulations to investigate soliton collisions. Sharma et al.⁵⁶ also used hybrid PIC simulations to study the propagation of solitons with varying amplitudes. They showed that larger amplitude solitons initially grow as they propagate, while smaller amplitude solitons propagate without significant changes from their initial state.⁵⁶ 

Furthermore, Qi et al. conducted another study using hybrid PIC simulations to demonstrate that as the soliton amplitude increases, reductive perturbation models fail to accurately predict its behavior.⁵⁷ They showed that when the amplitude surpasses a certain threshold, the wave becomes unstable and undergoes decay into a series of smaller amplitude waves.⁵⁷ 

In this work, we aim to study the propagation of Korteweg–de Vries (KdV) solitons launched in a one-dimensional (1D) plasma according to the solution of the KdV solution for the ion density, velocity, and potential using a fully kinetic plasma model. This study complements the works of Sharma et al.⁵⁶ and Qi et al.⁵⁷ who both modeled the electrons as a Boltzmann fluid. The validity of such hybrid models has not been verified for these types of studies. Their simulations ignored trapping effects and treated electrons as an isothermal background; therefore, they are unable to capture the nonlinear features and the deviation it causes. Furthermore, Kakad et al.⁵⁸ showed that when the phase velocity of the soliton is greater than the ion acoustic velocity, it is necessary to include kinetic effects (including electrons) in the model. This is the case for the solitons in the present work, underscoring the reason for using a fully kinetic model as opposed to a hybrid one.

Ott and Sudan⁵⁹ were one of the pioneers in looking at kinetic effects, namely, Landau damping of solitons. They derived a Korteweg–de Vries equation with a source term to model the lowest-order effects of resonant electrons, demonstrating the occurrence of linear Landau damping in solitons. However, they neglected trapped particle effects, which Meiss and Morrison⁶⁰ later showed to be potentially as significant as linear Landau damping terms. They predicted that nonlinear effects become important when the Landau time, $γL−1=[k2πmeTe/(8mi2)]−1/2$⁠, exceeds the electron bounce time, $ωbe−1=[me/(ek2ϕ)]1/2$⁠. This condition is met in this study; therefore, the role of these nonlinear trapping effects needs to be looked at.

Despite these theoretical advances, the time-dependent electron Landau damping problem remains largely unexplored using numerical methods such as PIC simulations. This gap in the literature also motivates our study, which aims to investigate the physics of nonlinear electron trapping and its effects on soliton properties, including the relationship between amplitude and speed. To our knowledge, this work represents the first comprehensive analysis of long-time soliton propagation using fully kinetic PIC simulations.

This paper is organized as follows: Sec. II provides a detailed description of the particle-in-cell (PIC) simulation employed in this study, including the simulation setup, initial conditions, and relevant parameters. Section III presents and discusses the results of our investigation, which is divided into several subsections. Subsection III A examines the full set of cold fluid equations and their ability to describe the soliton dynamics. Subsection III B focuses on the effects of electron trapping and is further divided into three sub-subsections. Sub-subsection III B 1 explores the impact of non-isothermal electrons on the soliton profile. Sub-subsection III B 2 introduces the modified KdV equation, which incorporates the nonlinearity arising from electron trapping. Sub-subsection III B 3 investigates the electron bounce frequency and its role in the soliton amplitude oscillations. Subsection III C analyzes the relationship between various soliton parameters, such as amplitude, width, and speed, and compares the findings with existing theoretical predictions. Subsection III D examines the effects of Landau damping on the soliton propagation. Finally, Sec. IV summarizes our conclusions, highlights the key findings of this study, and outlines potential directions for future research in the field of ion acoustic solitons.

The simulations were carried out using the open source PIC code called WarpX.^61,62 We chose to use WarpX's electrostatic solver due to its highly parallel and highly optimized nature, which allows it to run efficiently on GPUs. This was an attractive choice considering the large number of particles, the high resolution of both time and spatial scales, and the extended duration of the total run time.

A finite difference scheme is used to solve Poisson's equation. The electric potential, $ϕ$⁠, and the charge density, $ρ$⁠, are calculated at the nodes, while the electric field, E, is calculated at the cell centers. Both electrons and ions are treated as particles in our simulations. To accelerate ion dynamics, we use a fictitious electron-to-ion mass ratio (⁠ $me/mi$⁠) of 0.01. All simulations employ periodic boundary conditions with a domain length of $L=10$ m $=1184λD$⁠, where $λD$ is the Debye length defined as $λD=ε0kBTe/(n0qe2)$⁠. Here, $kB$ is the Boltzmann constant, $ε0$ is the vacuum permittivity, $qe$ is the elementary charge, $Te$ is the electron temperature, and $n0$ is the background plasma density. A longer domain length (⁠ $L=20$ m) is used in some simulations where long-time effects are analyzed. We set $Te=1500$ K and $n0=1011$ m⁻³, approximating typical Low Earth Orbit (LEO) space plasma conditions. Initially, the ion temperature, $Ti$⁠, is set to 0, representing a cold ion species consistent with the Korteweg–de Vries (KdV) formulation, but ion temperature is allowed to evolve with the simulation. A time step, $dt=0.01ωpe$ is used, where $ωpe=n0qe2/(ε0me)/(2π)$⁠. The spatial discretization is dx = 0.0024 m = 0.29  $λD$⁠. This dx was chosen to resolve for the narrow solitons.

In similar research studies, it is common practice to set the initial ion and electron densities equal to each other (⁠ $ni,0=ne,0$⁠).^56–58 This approach is acceptable because the PIC code iteratively adjusts the particle positions until the desired potential, given by Eq. (7), is achieved. However, by utilizing Eq. (9) to calculate the initial electron densities, the simulation can be optimized by starting with the target potential from the outset, thereby reducing the computational time required. Figure 1 shows an example of initial conditions of the normalized ion and electron number density, and electric potential.

The initial electron velocities are sampled from a Maxwellian distribution with an electron temperature (⁠ $Te$⁠) of 1500 K. The random motion of the electrons gives rise to thermal fluctuations, which manifest as numerical noise in the simulation. These fluctuations effectively act as numerical collisions between particles, leading to the numerical damping of the soliton. We conducted initial studies to investigate the sensitivity of this numerical damping to the number of particles per cell (ppc) used in the simulation. The results of these preliminary studies prompted the use of a high ppc value of $2.75×105$ to mitigate the effects of numerical damping and minimize disturbances to the soliton structure, an approach similar to that employed by Kakad et al.⁵⁸ A shape factor of the 3rd order was used to minimize the effect of noise. The simulations were run on NERSC's Perlmutter, utilizing 32 nodes with 4 GPUs per node.

The evolution of the initial solitons for different values of N is depicted in Fig. 2. In all cases, the soliton initially grows until it reaches a saturated state, as exemplified in Fig. 3 for N = 0.5. The particles in the simulation self-adjust until the soliton satisfies the full set of nonlinear equations, which is consistent with the results presented in previous studies.^56,57 Figure 4 demonstrates an approximately linear relationship between the starting and final amplitudes, similar to the findings of Sharma et al.⁵⁶ However, the change in amplitude is found to be higher than what is reported in their work. Furthermore, the growth rate of the soliton from its initial state to its final state depends on the initial amplitude, with smaller amplitude solitons (e.g., N = 0.1) growing to their saturated state much slower than larger amplitude solitons.

To further investigate the soliton propagation, a 2D FFT of the electric field data was performed, revealing the presence of plasma waves and ion acoustic waves in addition to the soliton propagation. This observation is in agreement with the findings of Kakad et al.⁵⁸ and Sharma et al.⁵⁶  Figure 5 presents an example of the power spectrum with N = 0.5.

In every simulation, it is observed that a tail is produced behind the soliton as predicted by Keener and McLaughlin,⁶³ Karpman and Maslove,⁶⁴ and Watanabe.⁶⁵ This tail grows in size until it eventually separates from the main soliton body and propagates at a lower amplitude and speed compared to the soliton. An example of this is shown in Fig. 6. The separation happens quicker for higher amplitude solitons.

Sagdeev's solution for $ϕ$ is obtained by solving Eq. (14) using the fourth-order Runge–Kutta (RK4) method. Figure 7 compares the results of solving (14) with the KdV solution and PIC simulation data after reaching saturation. The amplitude and position of the soliton were set to match the PIC data. For small amplitudes, Sagdeev's model and the KdV model show close agreement with each other; however, neither model accurately captures the steepening of the soliton potential (i.e., the width). As the amplitude increases, Sagdeev's model accurately describes the electric potential while the KdV model still cannot capture the steepening. The steepening of the soliton is a function of its nonlinearity. This suggests that at smaller amplitudes, there are nonlinearities not captured by either the Sagdeev or KdV model. As the amplitude grows, the standard higher order fluid nonlinearity begins to dominate, and the unmodeled nonlinearity has a diminished impact, allowing Sagdeev's model to accurately capture the electric potential. In Sec. III B 1, we propose that this unmodeled nonlinearity arises from electron trapping effects. To maintain consistency with previous studies,^56,57 we initialized our simulations using KdV solitons. However, we also discovered that initializing solitons using Sagdeev's equations resulted in the same saturated soliton state as the one obtained from the KdV solitons. For the sake of conciseness, these additional results are not presented in this work.

The resulting $ϕ$ can then be substituted into Eqs. (11) and (10) to calculate $ni$⁠. Additionally, M can be solved based on the amplitude of the electric potential by using (25). Figure 8 presents the ion density plots, revealing that the PIC soliton has evolved from the KdV profile to a new state. Here, the amplitude of the ion density is not fixed to the PIC result but is calculated self-consistently based on $ϕ$ and M. The KdV profile fails to capture the steepening and amplitude of this ion density, with the discrepancy increasing as the amplitude grows. At smaller amplitudes, Sagdeev's model accurately captures the amplitude, but not the width, consistent with the results in Fig. 7. However, as the amplitude increases, the opposite effect is observed—Sagdeev's model captures the steepening well, but not the amplitude. As previously stated, this suggests the presence of an unmodeled nonlinearity that affects the steepening at smaller amplitudes but has a diminishing effect on the soliton width at higher amplitudes. In terms of the height, this unmodeled nonlinearity appears to have an increasing effect on the ion density amplitude at higher soliton amplitudes. The ion density given by the PIC simulation is lower than what is expected in Sagdeev's model. The lower ion density is attributed to a non-Boltzmann distribution of electrons at a lower electron density, leading to a reduction in ion density to compensate for this effect. Section III B 1 on electron trapping provides more details on this phenomenon.

In deriving Eq. (14), an isothermal electron equation of state (i.e., the Boltzmann relation) was assumed. This assumption aligns with the hybrid PIC simulations carried out by Sharma et al.⁵⁶ and Qi et al.,⁵⁷ where electrons are treated as a Boltzmann fluid. However, Schamel⁶⁶ demonstrated that electron trapping within the soliton potential well leads to deviations from the Boltzmann relation for the electron density. This deviation introduces a new nonlinear term that is not accounted for in Sagdeev's method. Trapped particles are those that resonate with and move along with the potential well of the soliton. In phase space, these particles manifest as either a hump or a hole. Our simulations clearly exhibit electron trapping. Figure 9 illustrates the propagation of electron vortices in the x–v (phase) space, with each snapshot corresponding to a different time stage in the simulation for N = 0.5. As the soliton amplitude increases, the potential well deepens, enabling the trapping of more electrons. This trapping phenomenon creates a nonlinearity that causes the electron density to deviate from the Boltzmann relation. Figure 10 shows this deviation where the observed electron density is lower than the value predicted by the Boltzmann relation using the computed $ϕ$⁠.

Electron trapping effects manifest through the electron density. By forgoing the assumption of isothermal electrons and instead utilizing the numerical electron density, we can capture the nonlinearity that is not represented in Fig. 7. Figure 11 illustrates this result, where a modified version of Eq. (14), incorporating a general electron density, is solved using the electron density provided by the PIC simulation. This approach captures the width of the electric potential more accurately than in Fig. 7. The improved accuracy is due to the inclusion of the nonlinearity arising from electron trapping in the numerical electron density, which effectively models the steepening of the soliton. Solving Poisson's equation with this electron density and potential should also yield the correct ion density with the appropriate amplitude and width. However, it is important to note that this approach relies on the numerical result for electron density from the PIC simulation and does not employ an analytical model.

In Eq. (17), when $|b|≈O(N)$⁠, the strength of the trapping nonlinearity is on the same order as the usual nonlinearity. Furthermore, the nonlinearity due to trapping becomes dominant when $|b|≫O(N)$⁠.⁶⁷ Note that setting $β=1$ would simply yield the standard KdV solution. Figures 12 and 13 present the comparisons of using Eq. (17). These figures demonstrate that the $sech4$ solution given by the modified KdV equation provides a better fit for the data compared to the $sech2$ formulation of the KdV solution shown in Figs. 7 and 8. In Figs. 12 and 13, the soliton amplitudes were set to match the PIC result to directly compare the shape of the profile even though, strictly speaking, for the small amplitude limit, the amplitudes for the electric potential and ion densities are given by the same value. For smaller amplitudes, this model is able to accurately capture the soliton profile, while the error in the width increases slightly for larger amplitudes. This behavior is expected since the model is based on a small amplitude approximation. It is important to note that to achieve the best fit, the $β$ values used (shown in Table I) were different from those presented in Fig. 10. This implies that the “true” $β$ value does not produce the “best” solution. This discrepancy suggests that while the $sech4$ function is the best approximation to the soliton, the underlying partial differential equation (PDE) that describes it may be imprecise or incomplete. This could also be a limitation of using a KdV-type approximation, but it is interesting to note that the model still captures the shape of the soliton quite well. The determination of consistent $β$ values to describe all the soliton parameters is left for future work.

The KdV formulation increasingly deviates from the saturated soliton as the amplitude rises. While it is known that the KdV formulation begins to break down at higher amplitudes, the deviation presented in this paper is more significant than previously reported.^56–58 From an engineering perspective focused on designing sensors for soliton detection,^47,52,68 this deviation might seem relatively minor. However, in our effort to understand the fundamental physics of soliton propagation, stability, and damping, we propose that the nonlinearity due to electron trapping must be accounted for in the model. In summary, the modified KdV equation, which incorporates the effects of electron trapping through the $sech4$ solution, provides a better approximation to the soliton profile compared to the standard KdV equation. However, the discrepancy in the $β$ values used to obtain the best fit suggests that further refinements to the underlying PDE may be necessary to fully capture the soliton dynamics in the presence of electron trapping.

We saw electron trapping in the phase space shown in Fig. 9. As the soliton well reaches a saturated state, the regions of lower electron density are filled, and the phase space becomes mixed, remaining invariant for the remainder of the simulation. This phase space smearing effect at saturation was predicted by Meiss and Morrison,⁶⁰ who also anticipated that the soliton amplitude would oscillate at the electron bounce frequency (also known as the trapping frequency) $ωbe=[(ek2ϕ)/me]1/2$⁠. Here, we are using $k=2π/soliton width$⁠. To verify this prediction, we performed an FFT analysis on the soliton amplitude data presented in Fig. 3. The resulting frequency spectrum, shown in Fig. 14, reveals that the soliton amplitude indeed oscillates at a frequency close to the electron bounce frequency, $ωbe$⁠. This finding corroborates the theoretical predictions made by Meiss and Morrison, providing further insights into the dynamics of electron trapping and soliton propagation in the phase space.

The KdV equation predicts that normalized potential, ion density, and velocity should be equivalent in soliton solutions. However, our study shows results that deviates from this expectation for larger amplitudes. Table I presents the quantitative relationships between soliton starting amplitude, final ion density, final electric potential, and corresponding $β$ values. Here, we see discrepancies in the growth patterns of ion density and electric potential. While ion density exhibited substantial increases (but not as much as predicted by Sagdeev), the electric potential remained relatively constant, approximating its initial amplitude. This suggests a weaker coupling between these parameters than previously theorized and provides further evidence for the breakdown of the conventional KdV relationship under the conditions examined. Furthermore, the electric potential profile, like the ion density profile, demonstrates a better fit to a Schamel soliton model compared to the KdV soliton. One interesting thing to note here is the divergence in $β$ values that optimally fit the electric potential and ion density data for a given amplitude N. Table I delineates the distinct $β$ values for ion density and electric potential profiles.

The soliton speed remains relatively constant throughout the simulation, as evident in Fig. 5, where the slope of the bright yellow band (or the black dotted line) represents the soliton's velocity. Interestingly, the soliton speed reaches a saturated state much faster than the ion density and electric potential, almost immediately after the simulation commences. To compare the observed soliton velocity against theoretical predictions, we can use the quantities presented in Table I and perform three comparisons with (1) KdV, (2) Schamel, and (3) Sagdeev.

The comparison between the observed soliton speed from the PIC simulation and the theoretical predictions is presented in Fig. 15. For small-amplitude solitons, the KdV, the full Schamel, and Sagdeev models show relatively good agreement with the PIC soliton speed. However, as the soliton amplitude increases, this agreement diminishes. Among the three models, the Sagdeev model and the full Schamel model appears to provide the closest match to the PIC soliton speed and the PIC soliton always remains to lay in between these two models. As the amplitude increases, the trend seems to suggest that the soliton speed is better captured by the full Schamel model but further work is needed to verify this. Interestingly, when using the ion density to calculate the soliton speed, the theoretical models tend to overshoot the observed speed, while using the electric potential leads to an underestimation. Consequently, the PIC soliton speed is sandwiched between these two estimates in all cases. It is noteworthy that the mKdV Schamel formulation predicts a significantly higher speed than what is observed in the PIC simulation. This is particularly intriguing because the Schamel soliton provides a good fit to the saturated soliton shape, yet it does not accurately model the soliton's speed. This discrepancy is probably due to using inconsistent $β$ values and provides further evidence that the “best” fit $β$ is not the “true” $β$⁠. This mismatch could also suggest that the inertia created due to electron trapping may be greater than what has been theorized by Schamel,⁶⁶ resulting in a slower soliton speed. The theoretical models predict that the soliton speed increases with increasing soliton amplitude, and this trend is indeed observed in the simulation results. However, the dynamic evolution of the soliton speed as a function of the soliton amplitude is not captured in these findings. The results presented here focus on the final, saturated state of the soliton and do not provide insights into the transient behavior of the soliton speed during the growth phase. The soliton speed reaches saturation almost immediately after the simulation starts and therefore, to fully understand the relationship between soliton amplitude and speed, further investigation into the temporal dynamics of the soliton propagation is necessary. This would require a detailed analysis of the soliton speed at different stages of its evolution at much smaller timescales, from the initial formation to the final saturated state.

The investigation of soliton damping is of particular significance to the small debris research community. The concept of utilizing solitons for debris detection^47,52 relies heavily on whether the soliton persists long enough to be detected or if it damps away prematurely. Moreover, the damping rate can serve as an additional data point to help uniquely identify the characteristics of the debris, such as its size, material composition, charge, and velocity. Consequently, this section will concisely explore the effects of soliton damping.

Ion Landau damping can occur when a wave has a slow phase velocity that matches the thermal velocity of ions. The ion acoustic wave, in particular, is significantly influenced by Landau damping. If the electron temperature (⁠ $Te$⁠) is less than or equal to the ion temperature (⁠ $Ti$⁠), the phase velocity lies in the region where the velocity distribution function has a negative slope. Consequently, ion Landau damping occurs when $Te≤Ti$⁠. However, this condition is not met in our simulations. Table II presents the final ion temperature at the end of each simulation for various starting amplitudes. The ion temperature exhibits a slight increase from 0, with larger amplitudes resulting in a more pronounced increase. Nevertheless, the ion temperature remains substantially lower than the electron temperature. Thus, under the conditions investigated in our study, we conclude that ion Landau damping does not play a significant role.

Ions are not the only possible resonant particles. Ott and Sudan⁵⁹ were the first to study linear electron Landau damping of ion-acoustic solitons. They derived a modified Korteweg–de Vries (KdV) equation with a source term that models the lowest-order effects of resonant electrons. Their equation includes the lowest-order nonlinear terms and terms corresponding to the linear dispersion relation for ion-acoustic waves. The derivation is justified by a formal procedure in which these terms in the equation are of the same order. Ott and Sudan obtained a solution to this modified KdV equation, which predicts an amplitude that decays to zero over time. However, they neglected trapped-particle effects, which are of the same order as the linear Landau-damping terms included in their modified KdV equation.

To address this gap, Meiss and Morrion⁶⁰ derived a new time-dependent solution that incorporates effects associated with electron orbits from electron trapping. Their work predicts that nonlinear effects from trapped particles eventually balance out the damping from untrapped particles, leading to saturation at a finite amplitude. This is consistent with our observed results. They predicted that N ∼  $(me/mi)2$ is an effective threshold for the effect of electron linear Landau damping, which is met for all cases in this study. Above this threshold, the effects from electron Landau damping are negligible. They also showed that the soliton amplitude oscillates and decreases from its initial value before reaching a saturated amplitude, with this effect becoming more pronounced as $γL/ωbe$ increases. In our simulations thus far, this ratio has been approximately zero, and therefore, the dip and oscillation of the amplitude are not observed.

To isolate this effect, we performed another simulation with $me/mi=N=0.1$⁠, resulting in $γL/ωbe=0.18$⁠, where the nonlinear electron Landau damping should be more visible. The result of this simulation is shown in Fig. 16 and is compared against the analytical solution given by Meiss and Morrison.⁶⁰ In this figure, $ψ(t)$ represents the soliton ion density amplitude as a function of time. Since the ion mass is only ten times heavier than the electron mass, the noise from random electron fluctuations is amplified by the random ion motion. Therefore, the data were averaged to reduce noise. The final saturation amplitude differs from the theoretical prediction, but the overall behavior of the soliton amplitude, including the dip and oscillation, seems to be captured. It is also interesting to note that the amplitude experiences a sudden dip right after the simulation starts and then climbs back up before slowly dipping again. We believe this is a result of the particles in the simulation quickly readjusting themselves after the start of the simulation to fit the prescribed fields. An explanation of the differing final amplitudes will be left for future work.

As previously mentioned, for the data presented in Fig. 2, electron Landau damping effects are not observed. However, we wanted to investigate whether any effects could be detected if the simulation was run for a significantly longer time. As a test case, we ran the simulation for N = 0.5 for $γLt>5$ (⁠ $ωpit>325$⁠). Figure 17 shows a small damping of the amplitude on this time scale. We believe that this damping is due to numerical noise acting as collisions, which can damp the wave.

In this study, we found that the soliton's behavior is quite sensitive to simulation parameters such as the number of particles, cell size, time step, domain length, etc. Therefore, these parameters were carefully chosen through multiple iterations until convergence was achieved. However, there is always the possibility that some of these results may be influenced by numerical effects, and the authors acknowledge this possibility.

This comprehensive study, employing fully kinetic particle-in-cell (PIC) simulations, provides valuable insights into the behavior of ion acoustic solitary waves (solitons). The key findings reveal that solitons initialized according to the Korteweg–de Vries (KdV) solution undergo nonlinear evolution, reaching a saturated state at a higher amplitude that deviates from the KdV soliton for ion density and electric potential, as well as from the Boltzmann relation for electron density. This deviation intensifies with increasing soliton amplitude, accompanied by stronger nonlinear effects arising from electron trapping. The KdV model fails to capture the width of the soliton for all amplitudes and cannot model the different valued amplitudes for the ion density and electric potential. Solving the full set of nonlinear fluid equations (Sagdeev's model) allows for the approximation of the electric potential and ion density for small amplitudes, albeit with a slightly inaccurate width. As the amplitude grows, Sagdeev's model accurately captures the electric potential and the width of the ion density, but not its amplitude. This suggests that both the KdV model and Sagdeev's model fail to capture a nonlinearity that causes these discrepancies. We demonstrated that these discrepancies primarily arise from incorrectly assuming isothermal electrons with a Boltzmann description for the electron density. The deviation from the Boltzmann relation increases with increasing soliton amplitude and is caused by electron trapping effects. Electron trapping inside the soliton was observed through the phase space, and it was shown that this trapping induces an oscillation of the soliton amplitude at the electron bounce frequency. In correcting the electron density by using the numerical density in solving the full set of nonlinear fluid equations, we were able to model the soliton more accurately. Schamel derived an equation for the electron density under the influence of electron trapping, which includes a parameter $β$ to describe the additional nonlinearity from trapping. We found that there exists a $β$ value such that the electron density closely resembles the electron density from the PIC simulation. This $β$ decreases away from 1 as the soliton amplitude increases, suggesting increasing nonlinearity. Furthermore, Schamel derived a modified KdV (mKdV) equation, which is a small amplitude approximation of Sagdeev's model with his electron density model instead of a Boltzmann description. We discovered that there exists a $β$ value that enables the mKdV model to capture the soliton shape much better than the KdV model. This approximation works well for all amplitudes considered in this study, although it performs best for smaller amplitudes, as expected from a small amplitude approximation. However, like the KdV model, the mKdV model also assumes the same amplitude for electric potential and ion density, thus, failing to capture the saturated state. Interestingly, the higher amplitude ion density can still be approximated by the mKdV model by using a different $β$ value, suggesting that the final state of the soliton is best approximated by an analytical form of $sech4$⁠.

The soliton velocity was compared against predictions from the KdV, Schamel, and Sagdeev models. For small amplitudes, all models agree relatively well with the PIC soliton speed, but the agreement diminishes as amplitude increases. The Sagdeev and the full Schamel models provide the closest match, with the PIC speed lying between them. The mKdV Schamel formulation predicts significantly higher speeds than observed, despite the Schamel soliton fitting the saturated shape well. This discrepancy likely stems from using inconsistent $β$ values. Although the models predict increasing speed with amplitude, which is observed in the simulations, the dynamic evolution of the speed is not captured in these findings, as they focus on the saturated state. Further investigation into the temporal dynamics at smaller timescales is necessary to fully understand the amplitude-speed relationship.

While the differences between the models mentioned above might seem relatively minor from an engineering perspective, we believe that these discrepancies highlight the limitations of widely used models, such as the KdV model and Sagdeev's model, in capturing the fundamental physics of soliton formation, stability, and propagation. The primary reason for these limitations is that these models do not account for the nonlinearity arising from electron trapping. We conclude that Schamel's model, which incorporates the full set of nonlinear fluid equations, is a more suitable alternative to Sagdeev's model, and the modified-KdV (mKdV) equation is a better alternative to the standard KdV model. However, even Schamel's model has its limitations. These include the inability to accurately predict the $β$ parameter and the inconsistencies associated with requiring different $β$ values for modeling the electric potential and ion density profiles. Furthermore, the mKdV model fails to accurately capture the differing amplitudes between the electric potential and ion density, as well as the soliton speed. Additionally, Schamel's model does not include other potential nonlinearities, such as ion trapping, ion reflection, and radiation of the soliton. In light of these findings, we recommend that any future model development for soliton propagation should take electron trapping effects into account. By incorporating these nonlinear effects, more accurate and comprehensive models can be developed to better understand the complex dynamics of solitons in plasma environments.

Our study also examines Landau damping effects and concludes that both ion and electron Landau damping appear to be minimal under the simulated conditions. By using a smaller ion mass, we were able to better observe the effects of electron Landau damping as predicted by Meiss and Morrison. While we observed saturation at a smaller amplitude than the initial value and an oscillation of the amplitude, the final predicted amplitude did not match our simulation. Nonetheless, the possibility of numerical effects influencing the results is acknowledged, emphasizing the importance of careful parameter selection to mitigate numerical effects. The conditions used in these simulations were not conducive to Landau damping. Although we used realistic electron temperatures that mimic LEO (1500 K) conditions, using a cold ion species is not representative of a general plasma. Similarly, using a fictitious ion mass could also create artificial results. We anticipate that the use of realistic mass ratios in our simulations will lead to several effects: an increase in electron trapping, a decrease in electron Landau damping, and the bounce frequency remaining the same. Exploring the effects of Landau damping in more depth for a general plasma with realistic temperature ratios will be the focus of our next study.

In conclusion, this work provides a comprehensive, first-principles investigation into the propagation of KdV solitons using fully kinetic PIC simulations. The findings highlight the importance of nonlinear electron trapping effects in soliton dynamics and the limitations of conventional models such as KdV, Sagdeev's model and hybrid PIC simulations in capturing these phenomena. The insights gained from this study contribute to a deeper understanding of soliton behavior in plasma environments and lay the groundwork for further research into the fundamental physics of soliton propagation, shape, growth, stability, and damping.

The authors would like to express their sincere gratitude to Dr. Remi Lehe for his invaluable support in running simulations on WarpX and for providing insightful guidance during the data analysis process. His expertise and assistance have been instrumental in the successful completion of this study.

We would also like to extend our heartfelt appreciation to Dr. Sarvershar Sharma, Dr. Sudip Sengupta, Dr. Gurudas Ganguli, and Dr. Abhijit Sen for their thought-provoking discussions and invaluable insights. Their contributions have significantly enriched this work, and their knowledge and expertise have been essential in shaping the direction and quality of our research.

This work was supported by MIT Lincoln Laboratory in collaboration with IARPA's SINTRA program.

This research was supported in part by the Knight-Hennessy Fellowship program at Stanford University. This research was supported in part by the Exascale Computing Project (No. 17-SC-20-SC), a collaborative effort of the U.S. Department of Energy (DOE) Office of Science and the National Nuclear Security Administration; the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Exascale Computing Project under Contract No. DE-AC02-05CH11231.

This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility using NERSC Award No. NERSC DDR-ERCAP0030179.

The authors have no conflicts to disclose.

Ashwyn Sam: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead). Prabhat Kumar: Data curation (supporting); Methodology (supporting); Software (supporting); Visualization (supporting). Alex C. Fletcher: Conceptualization (equal); Formal analysis (supporting); Investigation (supporting); Methodology (supporting); Supervision (supporting); Validation (supporting); Visualization (supporting); Writing – review & editing (equal). Chris Crabtree: Conceptualization (equal); Formal analysis (supporting); Investigation (supporting); Methodology (supporting); Supervision (equal); Validation (supporting); Writing – review & editing (equal). Nicolas Lee: Funding acquisition (equal); Project administration (supporting); Writing – review & editing (supporting). Sigrid Elschot: Conceptualization (supporting); Funding acquisition (lead); Project administration (lead); Supervision (lead); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.