Kappa Distribution Function Effects on Landau Damping in Electrostatic Vlasov Simulation

The non-thermal high-energy electron effects on Langmuir wave-particle interaction are investigated using an initial value approach. A Vlasov-Poisson simulation is employed based on the splitting scheme by Cheng and Knorr (1976). The kappa distribution function is taken as an example of non-thermal electrons. The modification is manifested as an increase in the Landau damping rate and a decrease in the real frequency for a long wavelength limit. A part of the analyses using the modified plasma dispersion function (Summers and Thorne 1991) is reproduced for l = 2, 3 and 6. The dispersion relation from the initial value simulation and the plasma dispersion function compare favorably. (PACS numbers: 52.35.Fp, 52.35.Sb, 52.65.Ff).


INtrODUctION
Space plasma is far from being in a thermal equilibrium.Suprathermal electrons are often observed in space plasmas (Vasyliun 1968).The kappa distribution function (Leubner 2004) is one of the good examples of non-Maxwellian distribution functions.In one limit (Tsallis 1988) the kappa distribution function evolves toward Maxwellian.
To investigate wave-particle interaction, or Landau damping (Landau 1946;Jackson 1960) one needs to incorporate the velocity space dynamics using the Vlasov equation.One of the first pieces of work that made Vlasov equation numerical simulation available is the splitting scheme by Cheng and Knorr (1976) based on the method of characteristics.This method has become a standard method for the Vlasov type simulation.The method applies as long as the system is dissipation-less, in other words, if the phase volume of the system conserves.The Vlasov-type simulation in lower dimensional cases has advantages over Particle-in-Cell (PIC) simulation because the Vlasov simulation does not accompany statistical errors.The Vlasov simulation in lower dimension is suitable for investigating subtle effects such as a slight deviation in the equilibrium distribution function from Maxwellian.
One of our long term goals is to investigate the dynamics of Langmuir solitons (Zakharov 1972) which can then evolve into Langmuir turbulence (Wang et al. 1994(Wang et al. , 1995(Wang et al. , and 1996)).Using the splitting scheme (Cheng and Knorr 1976), the electrostatic Vlasov simulation has revealed plasma heating by Langmuir solitons (Lin et al. 1995).
The purpose of this paper is two-fold.The first purpose is to recapitulate the Cheng and Knorr (1976)'s method accurately for further advanced study of Langmuir solitons.In this work, the splitting scheme is revisited and the Cheng and Knorr (1976)'s results are reproduced.The second purpose is, as an initial exercise, to capture the non-Maxwellian distribution function effect on Landau damping.
This paper is organized as follows.Section 2 describes the basic computation model.Section 3, we start from verifying the simulation results with the free streaming case whose analytical solutions are known.The benchmark linear and nonlinear numerical simulation results are discussed in section 4. The non-Maxwellian kappa distribution function effects are discussed in section 5. Direct comparison of the simulation results with modified plasma dispersion function is discussed in section 6.We summarize this work in section 7.

MODEL EQUAtIONS AND NUMErIcAL MEtH-ODS
In this section, the model equation of the Vlasov-Poisson simulation is described.For the work transparency, we recapitulate Cheng and Knorr (1976) as precise as possible including the notations.A one-dimensional Vlasov-Poisson system in the MKS unit is given by and where the densities n j = n j (x, t) of each species (subscripts j = i for the ions and j = e for the electrons) are given by the distribution functions Here m j and q j are the mass and the charge of the species.The electric field is given by E and the vacuum permittivity is given by ε 0 .The coordinates in configuration and velocity spaces are given by x and υ, respectively.The equations are further normalized by the Debye length λ e = (ε 0 T e /n 0 m e ) 1/2 , the plasma frequency ω e = (n 0 e 2 /ε 0 m e ) 1/2 , and the electrostatic field is normalized by E = (eλ e E/T e ) which is equivalent to having the electrostatic energy being comparable to the electron thermal energy.Here, e is the unit charge and the electron temperature is given by T e .By employing the bars denoting the normalized values, we have x = x/λ e , t = ω e t, v = v/λ e ω e , and f = (λ e ω e /n 0 ) f, where n 0 is the equilibrium electron and ion densities.Note that the thermal velocity of the electrons is given by υ e = (T e /m e ) 1/2 = λ e ω e .
After the normalization, we obtain a Vlasov-Poisson system in the λ e (space) and ω e (time) scales, , , , , where the ion density is taken to be uniform [signified by unity in Eq. ( 5)] and the distribution function f is for the electrons.Equations ( 4) and (5) correspond to Eqs. (1a) and (1b) of Cheng and Knorr (1976).Hereafter we drop the bars.
The splitting scheme (Cheng and Knorr 1976) is based on the method of characteristics which is equivalent to the lowest order symplectic integrator (Ruth 1983).We evolve the distribution function by tracing the characteristic curves in the phase space.The method takes three steps which is given by and finally, whose kinetic energy part and the potential energy part are time advanced alternatively (the lowest order method corresponds to the well-known leap frog method employed frequently in PIC simulation).The superscript n stands for the time step.Note that the splitting scheme is not a finite difference scheme.In the splitting scheme, if the reference points along the characteristic curves "x -υΔt/2" or "υ + E(x)Δt" are exactly on the mesh points, the method is quite trivial.However, in general, the points of references are located in between mesh points of the x and v space.We thus need an interpolation technique to realize Eqs. ( 6), (7), and (8).As in Cheng and Knorr (1976), we use Fourier interpolation in the configuration space and linear interpolation in the velocity space (Watanabe 2005).The computational mesh we employed is exactly that of Fig. 1 in Cheng and Knorr (1976).Note that we do not have mesh points on the υ = 0 axis.
We need to solve the Poisson equation, Eq. ( 5) to make the simulation self-consistent.The Poisson equation is solved using the Fourier transform because we adopt a periodic configuration in x in this paper.All of the calculations in this paper employ periodic boundary conditions.The electric field is solved directly as in Eq. ( 5) without calculating the electrostatic potential.

FrEE StrEAMING cASE WItH E = 0
To begin with the numerical simulation, we start by verifying the solution of free streaming case [the Van Kampen mode, (Van Kampen 1955)] whose analytical solutions are known.Setting E = 0 in Eq. ( 4), we obtain , , , , Following Nicholson (1992), for example, the analytical solution can be given by setting , , ~exp For example, if we take an initial condition , , cos we obtain the analytical solution , , cos c os Shown in Fig. 1a is the solution of Vlasov equation in a free-streaming case.No force is acting on the electrons.The thick solid line is the initial condition of the distribution function at a fixed point x = 0.The time advanced distribution function by the numerical simulation is given by the black dots at t = 12.5 (time is normalized by the inverse of plasma frequency, ω e -1 ), which matches with the analytical solution given by the dash-dotted curves.Parameters employed in the simulation are the maximum cut-off velocity υ max = 10 υ e , and 0 ≤ x ≤ 4π.For the mesh points, 32 and 256 are taken in the x and v direction, respectively (although we did calculate, the regions y > 5.0 are not shown in the figure).Shown in Fig. 1b is the solution of Vlasov equation in the free-streaming case but with an initial condition given by , , The time advanced distribution function at t = 6.25 is given by the solid curve at x = π and by the dash-dotted curve at x = 3π.Note that, as demonstrated by the two solutions at x = π and x = 3π, the evolution of the distribution function exhibits point reflection across x = 2π.The distribution function streams in positive direction in υ > 0 while streams in negative direction in υ < 0, as a reminder.In both Figs.1a and b, A = 0.5 and k = 0.5 are taken.The calculation discussed in this section validates the interpolation scheme employed for the Vlasov equation.

bENcHMArK OF LINEAr AND NONLINEAr SIMULAtION rESULtS
In this section, linear and nonlinear simulation results are compared and benchmarked with those of Cheng and Knorr (1976).We take an initial condition of the form of Eq. ( 13) for both the linear and the nonlinear simulation.
The parameters employed in the linear simulation are exactly those of Cheng and Knorr (1976): the maximum cut-off velocity υ max = 4.0 υ e , and 0 ≤ x ≤ 4π, and the mesh points are 8 and 32 for x and υ, respectively.In Fig. 2, the Landau damping phase is shown in terms of the electric field strength E .As in Cheng and Knorr (1976) the recurrence  11) is taken as an initial condition.The thick solid line represents the initial distribution function.The time advanced distribution function is given by the black dots, while the analytical solution is given by the dash-dotted curve.(b) Time advanced distribution function at x = π (solid curve) and x = 3π (dash-dotted curve) when Eq. ( 13) is taken as an initial condition.The nonlinear simulation is shown in Fig. 3.In the simulation of Fig. 3, A = 0.5 and k = 0.5 are taken.Figure 3a shows the time evolution of the electric field at a fixed point x = π.Instead of amonotonic decrease the saturation of the amplitude can be seen after t = 20.Figure 3b shows the distribution function at t = 75 at a fixed point x = π.In Fig. 3b, we can see a local flattening of the distribution function in the vicinity of the phase velocity.The phase velocity of the Langmuir wave estimated using the linear theory is ω/k = 2.82.
Note that in the nonlinear simulation, as suggested in Cheng and Knorr (1976), the frequencies of all of the higher modes come into play at the later stage.As a result, resonance occurs at multiple locations in the velocity space and thus microscopic structures are generated [manifested as wrinkles in Cheng and Knorr (1976)] whose size can be comparable to the mesh size in the velocity space.To resolve all of the resonance one needs to employ an extremely high resolution in the velocity space.

