Nonlinear spectral analysis of ion acoustic solitons arising from a streaming charged object using the numerical inverse scattering transform

The Korteweg–De Vries (KdV) equation is a model of nonlinear collisionless waves in plasmas moving near the ion acoustic velocity.¹ It supports long term coherent structures, “solitons,” through balancing nonlinear and dispersive effects. Remarkably, these fully nonlinear waves are stable under collisons with other solitons,² which results in their existence for long time scales. When forcing is included in the KdV equations, two types of solitons can arise from a streaming disturbance:^3,4 pinned solitons which stay attached to the source and precursor solitons which propagate ahead of the source. These properties make solitons, especially precursor solitons, an attractive indirect sensing mechanism of generating sources.

The problem of detecting soliton sources naturally splits into two subproblems: automatically detecting the solitons and inferring properties of the source. This publication is concerned with the former. The latter problem needs to consider the soliton changes due to inhomogeneities of the plasma environments, damping rates, direction of propagation in all three dimensions, sensitivity of source shape to generated soliton amplitude, etc. Furthermore, the automatic detection method outlined in this paper cannot distinguish between precursor and pinned solitons. Some parts of the inference problem have been considered in Refs. 5 and 6. However, an inverse source model from measurements is still an open problem.

Spacecraft has observed electrostatic solitary waves (ESWs) at the Earth–Sun magnetopause,⁷ in front of and in an interplanetary shock,⁸ in the auoral regions,⁹ and other magnetic interface regions.¹⁰ Some of these ESW are modeled by the KdV equation and are direct observations of KdV-type solitons.^11,12 It is hypothesized that similar solitons could be generated by charged orbiting debris^5,6,13–15 as well. Detection of solitons could provide a mechanism to indirectly sense orbital debris in the ionosphere, including sub-cm debris that is too small to be detected using conventional (optical and radar based) methods.⁵ Such sub-centimeter debris poses a hazard to human presence in low Earth orbit on board space stations and manned space capsules due to high orbital speeds $v≈7 km/s$ and correspondingly high impact energies. There have been documented cases of Space Shuttle windows being cracked in orbit,¹⁶ requiring frequent replacement and large damage to solar panels on spacecraft.¹⁷ 

In spacecraft observations of ESW, solitons manifest as broadband data on a spectrogram even though they are a single coherent mode¹⁰ due to their large amplitude and correspondingly high nonlinearity. The solution of a partial differential equation is only a superposition of Fourier modes if it is linear in the solution. Using a traditional spectogram implicitly supposes that all coherent modes are Fourier modes. Previous attempts at dealing with coherent wave nonlinearity have relied upon statistical correlation between Fourier modes, termed bispectral analysis.^18,19 This method has been applied to sinusoidally and cnoidally driven KdV systems.²⁰ However, bispectral analysis still decomposes the soliton into many harmonics.

The homogeneous KdV equation, even though it is nonlinear, supports an exact spectra found through the inverse scattering transform (IST).^21,22 The spectra consists of a series of cnoidal waves. The IST reduces to the Fourier transform in the limit of small nonlinearity.²³ An IST filter could be used to screen for solitons from data in much the same way a spectrogram can isolate linear modes. In this paper, we use a numerical implementation of the IST as a spectral analysis tool to identify solitons from simulated data. This technique is demonstrated to be robust to noise.

The IST was created for the homogeneous KdV equation and has only recently been analytically extended to certain sources of the forced KdV (fKdV) equation.^24–27 Streaming localized sources produce solitary structures that are similar in form to the homogeneous solutions.³ For many of these sources, the IST decomposition is still a series of cnoidal waves. However, the mode's time evolution may be different from the homogeneous case. This paper demonstrates for the first time that the numerical IST can identify solitons in the forced case. Furthermore, the homogeneous IST has helped interpret the low frequency power spectrum and filter high frequency noise in a physically meaningful way in the nonlinear ocean waves community^28,29 because the KdV equation is also a model of surface water waves in the shallow depth limit. We believe that the method outlined in this paper can be used in the plasma physics community to filter and understand ion acoustic solitons, especially in the context of automatic detection by spacecraft and ion acoustic turbulence experiments, in a similar manner. This need has been remarked upon in hollow cathode experiments.³⁰ In this paper, we demonstrate that the IST is a better interpretable spectral technique than the Fourier transform for a simulated ion acoustic soliton system. We provide the source code and dataset used in this publication for the IST with the aim of it being used in future analysis (see Data Availability section).

