On the correct implementation of Fermi – Dirac statistics and electron trapping in nonlinear electrostatic plane wave propagation in collisionless plasmas

Strathprints is designed to allow users to access the research output of the University of Strathclyde. Unless otherwise explicitly stated on the manuscript, Copyright © and Moral Rights for the papers on this site are retained by the individual authors and/or other copyright owners. Please check the manuscript for details of any other licences that may have been applied. You may not engage in further distribution of the material for any profitmaking activities or any commercial gain. You may freely distribute both the url (https://strathprints.strath.ac.uk/) and the content of this paper for research or private study, educational, or not-for-profit purposes without prior permission or charge.

On the correct implementation of Fermi-Dirac statistics and electron trapping in nonlinear electrostatic plane wave propagation in collisionless plasmas

I. INTRODUCTION
Plasmas support a large number of nonlinear structures such as shocks, double layers, solitons, envelop solitons, and vortices.A particular class of nonlinear structures occurs due to particle trapping, where ions or electrons are trapped in the self-consistent potential wells in the plasma.Such structures arise naturally in streaming instabilities or other not necessarily large amplitude perturbations.3][14] At much smaller scales and in dense plasmas, the electron density is so high that quantum effects must be taken into account. 15,16The reason is that the inter-particle distance is small enough that the wave-functions of the electrons start to overlap, and Pauli's exclusion principle dictates that two electrons, which are fermions, cannot occupy the same quantum state.This changes the equilibrium distribution of the electrons from a Maxwell-Boltzmann distribution to a more general Fermi-Dirac (FD) distribution that takes into account electron degeneracy effects.In the limiting case of zero temperature, the equilibrium distribution for the Fermi gas is a flat-topped distribution in momentum space (or in velocity space for non-relativistic particles).At very small scales, there are also effects due to electron tunneling, and in this case, the dynamics of a single electron is described by the Schr€ odinger or Pauli equation, and statistical models are based on the Wigner equation.Some works have taken into account the first-order electron quantum tunneling effects in the description of electron holes, 17,18 but in this paper, we will neglect these effects.We deal with the general framework of structures within the Vlasov-Poisson (VP) description for which Bernstein, Greene, and Kruskal (BGK) 20 presented the first mathematically correct solution.Physically, however, another method, the pseudo-potential method 19 has proven preferential, as it is mathematically as general and correct but allows from the outset the selection of additional physically relevant solutions that are not possible to obtain within BGK theory. 20Hence, our model is semi-classical in that it takes into account the modification of the equilibrium electron distribution function to a Fermi-Dirac distribution, while we consider structures large enough that electron tunneling effects can be neglected.

II. CONSERVED QUANTITIES AND REDUCTION OF DIMENSIONALITY
For pedagogic reasons, we here briefly review some basic properties of the Vlasov-Poisson system which will be of importance in Sections III-VII.For brevity, we restrict the discussion to electrons, but the results are easily extended to include one or more ion species.The Vlasov equation describes the evolution of the electron distribution function f e in 3 spatial and 3 velocity dimensions, plus time, where e is the magnitude of the electron charge, and m e is the electron mass.The electrostatic potential / is obtained from Poisson's equation where n i is assumed to be an immobile and homogeneous neutralizing ion background density.Using that the Vlasov and Poisson equations are Galilei invariant, stationary waves in an inertial frame moving with the phase velocity of the wave can be obtained by setting the time-derivative to zero in Eq. ( 3), giving time-independent Vlasov equation Solutions of Eq. ( 3) are constant along particle trajectories in ðr; vÞ space, given by E ¼ m e v 2 =2 À e/ ¼ constant, where , and hence a solution can be expressed as f e ¼ f 0 ðEÞ, where f 0 ðEÞ is any function of one argument.(f 0 can also be a multi-valued, as discussed below.)The conserved energy E is essentially due to reversal and translational symmetries of the Vlasov equation in time.Depending on the shape of the potential /, there may exist additional symmetries leading to other conserved quantities.The most important case for our purposes is when / depends only one spatial variable, which we will take to be in the x-direction.In this case, there is a translational symmetry along the y and z directions, which leads to that the y and z components of the particle momentum (and velocity) are conserved (v y ¼ constant and v z ¼ constant along the particle's trajectory) and can be used to eliminate these quantities in E, so that a "one-dimensional" conserved energy can be found, E x ¼ m e v 2 x =2 À e/ ¼ constant.Hence, the solution of the time-independent Vlasov equation for this case can be written f e ¼ f 0 ðE x ; v y ; v z Þ, where f 0 is a function of three arguments.A last step to reduce the two velocity dimensions perpendicular to the x axis is to integrate the distribution function as This one-dimensional distribution function is a solution of the one-dimensional (1D) Vlasov equation obtained by integrating Eq. ( 3) over v y and v z , and where / is governed by the 1D Poisson equation If the potential instead has a rotational symmetry (spherical or cylindrical), an additional conserved quantity to E is the angular momentum with respect to the rotation axis, which can be used to construct classes of time-independent solutions of the Vlasov equation. 21Various extensions including a constant external magnetic field can be found in Refs.22-24 and relativistic effects in Refs. 25 and 26.