EFFEctS OF KAPPA DIStrIbUtION FUNc-tIONS
In this section we investigate the effects of high energy  The spatial distribution is given in the form of Eq. ( 13) for the initial condition, thus f (x, υ, 0) ~ [1 + A cos (kx)] f υ (υ) is given.We have normalized the kappa distribution function to satisfy d f 1 y y = 3 3 y -^h # , so that the effective number of electrons will be the same as in the Maxwellian case.Some of the notable features of the kappa distribution function are shown in Fig. 4. Figure 4a compares Maxwellian (black) and kappa distributions with l = 1 (red).One can see a large population high-energy tail in the l = 1 case.The functions are plotted in logarithmic scales for l = 1 (red) , l = 2 (green), l = 5 (blue), and Maxwellian (black) cases in Fig. 4b.
Figure 5 shows the linear damping with cases when l = 1 (red), l = 2 (green), and Maxwellian (black) are taken as initial distribution functions.The parameters employed are υ max = 10 υ e , and 0 ≤ x ≤ 4π with 32 and 256 mesh points in x and υ.As in Fig. 2 From the linear theory (Landau 1946;Jackson 1960) the Landau damping rate is proportional to the slope of the distribution function at the phase velocity of the wave [see, for example, Nicholson (1992)].Employing the measured real frequencies ω of Fig. 6b (and k = 0.5 for the wave number), the slope of the distribution function f 2 y y y^h at the phase velocity ω/k is estimated and shown in Fig. 7. Smaller l cases have more negative values (and thus larger damping rate) which supports the nature of Fig. 6a.For small values of k, these latter trends (damping rate increasing and the real frequency decreasing) are consistent with the analytical work (Chateau and Meyer-Vernet 1991;Summers and Thorne 1991;Thorne and Summers 1991).In the next section we conduct a direct comparison of the Vlasov simulation with the roots of the plasma dispersion relation by taking exactly the same distribution functions employed in Summers and Thorne (1991) and Thorne and Summers (1991).