For any fluid system where the nonlinear term is maintained and the dispersion relation includes a cubic term, the KdV equation will arise from a small frequency expansion. This has been found for a one-dimensional (1D) single species plasma system of cold ions and warm electrons in the limit of $ω≪vIA/λD$⁠, where ω is the frequency, $vIA=Te/mi$ is the ion acoustic velocity, T_e is the electron temperature in energy units, m_i is the ion mass, and λ_D is the Debye length. From here on, all quantities are nondimensionalized by parameters from the linearized dispersion relation as follows: $U→U/vIA, x→x/λD, t→tfp$⁠, and $ϕ→eϕ/Te$⁠, where $ϕ$ is the electric potential, f_p is the plasma frequency in nonangular units, and e is the electron charge. An important task in application of this method will be accurate estimation of these plasma parameters in order to obtain an accurate nondimensional quantity for filtering. However, a proper discussion of this problem is left for a future publication on a specific measurement technique. In this publication, estimates of plasma parameters come from the International Reference Ionosphere.³¹ The governing equations for perturbations along magnetic field lines in a two temperature plasma with $Te≫Ti$ are given by Eqs. (1)–(4), where n is the ion density, n_e is the electron density, and U is the ion bulk velocity,

The full derivation of the KdV equation is summarized well in the literature¹ or numerous textbooks.^32–34 All quantities are expanded to second order in a small parameter ϵ, where $ω2=ϵμ2$ and μ is the rescaled frequency. The zeroth order densities are all 1 and $U(0)=ϕ(0)=0$⁠,

It is found that the first-order expansion terms are linked: $ne(1)=n(1)=ϕ(1)=U(1)$⁠. Thus, densities and electric potential vary together. These first-order terms obey the same evolution equation, which is the KdV equation. This is given by Eq. (9), where $u=U(1)$ and t and x are scaled coordinates from the original ion and electron equations of motion. For the KdV equation governing the ion acoustic wave in a frame moving at the ion sound speed in an collisonless plasma of cold ions and warm electrons, α = 1 and $β=1/2$⁠,

The coefficients can be changed through scaling of the independent variables, so any equation of this general form is considered the KdV equation. In the following IST analysis, a new term $λ=α/(6β)$ is introduced to account for varying coefficients. λ is the nonlinearity parameter since it describes the ratio of nonlinear to dispersive terms. In order to perform IST analysis, one must assume a value of λ for the system. In the target application of soliton detection in the vicinity of Earth, a more realistic set of values for these coefficients are $α=1.0051$ and $β=0.4925$⁠. A full discussion of these values and a more general expression for α and β in terms of plasma parameters can be found in Sec. III A of Ref. 5. These coefficients are estimated from magnetospheric models of plasma composition. Unlike the case of the shallow water wave KdV equation where the coefficients depend strongly on the water height h, the ion acoustic KdV equation coefficients do not vary as strongly numerically.

When the nonlinear and dispersive terms exactly counteract one another, the wave is stationary. Forcing can be introduced by a free charge on the right hand side of Eq. (4). For the cases considered in this publication, S(x, t) will take the form of a moving $S(x−Vt)$ or stationary S(x) localized source,

When the same expansions are carried out with the source being of $O(ϵ2)$⁠, the result is now as Eq. (11). This is the forced KdV equation (fKdV),

As discussed at length in Ref. 5, the primary sources of damping will come from ion-neutral collisions and Landau damping. A full treatment of deriving damping times and estimated propagation distances can be found there. For centimeter and sub-centimeter sized debris induced precursor solitons are expected to propagate for tens of kilometers.

Equation (9) is exactly solvable through a transformation from the KdV to the time-independent Schrödinger equation (TISE).²¹ The TISE is interpreted as a scattering problem, where asymptotically planar waves $ψ∼ exp [±ikx]$ are modified by a localized scattering potential u(x) given by Eq. (12). Throughout this paper, we shall denote the scattering potential in the TISE as the TISE potential and the physical electrostatic potential as the electric potential. After nondimensionalization, these two quantities are actually the same through the KdV–TISE transformation, but they correspond to different uses in the two problems. The TISE is a functional eigenvalue problem in x at every t where $E=E(t)$ and $ψ=ψ(x,t)$⁠. If u(x, t) is a solution of the KdV equation, then E_t = 0. The eigenvalues of the associated $ψ(x,t)$ are the same at every time and all of the time dependence is in ψ,

