Evidence of a new class of cnoidal electron holes exhibiting intrinsic substructures, its impact on linear (and nonlinear) Vlasov theories and role in anomalous transport

Novel types of solitary electron holes or in more general terms a novel class of cnoidal electron holes propagating in collisionless, driven plasmas is found numerically and described analytically. These structural modes are distinguished by substructures both in the electron density in real space and in the trapped electron distribution function in velocity space and are due to a two-fold trapping scenario; a regular one, called β-trapping, and a weakly singular one, called γ-trapping. Although the latter already appears early in O(ϕ) in the ϕ -expansion of the density, where ϕ is the electrostatic potential, it is the regular β-trapping effect, appearing in O(ϕ3/2), that remains the main cause of these new class of self-consistent structures. Through a proper participation of both trapping mechanisms the shortcomings (inconsistencies) of undamped linear wave theories (Landau, van Kampen) are nonlinearly rectified becoming true solutions of the Vlasov-Poisson system at this kinetic level of plasma description only. They hence provide the right key to the strictly nonlinear, kinetic world of undamped, coherent, small amplitude structures in driven, collisionless plasmas. A cusp-like slope singularity of the otherwise continuous electron distribution fe along the contour of a separatrix, on the other hand, exposes thereby a limit for the reliability of a kinetic nonlinear Vlasov-Poisson description and a simulation of structure formation in intermittent turbulent plasmas. It points to a new, deeper settled physics caused by a local enhancement of collisionality near separatrices through the generation of correlations. It is suggested that multiple trapping has the same origin than chaos in the particle trajectories near resonance, namely the non-integrability of discrete particle dynamics.


Introduction
In 2019 the plasma community celebrated the 40th anniversary of the first experimental observation of solitary electron holes at Risø (Denmark) [1,2] and of its first concomitant Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. theoretical description [3]. Both works on coherent hole structures were stimulated by earlier computer simulations (see e.g. [4,5]) of the two-stream instability in which persistent holes have been found in long term runs. A development was subsequently set in motion that continues to this day in all three directions-experimental, numerical and theoretical. A huge number of papers revealing various aspects of hole formation especially with focuses on space physics have since then been published. Unfortunately, however, not all of them turn out reliable, as a not inconsiderable number of works do not stand up to a critical evaluation. As far as the theoretical description of undamped coherent structures is concerned, some of even the most respected authors prefer the traditional path and seek for solutions along the linear Vlasov equation, ignoring or being unaware of the invalidity of this approach. Other authors treat nonlinearity inadequately or deal with the shape of the structure only i.e. they don't consider its associated speed, the second ingredient of a complete theory. Often standing nonpropagating hole solutions are of interest only or no allowance is made of a drift between electrons and ions i.e. only a current-less plasma is considered.
In the present paper an attempt is made to find a way through this jungle of wave literature. Concentrating on situations in which undamped, propagating, coherent wave structures are excited and which, as a consequence, can be described theoretically by exact solutions of the governing equations (without any approximation such as linearization), we shall supply proof of the following statements: 1) coherent structures in a microscopic description require the full account of the nonlinear Vlasov equation, no matter how small, 2) non-uniqueness is a distinctive mark of hole solutions as several, distinguishable trapping scenarios can participate, 3) cusp-type slope singularities in phase space along separatrices inevitably occur in a kinetic description; they reflect the non-integrability seen in a deeper settled, discrete wave-particle interaction process and designate the ultimate limitation of a kinetic, mean-field, single particle Vlasov approach, 4) Vlasov-Poisson codes give up their relevance even earlier as they cannot resolve singularities, 5) the pseudo-potential method in the version of Schamel [6,7], provides the best theoretical access to these structures; it allows the mildest form of singularities, relates phase velocities with the trapping conditions in an unique manner, and provides the closest contact to the linear, but insufficient wave theories. The original BGK method [8] is less suited as it principally cannot deliver a phase velocity, the second ingredient of a complete theory (see later).
But let us now proceed to the actual implementation of the above.
The main progress made in the publication [3] with respect to an earlier derivation [9] was the flexibility in the construction of the trapped electron distribution-called for brevity: the β-trapping effect lateron-accounting for the numerically observed prominent depression of the electron distribution function f e , deeply in the trapped electron region.
With this in mind, we present in this paper new material that shows that the β-trapping effect alone cannot render the complex nonlinear dynamics of driven, collisionless plasmas in sufficient detail and generality. To capture some more subtle aspects we have to introduce a second trapping effect, the so-called 'γ-trapping' mechanism that affects mainly shallowly trapped particles. Our innovation in describing these two trapping scenarios provide the key to solve the unsatisfactory explanation of undamped coherent linear wave solutions, and through this development lifting them to properly describe the non-linear situation.
In this paper we take a first step up into a hitherto unexplored world of structure formation in evolutionary driven plasmas synonymously covered under the heading of intermittent plasma turbulence. The focus is on undamped coherent structures which are excited in the linearly subcritical, ion acoustic regime of current-driven plasmas where linear Landau theory prohibits long living structures [10,11]. This unexpected finding results from the coherency aspect of structure formation that, due to an absence of phase mixing, undermines Landau damping [12,13]. A coherent wave propagating with modest phase velocities is inevitably affected by particle trapping, a physical reality that is absent if the Vlasov equation is linearized first. A stationary, coherent structure, in which from the outset the γtrapping effect is neglected, generally belongs to the so-called 'privileged spectrum of cnoidal electron and ion holes'. This special class is represented and defined by modes that are characterized by regular (holomorph) trapped particle distributions. These modes are obtained within the kinetically upgraded version of the pseudo-potential method derived by Schamel [6,7], a method which may preferentially be used to obviate the unwelcome single-pole singularities that spoil linear Vlasov theories [14,15].