DIrEct cOMPArISON WItH MODIFIED PLASMA DISPErSION FUNctION
Using plasma dispersion function (Fried and Conte 1961) analogy for Maxwellian, Summers and Thorne [Summers and Thorne (1991); Thorne and Summers (1991)] extended their work to the kappa distribution function.We compare the damping rate and real frequencies in our simulation with their theoretical work for different values of к and different wave-numbers k.
We used exactly the same distribution function employed in Summers and Thorne [Summers and Thorne (1991); Thorne and Summers (1991)] to see the match between the two.We took initial distribution function in the form (see Appendix) Note the difference in the exponent part (" l -" dependence instead of " 1 l --") between Eqs. ( 16) and ( 14). Figure 8 plots the simulation results and numerical roots of the dispersion relation for l = 2, 3, and 6 [reproduced from the modified dispersion function of Summers and Thorne (1991) and Thorne and Summers (1991)].Figure 8a shows the damping rates versus the wave-numbers k. Figure 8b shows the real frequencies ω versus k.The roots from the (modified) dispersion relation are plotted as dashdotted curves.The black, red, and green dash-dotted curves are for l = 2, 3, and 6, respectively.
Those obtained from the linear numerical simulation are plotted as black circles.The damping rate and the real frequencies are obtained from electric field oscillation signals at a fixed point x = π.We performed the survey only up to k = 2.0.When the linear damping rate and real frequency magnitudes become comparable, the measurement for γ and ω becomes troublesome.The dispersion relation from the initial value simulation compares favorably with that from the plasma dispersion function.With smaller l values the real frequency decreases.Note that, however, with smaller к values the absolute value of the damping rates can also decrease for larger values of k, contrary to what we obtained  One of the advantages of the initial value approach is its application to nonlinear simulation.Our preliminary nonlinear simulation results employing a kappa distribution function as an initial condition are presented below.Figure 9 shows the electron distribution functions suggesting long time evolution up to t = 1500.The distribution functions are given at a fixed point x = π.Figure 9a is for a Maxwellian and Fig. 9b is for a kappa distribution function (l = 2).The dash-dotted curves are for t = 0 and the solid curves are for t = 1500.The integration of the distribution functions over the velocity space is the same for the two cases.The integration of the ion density over the configuration space is kept the same with the electron density (total numbers of ions and electrons in the system are the same).Figure 9 shows a relatively large cut-off velocity υ/υ the = 12.0 taken with a high resolution (1024 mesh points are taken in the velocity space).The l = 2 distribution function has a larger population at the high energy tail.In the simulation in Fig. 9, A = 0.5 and k = 0.5 are taken.Figure 9c shows the local expansion of Figs.9a and b near the resonant phase velocities.In the figure the perpendicular lines suggest the phase velocities.Note that the dash-dotted black line changes into solid black line (frequency down-shift for Maxwellian) while the dash-dotted red line changes into solid red line (frequency up-shift for kappa function).Over a very long time scale the normal mode frequency changes and the distribution func- tions evolve toward a similar equilibrium state.