In the TISE scattering problem, solving for $ψ(x,t)$ is not done directly for nontrivial TISE potentials. Rather, the problem is reformulated as a set of coefficients that relate the asymptotic solutions from the left and right of the scattering potential. These coefficients are functions of $E=k2$ and t. E is often called the energy level because in quantum mechanics E is the physical energy level of a wavefunction. There is no physical correspondence to energy here. One common set of coefficients consists of a—the transmission coefficient and b—the reflection coefficient. Under the condition that u(x, t) obeys the KdV equation and is localized (⁠$limx→±∞u=0$⁠), the time dependence of the scattering coefficients is explicitly given by Eqs. (13) and (14).³⁵ The process of finding the rooted eigenvalues k, and coefficients $a(k,0)$ and $b(k,0)$ for an initial value to the KdV equation $u(x,0)$ is known as forward scattering,

The scattering data and eigenvalues are sufficient to reconstruct the scattering potential through the Gelfand–Levitan–Marchenko equation. This process is called inverse scattering and is not used directly in this publication. Together the inverse scattering, forward scattering, and scattering data time integration is called the IST. This opens up the possibility of extending the numerical implementations of the IST to time integration of the KdV equation with sources and will be explored in the future work.

In the forced case, the time evolution of the scattering coefficients is not known. However, it has been shown that sources of certain analytical form are also solvable by the IST albeit with a different time evolution equation.^24–27,36

The Osborne algorithm²³ is used to compute the IST of a discrete signal with periodic boundary conditions. Periodic boundary conditions are natural for signal processing. The numerical implementation of the IST will be denoted the PIST to be consistent with other literature.³⁷ The term for ISTs for all nonlinear integrable equation is the nonlinear Fourier transform (NFT).³⁸ Note that some authors denote PIST as the KdV-nonlinear Fourier transform (KdV-NLFT).²⁹ PIST relies on a transfer matrix formulation of the TISE, defined in Eq. (15). The scattering matrix M is a complex matrix with unity determinant which takes the role of a and b. The reflection and transmission coefficients can be defined from M by $a=|M12|2$ and $b=|M11|2$⁠,

where ψ describes the relative amplitudes of the scattered waves. ψ contains the amplitude of both the left and right going waves evaluated at an arbitrary location. Both M and ψ will take different forms depending on a choice of basis functions which is not unique. In the simplest formulation, the choice of basis is planar waves. The full solution will be a superposition of $A exp (ikx)$ and $B exp (−ikx)$⁠. We, therefore, define $ψi$ as Eq. (16). For a constant scattering potential, Eq. (17) gives the transfer matrix, where $kn=k2−u$⁠. We use the subscript notation because this treatment of the scattering matrix is conducive to a piecewise constant representation of the signal, such as in the context of a discrete digital representation of a continuous signal. The T superscript indicates the transpose,

Similarly, for across a step in TISE potential M is given by Eq. (18), where k⁻ refers to the k_n before the step and k⁺ refers to the k_n after. Using these two basic solutions, the total M can be constructed sequentially, i.e., $Mx0→x2=Mx1→x2Mx0→x1$⁠,

As detailed in Refs. 23 and 39, the planar wave choice of basis for ψ, while being the most physically intuitive basis, is not the simplest one to deal with algebraically nor numerically. A computationally simpler basis is $[1,0]T$ and $[0,1]T$⁠.⁴⁰ Now the step transfer matrix is given by Eq. (19), where k_n has the same definition as above,

k_n is complex, but lies on either the real or imaginary axes. $sin (iRe(x))=isinh(Re(x))$ and $cos (iRe(x))=cosh(Re(x))$⁠, so the following substitutions can be made if $Im(kn)≠0$⁠:

This allows for the evaluation of the forward scattering problem without the need for complex algebra by checking the sign of $kn2$ at every loop iteration. M is now real.

