Head-on collision of nonlinear solitary solutions to Vlasov-Poisson equations

TL;DR

This study investigates the stability and collision dynamics of nonlinear solitary solutions to the Vlasov-Poisson equations in a plasma, revealing stability categories, instability mechanisms, and electron hole merging behaviors through kinetic simulations.

Contribution

It provides a comprehensive kinetic analysis of solitary solution stability, categorizes solutions by Mach number, and uncovers detailed electron hole interaction mechanisms.

Findings

Solutions are categorized as stable, semi-stable, and unstable based on Mach number.

Bipolar electron flux structures seed net-neutrality instability during collisions.

Electron holes can merge with varying degrees, exchanging phase space cores in non-integer rotations.

Abstract

Nonlinear solitary solutions to the Vlasov-Poisson set of equations are studied in order to investigate their stability by employing a fully-kinetic simulation approach. The study is carried out in the ion-acoustic regime for a collisionless, electrostatic and Maxwellian electron-ion plasma. The trapped population of electrons is modeled based on well-known Schamel distribution function. Head-on mutual collisions of nonlinear solutions are performed in order to examine their collisional stability. The findings include three major aspects: (I) These nonlinear solutions are found to be divided into three categories based on their Mach numbers, i.e. stable, semi-stable and unstable. Semi-stable solutions indicates a smooth transition from stable to unstable solutions for increasing Mach number. (II) The stability of solutions is traced back to a condition imposed on averaged velocities, i.e. net neutrality. It is shown that a bipolar structure is produced in the flux of electrons,early in the temporal evolution. This bipolar structure acts as the seed of the net-neutrality instability, which tips off the energy balance of nonlinear solution during collisions. As the Mach number increases, the amplitude of bipolar structure grows and results in a stronger instability. (III) It is established that during mutual collisions, a merging process of electron holes can happen to a variety of degrees, based on their velocity characteristics. Specifically, the number of rotations of electron holes around each other (in the merging phase) varies. Furthermore, it is observed that in case of a non-integer number of rotations, two electron holes exchange their phase space cores.

Figures

Equations

$$v\frac{\partial f_s(x,v)}{\partial x}+\frac{q_s}{m_s}E(x,t)$$

$$\frac{{\partial }^2\varphi (x)}{\partial x^2}$$

$$\int v{f}_{e}dv=\int v{f}_{i}dv.$$

$$f_s>0s=i,e.$$