III. GENERAL FORMULATION OF ELECTRON TRAPPING
Consider a collisionless plasma, in which the unperturbed homogeneous electron distribution function is f 0 ðṽÞ, with an equilibrium number density n 0 and an electron temperature T e .The plasma frequency x pe ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n 0 e 2 =m e 0 p , the thermal velocity v the ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k B T 0 =m e p , and the Debye screening length k De ¼ v the =x pe will be used for normalization of the dependent variables, here temporarily being denoted by tilde, into their dimensionless form through t ¼ x pe t; x ¼ x=k De ; v ¼ ṽ=v the ; / ¼ e /=k B T 0 , and f 0 ðvÞ ¼ ðv the =n 0 Þ f 0 ðṽÞ, which we will use below.
Consider next the underlying 1D Vlasov-Poisson (VP) system ( 5) and ( 6), which is perturbed by a coherent plane wave structure /ðx; tÞ ¼ /ðx À v 0 tÞ that propagates with velocity v 0 in the lab frame.It is then constant in time in a frame moving with the wave.In this wave frame, the normalized Vlasov equation hence becomes ½v@ x þ /ðxÞ@ v f e ðx; vÞ ¼ 0; (7)   where we, for brevity, have dropped the superscript "1D" on f e and have written v x ¼ v.For our purposes, we write the solution of Eq. ( 7) as where ¼ v 2 =2 À /ðxÞ is the single particle energy and hðxÞ is the Heaviside step function.The separatrix, which divides the phase space into a trapped and untrapped particle region, is given by ¼ 0. The unperturbed distribution function is thereby extended, as a consequence of /ðxÞ, to the inhomogeneous distribution f 0 ðnÞ, which represents the untrapped (or free) particles.The new argument n is a generalized velocity 27 and is defined by The solution involves a two-parametric (sgnðvÞ, ) constant of motion of (7), such that (8) solves (7).Note that n ¼ v when / ¼ 0 and that sgnðvÞ is a constant of motion only for the untrapped particles.Hence, for > 0, f 0 becomes multivalued with respect to and depends also on sgnðvÞ.In Eq. ( 8), f t ðnÞ represents the trapped electron distribution, which is restricted to the trapped region, À ffiffiffiffiffiffi 2/ p v ffiffiffiffiffiffi 2/ p , corresponding to 0. Here, we already have adopted a nonnegative /.The electron density is thus obtained as