Recall that $M=M(k)$⁠. It has been shown that for the periodic boundary conditions in the KdV equation, the associated TISE scattering problem can be simplified by using Bloch states.⁴¹ The result, given by Eq. (23), is that the solution to the KdV equation can be written as a finite linear superposition of hyperelliptic modes given by the zeros of the transfer matrix. $μj(x;x0,t)$ are hyperelliptic functions that contain all of the space and time dependence. E_j satisfy $(1/2)trM(k=Ej)±1=0$ and $M12(k=μj)=0$⁠. $N∈ℕ$ is the number of “degrees of freedom” in the problem,

The time evolution of the hyperelliptic functions is given by a coupled ODE system involving all of the μ_j and E_j at time t. It is omitted from this paper since it is not used. Each $μj(x,t)$ can be assigned a modulus $mj=(E2j+1−E2j)/(E2j+1−E2j−1)$ and an amplitude $aj=E2j+1−E2j$⁠. Each soliton will manifest itself as a separate μ_j in this decomposition with m_j = 1. When m_j is calculated numerically, if $mj≥0.999$ that mode is considered a soliton. This is an arbitrary choice based on numerical precision and benchmarking cases. All other modes are “radiation” modes. Thus the forward scattering transform reduces to finding the zeros of $trM(k)±1$ and $M12(k)$ with the initial data.

Given a discrete signal u_n which is of length N, the associated transfer matrix M for energy level $E=k2$ can be calculated with Algorithm 1. An intermediate scattering matrix at every piecewise step T is used to accumulate the signal. This is defined in Algorithm 1. The implementation of $M=TM$ uses Ballard's algorithm⁴² for improved numerical stability because $sinhx$ and $coshx$ both increase very quickly for small x. Matrix multiplication is done with CBLAS level 3 matrix multiplication subroutines. The M matrix calculation is implemented in C, while the overall zero finding and spectra construction passes are done in matlab. The matlab coder functionality is used to call the C subroutines.

Finding the zeros of M₁₂ and $trM±1$ is done with a four pass method outlined in Osborne.²³ The method relies on the fact that between each zero of M₁₁ $=Ej$⁠, there are one or two μ_i and between each E_i and $Ei+1$⁠, there is one zero of $tr(M)+1$ and one zero of $tr(M)−1$⁠. In the first pass, an accounting function $S(xn,E)=∑i=0n|sign(M11(xi+1,E))−sign(M11(xi,E))|/2$ is evaluated on a fixed resolution E grid from E_min to E_max. $S(xn,k)$ gives the number of sign changes in M₁₁ from x₁ to x_i which is assumed to be the same as the number of zeros of M₁₁ on a sufficiently fine E grid. E_min and E_max are chosen such that the graph of $M11(E)$ has no zeros left of E_min and E_max is in a region where $|M11|≤1$⁠. Osborne²³ discusses some theoretical properties analogous to Nyquist sampling theorems which constrain E_min and E_max. This application focuses on low frequency modes, so the sampling frequency is not rigorously considered. In the second pass, the exact zeros of M₁₁ are found using an optimization solver. The location of the sign change in S is used as an initial guess. The third pass finds the locations of μ_i through a bisection method root solver with $M11,i$ and $M11,i+1$ as the bounds. The sign of $M11,i$ and $M11,i+1$ is checked. If they are dissimilar, there is assumed to be one root, and the bisection solver is directly used. If the signs are the same, it is presumed that there are two roots. A third point of a dissimilar sign is searched for in-between the two bounds by first guessing the midpoint and then randomly guessing in that interval. Once a third point of dissimilar sign is found at E₃, two bisection searches between $[M11,i,E3]$ and $[E3,M11,i+1]$ are conducted. The fourth pass finds the zeros of $tr(M(E))+1$ and $tr(M(E))−1$ between each μ_i and $μi+1$ with two bisection searches. If $maxS=N−1$⁠, then there are N μ_i and $2N+1$ E_j. N is the number of modes in the hyperelliptic decomposition.

