Ultrafast electron holes in plasma phase space dynamics

Electron holes (EH) are localized modes in plasma kinetic theory which appear as vortices in phase space. Earlier research on EH is based on the Schamel distribution function (df). A novel df is proposed here, generalizing the original Schamel df in a recursive manner. Nonlinear solutions obtained by kinetic simulations are presented, with velocities twice the electron thermal speed. Using 1D-1V kinetic simulations, their propagation characteristics are traced and their stability is established by studying their long-time evolution and their behavior through mutual collisions.

www.nature.com/scientificreports/ structural characteristics (shape) and in the physical mechanism underlying their formation. (An interested reader is referred to the above references for details). The aim of this study is to characterize high-speed electron holes by establishing their occurrence in a kinetic framework, and by investigating their stability profile and probing their soliton-like features. For this purpose, a novel distribution function (df), the 'ELIN df ' , is introduced as a generalization of the Schamel df. The ELIN df adjusts the distribution function of the trapped population of electrons by relying on a dynamically varying parameter β so that its moments can fit a predetermined curve and all of the desired featured of the Schamel df are retained, such as consistency and smoothness in both spatial and velocity spaces inside the trapped region ( Fig. 1).
To show the stability of our nonlinear solutions, three series of simulations are reported. Firstly, by considering the long-time evolution of an initial condition we will confirm the stability of the solution's profile during propagation, thus establishing them as solitary waves. Then, two types of mutual collisions are reported, i.e. head-on collisions (with no overlapping in velocity space) and overtaking collisions (moving in parallel and with overlapping). The aforementioned phase shift through collisions has also been investigated, to corroborate the fact that electron holes behave as solitons.

Results
Long-term evolution. Figures 2 and 3 display the temporal evolution of EH1. The initial condition and the last step of temporal evolution can be compared and show that the overall shape of the electron hole (Fig. 2) and the corresponding potential or field profile (Fig. 3)stay unperturbed.
Head-on collision. Figure 4 depicts a head-on collision between EH1 and EH2. After the collision ( 0 < τ < 2 ), both solutions keep their shape and velocity compared to their initial state. Note that due to their large velocity, they are well-apart in the velocity direction, i.e. there is no overlapping, and hence their collision   Overtaking collision. Although the previous simulations demonstrate the stability of these EHs, the strongest test of the stability is their interaction via an overtaking collision when they overlap in the velocity direction.
In an overtaking simulation, we have used two EHs e.g. EH1 and EH3. Figure 5 presents the temporal evolution of electric field/potential around the collision time τ = 3.2 in a frame moving with M = 45 . Both EHs survive  www.nature.com/scientificreports/ the collision, and their respective velocity stays the same. Focusing on EH1, displacement can be witnessed after the collision. A phase shift can be measured by comparing EH profile with the red line, which is an extrapolation of an unperturbed path of this EH. This displacement is similar to the well-known effect of "phase shift" which observed to happen in mutual collisions of solitons [25][26][27][28][29] . We show in Fig. 6 the electron df during the overtaking collision, which demonstrates the considerable interaction between the EHs during the collision and their overlapping on velocity direction. Yet after the collision the EH1 is largely unperturbed, modulo the observed phase shift. Interestingly, data fitting has shown that the sech 2 curve form approximates the numerical data better than any other exponent, including the (expected, arguably) sech 4 form (see Eq. 39 in 14 ).

Discussion
In summary, we have provided a method to produce high-speed nonlinear solutions which move at a speed beyond the electron thermal speed. We showed that these electron holes are stable, retain their profile through collisions and remain so in the entire duration of the simulation. For mutual collisions with considerable overlap in the velocity direction, the EHs display a "phase shift" This phase shift represents a signature of soliton behavior and hence suggests that these EHs can be considered as solitons (at least approximately). This has been suggested for much lower-speed EHs before 5 but without the observed "phase shift" reported here. Simulation code. We have employed the Gkeyll simulation framework 33 to solve the Vlasov-Ampere system of equations [34][35][36] . Gkeyll discretizes the equations using the discontinuous Galerkin finite element method in space, with a strong stability-preserving Runge-Kutta method in time. We have adopted a piecewise cubic Serendipity Element space for the basis expansion 37 (further details can be found in Refs. 34,36 ). The Gkeyll method has been compared to the standard PIC method, where it was demonstrated that the effective phase space resolution of the method is very high, permitting detailed studies of df dynamics. Such high accuracy is of paramount importance for the resolution of EH dynamics in phase space 38 . Parameters. In our study, the temperature and mass ratio are