SUMMArY
This work revisited the splitting scheme and the simulation results were compared with Cheng and Knorr (1976).Based on the code validation, as an initial exercise, we discussed the non-Maxwellian distribution function effect by employing the kappa-distribution function.
The slope [the absolute value of f 2 y y ^h] of the distribution function at the phase velocity was estimated, which supports the nature of increasing damping rate at smaller l values.The simulation results compared favorably with the analyses based on the modified plasma dispersion function (Summers and Thorne 1991;Thorne and Summers 1991).The specific calculations that we demonstrated are for l = 2, 3, and 6, with the wave-numbers k = 0.5, 1.0 and 2.0.Our preliminary nonlinear simulation employing the kappa distribution function was presented.Here, the black dash-dotted line at ω/k = 2.82 is for Maxwellian at t = 0, the black solid line at ω /k = 2.66 is for Maxwellian at t = 1500, the red dash-dotted line at ω/k = 2.32 is for kappa function at t = 0, and the red solid line at ω/k = 2.47 is for kappa function at t = 1500, respectively.The subscripts in the figure, max and kap stand for the Maxwellian and the kappa distribution function, respectively.

Fig. 1 .
Fig. 1.Evolution of distribution functions in free streaming cases.(a) Distribution functions at a cross section x = 0 when Eq. (11) is taken as an initial condition.The thick solid line represents the initial distribution function.The time advanced distribution function is given by the black dots, while the analytical solution is given by the dash-dotted curve.(b) Time advanced distribution function at x = π (solid curve) and x = 3π (dash-dotted curve) when Eq. (13) is taken as an initial condition.
after t = 42.Note that only the Fourier component k = 0.5 is kept and the other modes are filtered out in the linear calculation.A mesh point x = π is chosen for a diagnostic point in Fig.2.The figure corresponds to Fig.3ofCheng and Knorr (1976) except that A = 0.01 instead of A = 0.5 is taken.The measured frequency and the damping rates are ω = 1.41 and γ = -0.155,respectively.