Figure 1 shows the intermediate calculations of the numerical IST of $u(x)=sech2(x/2)$ using the four pass process. This is an analytical solution to the homogeneous KdV equation corresponding to a single soliton moving at nondimensionalized speed c = 2. The variable speed in single soliton solution is discussed and given later by Eq. (25). The signal has 150 points from $[−12.5,87.5]$⁠. The fixed resolution E grid is 3000 points from $[−2,5]$⁠. The nonlinearity parameter λ = 1. The zero with $E≈0.5$ represents the soliton mode. All components of M vary from 10⁶⁰ to $−1027$ over a small width near $E≈−0.5$⁠. $−1027$ is a local minimum of this function, not a floating point representation of negative infinity. The zeros in that interval are not found to great precision because that function is steeper than double floating precision can resolve. However, this is acceptable for our purposes. In general, the leftmost zero is subject to this problem because all components of $M→∞$ as $E→−∞$⁠.

The solutions to the fKdV equation are solved using a pseudo-spectral fluid simulation which uses the Chan–Kerkhoven pseudospectral scheme⁴⁴ for time stepping the solution. The time evolution equation is given by Eq. (24), where k are the wavenumbers of the Fourier modes, $s=π/L, Δt$ is the time step, L is the length of the spatial domain, and $F$ represents a discrete Fourier transform. The time step $Δt$ is chosen to ensure a numerical stability condition

This has been demonstrated to be accurate to $unumerical− uanalytical≈10−4$ on known solutions to the fKdV equation.^5,6

The single soliton solution to the KdV equation is given by Eq. (25). c is the soliton speed and x₀ is the soliton initial location. For this numerical experiment, $x0=0$⁠, c = 2, λ = 1, and $x∈[−L/8,7L/8]$ with 150 points where L = 100. There are 200 time points with $t∈[0,5]$⁠. The numerical solution is visualized in Fig. 2(a). In Figs. 2(b)–2(d) k is the spectral wavenumber associated with the transforms. The PIST modes have a wavenumber associated with them. For the low modulus modes, this corresponds exactly to the wavenumber of equivalent Fourier mode. An interesting discussion of this can be found in Osborne and Bergamasco.⁴⁶ For the high modulus modes, this is related to the eigenvalue of the Schrödiner equation,

The PIST of this solution yields a single dominant mode with a modulus m > 0.999 indicating a soliton mode in Fig. 2(c). The PIST is taken at three times to indicate the spectrum is constant (i.e., $d/dt(E)=0$⁠) even though the soliton is spatially translating. $a1≈0.5$⁠, which corresponds to a u amplitude of $2a1≈1$⁠. From Eq. (23) with λ = 1 this a agrees with the analytical solution. The Fourier transform of the same signal is a constant broadband spectrum as seen in Fig. 2(b). The PIST is interpretable as a single solitary wave directly from the decomposition.

The initial condition $u(x,0)=12sech2(x)$ results in three solitons of differing heights as seen in Fig. 3(b). This spectral simulation has KdV parameters α = 6, β = 1, and λ = 1. From this numerical example forward, all data come from an fKdV simulation, in this case with zero forcing. The full spatial domain is $x∈[−300,300]$ discretized into 12 288 points. The time domain is $t∈[0,3]$⁠. The height of a soliton determines its speed, so these three solitons will spread out. PIST can identify these three solitons even when they are overlapping such as during the initial condition. Table I shows the first six modes of the PIST at t = 0. Three modes have a soliton modulus. The physical amplitude of these modes also corresponds to the amplitude of the resultant solitons in Figs. 3(a) and 3(b). Recall that $umax=2a$⁠.

Here, we demonstrate that PIST is robust to small noise motivated by ionospheric conditions. Solitary structures governed by the KdV equation will be a sech² pulse in $ϕ$⁠, not the electric field E. In one dimension $E=−(∂ϕ/∂x)$⁠, which would make a bidirectional electric field. However, we will assume that the strength of the nondimensional KdV signal $u∼ϕ∼|E|$ to discuss the relative strength of ion acoustic solitons to other plasma waves in the ionosphere. Table II shows the typical electric fields for several phenomena in the ionosphere.

The largest noise content will come from large magnitude whistler waves which are 30% as large as the ionospheric soliton. If the signal defined in Eq. (26) is superimposed on the single transiting soliton, then PIST still recovers that there is a single soliton. The first term represents a 5 kHz whistler wave with dimensionless frequency $ω′=(5 kHz)λD/vIA$⁠, where $λD=3.6$ cm and v_IA = 5.5 km/s. The second term is white noise with an amplitude one order of magnitude smaller than the whistler wave. This represents lower order noise signals,