IV. HALF-POWER EXPANSION OF THE ELECTRON DENSITY
In the following, we assume a small but finite amplitude of the perturbation such that 0 /ðxÞ w ( 1, where w is the small amplitude of the potential.The electron density in the presence of a structure is found by a velocity integration of Eq. ( 8).Utilizing the smallness of /, we get by a Taylor expansion of Eq. ( 10) around / ¼ 0, 27,28 where k 0 will later be a measure of the wave's wavenumber, and the coefficients are and where P denotes the principal value.In ( 12)-( 14), we made use of v ¼ n when / ¼ 0, which is valid for untrapped electrons.
Deviations from continuity and analyticity have recently been treated in Ref. 29.The coefficient c in Eq. ( 14) is given for completeness only and will not be used below.Trapping is manifested by the b/ 3=2 term in Eq. (11).Equation ( 11) is hence a half-power expansion in /, in which the / 1=2 term vanishes due to continuity at ¼ 0. 29 This provides the formal background under which a wide range of collisionless plasmas can be treated formally on the same footing.

V. THE MAXWELLIAN CASE
We start with the general degenerated state and perform subsequently the classical Maxwellian limit, as follows.For fermions in a degenerate state, the unperturbed single particle distribution is given by the Fermi-Dirac (FD) distribution, which in the three-dimensional (3D) velocity space reads in dimensional units 30

FFD
where l is the chemical potential and Ẽ ¼ ðṽ is the single particle energy in the 3D velocity space (in lab frame).In the limit T 0 !0, one obtains l ¼ ẽF ¼ ðh 2 =8m e Þð3n 0 =pÞ 2=3 , which is the Fermi energy and h is Planck's constant.
Returning to dimensionless quantities, we can hence write where the normalization constant is the normalized single particle energy in the 3D velocity space (in the lab frame).Assuming a unit equilibrium density, the relation between l and A is found by a 3D velocity integration using spherical coordinates in velocity space where is the polylogarithm function with index n.By a 2D velocity integration over ðv y ; v z Þ-space, using cylindrical coordinates, we get the reduced, 1D distribution in the lab frame The Maxwellian case is formally (for more details, see later) obtained in the limit l !À1, which by the small value expansion Li ðxÞ % x gives A ! ð2pÞ 3=2 e l , which together with lnð1 þ e Àv 2 x =2þl Þ ! e Àv 2 x =2þl give the Maxwellian distribution The free particle distribution function in the wave frame, introducing v ¼ v x À v 0 , then reads For the trapped particles, we use Note that continuity at the separatrix where n ¼ 0 for finite / is automatically guaranteed by Eqs. ( 21) and (22), and that (22) uses the lowest order terms of the Taylor expansion of e Àb j /¼0 .This latter exponential expression was used earlier for finite amplitude waves, e.g., in Refs.1-3 and 19, but the results are the same as long as w ( 1. Insertion of expressions ( 21) and ( 22) into Eqs.( 12) and ( 13), respectively, then yields which are familiar expressions found in many earlier papers of Schamel and coworkers.We note again that the b-dependence of f t and of b M is a crucial ingredient of getting suited trapped particle equilibria.If, as done by Landau and Lifshitz, 31 b is assumed zero, corresponding to a completely flat trapped particle distribution, then no solitary electron holes would exist theoretically, 1 to give an example of its limited application (besides its restriction to v 0 ¼ 0).The hollow nature of the trapped particle vortices found in phase space, being a general feature of all numerically observed phase space structures, demands the incorporation of a b < 0. The question is now how in a quantum statistical setting the corresponding structural equilibria turn out to be.

VI. THE FERMI-DIRAC CASE
For fermions in a degenerate state, we got ( 16) as the dimensionless Fermi-Dirac distribution F FD ðEÞ in 3D velocity space, from which the reduced 1D distribution f FD 0 ðv x Þ in the lab frame (19) was obtained.
The reduced free particle distribution function for !0 in the wave frame, introducing v ¼ v x À v 0 , then reads For the trapped particle distribution with < 0, we choose the same function as in the Maxwellian case but with a different factor that accounts for continuity at ¼ 0 (and at n ¼ 0) Equations ( 27) and ( 28) are the two desired functions needed to derive a and b for the density expression (11) in the FD case.By insertion of ( 27) into (12), we first get where we defined the Fermi-Dirac plasma dispersion function which is by definition real valued and whose derivative is given by The generalized complex plasma dispersion function Z l ðuÞ with a complex argument u would then be a one in which the principal value is replaced by the Landau contour.We note that Z l r ðuÞ !Z r ðuÞ and Z l r 0ðuÞ !Z 0 r ðuÞ in the classical limit l ( À1.In the case of a standing wave pattern, v 0 ¼ 0; a FD can be obtained by a direct integration and can be expressed in terms of a polylogarithm function as The function Li 1=2 ðxÞ is real-valued and negative for x < 0. On the other hand, for high phase velocities v 0 ) 1, we have the asymptotic expansion Figure 1 shows line-plots of a FD ðv 0 Þ given by ( 29) for different values of v 0 and l.While a FD ðv 0 Þ is a decreasing function of l for v 0 ¼ 0, the dependence of a FD ðv 0 Þ on l is more complicated for larger values of v 0 with either local maxima or minima (cf.Fig. 1(b)).
To obtain b FD , the two functions f FD 0 ðvÞ and f FD t ðvÞ have to be differentiated twice and, according to (13), the result has to be evaluated at v ¼ 0 and n ¼ 0, respectively.After some algebra, we get where we defined r ¼ e lÀv 2 0 =2 .For a wave with zero velocity with respect to the lab frame, v 0 ¼ 0, we hence get with a relative simple expression, which stands for the electron trapping nonlinearity in case of a non-propagating wave pattern in a degenerate plasma.The dependence of b FD ðv 0 Þ on v 0 , l, and b is shown in Fig. 2.

VII. CNOIDAL ELECTRON HOLES AND THEIR SPEED
To close the system, i.e., to find self-consistent solutions for both systems, we have to solve the second part, the Poisson equation, which in the immobile ion limit becomes In the second step, we have introduced the pseudo-potential Vð/Þ.Its knowledge allows by a quadrature to find the final shape of the potential structure /ðxÞ via the pseudo-energy where we without loss of generality have assumed Vð0Þ ¼ 0. Insertion of the electron density (11) into Eq.( 36) and a subsequent /-integration yield The speed of the potential structure is obtained from the nonlinear dispersion relation (NDR), VðwÞ ¼ 0, which becomes Eliminating a in Eq. ( 38) by Eq. ( 39) we get It represents a typical cnoidal wave solution, a one being represented by Jacobian elliptic functions such as cnðxÞ or snðxÞ.The first term in (40) represents a purely harmonic wave, and the second one is a solitary wave, provided that b < 0. The general two-parametric solution, which involves both terms, has been analyzed in Refs.32 and 33.It is a periodic function with an in general high content of Fourier harmonics.As said before, its phase velocity v 0 is determined by the NDR (39).
In case of a Maxwellian unperturbed plasma, where we have to replace (a; b) by (a M ; b M ), the NDR can be written as 6,27-29,32,33 where we used (24) in the last step to show coincidence with the previous result.In the more general degenerate electron state, we replace (a; b) by (a FD ; b FD ) in Eq. ( 39) to obtain the NDR with b FD ðv 0 Þ given by (34).
Figure 3 shows plots of the NDR (42), where we used x 0 ¼ k 0 v 0 , for b ¼ À1 and different values of l.Both the zero amplitude limit w ' 0 þ in Fig. 3(a) and the small amplitude case w ¼ 0:1 in Fig. 3 the typical thumb shape of the dispersion curves in the (x 0 , k 0 ) space.In general, the wave frequency increases with increasing values of l.The high-frequency Langmuir branch takes the frequency x 0 ! 1 as k 0 !0 as can also be shown by using the NDR (39) together with the asymptotic expansion (33) for a ¼ a FD together with Eq. ( 17) for A, and that b ¼ b FD !0 in the limit v 0 ! 1.
We note that the character of the dispersion relation changes and becomes "non-dispersive" with a lower cut-off for k 0 in the high frequency branch when the right hand side of the NDR is fixed rather than b and w (see Fig. 5 of Ref. 33).In case of a Maxwellian plasma, denoting the right hand side of (41) by B M , this cut-off is given by k0 :¼ ffiffiffiffiffiffi ffi B M p and the phase velocity v 0 diverges as v 2 0 $ 1 A fixed B M then corresponds to a b ¼ À 15 ffiffi p p 16 ffiffi ffi w p B M expðv 2 0 =2Þ, which assumes extremely large negative values.In such a situation, the trapped electron range is more or less empty.

VIII. CONCLUSION
In conclusion, we have presented a general framework of weak coherent structures in partially Fermi-Dirac degenerate plasmas.These coherent structures are necessarily determined by the electron trapping nonlinearity and hence are beyond any realm of linear Vlasov theory.This implies that even in the infinitesimal amplitude limit, linear Landau and/or van Kampen theory fail as possible descriptive approaches to these structures. 29,33The spatial profiles of the potential /ðxÞ in the small amplitude limit are of cnoidal hole character, similar is in the classical case, and are given by the already known expressions (see Ref. 33 and references therein).The NDR on the other hand is for given b; w again of thumb shape type but appears to be rather sensitive to the electron degeneracy characterized by the normalized chemical potential l.An important aspect of the nonlinear theory in the small but finite amplitude limit is that the electron density is described by half-power expansions of the electrostatic potential.Different results 34 not giving rise to half-power expansions of the potential, which we find questionable, are found by not following the recipe of first defining trapped and free populations of electrons before calculating the electron density and by an incorrect reduction in the distribution function from 3D to 1D.
Finally, a mention on the validity of the present approach is necessary.We note that hot and non-degenerate plasmas are characterized by k n ) k dB , where k n ¼ ð4pn 0 =3Þ À1=3 is the mean particle distance and k dB ¼ h=ðm e v the Þ is the thermal de Broglie wavelength with the reduced Planck's constant h ¼ h=ð2pÞ, while cold and degenerate plasmas are characterized by k n Շ k dB .In the latter limit, the equilibrium distribution changes from a Maxwell-Boltzmann distribution to a Fermi-Dirac distribution.Moreover, the collisionless nature of the plasma requires that the mean kinetic energy of the particles is much larger than the mean potential energy due to Coulomb interactions between the particles.For a Maxwell-Boltzmann distributed plasma, this leads to the requirement 3k B T=2 > e 2 k n =ð4p 0 Þ, or k n < k De , while in the cold limit of Fermi-Dirac distributed electrons, we have 3E Fe =5 > e 2 k n =ð4p 0 Þ, which is fulfilled for sufficiently high densities n 0 տ a À3 0 where a 0 ¼ 4p 0 h 2 =ðm e e 2 Þ % 5:3 Â10 À11 m is the Bohr radius, giving n 0 տ 6:3 Â 10 30 m À3 .This is larger than typical densities of conduction electrons in metals, $5 Â 10 28 m À3 , but this restriction may be relaxed due to Pauli blocking reducing the collisionality of the electrons. 35,36The use of the Vlasov approach against the more general Wigner framework of quasi-distributions 6 necessitates the semi-classical limit DxDv ) h=m e , where Dx is the spatial length-scale of the structure and Dv is the velocity scale.As an example, if Dx $ 1=k 0 is the typical length-scale of a periodic wave-train and the velocity scale in an ultracold plasma is the Fermi speed Dv $ v Fe ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2E F =m e p , then the Wigner approach has to be used when the wavelength k ¼ 2p=k 0 is small enough to be comparable to the interparticle distance k n .If the velocity scale is instead given by the thermal speed v the in a warm, Maxwellian plasma, then the Wigner approach would only be necessary at extremely short wavelengths when k is comparable to k dB ( k n .The Maxwellian limit l !À1 is only a formal one, since relativistic effects have to be taken into account when k B T becomes comparable to the electron rest mass energy m e c 2 , which occurs at temperatures higher than $10 9 K, and in which case relativistic effects have to be taken into account in the treatment of electron holes. 25A similar constraint applies to very high densities when the electron Fermi energy E F becomes comparable to m e c 2 , which happens when the k n տ k C , where k C ¼ h=ðm e cÞ is the Compton wavelength; this happens at densities n 0 տ 10 34 m À3 that exist in dense astrophysical bodies such as white dwarf stars. 37,38lanned future extensions of the theory are to take the above mentioned relativistic effects into account.A detailed examination of the NDR (42) is also delegated to a forthcoming paper where analytic expressions will be presented for the l dependence of the Langmuir branch, x 0 $ 1; k 0 ( 1, and the slow electron acoustic branch x 0 / k 0 ; k 0 ( 1.As Fig. 3 shows, the slow branch is especially sensitive to l.Another extension is the replacement of the phase velocity v 0 by ṽD :¼ v D À v 0 , where v D is the drift velocity, in the above expressions to handle a current in the electronic system. 6,32,39,40 a) hans.schamel@uni-bayreuth.deb) bengt.eliasson@strath.ac.uk 1070-664X/2016/23(5)/052114/7 V C Author(s) 2016.23, 052114-1 PHYSICS OF PLASMAS 23, 052114 (2016)

FIG. 2 .
FIG.1.The linear coefficient a FD for different values of v 0 and l.Note that a FD ! a M in the classical limit l ( À1.
under the Reuse of AIP Publishing content is subject to the terms at: https://publishing.aip.org/authors/rights-and-permissions.Downloaded to IP: 130.159.82.198On: Thu, 02 Jun assumption that f e ðx; vÞ is continuous across the separatrix and that f t ðnÞ is an analytic function in n that ) ) Reuse of AIP Publishing content is subject to the terms at: https://publishing.aip.org/authors/rights-and-permissions.Downloaded to IP: 130.159.82.198On: Thu, 02 Jun