Numerical simulation-Part 1
A typical evolution triggered by a tiny seed-like initial perturbation in the electron distribution function f e is shown by figure 3 of [16] which is reproduced here in figure 1 for convenience.
The evolutionary process is represented by six snapshots of the electron phase space, left column, and of the electron (red) and ion (blue) density perturbation, right column. The selected plasma parameters are: relatively cold ions, q = T T 10 e i ≔ , a subcritical electron drift, v D =0.01 being below the critical drift v D *=0.053, and an initial position of the localized perturbation at x=15, v=0.01. We see that after a temporal transition phase that involves creation of filaments and acceleration a rather robust solitary electron hole emerges asymptotically that moves with a phase speed of approximately v 0 =0.035. It is worth noticing that this phase velocity, which is normalized by electron thermal velocity, is supersonic with respect to the ion acoustic speed given for a hydrogen plasma by = = c m m 0.023 s e i . It is located microscopically at the high velocity tail of the electron distribution function, properties, especially the acceleration and jump of velocity, which have been explained in some detail in terms of the ionic shielding response function and of the negative energy property of the hole structure in [16].
Our attention is directed here on a different aspect, the central dip in n e which lasts for all times. It is seen in all of our simulations and, as we will show, cannot be explained by existing theories, inclusively the (privileged) cnoidal hole theory of Schamel and coworkers [3,[17][18][19][20][21][22][23][24][25].
In the following we present an extension of this theory that can be attributed to this phenomenon and explain its appearance in theoretical terms.

Derivation of nonlinear dispersion relation and pseudo-potential
With reference to [21] and [24] we start from exact solutions of the time-independent electron Vlasov equation, In this equation θ(x) represents the Heavyside step function, where v D is the electron drift velocity, v 0 is the phase velocity in electron space; k 0 is related with the wave length of the structure. Notice that f e (x, v) is a function of two constants of motion, ò and σ, which is a necessary requisite for a propagating wave solution. It consists of two parts, the contribution of untrapped particles, ò>0 [26], and the one of trapped particles, ò0. Trapping is controled firstly by the parameter β, which yields a notch in the electron distribution in the resonant trapped particle region when β<0. A further trapping term has secondly been added by the term g - . It independently contributes to trapping and introduces a weak singularity of the slope of f e (x, v) at the separatrix ò=0. It will be the source of our extended view.
We emphasize that a similar  term has already been inaugurated by the free particle distribution function. In other words, f et with the - term included belongs to the same mathematical class of admitted functions.
Note that f e (x, v) is continuous across the separatrix and it is assumed that 0f(x)ψ=1. Note also that the distribution ansatz (1) is already perturbative in character utilizing the smallness of ψ and ò, respectively, in f et . By this we get a two parametric solution with the best possible connection to common linear wave theories. Other choices of f et , involving e.g. higher powers of ò or an exponential of ò, are possible but make sense only when the amplitudes are finite, ψO(1). Such structures are often found detached from small amplitude wave spectra as they typically don't have a y  + 0 limit anymore. An example is the strong double layer solution of [27].
Other extensions may involve distributions that are discontinuous at the separatrix, as analysed by [24,25,28]. This corresponds to a line discontinuity in the terminology of [29] being defined as a discontinuity along a 1D manifold (separatrix) embedded in a 2D surface (phase space). These solutions not only yield a different potential (a solitary-like cos 4 (x) solution, see appendix of [24]) but require a tiny overpopulation of the trapped particle component, being, however, not suggested by our simulation.
Our continuous distributions (1) exhibit a different type of line singularity at separatrix, namely a cusp-type slope singularity of f e introduced by  | | that we may call a linecusp singularity.
A further class of different potential profiles is obtained by the BGK method [30][31][32][33][34] where f(x) is posed in real space rather than the general structure of f et in phase space, as done here. A member of this class, however, has a certain disadvantage as it cannot be associated with a definite phase velocity [24,25]. Moreover, it has a stronger slope singularity of f et at the separatrix ( ¶~- . We shall come back to the presence and consequences of singularities at the end of the paper. Velocity integration of (1) gives the electron density [6,7,20,24,25] f y g p f )˜represents the β -trapping effect and appears one order later in the f expansion of the density than the γtrapping effect which already becomes manifest in linear order.
In a similar way the ion density can be found as (see (20a) and (21) of [21]) where a corresponding ion γ-trapping term has already been ignored. Also B i has been neglected since we are mainly interested in the electron nonlinear effects only. Note that u 0 =μ v 0 is the phase velocity in the ion phase space and it where in the last step the pseudo-potential f ( ) has been introduced, we arrive at represents the nonlinear dispersion relation (NDR) and determines the phase velocity. The pseudo-potential (5) is presented in its canonical form that is obtained by a subtraction of y =  0 ( ) through which the γtrapping effect drops out. f ( ) stands for potential structure f(x) which is given via a quadrature of the pseudo-energy, . Both are the two necessary ingredients of a complete nonlinear wave theory.
(5) stands for a wide class of cnoidal hole modes representing a continuous wave spectrum, the most prominent member of which being the solitary electron hole mode [3,9]. It is given in the limit k 0 =0 by f Note that this continuous wave spectrum remains unaffected by the γtrapping effect.

Preliminary view of trapping in subcritical plasmas
We return to our plasma simulation and utilize its special

( )
A further simplification is achieved for a solitary wave given through k 0 =0 by which we get: The NDR (4) is in this limit given by These are the two equations that allow us to enter into a new world of plasma structures.
We first see that through the bracket term in (7)  This condition tells us that γ has a decisive influence on the existence of hole solutions exhibiting a central density dip.
As an example we take our observed values for θ=10 for a hydrogen plasma: v 0 =0.035, u 0 =4.74, , in which case the NDR (8) becomes B e =0.48+0.886γ. Putting this B e into our dip constraint (9) we get 0.13<γ. We thus learn that within the privileged hole spectrum, given by γ=0, a local density dip in x-space is not available.
If we however take γ=1, as a first demonstration for the effect of a nonvanishing γ, and ψ=0.04, we get B e = 1. 2 where y = 2 +0.2828, shows that shoulders appear near the separatrix for shallowly trapped particles.
Or said in different words, in v-space the distribution f et exposes maxima near the separatrix stemming from the γ-trapping effect. For more deeply trapped particles, on the other hand, the β-trapping effect remains the striking protagonist causing the 'deep' excavation. Although it appears one order later in the f expansion of n e and may hence be considered as a higher order, ignorable effect it is the main reason for the formation of these structures f(x). It hence indeed remains essential for the existence of this 'extended privileged class of cnoidal hole solutions', as we may call it, with its members showing a delicate substructure both in real space and phase space. We finally mention that shoulders obviously represent a general characteristic of hole solutions as they have been seen in finite amplitude BGK solutions, as well [32]. Their possible physical origin and consequences will be discussed in section 7 later.

Nonlinear revaluation of Landau and van Kampen solutions
As a by-product of not less fundamental importance, we address an old issue of Vlasov-Poisson waves, the van Kampen approach. Treating undamped harmonic waves by the truncated, linear Vlasov equation van Kampen arrives at his famous continuous wave spectrum given by ω=ω(k, λ) where λ describes the 'non-dispersion' part of the continuum being introduced by a delta-function contribution of f e at resonant velocity [21,35]. Unfortunately, he misses to solve also the full Vlasov equation due to the well-known inconsistency, namely the failure of ¶ < | at resonant velocity, a fact being still overlooked by a majority of plasma physicists.
Within our approach, however, we can upgrade this solution by finding an equilvalent nonlinear description of  this spectrum which is valid for the untruncated, full Vlasov-Poisson system without any reservation. If we set B e =0, v D =0 and k 0 =k into (4), (5) and use immobile ions (formally represented by θ=0, see (3)) we obtain which is exactly van Kampen's final set of equations. (We note in parenthesis that the correctness of this set obtained in a linear fashion does not imply that the linear route chosen for its derivation must be correct as well. Linearization is a prohibited procedure violating intrinsic consistency [15,36]. Only the nonlinear route can claim universal validity!) This set describes the continuous spectrum in which −λ represents the factor in front of the δ-function ( [21,35]) and k is the wave number of the monochromatic wave. ). Note that high phase velocities demand a negative λ which is consistent with γ being positive. The constraint B e =0 is either satisfied in the fluid limit, which typically yields a negative β (see also [37]). We conclude that by both trapping mechanisms, the β-one being constraint by B e =0, the van Kampen wave solution can be uplifted on a relieble nonlinear level becoming only then a true solution of the VP system.
Similarly, we can upgrade Landau's zero-damped harmonic wave solution too, satisfying the dispersion part of the dispersion relation, the thumb curve, by setting B e =0, k 0 =k and γ=0 [15,36] . In this case it is the β-trapping effect alone that accounts for this upgraded, trustworthy solution of the full VP system.
Generally speaking, both trapping effects, the γas well as the β-effect, are needed to complete the undamped coherent wave theory being hence its elementary ingredients.
Hence, linear Vlasov theories, although macroscopically indistinguishable from their nonlinear counterpart, loose their validity if treated more deeply in the microscopic domain. This holds at least for stationary waves but of course has to be extended to transient processes, as well, whenever coherent structures interact resonantly with particles.

Stability and singularity of hole structures
Furthermore, let us regard the stability of a nonlinearly upgraded van Kampen mode, the corresponding Landau one having already been considered in [38]. It is easily seen by a slight modification of the Landau-result [38] that such modes in current-carrying plasmas are unconditionally marginal stable independent of the chosen v D , whether it is above or below the critical one v D *, again in conflict with standard linear Landau theory.
With this we have found an intrinsically consistent, selfconsistent description of our numerical findings and have thereby lifted up linear Vlasov theories on their appropriate nonlinear level.
A remark on the singularity of the involved distribution function f e in (1) is called for. A N-body wave particle description of this phenomenon (see e.g. [39][40][41][42][43][44] ) reveals a chaotic region around the separatrices as the particle trajectories behave erratic, ergodic rather than controlled, deterministic there. This effect is discrete in nature and can of course not be resolved by a kinetic, mean-field approach with its smoothed distribution functions. An overstrained formalism therefore responds this way by forming line singularities in this restricted area of phase space. The mathematical singularity hence has its physical origin in the non-integrability of particle motion in the region located near the separatrix.

Numerical simulation-Part 2
Having presented these findings we can now come to a closer inspection and interpretation of our results. The progress in the numerical evolution in a refined examination is seen by the snapshot of f e at the hole center for t=155 in figure 4 and for t=235 in figure 5.
We see that stationarity is not fully achieved as the trapped distribution function still shows non-symmetry, the former more than the latter. What is already essentially established are the phase velocity v 0 =0.035 and the wave amplitude ψ=5.2 * 10 −5 . Especially the region near the separatrices is still exposed to variations, showing that separatrix crossing processes (of the characteristics) still need time to be finalized. Also the depth of the dip is a little bit less pronounced in the later stage, an indication that the filling of the deeply trapped region is still going on. An indication of shoulders is seen but not in a systematic, symmetric manner. This shows that the code has to struggle in the vicinity of the separatrices to meet the requirements of the contiuous Vlasov-Poisson (VP) system.
We don't see a full agreement between the numerical and analytical results, in the shoulders as well as in the depth of the trapped region. And, as we would like to see, this should hold true especially if we give the numerics more time to approach the asymptotic regime of stationarity.
The question now is: why does the code not achieve an identical result?
The reason lies, of course, in the discretization that prevents the code from resolving the slope singularity of the VP solution. The code is overtaxed to approach the continuous exact result and is hence not in the position to resolve the details of the analytic solution.
The deviation between numerics and analysis, as seen in the figures, has hence its origin in the overstraining of the code. The numerics rather than the hole solution is responsible for their differences in the resonant particle region.
The true physics, however, isn't represented by the analysis either as it does not include the chaotic dynamics of trajectories (characteristics) in the resonant region, as pointed out already.
Our claim however is that singularity and non-integrability of the resonant wave-particle interaction process are twinned and that the singularity in the continuous, mean-field, single particle VP system signalizes the presence of a new physics that can however only be met on a deeper level of plasma description, such as e.g. through the participation of the two-body correlation function in a BBGKY hierarchy [11].
The singularity of the slope of f e at the separatrix gives a hint that the two-body correlation function g e , which is usually of O(1/n λ D 3 ) and hence very small, where λ D is the Debye length, is dramatically increased in this region as g e , which appears as a driving term in the single particle equation of the hierarchy, becomes in lowest order proportional to ∂ v f e ∂ v f e and is hence at least in this region no longer negligible. The same quality probably holds also true for the ternary correlation function as it is driven by a term ∂ v f e ∂ v f e ∂ v f e and so forth for the higher correlations as well. Forthcoming investigations may clarify this issue.
Numerically, one can argue that the discretization procedure, by its smoothing property, anticipates this impact from the deeper lying 'chaotic/collisional' physics but one should be careful and not intermix both effects since they have essentially nothing to do with each other.
Even a PIC code cannot offer a way out as it still employs a discrete grid and would probably need a tremendously big number of particles in order that a macroparticle, representing a cloud of individual particles, performs the same 'chaotic/collisional' dynamics than the particles itself.
The numerically best representation of kinetic VP plasmas is presumably a Fourier-Hermite decomposition of f e because it can shift its limits of reliability furthest out [20,22]. Simply by increasing the number of Fourier and Hermite function modes one can get the desired resolution in space and velocity. One should, however, keep in mind that in order to represent the nonlinear structure correctly one would indeed need an infinite number of them.
We are now retrieving the determination of ψ and B e , used in our data, by presenting the best fit between our ana-   electron hole is superimposed by an electron plasma oscillation (EPO) which oscillates in time with Langmuir frequency and is situated near v D =0.01. The EPO is also seen in figure 5 as a small perturbation in f e , noting that we have used in figures 5 and 7 an instant where the EPO is not completely absent, i.e. at zero level. The simultaneous excitation of several members of the zoo of the (extended) privileged cnoidal hole spectrum, initiated by our initial perturbation, can thus be exposed.

Summary and conclusions
In summary, the existence of a new class of electron holes has been shown numerically by our simulations of current-driven plasmas in the ion acoustic, linearly subcritical regime and proven analytically by an extension of the standard (privileged) cnoidal hole theory of Schamel. The main characteristics of a typical hole mode is the development of a substructure in the perturbed electron density in terms of an additional local depression at its center being superimposed on the clock-like profile. It exists at all available phase velocities especially at supersonic speeds where standard hole theory prohibits a central dip. Responsible for this new phenomenon with its measurable macroscopic signature is an additional weak slope singularity of the underlying trapped electron distribution function at the border line in electron phase space that separates trapped from untrapped particles. The latter, called γ-trapping effect, contributes constructively to the main trapping effect (β-trapping), that is mainly responsible for the potential structure, causing a shoulder in the trapped particle region near the separatrix. Both trapping mechanism turn out to be the key to lift undamped linear wave solutions (Landau, van Kampen) on their appropriate nonlinear level.
In contrast to linear wave theory, these upgraded structures in the harmonic wave limit turn out marginally stable in subcritical regions of current-carrying plasmas independent of drift and phase velocities.
Our analysis has opened a window to a new physics beyond the kinetic description in which a local enhancement of the otherwise negligible collision operator near separatrices is involved. This not only affects kinetic plasma theory by introducing collisionality/chaos near separatices of structures in intermittent turbulent plasmas but also numerical (Vlasov) codes that cease to be trustworthy even earlier. The latter can at best provide the strength and the phase velocity of a structure but cannot record the intrinsic details of the realistic wave-particle interaction physics in its resonant domain and hence in general terms the physics associated with anomalous resistivity and transport [20,22].