Under this noise perturbation $u→u+unoise$⁠, PIST decomposes the signal into a soliton and many low amplitude modes as seen in Figs. 4 and 5. The question may arise, is there an effective signal to noise ratio dependent upon the signal discretization or fundamental to the IST? This is harder to establish due to the nonlinear nature of the governing equation. If one adds a significantly large noise perturbation, the noise perturbation itself will quickly break into a train of solitons. An example of this is the original Zabusky–Kruskal problem² where a large amplitude cosine wave was observed to break into many solitons. This analysis simply demonstrates that the method is robust to small perturbations which have not steepened enough for the nonlinear term to be dominant.

When localized forcing is added to the KdV equation, this models a moving debris object with a nonzero surface potential. Due to plasma sheath conditions, it is very likely that a debris object immersed in plasma will have a nonzero surface potential. For a spherical object in a Maxwellian plasma, the surface potential is determined by the relative length scales of the sphere and the Debye length of the plasma it is immersed in Ref. 50. We apply PIST to the plasma near the debris to search for automatic detection of this signature. Through our numerical experiments, we have found for the first time that the PIST decomposition finds the correct number of solitons on data from the fKdV equation with localized forcing.

Wu showed³ that a sech² choice of forcing, given by Eq. (27), with an initial condition of $u0(x)=a sech2(Kx)$ results in a stable nonmoving soliton solution where $u(x,t)=u0(x)$⁠. In the case of a zero initial condition, this forcing results in the growth of a soliton at x = 0 and then subsequent detachment once it reaches a critical amplitude. The growth and detachment cycle results in a periodic generation of solitons from the source.

The PIST is capable of detecting a soliton, even in a forced setting. In Fig. 6, the fKdV equation with a forcing term given by Eq. (27) with $b=−0.1250$ and K = 0.6124,

The discretization parameters were $x∈[−250,250]$ with 2048 points, $t∈[0,2]$ with $Δt=0.001$⁠. $α=−3/2$ and $β=−1/6$ so $λ=3/2$⁠. The fixed E grid is set from $E∈[−1.4,2]$ with $ΔE=6.8×10−5$⁠. Figure 6 shows the simulation which is spatially clipped near the nonzero signal. The generation and detachment pattern happens every $Δt≈45$⁠. The solitons then propagate in the –x direction. For analysis, truncated spatial cuts are taken at noteworthy times. t = 20, 55, and 90 are times when the precursor soliton initially detaches from the forcing term [Figs. 7(a), 7(c), and 7(d)]. At t = 40, the first soliton has already detached and the second one is forming [Fig. 7(b)]. At t = 100, a fourth soliton just starts to form and the first three are cleanly detached [Fig. 7(e)]. Figures 8(a)–8(e) show the PIST mode amplitudes and moduli. The moduli m are near m = 1, so $σ(m)$ is plotted where $σ(m)=log (1−9(m/9))/ log (1/10)$⁠. This defines a logarithmic scale which counts the number of post-decimal point nines. For example: $σ(1)=∞, σ(0.9)=1, σ(0.999)=3$⁠, and $σ(0)=0$⁠. Under the definition that a soliton is defined as any mode with m > 0.999, then $σ(m)>3$ is a soliton. This is plotted as horizontal line in Figs. 8(a)–8(e).

Figures 8(a) and 8(b) show that at t = 20 and t = 40, respectively, there are one and two solitons seen in the system by PIST. This matches the visual observation in Figs. 7(a) and 7(b). The amplitudes for the single soliton in PIST (⁠$2×0.75=1.5$⁠) is close to the amplitude of the observed signal. What is observed may be partially the superposition of the other modes of the system and the single soliton; therefore, this is good agreement. In the two soliton case [Fig. 7(c)], it is not clear that the amplitudes of the observed solitons match the PIST readings. However, the non-solitonic modes are large amplitude during this cut, as seen by the large negative of value of the soliton. A KdV soliton, by definition, must be single signed. At t = 55, there are five solitons in the decomposition, even though there are only two in the simulation. Most of the energy of the system is in the first two solitons. Similarly, at t = 90, there are six solitons even though the system only has three. The energy is mostly contained in the first three modes. After the third soliton detaches, the decomposition at t = 100 shows four solitons as expected by visually inspecting Fig. 7(e). At t = 90, the solitons are nonlinearly interacting. During the time where the PIST number of solitons does not match the simulation, the solitons lose coherence because they are interacting with each other. This suggests that in the forced setting, PIST can identify solitons as long as they are not interacting with each other.