$$f_{s}(v)=\left\{\begin{array}[]{lr}A\ exp\Big{[}-\big{(}\sqrt{\frac{\xi_{s}}{2}}v_{0}+\sqrt{\varepsilon(v)}\big{)}^{2}\Big{]}&\textrm{if}\left\{\begin{array}[]{lr}v<v_{0}-\sqrt{\frac{2\varepsilon_{\phi}}{m_{s}}}\\ v>v_{0}+\sqrt{\frac{2\varepsilon_{\phi}}{m_{s}}}\end{array}\right.\\ A\ exp\Big{[}-\big{(}\frac{\xi_{s}}{2}v_{0}^{2}+\beta_{s}\varepsilon(v)\big{)}\Big{]}&\textrm{if}\left\{\begin{array}[]{lr}v>v_{0}-\sqrt{\frac{2\varepsilon_{\phi}}{m_{s}}}\\ v<v_{0}+\sqrt{\frac{2\varepsilon_{\phi}}{m_{s}}}\end{array}\right.\end{array}\right.$$

$$n_s(\varphi )=n_{0s}\int f_s(\varphi )\text{d}v.$$

$$S(\varphi )={\int }_0^{\varphi }\sum _s{n}_s(\varphi )\text{d}\varphi -1+\text{const}s=i,e$$

$$\frac{1}{2}\Big{(}\frac{d\phi}{dx}\Big{)}^{2}+S(\phi)=0.$$

$$\begin{gathered} \frac{\partial f_s(x,v,t)}{\partial t}+v\frac{\partial f_s(x,v,t)}{\partial x} \\ +\frac{q_s}{m_s}E(x,t)\frac{\partial f_s(x,v,t)}{\partial v}=0,s=i,e. \end{gathered}$$

$$\frac{{\partial }^2\varphi (x,t)}{\partial x^2}=\sum q_s{n}_s(x,t),s=i,e.$$

$$n_s(x,t)=n_{0s}{N}_s(x,t),$$

$$N_s(x,t)=\int f_s(x,v,t)dv.$$

$$\sum _s{n}_{0s}=\sum _s{N}_{0s}{q}_s=0.$$

$$\frac{\text{d}{ϵ}_{\varphi }}{\text{d}t}=E\sum _s{q}_s\int v{f}_s\text{d}v,s=i,e$$

$$(\overline{n_e{v}_e}=)\int v{f}_e\text{d}v=\int v{f}_i\text{d}v(=\overline{n_i{v}_i})$$

Full text

I Introduction

Electron holes have been proven to be ubiquitous in the Earth’s magnetosphere via multiple satellite observations Bale et al. (1998); Franz, Kintner, and Pickett (1998); Matsumoto et al. (1994); Ergun et al. (1998); Deng et al. (2006); Kojima et al. (1997). These space-borne experiments display the existence of electrostatic solitary waves (ESWs) which implies holes in phase space. In 1956, Bernstein, Greene, and KruskalBernstein, Greene, and Kruskal (1957) derived the solution of the time-independent Vlasov-Poisson equations and proved the existence of general solutions for nonlinear plasma waves in the electrostatic regime, which is known as BGK modes. BGK modes include localized (solitary), shock and double layers solutions. Their study established that trapped population diversifies nonlinear solutions, way beyond the limits of fluid theory predictions. A broad class of nonlinear solutions, BGK modes, can be constructed based on the BGK approach. This method, basically, divides the distribution function into two parts, namely free and trapped populations. The existence of trapped population makes the solution fundamentally nonlinear, with no linear counterparts. Hence, even in small-amplitude limit BGK solutions are not reduced to linear solutions Bohm and Gross (1949); Schamel (2000, 2012). In the years after BGK’s seminal work, it was suggested that two more conditions had to be added to the Vlasov-Poisson equations in order to have physically meaningful solutions, i.e. net-neutrality conditionSmith (1970b) and positiveness of distribution functionSchamel (1971). These two conditions provide critical limitations which make it easier to pinpoint the physically possible solutions, among the infinite number of BGK modes. Hence, the set of equations to be solved simultaneously can be listed as:

[TABLE]

In general (when the fully time-dependent system is considered), the Vlasov equation for each species imply charge conservation, which for our case of one-dimensional spatial dependence reads $\partial\rho_{c}/\partial t+$ $\partial j/\partial x=0$, where $\rho_{c}=\sum q_{s}\int f_{s}dv$ and $j=\sum q_{s}\overline{n_{s}v_{s}}$ which basically a sum over flux of particles ($\overline{n_{s}v_{s}}=\int vf_{s}dv$). For time-independent solutions this implies $j=\mathrm{const}$. However, considering cases where the constant current deviates largely from zero is unphysical, as this would imply a substantial magnetic field generation, invalidating an analysis based purely on the Vlasov-Poisson system. The most strict choice is to limit ourselves to $j=0$, which for a hydrogen plasma ($q_{e}=-q_{e}$) reduces to

[TABLE]

We do not need net neutrality to hold exactly, as a small deviation only implies a weak magnetic field that can be neglected, but Eq. (3) must apply at last approximately.

Finally, we cannot allow for the distribution function to be negative , i.e. we must have

[TABLE]

Thus, the full set includes the time-independent Vlasov equations for each species (1), Poisson equation (2), the net-neutrality condition (3) and the positiveness of distribution function (4). Since we are focusing on the solitary (localized) solutions here, it is assumed that the electric potential approaches zero far enough from the peak of the solution, $(\phi\longrightarrow 0$ when $x\longrightarrow\infty)$.

Solving the above equations can be translated into finding two functions, according to:

Distribution functions of each species $(f_{i},f_{e})$: Each should be a function of the total energy (kinetic and potential energy) and has to be always positive. This applies to both free and trapped populations. It also should be continuous in the whole phase space domain. 2. 2.

Sagdeev pseudo-potential function $(S)$: It originates from the Poisson’s equation.

The Sagdeev pseudo-potential approachDavis, Lust, and Schlüter (1958); Sagdeev (1966) (also known as the phase plane approachInfeld and Rowlands (2000)) transforms the Poisson’s equation into an equation of the same form as for a particle moving in a potential.

Additionally, the flux of each species $(\overline{n_{i}v_{i}},\overline{n_{e}v_{e}})$ should almost cancel out, in accordance with the net-neutrality condition (3).

The BGK method approaches the problem by guessing a solution for $\phi(x)$ and then finding the desired trapped distribution function. The BGK approach has been used to find shock solutions Montgomery and Joyce (1969); Smith (1970a) as well as solitary waves Turikov (1984); Muschietti et al. (1999). The BGK approach predicts the width and potential range that supports the nonlinear solitary solution.

In contrast, Schamel Schamel (1972b) suggested a known trapped distribution function as the starting point, which results in the electric potential profile ($\phi(x)$). Specifically, SchamelSchamel (1971, 1972a) proposed a trapped distribution function which is positive everywhere, and continuous in both velocity and spatial directions:

[TABLE]

Here $A=\sqrt{\frac{\xi_{s}}{2\pi}}n_{0s}$ and $\xi_{s}=\frac{m_{s}}{T_{s}}$ are the amplitude and normalization factor respectively. Moreover, $\varepsilon(v)=\frac{\xi_{s}}{2}(v-v_{0})^{2}+\phi\frac{q_{s}}{T_{s}}$ represents the (normalized) energy of particles. The shape of the trapped population is determined by the trapping parameter ($\beta$). It can take three forms, e.g. hollow ($\beta<0$), flat ($\beta=0$) or hump ($\beta>0$). Other types of distribution function has been proposed to model trapped population. For example, SmithSmith (1970b) suggested monotonic transitions and distribution functions. This selection leads to monotonic particle densities. More recent studies for space plasmas have treated the trapped distribution function by considering a kappa distribution for the ambient plasma.Aravindakshan, Kakad, and Kakad (2018b, a).

Kinetic simulation studies of BGK modes (which may be generated for example by the Schamel approach or the BGK approach), can be divided into two main genres, i.e. production studies and propagation studies. There is a large body of numerical studies about different methods of exciting solitary BGK modes (referred to here as production studies). Interested readers might follow this topic in review papers such as Ref.[8] and [16] and the references within.

Both for production studies and stability studies, ions are typically assumed to be immobile (with infinite mass). Immobile ions result in merging (coalescence) behavior of electron holes Rasmussen (1982); Fijalkow and Nocera (2003); Ghizzo et al. (1988); Lynov et al. (1980); Shoucri (2017). When the ion dynamics is included in the simulation box, a splitting process of the merged hole is reported, which suggests that the BGK modes can be considered as solitons of the Vlasov-Poisson set of equations Hosseini Jenab and Spanier (2017); Zhou and Hutchinson (2018); Kakad, Kakad, and Omura (2017).

When it comes to the propagation studies of BGK modes, two questions is in focus, the nonlinear dispersion relation (NDR) and the stability of the solutions. The NDR refers to a relationship between the key aspects of the BGK modes, i.e. width, amplitude and velocity. Turikov Turikov (1984) used the BGK approach to produce electron holes (nonlinear solitary waves) in the simulation box, by assuming two types of bell shaped functions for the electrical potential. They found that the width of the solitary wave increases with increasing amplitude, unlike the KdV soliton. Concerning the stability of solitary waves, they reported that for Mach numbers bigger than two $(M>2)$ the electron holes decay during several plasma periods, and this time decreases with increasing Mach number. Other simulation studies were mostly concerned with NDR rather than stability Muschietti et al. (1999); Lynov et al. (1985).

In this study, we adopt the Schamel approach in order to produce BGK modes. Hence, the shape of the electric potential is determined self-consistently through the Sagdeev pseudo-potential equation. Amplitude, width depends on the trapping parameter and velocity of nonlinear solution. In our numerical approach both electrons and ions are described by the Vlasov equation.

The main point of inquiry for this study can be formulated as follows:

Are nonlinear solutions (produced by the Schamel approach) stable against collisions?

We answer this question by temporal simulations of nonlinear solutions with different velocities. In other words, this study explores the effect of velocity (as a major factor in the NDR) on the collisional stability of nonlinear solutions.

Our simulation results indicate that the nonlinear solution (provided by the Schamel distribution function) can be divided into three categories, e.g. stable, semi-stable and unstable solutions. Stable solutions are clustered slightly above the ion-thermal velocity (Mach number$>1$). As the velocity (Mach number) increases, the solutions transform from stable to unstable in a smooth transition. Hence, there exist solutions which are semi-stable. Some fluid simulations have suggested the instability of Sagdeev solutions. However our study stands as a the first attempt to explore this in the fully-kinetic regimeLotekar, Kakad, and Kakad (2017).

Our study reveals that net-neutrality violation is a cause of the instability of the nonlinear solutions in our study. In all the cases of simulations, a bipolar structure in flux of electrons is observed to form early in the simulation, due to an initial violation of equation (3). A larger deviation from (3) results in a stronger bipolar structure and a more rapid instability. The study of nonlinear solutions to the Vlasov-Poisson system shows that during collisions, the electron holes merge and exchange trapped populations. However, the merged holes break up due to the presence of the ion dynamic. Our study shows that the details of the merging and breaking up process varies for different electron hole set-ups. We focused on the number of rotations of the core of the electron holes during collisions, and the exchange of the trapped population, in order to parameterize different behaviors. It was observed that the number of rotations grows with the relative velocity of the electron holes. In the case when each of the electron hole core makes half turns during a collision (e.g. if the number of turns are $2.5$ or $3.5$), the electron holes exchange their cores with each other.

II Basic Equations and the Numerical Schemes

The normalization of the equations and variables are carried out with respect to the ion parameters. The spatial direction (length) and time are normalized by the ion Debye length ($\lambda_{Di}=\sqrt{\frac{\varepsilon_{0}K_{B}T_{i}}{n_{i0}e^{2}}}$) and the ion plasma frequency (\omega_{pi}={\big{(}\frac{n_{i0}e^{2}}{m_{i}\varepsilon_{0}}\big{)}^{\frac{1}{2}}}), respectively. The velocity, energy, and potential are rescaled by the ion thermal velocity ($v_{th_{i}}=\sqrt{\frac{K_{B}T_{i}}{m_{i}}}$), ion thermal energy ($K_{B}T_{i}$) and ($\frac{K_{B}T_{i}}{e}$ ), respectively. The mass and charge of the species are normalized by the ion mass ($m_{i}$) and the elementary charge ($e$).

II.1 Initialize the nonlinear solutions

As mentioned in the introduction, we are following the approach proposed by Schamel Schamel (1972b). Consequently, we start with choosing the distribution function of the trapped population, i.e. setting the value of the trapping parameter ($\beta$). Note that the Schamel distribution function guarantees positiveness and continuity of the distribution function, unconditionally for any value of the trapping parameter ($\beta$) in the whole phase space.

In the next step, one can determine the dependence of the distribution function on the electric potential ($f(\phi)$), as implied by the Schamel distribution function. Next, the number density dependency on the electric potential $n(\phi)$ can be determined by integration over the velocity:

[TABLE]

Next, the Sagdeev pseudo-Potential function $S(\phi)$ should be established by integration over $\phi$:

[TABLE]

Finally, the potential profile in the spatial direction $\phi(x)$ can constructed from the Sagdeev pseudo-potential $S(\phi)$ by integrating:

[TABLE]

The constructed electric potential and the Schamel distribution function associated with it, is used then to initialize the Vlasov-Poisson simulations, i.e. following the temporal evolution of the distribution functions. The details of this method can be found in Ref.[15].

II.2 The fully kinetic Vlasov-Poisson method

In order to perform fully kinetic simulations, both species are studied numerically using the Vlasov equations:

[TABLE]

The Poisson’s equation provides the force (electric) field:

[TABLE]

where $s=i,e$ represents the corresponding species. The Vlasov and Poisson equations are coupled by density integrations for each species to form a closed set of equations:

[TABLE]

[TABLE]

In which $N_{s}$ stands for the number density. $n_{0s}(=N_{0s}q_{s})$ is the (normalized) unperturbed value of the charge density. The background charge density is taken to be zero, i.e.

[TABLE]

Here, the Vlasov-Hybrid Simulation (VHS) method is adopted, which is closely related to the Particle-In-Cell (PIC) simulation method. In VHS, the distribution function in phase space of each species is sampled by enough phase points, i.e. all the details of the trapped population are covered Nunn (1993); Kazeminezhad, Kuhn, and Tavakoli (2003); Abbasi, Jenab, and Pajouh (2011); Kilian, Schreiner, and Spanier (2018). At each time step, the arrangement of phase points in phase space produces the associated distribution function for that time step. By knowing the distribution function for each step, all the useful quantities at the fluid level, e.g. the momentum density of the distribution function, can be easily deduced. Quantities such as kinetic energy density and entropy density are used for monitoring the conserved quantities, i.e. the total energy and entropy. Finally, the number density (and thereby charge density) is used to calculate the electric field through Poisson’s equation, which in turn is fed back to the Vlasov equation to push the phase space for the next time step. Note that the initial value of the distribution function associated with each of the phase points stays intact during the simulation. This property of the simulation method guarantees the positiveness of the distribution function under any circumstances during the temporal evolution Kazeminezhad, Kuhn, and Tavakoli (2003).

II.3 Variables and parameters

A number of parameters remain fixed through all of our simulations. Specifically this include: the mass ratio $\frac{m_{i}}{m_{e}}=1836$, $\Delta x=1.0$ (the grid size on the spatial direction) and $\Delta v_{x}=0.0025$ ($\Delta v_{x}=0.0025/v_{th_{e}}$) which is the grid size in the velocity direction for ions (electrons). Periodic boundary conditions are adopted in the spatial direction. Furthermore, the input parameters for each set of the simulations consist of two variables, the trapping parameter ($\beta$) and the velocity of the nonlinear solution ($v_{p}$). Note that by our choice of normalization, the ion-acoustic velocity, the electron plasma frequency and the electron thermal velocity are $v_{C}=\sqrt{1+\frac{T_{e}}{T_{i}}}$, $\omega_{pe}=\sqrt{\frac{m_{i}}{m_{e}}}$ and $v_{th_{e}}=\sqrt{\frac{m_{i}}{m_{e}}\frac{T_{e}}{T_{i}}}$; respectively. In what follows, we will express the velocity of solitary waves using the Mach number: $M=\frac{v_{p}}{v_{C}}$, in order to follow the notation of many theoretical papers.

The length of the simulation box is chosen as $L=1024$. The number of grid points in the simulation grid is $(S_{x},S_{v})=(1024,4800)$ for each species. Initially there are 16 phase points per cell in phase space, which are chosen randomly Abbasi, Jenab, and Pajouh (2011). Hence, the total number of phase points for each species is $\Sigma=78\times 10^{6}$. The velocity cut-off is (-6$V_{ths}$ ; 6$V_{ths}$): for electrons ($-6V_{th_{e}},6V_{th_{e}}$) with $V_{th_{e}}=\sqrt{T_{e}/M_{e}}$ and ions ($-6,6$).

The Vlasov equation, as a collisionless Boltzmann equation, should conserve entropy and other forms of Casimir invariants Elskens, Escande, and Doveil (2014) as well as the total energy. These conservation quantities are used as a yardstick for the simulation results. All simulation results reported here have a deviation strictly less than one percent in total energy and entropy.

III Results and Discussion

III.1 Three categories of the solutions:

In order to study the stability of the nonlinear solutions, we have made a simulation for each combination of $\beta$ (trapping parameter) and $M$ (Mach number). In each of these simulations, two identical nonlinear solutions (with opposite velocity, i.e. $\pm M$) are positioned on each side of the simulation box in the spatial direction. The temporal evolution of these simulations result in symmetric mutual head-on collisions. Since periodic boundary conditions are adopted, collisions happen repeatedly as long as the simulation is executed.

Fig. 1 presents the results regarding the electron-ion energy exchange for $(M,\beta)=$ $(1.1,0)$, $(1.25,-1.25)$ $(1.5,-2.5)$ and $(1.75,-2.75)$. The ion-acoustic speed for these simulations is $v_{C}=8$. The exchange of kinetic energy between electron and ions provides a reasonable estimator for the stability of the solutions. For stable solutions, the exchange of energy is perfectly symmetrical in time. Electrons gain kinetic energy from ions during the first half of the collision, and opposite flow of energy happens in the second half. After the collision, if the solutions are stable, both kinetic energies return to the same level as before the collision.

In case of $(M,\beta)=$ $(1.5,-2.5)$, the deviation from the initial value starts to diverge after each collision. This accumulates, and after a few number of collisions the solution can not recover its original level of energy balance even approximately, and the initial structure is lost. Such solutions are considered as semi-stable solutions. As for stable solutions, they can recover their original balance of energy after numerous collisions, for example the case of $(M,\beta)=$ $(1.1,0)$ and $(1.25,-1.25)$. The case of $(M,\beta)=$ $(2,-3)$ represents an unstable solution. After just one collision, the energy balance is destroyed and two electron holes merge with each other and they lose their initial characteristics.

Therefore, we have distinguished three categories of solutions which can be considered as a continuum, starting from stable solutions up to completely unstable solutions. On one side of this continuum, collisionally stable solutions are situated slightly above the ion-acoustic speed and survive multiple mutual head-on collisions. As Mach number increases, solutions start to get more unstable, i.e. they can survive a smaller number of collisions before destruction. On the other side of the continuum, for large Mach numbers, solutions are so unstable that they do not survive even one collision.

The conservation of energy and entropy of all these simulations were closely monitored and the deviation from exact conservation stood below $0.1\%$. The deviation of a variable ($X$) is defined here as the deviation from its initial value: $\Delta X=X(t)-X(0)$. In case of unstable nonlinear solutions, each collision creates some ion-acoustic and Langmuir noise in the simulation box. The Langmuir and the ion-acoustic modes are identified by using the Fourier transform and by checking their respective dispersion relation.Hosseini Jenab and Spanier (2016). When high frequency noise are produced, with a wavelength smaller than the grid spacing in the spatial direction, the small scale structures can no longer be tracked in the simulation and this causes a small deviation in the energy conservation. Electron (due to their light weight) reaches these level faster than the ions, and hence, one can observe an artificial increase in their energy level for unstable solutions. The more unstable the solution, the more prominent this effect appears in the results. However, the kinetic energy of electrons is in order of $10^{5}$ and even in the most unstable results reported here the deviation stands around $10^{2}$ which is around $0.1\%$.

Fig.2 explores the temporal evolution of the kinetic energy density (KED) of ions and electrons during head-on mutual collision. The figure shows results for the second collision between two nonlinear solutions with $M=$±$1.1$. As the two nonlinear solutions collide, the electron KED of the two solutions add up by following the superposition principle. However, during the collision, the ion KED stays the same as before or after the collision. This signals the duality that exists in the nature of these nonlinear solutions. They can be considered as “Coupled electron Holes with ion-acoustic Solitons (CHS)”Saeki and Genma (1998); Zhou and Hutchinson (2018). This theoretical framework is established when considering the electron- and ion-dynamic separately. On one hand, in the ion-acoustic regime and in the absence of ions’ motions, electron holes are considered as vortexes in phase space with a nonlinear structure prone to merging. When they merge, the electron KED adds up with little or no signs of a nonlinear interaction. On the other hands, considering ions without independent electron-dynamics (electrons act as Boltzmannian fluid) one can produce solitons, described by the KdV equation (for a weakly nonlinear amplitude). During the collision of two solitons, due to the nonlinearity, the superposition principle does not apply, and hence, it is not a surprise that the ion KED of two solitons do not add up.

Detailed analysis, however, revealed that during a short time in the middle of the merging stage and before the start of splitting, the energy of the electron holes goes a bit above the sum of the kinetic energy of the two electron holes. This is what can be observed in the Fig. 1 as an symmetric exchange of energy between electrons and ions. As the two ion-acoustic solitons collide, the extra energy of ions is converted to electrons.

III.2 net-neutrality instability:

Inspired by the symmetric energy exchange between electron and ions during collisions (in case of stable solutions), we have studied the balance of energy between electrons and ions, and the relation to the stability properties. For stable solutions, the potential energy (directly related to the electric potential profile) should stay stable. The changes in the potential energy is directly related to the flux of species:

[TABLE]

In order to have a stable solution, there should be no flow of potential energy ($\frac{\,\textrm{d}{\epsilon}_{\phi}}{\,\textrm{d}{t}}=0$). This condition (in its strongest form) results in

[TABLE]

which is called “net neutrality” condition in some of the theoretical studies of Vlasov-Poisson equation Schamel (1972b); Smith (1970b); Montgomery and Joyce (1969). However, we found numerically that if one wants to adhere to the strong form of this condition, the Sagdeev pseudo-potential will not converge. For the case of $\overline{n_{e}v_{e}}=\overline{n_{i}v_{i}}$, $S(\phi)$ appears as a straight line with a small slope. Hence, for the Sagdeev pseudo-potential to converge, a small deviation in the net neutrality should exist. The amplitude of this deviation depends on the details of the shape of the electric field and the potential. However, the strongest form of deviation from the net-neutrality condition happens when the electron average velocity has a bipolar structure, with opposite signs on the frontal side and the backside of the pulse.

Fig. 3 represents the profile of the flux of the electrons and ions early in the simulations, presented in Fig. 1. Here the nonlinear solution has a positive velocity and hence moves to the right. One can observe that the flux of electrons forms a bipolar structure early in the simulation ($\tau\simeq 5$). As the Mach number increases from $M=1.1$ to $M=2$, so is the amplitude of the bipolar structure and the instability of the solutions (compare this with Fig. 1). Note that the trapping parameter is selected in order to have the same profile for the flux of ions.

III.3 Internal dynamics of electron holes during collisions:

Fig. 4 shows the temporal evolution of the electron holes during collisions for different cases of relative speed. During the merging phase of a collision, two electron holes rotate around the collective center of mass. Our results show that the number of rotations depends on the relative speed of the electron holes. The expectation would be that as the relative speed increases, the time span of collision decreases and hence a fewer number of rotations can take place. However, it can be observed that as the relative speed increases, the number of rotations grow. In other words, electron holes spin faster when their relative speed is higher.

As the Mach number increases, the instability grows and hence the electron hole after collisions show a larger deviation from the stable shape. In cases of $(M=1.5,\beta=-2.5)$ and $(M=1.75,\beta=-2.75)$, it can be seen that the electron holes lose a rather large chunk of their trapped population after the collision. The instability originates from the bipolar structure in electron flux, that pushes the trapped population in opposite direction on the two ends of the electron hole. Hence, the electron hole loses parts of the trapped population on both ends of its domain. At the front of the electron hole (right/left side if the electron hole propagate to right/left), the trapped population is pushed upwards and leave the electron hole as a stretched arms in front of it. While on the rear side, the trapped population leaks out in the form of smaller electron holes leaving the trapped area.

The number of rotations during the merging phase plays a crucial role in the internal structures of the electron holes after the collision as well as during the collision process itself. Fig. 5 shows that in case of a non-integer number of rotations, the two electron holes exchange the cores with each other. Fig. 6 presents the trace of cores for the case of $M=1.1,\beta=0$. During collision, each core makes $2.5$ rotations and then they change their cores. Phase shift during collision of solitons has been predicted from a few theoretical studies of two-soliton solutions in KdV regime. Phase shift refers to the discontinuity in the time-line of soliton’s propagation before and after collision. This discontinuity can be observed in Fig. 6, as the time line shifts up during collision. In other words, on a kinetic level, the phase shift happens due to the rotations during the collision and the time lag that it produces in the time-line of propagations.

IV Conclusions

In this work, nonlinear solutions for the Vlasov-Poisson set of equations were studied using a fully kinetic simulation approach. The solutions has been provided by modeling the trapped population based on the Schamel approach. Simulations are carried out in a plasma consisting of electrons and ions with Maxwellian distribution functions. The main point of inquiry focused on the collisional stability of these solutions and its dependency on the propagation velocity (Mach number) as a major parameter. The Mach number plays a crucial role in the nonlinear dispersion relation (NDR) of the nonlinear solutions.

The kinetic energy of the plasma constituents are presented for different cases of stability. The exchange of energy between electron and ions provides a reliable measure of the stability. It is found that the solutions can be divided into three categories. (I) Stable solutions which are basically solitons of the Vlasov-Poisson equations. These are the solutions which can survive multiple mutual head-on collisions (at least nine successive collisions are considered in our simulations). (II) Semi-stable solutions which can hold their shapes for a few collisions. (III) Unstable solutions which get destroyed during the first collision. It is demonstrated here that the solutions smoothly transform from stable to semi-stable and then to unstable forms by increasing the propagation velocity, i.e. increasing the Mach number from 1.

Details of the exchange of energy between electrons and ions present two different natures of dynamics. Electron dynamics follows a linear superposition of kinetic energy during the collisions. However, ions do not follow a linear superposition principle. The combined dynamics result in an exchange pattern of energy between electron and ions. This exchange pattern turns out to be a reliable criterion to determine the stability of the solutions.

Our investigation shows that there exist a violation of the net-neutrality for the nonlinear solutions studied here. Violation of the net-neutrality results in a break-up of the energy balance of the nonlinear solutions. This violation appears as a bipolar structure in the flux of electrons. By increasing the Mach number of the solutions, the strength of this bipolar structure grows in comparison to the flux of ions. Therefore, it causes a stronger break-up of the energy balance sooner in the nonlinear solutions. Hence, the stability of the solutions decreases.

Furthermore, we have presented details of the electron phase space dynamics during collisions. These details show that the process of collision is built up of two major steps, namely merging and splitting. During the merging process electron holes rotate around their collective center of mass. The number of rotations depends on their relative speed. As the relative speed increases, the time span of this phase of collision decreases. However, the number of rotations increases simultaneously. Hence, different internal dynamics are possible. As a special case, it is shown that the collision of electron holes can be approximately considered as an exchange of their cores with each other. Furthermore, simulation results demonstrate a phase shift in the hole’s propagation time-line before and after collision. The phase shift can be interpreted as a time lag during collision in which the rotations of holes around their collective center of mass takes place.

Unstable solutions can be considered as a warning sign of over-interpreting the meaning of Sagdeev solutions. Our simulation show that despite that the solutions originate from the Sagdeev approach, not all of them are solitons or solitary waves, i.e. stable against collisions. Hence, the stability of the Sagdeev solution has to be considered as a non-granted problem in itself when the Sagdeev pseudo-potential is employed. There are many studies in multi-species plasmas which basically are attempts to find Sagdeev solutions to the set of fluid or kinetic equations of the system. It seems that this approach has taken for granted that the Sagdeev solutions can be considered as stable solutions.

Our simulation results shows that balance of energy and momentum should be considered as an important factor in the stability studies. This study also hints that in constructing BGK solutions one should pay attention to the flux profile. In order to construct BGK modes which are stable against mutual collisions, initial conditions should be adjusted in order to avoid the formation of bipolar structures in the electron flux, or at least to decrease its amplitude.

Acknowledgements.

G. Brodin and S. M. Hosseini Jenab would like to acknowledgment financial support by the Swedish Research Council, grant number 2016-03806. This work is supported by the Knut and Alice Wallenberg Foundation through the PLIONA project.