Fig. 2 .
Fig. 2. Linear simulation results.Damping of the electric field strength E is shown.Here, υ max = 4 υ e , and 0 ≤ x ≤ 4π.The mesh points are 8 and 32 for x and υ, respectively.
Figure5shows the linear damping with cases when l = 1 (red), l = 2 (green), and Maxwellian (black) are taken as initial distribution functions.The parameters employed are υ max = 10 υ e , and 0 ≤ x ≤ 4π with 32 and 256 mesh points in x and υ.As in Fig.2, we have taken A = 0.01 and k = 0.5.With the kappa distribution function the Landau damping rate increases and the real frequency decreases (smaller l has larger effects).The variation in the linear damping rate γ and the real frequencies ω are summarized in Figs.6a and b as a function of l.The value l varies from 1.0 to 10.0.The values plotted in both Figs.6a and b are normalized by that of Maxwellian (thus in the figures, 1 max c c =-and 1 max ~~= for a Maxwellian initial condition).From the linear theory(Landau 1946; Jackson 1960) the Landau damping rate is proportional to the slope of the distribution function at the phase velocity of the wave [see, for example,Nicholson (1992)].Employing the measured real frequencies ω of Fig.6b(and k = 0.5 for the wave number), the slope of the distribution function f 2 y y y^h at the phase velocity ω/k is estimated and shown in Fig.7.Smaller l cases have more negative values (and thus larger damping rate) which supports the nature of Fig.6a.For small values of k, these latter trends (damping rate increasing and the real frequency decreasing) are consistent with the analytical work(Chateau and Meyer-Vernet 1991;Summers and Thorne 1991;Thorne and Summers 1991).In the next section we conduct a direct comparison of the Vlasov simulation with the roots of the plasma dispersion relation by taking exactly the same distribution functions employed inSummers and Thorne (1991) andThorne and Summers (1991).

Fig. 5 .
Fig. 5. Linear damping of electric field strength E comparing the cases when Maxwellian (black) and l = 1 (red) and l = 2 distribution functions are taken as initial conditions.

Fig. 6 .
Fig. 6.Effect of kappa function summarized: (a) damping rate normalized by the absolute value of damping rate with Maxwellian; (b) real frequency normalized by the value of frequency with Maxwellian.

Fig. 7 .
Fig. 7.The absolute value of f 2 y y ^h at the phase velocity.For the phase velocity ω/k, the values of ω are employed from Fig. 6b.The plotted values are normalized by that of the Maxwellian.

Fig. 8 .
Fig. 8.Comparison of the numerical roots of the dispersion relation and the simulation results: (a) the damping rates γ versus the wave vector k; (b) the real frequencies ω versus k.The roots from the (modified) dispersion relation are plotted as dash-dotted curves.The black, red, and green dashdotted curves are for l = 2, 3, and 6, respectively.Black circles are obtained from linear numerical simulation.Numerical simulation is done for k = 0.5, 1.0, and 2.0 for l = 2, 3, and 6.The simulation results are given at a fixed point x = π.

Acknowledgements
A part of this work is supported by National Cheng Kung University Top University Project and a part by the National Science Council of Taiwan, NSC 100-2112-M-006-021-MY3.One of the authors Yasutaro Nishimura would like to thank Professor Yasushi Nishida for discussions.The authors thank one of the Referees for suggesting several pieces of work on the kappa distribution function analysis.A part of computation is perfomed on ALPS Cluster at National Center for High-Performance Computing (NCHC), Taiwan.

Fig. 9 .
Fig. 9. Electron distribution functions suggesting long time evolution up to t = 1500.The distribution functions are given at a fixed point x = π.(a) For a Maxwellian, (b) for a kappa distribution function (l = 2).The dash-dotted curves are for t = 0 and the solid curves are for t = 1500.(c) Local expansion of (a) and (b) near the resonant phase velocities.Here, the black dash-dotted line at ω/k = 2.82 is for Maxwellian at t = 0, the black solid line at ω /k = 2.66 is for Maxwellian at t = 1500, the red dash-dotted line at ω/k = 2.32 is for kappa function at t = 0, and the red solid line at ω/k = 2.47 is for kappa function at t = 1500, respectively.The subscripts in the figure, max and kap stand for the Maxwellian and the kappa distribution function, respectively.