Iterative method to find stable solutions. Our method follows the BGK method and starts by adopting
an arbitrary function for the electrostatic potential φ(x) and by choosing the value of the electron hole speed (v EH ) . We then use the ELIN df to produce the electron distribution function. Given that the potential profile provides the charge density (ρ(x)) , and using the Schamel df for the ions to obtain n i (x) , we then use the total charge density (profile) n e (x) = ρ(x) − n i (x) as a "guiding equation" for the ELIN df and thus construct the electron hole. We have adopted, to start with, the simplest form of potential profile suggested for electron holes i.e. φ = A sech p (x/L) in which p = 2 and A and L are the EH amplitude and length, respectively. The amplitude and length (values) are chosen randomly; however the system will damp/break the forced profile if it is not close-enough to a self-consistent nonlinear solution. The resulting electron hole may have different size and velocity, but with an iterative process, one can find the combination of {A, L} for which the solution will be stable enough for a specific ( Elin DF method to construct electron holes. In order to explain our novel distribution function approach, firstly we need to represent the Schamel distribution function in energy-dependent format. Here we briefly discuss this, more details can be found in the reference 19 . Schamel approach devides the distribution function into two parts, namely free and trapped particles which are separated by a separatrix. Focusing on the free particles, the following steps are taken to determine their distribution function ( f f ), assuming a pulse moving with a velocity ( v EH ) in the laboratory frame: 1. the shifted kinetic energy is found in the co-moving frame: T v ′2 and v ′ = v − v EH is the velocity in the co-moving frame. 2. the shifted kinetic energy is calculated in the laboratory frame: Free particles fulfill the condition ε ′ K > ε φ . Note that, in order to calculate the df at point v, we use the df at the point v sh , which can be written as f = D g (ε K sh ) in energy format. Here, v sh presents the velocity of particles before their interaction with the potential profile. By D g we denote a general distribution function satisfying the Vlasov equation, i.e. in principle any function depending on the constant(s) of motion. Here, the energy is used to construct a valid function. Well-known examples of D g are the Maxwell-Boltzmann df, the κ df 41-44 and the Cairns 45 distribution function(s).
In other words we trace the characteristics of the particle back in phase space. Then we use the value of df at v sh as the value of df for v since the df stays constant on the characteristics of Vlasov equation 46 .
The distribution function of trapped particles ( f t ) which are subject to the trapping condition ( ε ′ K < ε φ ) can be achieved by following the steps below: 1. the shifted kinetic energy is found in the co-moving frame: ε ′ K sh = |ε ′ K − ε φ | , using a Maxwellian df on top of this kinetic energy with a coefficient β , will provide the shape of trapped distribution function: In order to have continuity between trapped and free df where they meet in the velocity direction, f shape is multiplied by f base = D g (ε S ) . Hence f t = f base × f shape .
Here, D g (ε S ) stands for the distribution function at the separatrix where ε ′ K = ε φ and works as a constant value which can increase or decrese the f t , in order to adjust it with the free distribution function. The second component, f shape is velocity-dependent and is controlled by β . It may appear in three qualitative shapes, i.e. flat, a bump or a hollow curve, if β = 0 , β > 0 or β < 0 , respectively (see Fig. 1).
Hence, the total form of the Schamel distribution function 12 can be written in terms of the energy as: f = af (ε K ) in which a is a normalization constant and www.nature.com/scientificreports/ One can understand the Schamel df as carving up a given general distribution function ( D g ) around a particular velocity (hole velocity) and inserting a Maxwellian df with arbitrary temperature inside the hole ( f t ).
In the above representation of the Schamel df we used the analytical form of φ(x) . However, one can equally use the discretizied form of φ(x) (by deviding it into a n intervals of �φ ). Schamel df can then be retrieved by n → ∞ ( �φ → 0 ). In terms of simulation approach, these two methods are equal since even when using the analytical approach, one had to use discretization for φ(x) and there is limit on how small �φ can get.
In other words, to generate distribution function ( f f and f t ) for each interval, we only need update the value of f base in our approach and repeat the process. This results in multiple carvings, each based on the previous distribution function and it recursively progresses.
We assume φ = A sech 2 (x/L) as crude approach (stablished by trial and error in the beginning), and then we dicretize the first half of φ(x) into n intervals in the following form ( �φ = A n ): The second half will be the same as the first half except for a simple inversion. Hence we just build the first half of df and the second half is just simple inverted copy of it.
In this approach β can be changed for each interval, and this add a new degree of freedom to the Schamel df. We call this ELIN (rEcursiveLy extendable distribution for a trapped populatIoN) distribution function. The distribution function for each interval can be presented by the following equation. In which the D g (in Schamel df) is replaced by the distribution function f i−1 of previous interval and each interval has its own β i : in which f 0 is the initial unperturbed df (here assuming Maxwellian df, i.e. f 0 = D m ). β i can change arbitrarily in order for moments of df to fit a "guiding equation" (here, the equation for the electron density). To obtain a smooth distribution function in the x direction, one can increase n until the numerically-desired level of smoothness is achieved. An example of the ELIN df profile is presented at Fig. 1 which shows 10 successive (carving) iterations with β approaching zero from below (negative side). Note that since β originates from a continuous guiding equation, hence their successive values follow a pattern and are not randomly chosen.
To conclude, we have introduced a new method for constructing electron holes within a kinetic framework, which relies on a successive multi-step extension of the Schamel df (here represented in energy-dependent form), i.e. the ELIN df method. The ELIN df adopts a continuously varying value for β , in contrast to the Schamel df where β is a constant. This extension provides an infinite number of parameters for the ELIN df, which enables it to construct an electron hole for any given bell-shaped potential profile. In our computational approach, the number of free parameters in the ELIN df is finite and equals the number of intervals (n). We have adopted an iterative method (inspired by Newton's iterative scheme), built on top of the ELIN df method, to find the stable solutions. Starting from an initial guess, in each iteration of this method firstly we use the ELIN df to build an electron hole and then utilize the Vlasov-Poisson simulation method to follow the temporal evolution of the electron hole for a short time. We use the potential profile at the end of each iteration as an input for the next round of iteration. After a few iterations, the initial and final potential profiles are close enough for this to be considered as a stable configuration, for closure. Then, one can move on to longer-time numerical experiments, to investigate the long-time evolution of these localized structures and their behavior through mutual collisions. As a representative set, three stable solutions (i.e. EH1, EH2 and EH3; see above) have been reported in detail. (