In dusty plasmas, Gaussian forcing is used to approximate the surface charge on debris objects⁵ such as in Eq. (28) where $Φsurf$ is a surface potential determined by sheath theory. a is a characteristic size of the debris and V is the relative speed in the plasma medium. This is discussed in Truitt and Hartzell⁶ 

Figure 9 shows a simulation of such a forcing term centered at x = 0. Qualitatively, it has the same behavior of growing and detaching solitons from the forcing source as was shown in the prior examples.

The simulation has the following physical parameters: $α=1.0051, β=0.4925$⁠, a = 0.2783, $λD=3.6$ cm, v_IA = 6.1 km/s, and V = 1.1240. a and V are nondimensional quantities normalized by λ_D and v_IA, respectively. Relevant numerical parameters are $tmax=200, Δt=0.001$⁠, L = 300, $N=512×15=7680, E∈[−1.5,3]$⁠, and $ΔE=1.5×10−3$⁠. These quantities correspond to a 1 cm sized debris object in a 2000 km altitude circular Earth orbit. The variation in α and β is to account for the observed κ value of electrons in Earth magnetospheric models⁵ at that location.⁵ The data used by the PIST are clipped from x = 20.625 to x = 59.6875 to reduce the PIST run time and isolate the precursor solitons.

It can be seen that in each of these three cuts, the correct number of solitons is identified. At t = 70, a single soliton mode is detected in Fig. 10(a). At t = 130, two solitons are detected in Fig. 10(b). Finally, at t = 180, three solitons are detected in Fig. 10(c). It was observed that if the signal is not clipped close to the source and the wake behind the moving debris is part of the signal, the PIST decomposition becomes less accurate.

When two solitons are present in the frame, the PIST does not recreate the correct amplitude as seen in Fig. 10(b). One soliton mode has the correct amplitude while the other one is significantly lower than the truth. However, when the second soliton is clipped out of the frame, the IST reports the remaining soliton amplitude correctly. An example of this is shown in Fig. 11. This phenomena also occur in the unforced setting. It shall be investigated in future work.

In the use case where the IST would be used to post-process spacecraft data, multiple domain windows could be used to analyze the same data. Since the full IST identifies the correct number of solitons, the full domain filtering would still be admissible as an automated soliton detector.

This study applied a numerical nonlinear spectral analysis tool, the inverse scattering transform, to the forced and unforced Korteweg–De Vries equation in the context of an ion acoustic wave generated from a debris object. For the unforced KdV equation, the PIST provides a more useful decompositon than a Fourier transform because it correctly identifies nonlinear coherent modes as a single mode instead of broadband noise. It can also identify multiple colocated modes which visually look similar. This decomposition is shown to be stable to small perturbations in the signal. For some localized forcing functions, the PIST was found to identify the correct number of solitons, even though it is derived for homogeneous solutions to the KdV equation. Specifically, we have shown that the IST correctly identifies solitons predicted to be produced by charged orbital debris.

This material is based upon the work supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Department of Energy Computational Science Graduate Fellowship under Award No. DE-SC0021110, the Flagship Fellowship Program at the University of Maryland Graduate School, College Park, and the A. James and Alice B. Clark Foundation through the Clark Doctoral Fellows Program. The first author would like to thank Dr. James Drake and Dr. Thomas Antonsen for their class notes and discussion on the KdV equation as well as Dr. Raymond Sedwick for reviewing the manuscript. The first author also thanks the referees for their insight and comments.

This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

The authors have no conflicts to disclose.

Ian DesJardin: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). Christine M Hartzell: Conceptualization (equal); Project administration (equal); Supervision (equal); Writing – review & editing (equal).

The PIST code can be found at https://github.com/Planetary-Surfaces-and-Spacecraft-Lab/KVIST and also on Zenodo.⁵¹ The simulation code¹⁵ can be found at the following DOI: https://doi.org/10.13016/for0-xjyd. The dataset used by KVIST can be found separately on Zenodo.⁴⁵ The PIST code, simulation code, and simulation dataset can also be provided by the corresponding author upon reasonable request.