# Time evolutions of information entropies in a one-dimensional Vlasov-Poisson
system

K. Maekaku<EMAIL_ADDRESS>Graduate School of Frontier
Science, The University of Tokyo, Kashiwa 277-856, Japan  H. Sugama  National
Institute for Fusion Science, Toki 509-5292, Japan  Graduate School of
Frontier Science, The University of Tokyo, Kashiwa 277-856, Japan  T.-H.
Watanabe  Department of Physics, Nagoya University, Nagoya 464-8602, Japan

###### Abstract

A one-dimensional Vlasov-Poisson system is considered to elucidate how the
information entropies of the probability distribution functions of the
electron position and velocity variables evolve in the Landau damping process.
Considering the initial condition given by the Maxwellian velocity
distribution with the spatial density perturbation in the form of the cosine
function of the position, we derive linear and quasilinear analytical
solutions that accurately describe both early and late time behaviors of the
distribution function and the electric field. The validity of these solutions
is confirmed by comparison with numerical simulations based on contour
dynamics. Using the quasilinear analytical solution, the time evolutions of
the velocity distribution function and its kurtosis indicating deviation from
the Gaussian distribution are evaluated with the accuracy of the squared
perturbation amplitude. We also determine the time evolutions of the
information entropies of the electron position and velocity variables and
their mutual information. We further consider Coulomb collisions which relax
the state in the late-time limit in the collisionless process to the thermal
equilibrium state. In this collisional relaxation process, the mutual
information of the position and velocity variables decreases to zero while the
total information entropy of the phase-space distribution function increases
by the decrease in the mutual information and demonstrates the validity of
Boltzmann’s H-theorem.

††preprint: AIP/123-QED

##  I Introduction

Landau damping is one of the most intriguing physical processes involving
collective interactions between waves and particles in collisionless plasmas.
Landau  Despite the time reversibility of the Vlasov-Poisson equations, their
solution shows that, as time progresses, plasma oscillations are Landau-damped
and the electric field energy is converted into the kinetic energy of
particles. The effects of Landau damping are universally observed in various
plasma wave-particle resonance phenomena, such as drift waves and geodesic
acoustic mode (GAM) oscillations in slab and toroidal magnetic field
configurations in space and fusion plasmas, and have been the subject of
extensive theoretical research. Case ; VanKampen ; Nicholson ; Hazeltine-
Waelbroeck ; Prigogine ; Zocco ; FFChen ; Hammett ; Villani ; Zonca ;
Sugama2006a ; Sugama2006b ; Smolyakov

The most commonly used theoretical model to study Landau damping is the one-
dimensional Vlasov-Poisson system which consists of electrons with ions
treated as uniform background positive charge with infinite mass. The solution
of the linear initial value problem derived by Case and Van Kampen for this
system is well known. Case ; VanKampen ; Nicholson ; Hazeltine-Waelbroeck  The
behavior of the late-time solution exhibiting Landau damping can be well
approximated by using only eigenfrequencies with the slowest damping rate. In
this study, the effect of an infinite number of complex eigenfrequencies is
included to represent the early-time behavior of the electric field.

When treating the position $x$ and velocity $v$ of electrons as random
variables, denoted as $X$ and $V$ respectively, the information entropy
Information  $S_{p}(X,V)\equiv-\int\int p(x,v,t)\log p(x,v,t)\,dx\,dv$ derived
from the joint probability density function $p(x,v,t)$ of $(X,V)$ is known to
be one of the Casimir invariants in the Vlasov-Poisson system, Casimir ;
Maekaku  and it does not depend on time $t$. On the other hand, the entropies
$S_{p}(X)\equiv-\int p_{X}(x,t)\log p_{X}(x,t)\,dx$ and $S_{p}(V)\equiv-\int
p_{V}(v,t)\log p_{V}(v,t)\,dv$ determined from the marginal probability
density functions $p_{X}(x,t)=\int p(x,v)\,dv$ and $p_{V}(v,t)=\int
p(x,v)\,dx$ of $X$ and $V$ respectively, are allowed to vary with time. In the
present work, we derive novel analytical expressions that accurately represent
the early-time behaviors of the linear solutions for the electric field and
the distribution function by series expansions in time and velocity variables.
Combining these early-time expressions with the late-time approximate
solutions, we accurately determine the time evolution of the spatially
averaged velocity distribution function of electrons obtained from quasilinear
theory. Drummond ; Vedenov  Using these linear and quasilinear analytical
solutions, we clarify the time evolution of the information entropies
$S_{p}(X)$, $S_{p}(V)$, and the mutual information $I(X,V)\equiv
S_{p}(X)+S_{p}(V)-S_{p}(X,V)$ associated with the Landau damping.

We note that the quasilinear solution for the background velocity distribution
function at time $t$ is obtained by integrating the product of the electric
field and the linear solution of the perturbed distribution function over time
from the initial time to $t$. Therefore, using only the linear solutions
represented by a few complex eigenfrequencies is insufficient. We emphasize
that, to describe the time evolutions of the entropies and the mutual
information, it is necessary to use accurate expressions of the linear
solutions near the initial time as derived in this study.

The validity of the obtained linear and quasilinear solutions is confirmed by
comparison with simulation results based on contour dynamics. Zabusky ; Sato
We also obtain the velocity distribution $p_{V}$ in the limit of
$t\rightarrow+\infty$ and show how it deviates from the Gaussian distribution.
In addition, we consider the effects of Coulomb collisions, which relax the
distribution function to the thermal equilibrium state, decrease the mutual
information $I(X,V,t)$ to zero, increase the entropy $S_{p}(X,V)$, and thus
validate Boltzmann’s H-theorem. H-theorem

The rest of this paper is organized as follows. In Sec. II, the basic
equations of the Vlasov-Poisson system are presented, and the solution to its
linear initial value problem is provided using the Laplace transform. The
validity of the linear analytical solutions for the electric field and the
perturbed distribution function is confirmed by comparison with contour
dynamics simulations with small initial perturbation amplitudes. In Sec. III,
using approximate expressions for complex eigenfrequencies with large absolute
values, an approximate integral formula for the electric field incorporating
the effects of an infinite number of complex eigenfrequencies is derived. In
addition, the linear analytical solutions for the electric field and the
distribution function near the initial time are expressed as series expansions
in time and velocity variables. In Sec. IV, using the results of Secs. II and
III, a quasilinear analytical solution describing the time evolution of the
background velocity distribution function of electrons is derived, and its
validity is confirmed by contour dynamics simulations as well. By using this
quasilinear analytical solution, the time evolution of the electron kinetic
energy increases due to Landau damping, the velocity distribution function in
the limit of $t\rightarrow+\infty$, and physical quantities, such as the
kurtosis representing a deviation from the Gaussian distribution, are
accurately determined. In Sec. V, the time evolution of the information
entropies of the position and velocity variables of electrons and the mutual
information is determined with an accuracy of the order of the squared
perturbation amplitude. In Sec. VI, the changes in the entropies and the
mutual information from the collisionless process to the thermal equilibrium
state due to Coulomb collisions are evaluated. Finally, conclusions and
discussion are given in Sec. VII.

##  II Linear Analysis of Vlasov-Poisson system

###  II.1 Vlasov-Poisson system

Under the assumption that there is no magnetic field, the Vlasov equation for
electrons is given by

$\frac{\partial f({\bf x},{\bf v},t)}{\partial t}+{\bf v}\cdot\frac{\partial
f({\bf x},{\bf v},t)}{\partial{\bf x}}-\frac{e}{m}{\bf E}({\bf
x},t)\cdot\frac{\partial f({\bf x},{\bf v},t)}{\partial{\bf v}}=0,$ (1)

where $-e$ and $m$ are the electron charge and mass, respectively, and
$E(x,t)$ represents the electric field. The distribution function $f({\bf
x},{\bf v},t)$ of electrons is defined such that $f({\bf x},{\bf
v},t)d^{3}xd^{3}v$ represents the number of electrons in the phase space
volume element $d^{3}x\;d^{3}v$ around $({\bf x},{\bf v})$ at time $t$. We
ignore the motion of ions that have the uniform density $n_{0}$ and infinite
mass. Then, the electric field ${\bf E}({\bf x},t)$ is determined by Poisson’s
equation,

$\nabla\cdot{\bf E}({\bf x},t)=4\pi e\left(n_{0}-\int d^{3}v\;f({\bf x},{\bf
v},t)\right).$ (2)

We now assume that $f({\bf x},{\bf v},t)$ and ${\bf E}({\bf x},t)$ are
independent of $y$ and $z$, and that $E_{y}=E_{z}=0$. Integration of $f({\bf
x},{\bf v},t)$ with respect to $v_{y}$ and $v_{z}$ is done to define

$f(x,v_{x},t)\equiv\int_{-\infty}^{+\infty}dv_{y}\int_{-\infty}^{+\infty}dv_{z}\;f({\bf
x},{\bf v},t).$ (3)

Then, from Eq. (1), we obtain the Vlasov equation in the $(x,v)$ phase space
space as

$\frac{\partial f(x,v,t)}{\partial t}+v\frac{\partial f(x,v,t)}{\partial
x}-\frac{e}{m}E(x,t)\frac{\partial f(x,v,t)}{\partial v}=0,$ (4)

where $v=v_{x}$ and $E(x,t)=E_{x}(x,t)$ are used. Poisson’s equation in Eq.
(2) is rewritten as

$\frac{\partial E(x,t)}{\partial x}=4\pi
e\left(n_{0}-\int_{-\infty}^{+\infty}dv\;f(x,v,t)\right).$ (5)

Hereafter, we consider the structure of the distribution function $f(x,v,t)$
on the two-dimensional $(x,v)$ phase space instead of the six-dimensional
$({\bf x},{\bf v})$ space. The number of electrons in the area element
$dx\;dv$ about the point $(x,v)$ in the phase space at time $t$ is given by
$f(x,v,t)dx\;dv$.

###  II.2 Linearization

We now write the distribution function $f(x,v,t)$ as the sum of the
equilibrium and perturbation parts,

$f(x,v,t)=f_{0}(v)+f_{1}(x,v,t).$ (6)

Here the equilibrium distribution function is given by $f_{0}(v)$ which
satisfies

$\int d^{3}v\;f_{0}(v)=n_{0}.$ (7)

The perturbed distribution function $f_{1}(x,v,t)$ and the electric field
$E(x,t)$ are assumed to be given by

$\displaystyle f_{1}(x,v,t)$ $\displaystyle=$
$\displaystyle\mbox{Re}[f_{1}(k,v,t)\exp(i\;k\;x)],$ $\displaystyle E(x,t)$
$\displaystyle=$ $\displaystyle\mbox{Re}[E(k,t)\exp(i\;k\;x)],$ (8)

where the wavenumber in the $x$-direction is given by $k>0$. Substituting Eq.
(6) into Eq. (4) and using Eq. (II.2), we obtain the linearized Vlasov
equation,

$\frac{\partial f_{1}(k,v,t)}{\partial
t}+i\;k\;v\;f_{1}(k,v,t)-\frac{e}{m}E(k,t)\frac{\partial f_{0}(v)}{\partial
v}=0,$ (9)

where the nonlinear term

$-\frac{e}{m}E(x,t)\frac{\partial f_{1}(x,v,t)}{\partial v}$ (10)

is neglected as a small term. Using Eqs. (6) and (7), Poisson’s equation in
Eq. (5) is rewritten as

$i\>k\;E(k,t)=-4\pi e\int_{-\infty}^{+\infty}dv\;f_{1}(k,v,t).$ (11)

###  II.3 Laplace transform

We now use the Laplace transform to solve the linearized Vlasov-Poisson
equations given by Eq. (9) and (11). The Laplace transforms of $f_{1}(k,v,t)$
and $E(k,t)$ are denoted by $f_{1}(k,v,\omega)$ and $E(k,\omega)$,
respectively. The inverse Laplace transform gives

$f_{1}(k,v,t)=\int_{L_{f}}\frac{d\omega}{2\pi}f_{1}(k,v,\omega)e^{-i\omega
t},$ (12)

and

$E(k,t)=\int_{L_{E}}\frac{d\omega}{2\pi}E(k,\omega)e^{-i\omega t},$ (13)

where the Laplace contours $L_{f}$ and $L_{E}$ need to pass above all poles of
$f_{1}(k,v,\omega)$ and $E(k,\omega)$ on the complex $\omega$-plane,
respectively. Now, the Vlasov equation and Poisson’s equation in Eqs. (9) and
(11) are represented in terms of $f_{1}(k,v,\omega)$ and $E(k,\omega)$ by

$(-i\;\omega+i\;k\;v)f_{1}(k,v,\omega)-\frac{e}{m}E(k,\omega)\frac{\partial
f_{0}(v)}{\partial v}=f_{1}(k,v,t=0),$ (14)

and

$i\;k\;E(k,\omega)=-4\pi e\int_{-\infty}^{+\infty}dv\;f_{1}(k,v,\omega),$ (15)

respectively.

Solving Eqs. (14) and (15) for the perturbed distribution function
$f_{1}(k,v,\omega)$ and the electric field $E(k,\omega)$, we obtain

$f_{1}(k,v,\omega)=\frac{(e/m)E(k,\omega)\partial f_{0}(v)/\partial
v+f_{1}(k,v,t=0)}{-i\;\omega+i\;k\;v},$ (16)

and

$E(k,\omega)=\frac{4\pi
e}{k^{2}\epsilon(k,\omega)}\int_{-\infty}^{+\infty}dv\;\frac{f_{1}(k,v,t=0)}{v-\omega/k},$
(17)

where the dielectric function $\epsilon(k,\omega)$ is defined by

$\epsilon(k,\omega)=1-\frac{\omega_{p}^{2}}{n_{0}k^{2}}\int_{-\infty}^{+\infty}dv\;\frac{\partial
f_{0}(v)/\partial v}{v-\omega/k}.$ (18)

Here, the plasma frequency is defined by $\omega_{p}\equiv(4\pi
n_{0}e^{2}/m)^{1/2}$. Now, we define

$F(k,\omega)\equiv\int_{-\infty}^{+\infty}dv\;\frac{f_{1}(k,v,t=0)}{v-\omega/k},$
(19)

to write $E(k,t)$ and $f_{1}(k,v,t)$ for $t>0$ as

$E(k,t)=\frac{4\pi e}{k^{2}}\int_{L_{E}}\frac{d\omega}{2\pi}e^{-i\omega
t}\frac{F(k,\omega)}{\epsilon(k,\omega)},$ (20)

and

$\displaystyle f_{1}(k,v,t)$ $\displaystyle=$
$\displaystyle\int_{L_{f}}\frac{d\omega}{2\pi}\frac{e^{-i\omega
t}}{-i(\omega-k\,v)}\times$ (21)
$\displaystyle\left[\frac{\omega_{p}^{2}}{n_{0}k^{2}}\frac{F(k,\omega)}{\epsilon(k,\omega)}\frac{\partial
f_{0}(v)}{\partial v}+f_{1}(k,v,t=0)\right],$

respectively. Let us define complex-valued eigenfrequencies
$\\{\omega_{\mu}\\}$ as zeros of $\epsilon(k,\omega)$,

$\epsilon(k,\omega_{\mu})=0.$ (22)

Then, we can write

$E(k,t)=-i\frac{4\pi
e}{k^{2}}\sum_{\mu}e^{-i\omega_{\mu}t}\frac{F(k,\omega_{\mu})}{\partial_{\omega}\epsilon(k,\omega_{\mu})},$
(23)

and

$\displaystyle f_{1}(k,v,t)$ $\displaystyle=$ $\displaystyle
e^{-ik\,v\,t}\left[f_{1}(k,v,t=0)+\frac{\omega_{p}^{2}}{n_{0}k^{2}}\frac{F(k,k\,v)}{\epsilon(k,k\,v)}\frac{\partial
f_{0}(v)}{\partial v}\right]$ (24)
$\displaystyle\mbox{}+\frac{\omega_{p}^{2}}{n_{0}k^{2}}\frac{\partial
f_{0}(v)}{\partial
v}\sum_{\mu}\frac{e^{-i\omega_{\mu}t}}{(\omega_{\mu}-k\,v)}\frac{F(k,\omega_{\mu})}{\partial_{\omega}\epsilon(k,\omega_{\mu})},$

where the derivative with respect to the complex-valued frequency $\omega$ is
represented by $\partial_{\omega}\equiv\partial/\partial\omega$.

###  II.4 Conditions of equilibrium and initial perturbation

Hereafter, we assume that the equilibrium function $f_{0}(v)$ is given by the
Maxwellian

$f_{0}(v)=n_{0}\sqrt{\frac{m}{2\pi
T}}\exp\left(-\frac{mv^{2}}{2T}\right)=\frac{n_{0}}{\sqrt{\pi}v_{T}}\exp\left(-\frac{v^{2}}{v_{T}^{2}}\right),$
(25)

where $T$ represents the equilibrium temperature and

$v_{T}\equiv\sqrt{2}v_{t}\equiv\sqrt{2T/m}$ (26)

is used. Then, the dielectric function $\epsilon(k,\omega)$ in Eq. (18) is
expressed by

$\epsilon(k,\omega)=1+\frac{1}{k^{2}\lambda_{D}^{2}}\left[1+\frac{\omega}{kv_{T}}Z\left(\frac{\omega}{kv_{T}}\right)\right],$
(27)

where the Debye length $\lambda_{D}$ is defined by

$\lambda_{D}\equiv\frac{v_{t}}{\omega_{p}}\equiv\sqrt{\frac{T}{4\pi
n_{0}e^{2}}}$ (28)

and the plasma dispersion function $Z(\zeta)$ is defined by Eq. (162) in
Appendix A.

Using Eq. (27), the dispersion relation in Eq. (22) is rewritten as

$\epsilon(k,\omega_{\mu})=1+\kappa^{-2}\left[1+\zeta_{\mu}Z(\zeta_{\mu})\right]=0,$
(29)

where $\zeta_{\mu}=\omega_{\mu}/kv_{T}$ and $\kappa\equiv k\lambda_{D}$.
Figure 1 shows the distribution of complex eigenfrequencies
$\omega_{\mu}=\omega_{\mu r}+i\gamma_{\mu}$ in the complex plane, calculated
from Eq. (29) for $k\lambda_{D}=1/2$. There are an infinite number of complex
eigenfrequencies. The complex eigenfrequency with $\omega_{r}>0$ and the
smallest $|\gamma|$ is given by
$\omega_{\mbox{\o}}=(1.41566-0.153359i)\omega_{p}$. From Eqs. (29) and (166),
we find

$[\epsilon(k,-\omega^{*})]^{*}=\epsilon(k,\omega),$ (30)

and

$\epsilon(k,\omega_{\mu})=0\hskip 14.22636pt\Rightarrow\hskip
14.22636pt\epsilon(k,-\omega_{\mu}^{*})=0,$ (31)

where ∗ represents the complex conjugate. Thus, if $\omega_{\mu}=\omega_{\mu
r}+i\gamma_{\mu}$ is a zero of the dielectric function,
$-\omega_{\mu}^{*}=-\omega_{\mu r}+i\gamma_{\mu}$ is so, too. Therefore,
complex eigenfrequencies are distributed in pairs symmetric with respect to
the ${\rm Im}(\omega)$-axis as seen in Fig. 1.

Figure 1: Complex eigenvalues of frequency $\omega=\omega_{r}+i\gamma$
($k\lambda_{D}=1/2$). The horizontal and vertical axes represent the real part
$\omega_{r}$ and the imaginary part (or the growth rate) $\gamma$,
respectively. The complex eigenfrequency with $\omega_{r}>0$ and the smallest
$|\gamma|$ is given by $\omega_{\mbox{\o}}=(1.41566-0.153359i)\omega_{p}$.
There are an infinite number of eigenfrequencies other than
$\omega_{\mbox{\o}}$.

The derivative of $\epsilon(k,\omega)$ with respect to $\omega$ is given by

$\frac{\partial\epsilon(k,\omega)}{\partial\omega}=\frac{1}{kv_{T}}\frac{\partial\epsilon}{\partial\zeta}=\frac{1}{kv_{T}}\frac{1}{k^{2}\lambda_{D}^{2}}\left[-2\zeta+(1-2\zeta^{2})Z(\zeta)\right]$
(32)

where Eqs. (27), (167), and $\zeta=\omega/kv_{T}$ are used. Then, from Eq.
(29) and (32), we find

$\displaystyle\partial_{\omega}\epsilon(k,\omega_{\mu})$ $\displaystyle=$
$\displaystyle\frac{1}{\omega_{\mu}}\left[-\left(1+\frac{1}{k^{2}\lambda_{D}^{2}}\right)+2\zeta_{\mu}^{2}\right]$
$\displaystyle=$
$\displaystyle\frac{1}{k^{2}\lambda_{D}^{2}\omega_{\mu}}\left[-1-k^{2}\lambda_{D}^{2}+\frac{\omega_{\mu}^{2}}{\omega_{p}^{2}}\right]$

We now impose the initial condition given by

$f_{1}(k,v,t=0)=\alpha f_{0}(v),$ (33)

where $f_{0}(v)$ is the Maxwellian equilibrium distribution function given by
Eq. (25) and $\alpha$ is a small constant. Then, from Poisson’s equation, we
have

$E(k,t=0)=i\frac{4\pi n_{0}e}{k}\alpha=i\frac{m}{e}\frac{\omega_{p}}{k}\alpha$
(34)

Using the initial condition in Eq. (33), the plasma dispersion function in Eq.
(162), and Eq. (19), we obtain

$F(k,\omega)\equiv\int_{-\infty}^{+\infty}dv\;\frac{f_{1}(k,v,t=0)}{v-\omega/k}=\alpha\frac{n_{0}}{v_{T}}Z\left(\frac{\omega}{kv_{T}}\right),$
(35)

and

$F(k,kv)\equiv\int_{-\infty}^{+\infty}dv^{\prime}\;\frac{f_{1}(k,v^{\prime},t=0)}{v^{\prime}-v}=\alpha\frac{n_{0}}{v_{T}}Z\left(\frac{v}{v_{T}}\right),$
(36)

which are used to evaluate $E(k,t)$ and $f_{1}(k,v,t)$ given in Eqs. (23) and
(24), respectively.

###  II.5 Solution of initial value problem

We now substitute Eqs. (II.4), (35), and (36) into Eqs. (23) and (24) to
express the electric field and the perturbed distribution function at time
$t>0$ as

$\displaystyle E(k,t)$ $\displaystyle=$ $\displaystyle\frac{4\pi
e}{k^{2}}\int_{L_{E}}\frac{d\omega}{2\pi}e^{-i\omega
t}\frac{F(k,\omega)}{\epsilon(k,\omega)}$ (37) $\displaystyle=$
$\displaystyle\alpha\frac{4\pi
en_{0}}{k}\int_{-\infty}^{+\infty}\frac{d\zeta}{2\pi}\;\frac{Z(\zeta)e^{-i\zeta\tau}}{1+\kappa^{-2}[1+\zeta
Z(\zeta)]}$ $\displaystyle=$ $\displaystyle-i\frac{4\pi
e}{k^{2}}\sum_{\mu}e^{-i\omega_{\mu}t}\frac{F(k,\omega_{\mu})}{\partial_{\omega}\epsilon(k,\omega_{\mu})}$
$\displaystyle=$
$\displaystyle-i\;\alpha\frac{m\omega_{p}^{2}}{ek}k^{2}\lambda_{D}^{2}(1+k^{2}\lambda_{D}^{2})$
$\displaystyle\times\sum_{\mu}\frac{e^{-i\omega_{\mu}t}}{1+k^{2}\lambda_{D}^{2}-\omega_{\mu}^{2}/\omega_{p}^{2}}$
$\displaystyle=$
$\displaystyle-i\;\alpha\frac{m\omega_{p}^{2}}{ek}k^{2}\lambda_{D}^{2}(1+k^{2}\lambda_{D}^{2})$
$\displaystyle\times\sum_{\mu,\;{\rm Re}\omega_{\mu}>0}2{\rm
Re}\left[\frac{e^{-i\omega_{\mu}t}}{1+k^{2}\lambda_{D}^{2}-\omega_{\mu}^{2}/\omega_{p}^{2}}\right],$

and

$\displaystyle f_{1}(k,v,t)=$ (38)
$\displaystyle\int_{-\infty}^{+\infty}\frac{d\omega}{2\pi}\frac{e^{-i\omega
t}}{-i(\omega-k\,v)}\Biggl{[}f_{1}(k,v,t=0)$
$\displaystyle+\frac{\omega_{p}^{2}}{n_{0}k^{2}}\frac{F(k,\omega)}{\epsilon(k,\omega)}\frac{\partial
f_{0}(v)}{\partial v}\Biggr{]}$ $\displaystyle=\alpha
f_{0}(v)\Biggl{[}\exp\left(-i\frac{v}{v_{T}}\tau\right)$ $\displaystyle-$
$\displaystyle
i\kappa^{-2}\frac{v}{v_{T}}\int_{-\infty}^{+\infty}\frac{d\zeta}{2\pi}\;\frac{Z(\zeta)e^{-i\zeta\tau}}{(\zeta-v/v_{T})\\{1+\kappa^{-2}[1+\zeta
Z(\zeta)]\\}}\Biggr{]}$
$\displaystyle=e^{-ik\,v\,t}\Biggl{[}f_{1}(k,v,t=0)+\frac{\omega_{p}^{2}}{n_{0}k^{2}}\frac{F(k,k\,v)}{\epsilon(k,k\,v)}\frac{\partial
f_{0}(v)}{\partial v}\Biggr{]}$ $\displaystyle+$
$\displaystyle\frac{\omega_{p}^{2}}{n_{0}k^{2}}\frac{\partial
f_{0}(v)}{\partial
v}\sum_{\mu}\frac{e^{-i\omega_{\mu}t}}{(\omega_{\mu}-k\,v)}\frac{F(k,\omega_{\mu})}{\partial_{\omega}\epsilon(k,\omega_{\mu})}$
$\displaystyle=\alpha
f_{0}(v)\Bigg{[}e^{-ik\,v\,t}\left\\{1-\frac{(v/v_{T})Z(v/v_{T})}{k^{2}\lambda_{D}^{2}+[1+(v/v_{T})Z(v/v_{T})]}\right\\}$
$\displaystyle-$ $\displaystyle
k\;v\;(1+k^{2}\lambda_{D}^{2})\sum_{\mu}\frac{e^{-i\omega_{\mu}t}}{(\omega_{\mu}-k\,v)(1+k^{2}\lambda_{D}^{2}-\omega_{\mu}^{2}/\omega_{p}^{2})}\Bigg{]},$

respectively, where the normalized time $\tau\equiv kv_{T}t$ and the residue
theorem are used.

Figure 2 shows $E(k,t)$ calculated as a function of time using Eq. (37) for
$k\lambda_{D}=1/2$. The summation over $\mu$ in Eq. (37) is done by including
a finite number of pairs of the complex eigenfrequencies
$(\omega_{\mu},-\omega_{\mu}^{*})$ from
$(\omega_{\mbox{\o}},-\omega_{\mbox{\o}}^{*})$ to the pair with the $N$th
smallest decay rate $|\gamma_{\mu}|$. The results obtained for the cases of
$N=2$ and $N=100$ are shown in Fig. 2 where the result from the contour
dynamics (CD) simulation for the case of $\alpha=0.01$ are also plotted by the
red dashed curve. The CD method and the simulation conditions used in this
paper are explained in Appendix C. A good agreement among the three results is
confirmed except near $t=0$. When using the larger number $N$ of the pairs of
complex eigenfrequencies, the better agreement with the CD simulation result
is confirmed. Near $t=0$, even in the large $N$ limit, the sum over $\mu$ does
not converge uniformly, and the number of pairs $N$ required for a good
convergence increases to infinity as $t$ approaches 0. In Sec. III, we derive
an approximate expression for $E(k,t)$ near $t=0$ that includes the effect of
an infinite number of complex eigenfrequencies $\\{\omega_{\mu}\\}$.

Figure 2: Time evolution of $E(k,t)/E(k,0)$ for $k\lambda_{D}=1/2$. The red
dashed curve represents the result from the CD simulation using $\alpha=0.01$.
The solid curves are obtained from Eq. (37) where the summation over $\mu$ is
done by including a finite number of pairs of the complex eigenfrequencies
$(\omega_{\mu},-\omega_{\mu}^{*})$ from
$(\omega_{\mbox{\o}},-\omega_{\mbox{\o}}^{*})$ to the pair with the $N$th
smallest decay rate $|\gamma_{\mu}|$. Here, the case for $N=2$ and 100 are
shown.

Similar to the case of $E(k,t)$, the effect of a larger number of the complex
eigenfrequencies $\omega_{\mu}$ on $f_{1}(k,v,t)$ becomes more significant
near $t=0$. To accurately evaluate $f_{1}(k,v,t)$ from Eq. (38), we need to
include a larger number $N$ of complex eigenfrequency pairs as
$t\rightarrow+0$. However, we should note that the denominator of each term in
the sum over $\mu$ in the analytical solution for $E(k,t)$ [Eq. (37)] is a
quadratic function of $\omega_{\mu}$ while in the analytical solution for
$f_{1}(k,v,t)$ [Eq. (38)], it is a cubic function of $\omega_{\mu}$.
Therefore, $f_{1}(k,v,t)$ converges faster than $E(k,t)$ as the number $N$ of
included eigenfrequency pairs $\omega_{\mu}$ increases, and $f_{1}(k,v,t)$
shows a uniform convergence even near $t=0$. Figure 3 shows contour plots of
$f_{1}(x,v,t)/(\alpha^{2}f_{0}(v))$ in the $(v,t)$-plane for
$\omega_{p}t=0,0.1,1,10$. For $t>0$, the plots show
$f_{1}(x,v,t)/(\alpha^{2}f_{0}(v))$ calculated using 2 pairs (upper row) and
100 pairs (middle row) of complex eigenfrequencies from Eq. (38) and
$f_{1}(x,v,t)/(\alpha^{2}f_{0}(v))$ from CD simulation results for
$\alpha=0.01$ (lower row). The results using 2 pairs of complex
eigenfrequencies agree with the results of CD simulation for $\omega_{p}t\geq
1$. The results using 1000 pairs show a better agreement with CD simulation
results for all $t>0$.

Figure 3:  Contours of $f_{1}(x,v,t)/(\alpha^{2}f_{0}(v))$ in the
$(x,v)$-space for the case of $k\lambda_{D}=1/2$. The results obtained from
Eq. (38) using 2 pairs (upper row) and 100 pairs (middle row) of complex
eigenfrequencies and those from the CD simulation using $\alpha=0.01$ (lower
row for $t>0$) are shown.

##  III Approximate representation of linear solutions in early time

###  III.1 Approximate expression for complex eigenfrequencies with large
absolute values

We here recall that the complex-valued eigenfrequency $\omega_{\mu}$ is
determined by Eq. (29) using the plasma dispersion function $Z(\zeta)$. When
$|\zeta|\gg 1$, $|{\rm Re}\zeta|>|{\rm Im}\zeta|$, and ${\rm Im}\zeta<0$, we
have

$1+\zeta Z(\zeta)=i\sqrt{\pi}\zeta e^{-\zeta^{2}}+{\cal O}(\zeta^{-2}).$ (39)

Therefore, it is found from Eqs. (29) and (39) that, when $|\zeta_{\mu}|\gg 1$
and ${\rm Im}\zeta_{\mu}<0$, $\zeta_{\mu}$ satisfies

$i\sqrt{\pi}\zeta_{\mu}e^{-\zeta_{\mu}^{2}}=-\kappa^{2}+{\cal
O}(\zeta_{\mu}^{-2}).$ (40)

Using Eq. (40) and $\zeta_{\mu}\equiv|\zeta_{\mu}|e^{i\theta_{\mu}}$, we
obtain

$\zeta_{\mu}\simeq
i\frac{\kappa^{2}}{\sqrt{\pi}}e^{\zeta_{\mu}^{2}}=i\frac{\kappa^{2}}{\sqrt{\pi}}\exp[|\zeta_{\mu}|^{2}(\cos
2\theta_{\mu}+i\sin 2\theta_{\mu})],$ (41)

from which we get

$|\zeta_{\mu}|\simeq\frac{\kappa^{2}}{\sqrt{\pi}}\exp(|\zeta_{\mu}|^{2}\cos
2\theta_{\mu}),$ (42)

and

$e^{i\theta_{\mu}}\simeq i\exp(i|\zeta_{\mu}|^{2}\sin 2\theta_{\mu}).$ (43)

Therefore, $\theta_{\mu}$ satisfies

$\theta_{\mu}\simeq\frac{\pi}{2}+|\zeta_{\mu}|^{2}\sin 2\theta_{\mu}+2\pi n,$
(44)

where $n$ is an integer. Then, we define $\delta_{\mu}$ by

$2\theta_{\mu}=-\frac{\pi}{2}+\delta_{\mu}$ (45)

to have

$|\zeta_{\mu}|\simeq\frac{\kappa^{2}}{\sqrt{\pi}}\exp(|\zeta_{\mu}|^{2}\sin\delta_{\mu}),$
(46)

and

$\log\left(\frac{\sqrt{\pi}}{\kappa^{2}}|\zeta_{\mu}|\right)\simeq|\zeta_{\mu}|^{2}\sin\delta_{\mu}\simeq|\zeta_{\mu}|^{2}\delta_{\mu},$
(47)

which leads to

$\delta_{\mu}\simeq\frac{1}{|\zeta_{\mu}|^{2}}\log\left(\frac{\sqrt{\pi}}{\kappa^{2}}|\zeta_{\mu}|\right),$
(48)

where $0<\delta_{\mu}\ll 1$ is assumed. Substituting Eq. (45) into Eq. (44),
we have

$\displaystyle-\frac{\pi}{4}+\frac{\delta_{\mu}}{2}$ $\displaystyle\simeq$
$\displaystyle\frac{\pi}{2}-|\zeta_{\mu}|^{2}\cos\delta_{\mu}+2\pi n$ (49)
$\displaystyle\simeq$
$\displaystyle\frac{\pi}{2}-|\zeta_{\mu}|^{2}\left(1-\frac{\delta_{\mu}^{2}}{2}\right)+2\pi
n,$

which yields

$|\zeta_{\mu}|^{2}\simeq\frac{3\pi}{4}-\frac{\delta_{\mu}}{2}+2\pi n.$ (50)

Thus, the approximate value of $\zeta_{\mu}\equiv|\zeta_{\mu}|e^{i\theta}$ for
$|\zeta_{\mu}|\gg 1$ is given by

$|\zeta_{\mu}|\simeq\sqrt{\left(2n+\frac{3}{4}\right)\pi}$ (51)

with

$\theta_{\mu}=-\frac{\pi}{4}+\frac{\delta_{\mu}}{2}\simeq-\frac{\pi}{4}+\frac{1}{2|\zeta_{\mu}|^{2}}\log\left(\frac{\sqrt{\pi}}{\kappa^{2}}|\zeta_{\mu}|\right)$
(52)

where the integer $n$ needs to be large, $n\gg 1$, to satisfy the condition
$|\zeta_{\mu}|\gg 1$. Recall that when $\omega_{\mu}$ is a complex-valued
eigenfrequency, $-\omega_{\mu}^{*}$ is so, too. Then, the approximate values
of complex-valued eigenfrequency with large absolute values are denoted by
$\omega_{n}$ and $-\omega_{n}^{*}$ which are given by

$\frac{\omega_{n}}{kv_{T}}\equiv\zeta_{n}\equiv\sqrt{\left(2n+\frac{3}{4}\right)\pi}\;\;e^{i\theta_{n}},$
(53)
$\frac{-\omega_{n}^{*}}{kv_{T}}\equiv-\zeta_{n}^{*}\equiv\sqrt{\left(2n+\frac{3}{4}\right)\pi}\;\;e^{i(\pi-\theta_{n})}$
(54)

and

$\theta_{n}\equiv-\frac{\pi}{4}+\frac{1}{2|\zeta_{n}|^{2}}\log\left(\frac{\sqrt{\pi}}{\kappa^{2}}|\zeta_{n}|\right)$
(55)

where the integer $n$ is used instead of $\mu$ to specify different
approximate eigenfrequencies. The approximation is good for $n\gg 1$.

In Fig. 4, the exact complex eigenfrequencies and the approximate
eigenfrequencies for $k\lambda_{D}=1/2$ are plotted as red and blue dots,
respectively, in the region where ${\rm Re}(\omega)>0$. Here, Eqs. (53) and
(55) with $n=0,1,2,\cdots$ are used to evaluate the approximate
eigenfrequencies. It is observed that the approximate values approach the
exact values as the magnitude of the complex eigenfrequencies increases.

Figure 4: Distributions of the exact values (red dots) and the approximate
values (blue dots) of the complex eigenfrequencies for $k\lambda_{D}=1/2$ in
the region where $\omega_{r}>0$. Here, Eqs. (53) and (55) with
$n=0,1,2,\cdots$ are used to evaluate the approximate eigenfrequencies.

###  III.2  Approximate expression for effects of an infinite number of
complex-valued eigenfrequencies

We here consider an infinite series, $\sum_{n=N}^{\infty}{\rm
Re}[e^{i\zeta_{n}\tau}/(C-\zeta_{n}^{2})]$, where $C$ is a positive constant
and $\tau\equiv kv_{T}t$. This form of the infinite series is included in Eq.
(37) for $E(k,t)$. When $n\gg 1$ and $|\zeta_{n}|\gg C$, we use Eqs. (53) and
(55) to make the following approximations,

$\displaystyle{\rm Re}\left[\frac{e^{i\zeta_{n}\tau}}{C-\zeta_{n}^{2}}\right]$
(56) $\displaystyle\simeq$
$\displaystyle-\frac{1}{|\zeta_{n}|^{2}}\exp(|\zeta_{n}|\tau\sin\theta_{n})\cos(2\theta_{n}+|\zeta_{n}|\tau\cos\theta_{n})$
$\displaystyle\simeq$
$\displaystyle-\frac{1}{|\zeta_{n}|^{2}}\exp\left(-\frac{|\zeta_{n}|\tau}{\sqrt{2}}\right)\sin\left(\delta_{n}+\frac{|\zeta_{n}|\tau}{\sqrt{2}}\right),$

where

$|\zeta_{n}|\equiv\sqrt{\left(2n+\frac{3}{4}\right)\pi},\hskip
14.22636pt\delta_{n}\equiv
2\theta_{n}+\frac{\pi}{2}\equiv\frac{1}{|\zeta_{n}|^{2}}\log\left(\frac{\sqrt{\pi}}{\kappa^{2}}|\zeta_{n}|\right).$
(57)

For $\tau=0$, using Eq. (56) and $|\delta_{n}|\ll 1$ gives

${\rm
Re}\left[\frac{1}{C-\zeta_{n}^{2}}\right]\simeq-\frac{\delta_{n}}{|\zeta_{n}|^{2}}.$
(58)

For $\tau>0$, we can take $N\gg 1$ such that
$|\zeta_{N}|\tau/\sqrt{2}\gg\delta_{N}>0$. Then, Eq. (56) is rewritten as

${\rm
Re}\left[\frac{e^{i\zeta_{n}\tau}}{C-\zeta_{n}^{2}}\right]\simeq-\frac{1}{|\zeta_{n}|^{2}}\exp\left(-\frac{|\zeta_{n}|\tau}{\sqrt{2}}\right)\sin\left(\frac{|\zeta_{n}|\tau}{\sqrt{2}}\right)$
(59)

where $n\geq N$ and $|\zeta_{n}|^{2}\propto n$. We find that, for a large but
fixed integer $N$,

$\displaystyle S_{N}^{N^{\prime}}(\tau)$ $\displaystyle\equiv$
$\displaystyle\sum_{n=N}^{N^{\prime}}{\rm
Re}\left[\frac{e^{i\zeta_{n}\tau}}{C-\zeta_{n}^{2}}\right]$ (60)
$\displaystyle\simeq$
$\displaystyle-\sum_{n=N}^{N^{\prime}}\frac{1}{|\zeta_{n}|^{2}}\exp\left(-\frac{|\zeta_{n}|\tau}{\sqrt{2}}\right)\sin\left(\frac{|\zeta_{n}|\tau}{\sqrt{2}}\right)$

does not uniformly converge to
$S_{N}(\tau)\equiv\lim_{N^{\prime}\rightarrow\infty}S_{N}(\tau)$ because
$|S_{N}^{N^{\prime}}(\tau)-S_{N}(\tau)|<\epsilon$ requires $N^{\prime}\geq
N^{\prime}(\epsilon,\tau)$ where
$N^{\prime}(\epsilon,\tau)\propto|\zeta_{N^{\prime}(\epsilon,\tau)}|^{2}\propto\delta(\epsilon)/\tau^{2}$
which diverges as $\tau\rightarrow+0$. Therefore, the number of the complex-
valued eigenfrequencies, which are necessary to be included for accurately
evaluating $E(k,t)$, diverges to infinity as $t\rightarrow+0$. On the other
hand, when $v$ is fixed, the series expansion for representing $f_{1}(k,v,t)$
uniformly converges. It is because the extra factor $1/(\omega_{n}-kv)$
included in the series expansion accelerates the convergence as
$n\rightarrow\infty$. Now, we use

$\displaystyle 2\pi$ $\displaystyle=$
$\displaystyle|\zeta_{n+1}|^{2}-|\zeta_{n}|^{2}=(|\zeta_{n+1}|+|\zeta_{n}|)(|\zeta_{n+1}|-|\zeta_{n}|)$
(61) $\displaystyle\simeq$ $\displaystyle
2|\zeta_{n}|(|\zeta_{n+1}|-|\zeta_{n}|)$

to write

$\displaystyle\sum_{n=N}^{\infty}{\rm
Re}\left[\frac{e^{i\zeta_{n}\tau}}{C-\zeta_{n}^{2}}\right]\simeq-\sum_{n=N}^{\infty}\frac{1}{|\zeta_{n}|^{2}}\exp\left(-\frac{|\zeta_{n}|\tau}{\sqrt{2}}\right)\sin\left(\frac{|\zeta_{n}|\tau}{\sqrt{2}}\right)$
(62) $\displaystyle\simeq$
$\displaystyle-\frac{1}{\pi}\sum_{n=N}^{\infty}\frac{(|\zeta_{n+1}|-|\zeta_{n}|)}{|\zeta_{n}|}\exp\left(-\frac{|\zeta_{n}|\tau}{\sqrt{2}}\right)\sin\left(\frac{|\zeta_{n}|\tau}{\sqrt{2}}\right)$
$\displaystyle\simeq$
$\displaystyle-\frac{1}{\pi}\int_{|\zeta_{N}|}^{+\infty}\frac{d|\zeta|}{|\zeta|}\exp\left(-\frac{|\zeta|\tau}{\sqrt{2}}\right)\sin\left(\frac{|\zeta|\tau}{\sqrt{2}}\right)$
$\displaystyle=$
$\displaystyle-\frac{1}{\pi}\int_{|\zeta_{N}|\tau/\sqrt{2}}^{+\infty}dx\frac{e^{-x}\sin
x}{x}$ $\displaystyle=$
$\displaystyle-\frac{1}{4}+\frac{1}{\pi}\int_{0}^{|\zeta_{N}|\tau/\sqrt{2}}dx\frac{e^{-x}\sin
x}{x},$

where

$\int_{0}^{+\infty}dx\frac{e^{-x}\sin x}{x}=\frac{\pi}{4}$ (63)

is used.

Now, using Eqs. (37), (62), and

${\rm
Re}\biggl{[}\frac{e^{-i\omega_{\mu}t}}{1+\kappa^{2}-\omega_{\mu}^{2}/\omega_{p}^{2}}\biggr{]}=\frac{1}{2\kappa^{2}}{\rm
Re}\biggl{[}\frac{e^{-i\zeta_{\mu}\tau}}{(1+\kappa^{2})/2\kappa^{2}-\zeta_{\mu}^{2}}\biggr{]},$
(64)

we can express $E(k,t)$ as

$\displaystyle E(k,t)$ (65) $\displaystyle=$
$\displaystyle-i\;\alpha\frac{m\omega_{p}^{2}}{ek}(1+\kappa^{2})\sum_{{\rm
Re}\zeta_{\mu}>0}{\rm
Re}\biggl{[}\frac{e^{-i\zeta_{\mu}\tau}}{(1+\kappa^{2})/2\kappa^{2}-\zeta_{\mu}^{2}}\biggr{]}$
$\displaystyle\simeq$
$\displaystyle-i\;\alpha\frac{m\omega_{p}^{2}}{ek}(1+\kappa^{2})$
$\displaystyle\mbox{}\times\biggl{(}\sum_{0<{\rm Re}\zeta_{\mu}<{\rm
Re}\zeta_{N}}{\rm
Re}\biggl{[}\frac{e^{-i\zeta_{\mu}\tau}}{(1+\kappa^{2})/2\kappa^{2}-\zeta_{\mu}^{2}}\biggr{]}$
$\displaystyle-\frac{1}{\pi}\int_{|\zeta_{N}|\tau/\sqrt{2}}^{+\infty}dx\frac{e^{-x}\sin
x}{x}\biggr{)}$ $\displaystyle\equiv$ $\displaystyle E_{N}(k,t),$

where $\zeta_{n}\equiv\omega_{n}/kv_{T}$ and $\zeta_{n}\tau\equiv\omega_{n}t$.
Recall that we need to choose $N\gg 1$ such that
$|\omega_{N}|t/\sqrt{2}\gg\delta_{N}>0$. However, when $N\gg 1$,
$E(k,t=0)\simeq E_{N}(k,t=0)=i(4\pi
n_{0}e/k)\alpha=i(m/e)(\omega_{p}/k)\alpha$ also holds (as shown later) and
$E(k,t)\simeq E_{N}(k,t)$ for any $t\geq 0$.

From the residue theorem, we obtain

$\displaystyle\frac{1}{2\pi}\int_{-\infty}^{+\infty}d\zeta\frac{Z(\zeta)e^{-i\zeta\tau}}{1+\kappa^{-2}[1+\zeta
Z(\zeta)]}$ (66) $\displaystyle=$
$\displaystyle-i\sum_{\mu}\frac{\kappa^{2}(1+\kappa^{2})e^{-i\omega_{\mu}t}}{1+\kappa^{2}-\omega_{\mu}^{2}/\omega_{p}^{2}}$
$\displaystyle=$ $\displaystyle-2i\sum_{{\rm Re}\omega_{\mu}>0}{\rm
Re}\biggl{[}\frac{\kappa^{2}(1+\kappa^{2})e^{-i\omega_{\mu}t}}{1+\kappa^{2}-\omega_{\mu}^{2}/\omega_{p}^{2}}\biggr{]}$

where $\tau>0$. Using the asymptotic expansion of $Z(\zeta)$ for $|\zeta|\gg
1$ and taking the limit for $\tau\rightarrow+0$, we find

$\displaystyle\lim_{\tau\rightarrow+0}\frac{1}{2\pi}\int_{-\infty}^{+\infty}d\zeta\frac{Z(\zeta)e^{-i\zeta\tau}}{1+\kappa^{-2}[1+\zeta
Z(\zeta)]}$ (67) $\displaystyle=$
$\displaystyle-2i\lim_{\tau\rightarrow+0}\sum_{{\rm Re}\omega_{\mu}>0}{\rm
Re}\biggl{[}\frac{\kappa^{2}(1+\kappa^{2})e^{-i\omega_{\mu}t}}{1+\kappa^{2}-\omega_{\mu}^{2}/\omega_{p}^{2}}\biggr{]}=i.$

However, for $\tau=0$, we obtain

$\displaystyle\frac{1}{2\pi}\int_{C}d\zeta\frac{Z(\zeta)}{1+\kappa^{-2}[1+\zeta
Z(\zeta)]}$ (68) $\displaystyle=$ $\displaystyle-2i\sum_{{\rm
Re}\omega_{\mu}>0}{\rm
Re}\biggl{[}\frac{\kappa^{2}(1+\kappa^{2})}{1+\kappa^{2}-\omega_{\mu}^{2}/\omega_{p}^{2}}\biggr{]}=\frac{i}{4}(3-\kappa^{2}),$

which is derived using the residue theorem and noting that, for
$|\zeta|\rightarrow+\infty$,

$\displaystyle Z(\zeta)\sim-\frac{1}{\zeta}$ $\displaystyle\mbox{for}\hskip
8.53581pt-\frac{\pi}{4}<\arg\zeta<\frac{5\pi}{4},$
$\displaystyle|Z(\zeta)|>>1$ $\displaystyle\mbox{for}\hskip
8.53581pt-\frac{3\pi}{4}<\arg\zeta<-\frac{\pi}{4}.$

In Eq. (68), the integral $\int_{C}d\zeta$ is done along the circle
$C:\zeta=Re^{i\theta}$ in the complex $\zeta$-plane where the radius of $C$ is
given by $R\gg 1$ and the orientation of $C$ is taken clockwise by varying the
argument $\theta$ of $\zeta$ from $\pi$ to $-\pi$. Comparing Eqs. (67) and
(68) shows that changing the order of operations $\lim_{\tau\rightarrow+0}$
and $\sum_{\mu}$ results in different values. Thus, when any large but finite
fixed number of eigenfrequencies are included to evaluate $E(k,t)$, the limit
value for $\tau\rightarrow+0$ does not converge to the correct initial value
but to another value, the ratio of which to the correct one is estimated from
Eqs. (67) and (68) as $\frac{1}{4}(3-\kappa^{2})$. We use Eq. (68) to confirm
that the approximate solution $E_{N}(k,t)$ given in Eq. (65) can correctly
approach to the correct initial value as $t\rightarrow+0$,

$\displaystyle\lim_{N\rightarrow\infty}E_{N}(k,t=0)$ $\displaystyle=$
$\displaystyle-i\;\alpha\frac{m\omega_{p}^{2}}{ek}(1+\kappa^{2})$ (69)
$\displaystyle\mbox{}\times\biggl{(}\sum_{{\rm Re}\zeta_{\mu}>0}{\rm
Re}\biggl{[}\frac{1}{(1+\kappa^{2})/2\kappa^{2}-\zeta_{\mu}^{2}}\biggr{]}$
$\displaystyle-\frac{1}{\pi}\int_{0}^{+\infty}dx\frac{e^{-x}\sin
x}{x}\biggr{)}$ $\displaystyle=$ $\displaystyle
i\alpha\frac{m\omega_{p}^{2}}{ek}\biggl{[}\frac{1}{4}(3-\kappa^{2})+\frac{1}{4}(1+\kappa^{2})\biggr{]}$
$\displaystyle=$ $\displaystyle
i\alpha\frac{m\omega_{p}^{2}}{ek}=i\alpha\frac{4\pi n_{0}e}{k}$
$\displaystyle=$ $\displaystyle E(k,t=0)$

Now, we express the approximate solution $E_{N}(k,t)$ in Eq.(65) by the sum of
two parts as

$E_{N}(k,t)=E_{<N}(k,t)+E_{>N}(k,t),$ (70)

where

$\displaystyle E_{<N}(k,t)$ $\displaystyle\equiv$
$\displaystyle-i\;\alpha\frac{m\omega_{p}^{2}}{ek}(1+\kappa^{2})$
$\displaystyle\mbox{}\times\sum_{0<{\rm Re}\zeta_{\mu}<{\rm Re}\zeta_{N}}{\rm
Re}\biggl{[}\frac{e^{-i\zeta_{\mu}\tau}}{(1+\kappa^{2})/2\kappa^{2}-\zeta_{\mu}^{2}}\biggr{]},$
$\displaystyle E_{>N}(k,t)$ $\displaystyle\equiv$ $\displaystyle
i\;\alpha\frac{m\omega_{p}^{2}}{ek}(1+\kappa^{2})$ (71)
$\displaystyle\mbox{}\times\frac{1}{\pi}\int_{|\zeta_{N}|\tau/\sqrt{2}}^{+\infty}dx\frac{e^{-x}\sin
x}{x}.$

The top panel of Fig. 5 shows $E_{<N}(k,t)$, $E_{>N}(k,t)$, and
$E_{N}(k,t)=E_{<N}(k,t)+E_{>N}(k,t)$, for $N=100$. There, the approximate
solution $E_{N}(k,t)$ for $N=100$ and the CD simulation results for
$\alpha=0.01$ are plotted by blue solid and red dashed lines, respectively.
Since the effects of complex eigenfrequencies $\omega_{\mu}$ with large
absolute values included in $E_{>N}(k,t)$ become more significant as $t$
approaches 0, $E_{<N}(k,t)$ deviates from the CD simulation results for small
$t$ (specifically $t<0.2$ in Fig. 5). The diamonds plotted in the vertical
axes of the top and bottom panels in Fig. 5 correspond to the value given by
$\lim_{N\rightarrow\infty}E_{<N}(k,0)/E(k,0)=(3-\kappa^{2})/4=11/16=0.6875$
for $\kappa\equiv k\lambda_{D}=1/2$. However, the approximate solution agrees
with the CD simulation result by adding $E_{>N}(k,t)$. The bottom panel shows
$E_{<N}(k,t)$, $E_{>N}(k,t)$, and $E_{N}(k,t)=E_{<N}(k,t)+E_{>N}(k,t)$ for
$N=1$. Even in this case of $N=1$, the difference between $E_{N}(k,t)$ and the
CD simulation result near $t=0$ is as small as about 5%.

Figure 5:  Time evolutions of the normalized electric field $E(k,t)/E(k,0)$
for $k\lambda_{D}=1/2$. The top panel shows $E_{<N}(k,t)$ (green line),
$E_{>N}(k,t)$ (orange line), and $E_{N}(k,t)=E_{<N}(k,t)+E_{>N}(k,t)$ (blue
line) for $N=100$. The bottom panel shows $E_{<N}(k,t)$ (green line),
$E_{>N}(k,t)$ (orange line), and $E_{N}(k,t)=E_{<N}(k,t)+E_{>N}(k,t)$ (blue
line) for $N=1$. The red dotted lines in both panels represent the CD
simulation results for $\alpha=0.01$. The diamonds plotted in the vertical
axes of the top and bottom panels correspond to the value given by
$\lim_{N\rightarrow\infty}E_{<N}(k,0)/E(k,0)=(3-\kappa^{2})/4=11/16=0.6875$
for $\kappa\equiv k\lambda_{D}=1/2$.

As described in Appendix B, the Vlasov-Poisson system is symmetric under the
time reversal transformation whether the linear approximation is made or not.
Especially, under the initial condition made in this section, the distribution
function at $t=0$ is an even function of $v$, from which we find that the
distribution function is an even function of time $t$ as explained in Appendix
B. Accordingly, $f_{1}(k,-v,-t)=f_{1}(k,v,t)$ and $E(k,-t)=E(k,t)$ hold for
any $t$ and $v$. Therefore, $E(k,t)$ and $f_{1}(k,v,t)$ for $t<0$ are
immediately given from Eqs. (37) and (38) with making use of the time
reversal.

###  III.3 Approximate solution in early time

We here use the asymptotic expansion of $Z(\zeta)$ for $|\zeta|\gg 1$ given in
Eq. (169) to derive

$\displaystyle 1+\kappa^{-2}[1+\zeta Z(\zeta)]$ (72) $\displaystyle=$
$\displaystyle
1-\kappa^{-2}\sum_{n=1}^{N-1}\frac{\Gamma(n+1/2)}{\sqrt{\pi}}\frac{1}{\zeta^{2n}}$
$\displaystyle+$ $\displaystyle\kappa^{2}\frac{\Gamma(N+1/2)}{\sqrt{\pi}}\zeta
e^{-\zeta^{2}}\int_{i\infty}^{\zeta}\frac{e^{s^{2}}}{s^{2N}}ds.$

Then, for $|\zeta|\gg 1$ and $-\pi/4<\arg\zeta<5\pi/4$, we obtain

$\displaystyle\frac{Z(\zeta)}{1+\kappa^{-2}[1+\zeta Z(\zeta)]}$ (73)
$\displaystyle=$
$\displaystyle-\frac{1}{\zeta}\left[1+\sum_{n=1}^{N}e_{n}(\kappa^{2})\zeta^{-2n}+{\cal
O}(\zeta^{-2N-2})\right],$

where $e_{n}(\kappa^{2})$ $(n=1,2,\cdots)$ are recursively defined by

$e_{n}(\kappa^{2})\equiv\frac{\Gamma(n+1/2)}{\Gamma(1/2)}+\kappa^{-2}\sum_{j=0}^{n-1}\frac{\Gamma(n-j+1/2)}{\Gamma(1/2)}e_{j}(\kappa^{2}),$
(74)

with $e_{0}(\kappa^{2})\equiv 1$. Thus, we find

$\displaystyle e_{1}(\kappa^{2})$ $\displaystyle=$
$\displaystyle\frac{1}{2}+\frac{1}{2}\kappa^{-2},$ $\displaystyle
e_{2}(\kappa^{2})$ $\displaystyle\equiv$
$\displaystyle\frac{3}{4}+\kappa^{-2}\left(\frac{3}{4}+\frac{1}{2}e_{1}(\kappa^{2})\right)=\frac{3}{4}+\kappa^{-2}+\frac{1}{4}\kappa^{-4},$
$\displaystyle\cdots.$ (75)

We now use Eq. (73) to get

$\displaystyle\int_{C}\frac{d\zeta}{2\pi}\frac{Z(\zeta)e^{-i\zeta\tau}}{1+\kappa^{-2}[1+\zeta
Z(\zeta)]}$ (76) $\displaystyle=$
$\displaystyle-\int_{C}\frac{d\zeta}{2\pi\zeta}\biggl{[}\sum_{n=0}^{N}\frac{e_{n}(\kappa^{2})}{\zeta^{2n}}+{\cal
O}(\zeta^{-2N-2})\biggr{]}e^{-i\zeta\tau}$ $\displaystyle=$ $\displaystyle
i\sum_{n=0}^{N}\frac{(-1)^{n}}{(2n)!}e_{n}(\kappa^{2})\tau^{2n}+{\cal
O}(\tau^{2N+2}),$

where the integral $\int_{C}d\zeta$ is done along the circle
$C:\zeta=Re^{i\theta}$ in the complex $\zeta$-plane. Here, the radius of $C$
is given by $R\gg 1$ and the orientation of $C$ is taken clockwise by varying
the argument $\theta$ of $\zeta$ from $\pi$ to $-\pi$. Equation (37) is
rewritten here for $\tau\equiv kv_{T}t>0$ as

$\displaystyle E(k,t)$ $\displaystyle=$ $\displaystyle\alpha\frac{4\pi
en_{0}}{k}\int_{C}\frac{d\zeta}{2\pi}\;\frac{Z(\zeta)e^{-i\zeta\tau}}{1+\kappa^{-2}[1+\zeta
Z(\zeta)]},$ (77)

into which Eq. (76) is substituted to derive the Taylor expansion of $E(k,t)$
about $\tau=kv_{T}t=0$ is written as

$E(k,t)=i\alpha\frac{4\pi
en_{0}}{k}\sum_{n=0}^{N}\frac{(-1)^{n}}{(2n)!}e_{n}(\kappa^{2})\tau^{2n}+{\cal
O}(\tau^{2N+2}).$ (78)

Next, Eq. (38) for $f_{1}(k,v,t)$ is rewritten here as

$\displaystyle f_{1}(k,v,t)=$ $\displaystyle\alpha
f_{0}(v)\Biggl{[}\exp\left(-i\frac{v}{v_{T}}\tau\right)-i\kappa^{-2}\frac{v}{v_{T}}\int_{-\infty}^{+\infty}$
(79)
$\displaystyle\frac{d\zeta}{2\pi}\;\frac{Z(\zeta)e^{-i\zeta\tau}}{(\zeta-v/v_{T})\\{1+\kappa^{-2}[1+\zeta
Z(\zeta)]\\}}\Biggr{]}.$

Then, we use

$\frac{1}{\zeta-v/v_{T}}=\frac{1}{\zeta}\left[1+\sum_{n=1}^{\infty}\left(\frac{v}{v_{T}}\right)^{n}\zeta^{-n}\right],$
(80)

and Eq. (76) to obtain

$\displaystyle\frac{Z(\zeta)}{(\zeta-v/v_{T})\\{1+\kappa^{-2}[1+\zeta
Z(\zeta)]\\}}$ (81) $\displaystyle=$
$\displaystyle-\frac{1}{\zeta^{2}}\left[1+\sum_{n=1}^{N}d_{n}(\kappa^{2},v/v_{T})\zeta^{-n}+{\cal
O}(\zeta^{-N-1})\right],$

where $d_{n}(\kappa^{2},v/v_{T})$ $(n=0,1,2,\cdots)$ are defined by

$d_{n}(\kappa^{2},v/v_{T})\equiv\sum_{j=0}^{[n/2]}(v/v_{T})^{n-2j}e_{j}(\kappa^{2}),$
(82)

and $[n/2]$ denotes the greatest integer less than or equal to $n/2$. From Eq.
(82), we have

$\displaystyle d_{0}(\kappa^{2},v/v_{T})$ $\displaystyle=$ $\displaystyle 1$
$\displaystyle d_{1}(\kappa^{2},v/v_{T})$ $\displaystyle=$ $\displaystyle
v/v_{T}$ $\displaystyle d_{2}(\kappa^{2},v/v_{T})$ $\displaystyle=$
$\displaystyle(v/v_{T})^{2}+e_{1}(\kappa^{2})$ $\displaystyle
d_{3}(\kappa^{2},v/v_{T})$ $\displaystyle=$
$\displaystyle(v/v_{T})^{3}+(v/v_{T})\,e_{1}(\kappa^{2})$ $\displaystyle
d_{4}(\kappa^{2},v/v_{T})$ $\displaystyle=$
$\displaystyle(v/v_{T})^{4}+(v/v_{T})^{2}\,e_{1}(\kappa^{2})+e_{2}(\kappa^{2})$
$\displaystyle d_{5}(\kappa^{2},v/v_{T})$ $\displaystyle=$
$\displaystyle(v/v_{T})^{5}+(v/v_{T})^{3}\,e_{1}(\kappa^{2})+(v/v_{T})\,e_{2}(\kappa^{2})$
$\displaystyle\cdots.$ (83)

Integrating Eq. (81) with respect to $\zeta$ yields

$\displaystyle\int_{-\infty}^{+\infty}\frac{d\zeta}{2\pi}\;\frac{Z(\zeta)e^{-i\zeta\tau}}{(\zeta-v/v_{T})\\{1+\kappa^{-2}[1+\zeta
Z(\zeta)]\\}}$ (84) $\displaystyle=$
$\displaystyle-\int_{C}\frac{d\zeta}{2\pi}\;\frac{e^{-i\zeta\tau}}{\zeta^{2}}\biggl{[}1+\sum_{n=1}^{N}d_{n}(\kappa^{2},v/v_{T})\zeta^{-n}+{\cal
O}(\zeta^{-N-1})\biggr{]}$ $\displaystyle=$ $\displaystyle
i\sum_{n=0}^{N}\frac{(-i\tau)^{n+1}}{(n+1)!}d_{n}(\kappa^{2},v/v_{T})+{\cal
O}(\tau^{N+2}).$

Combining Eqs. (79) and (84), the Taylor expansion of $f_{1}(k,v,t)$ about
$\tau=kv_{T}t=0$ is derived as

$\displaystyle f_{1}(k,v,t)$ $\displaystyle=$ $\displaystyle\alpha
f_{0}(v)\biggl{[}\exp\left(-i\frac{v}{v_{T}}\tau\right)$ (85)
$\displaystyle+\kappa^{-2}\frac{v}{v_{T}}\sum_{n=1}^{N}\frac{(-i\tau)^{n}}{n!}d_{n-1}(\kappa^{2},v/v_{T})+{\cal
O}(\tau^{N+1})\biggr{]}$ $\displaystyle=$ $\displaystyle\alpha
f_{0}(v)\biggl{[}1+\sum_{n=1}^{N}\frac{(-i\tau)^{n}}{n!}\biggl{\\{}\left(\frac{v}{v_{T}}\right)^{n}$
$\displaystyle+\kappa^{-2}\frac{v}{v_{T}}d_{n-1}(\kappa^{2},v/v_{T})\biggr{\\}}+{\cal
O}(\tau^{N+1})\biggr{]}$ $\displaystyle=$ $\displaystyle\alpha
f_{0}(v)\biggl{[}1+\sum_{n=1}^{N}\frac{(-i\tau)^{n}}{n!}f_{n}(\kappa^{2},v/v_{T})$
$\displaystyle+{\cal O}(\tau^{N+1})\biggr{]},$

where $f_{n}(\kappa^{2},v/v_{T})$ $(n=1,2,\cdots)$ are defined by

$f_{n}(\kappa^{2},v/v_{T})\equiv\left(\frac{v}{v_{T}}\right)^{n}+\kappa^{-2}\frac{v}{v_{T}}d_{n-1}(\kappa^{2},v/v_{T}).$
(86)

We now compare the analytical solutions obtained by the series expansions in
Eqs. (78) and (85) with the CD simulation results for $\alpha=0.01$.

Figure 6: Time evolution of $E(k,t)/E(k,0)$ for $k\lambda_{D}=1/2$ The red
dashed line represents the CD simulation result for $\alpha=0.01$. Results
obtained from Eq. (78) including terms up to the orders of $\tau^{2}$,
$\tau^{6}$, and $\tau^{12}$ are shown by the blue, orange, and green lines,
respectively.

Figure 6 shows the analytical solutions of Eq. (78) for $k\lambda_{D}=1/2$
including terms up to orders of $\tau^{2},\tau^{6}$, and $\tau^{1}2$. We can
confirm that the discrepancy from the simulation results decreases as the
number of the included terms increases. The calculation up to $\tau^{6}$ fits
well with the CD simulation results when $\omega_{p}t<\sqrt{2}$ (which
corresponds to $\tau\equiv kv_{T}t\equiv\sqrt{2}(k\lambda_{D})\omega_{p}t<1$).
Figure 7 shows the plot of $f_{1}(x,v,t)/f_{0}(v)$ for $k\lambda_{D}=1/2$
calculated using Eq. (85) with terms up to the order of $\tau^{13}$ included.
It shows a good agreement with the CD simulation results shown in the bottom
row of Fig. 3 for $\omega_{p}t\leq 1$.

Figure 7: Contours of $f_{1}(x,v,t)/(\alpha^{2}f_{0}(v))$ in the $(x,v)$-space
for $k\lambda_{D}=1/2$ calculated using Eq. (85) with terms up to the order of
$\tau^{13}$ included. It shows a good agreement with the CD simulation results
shown in the bottom row of Fig. 3 for $\omega_{p}t\leq 1$.

##  IV  Analytical solution based on quasilinear theory and its verification
by numerical simulation

###  IV.1 Spatially averaged distribution function

In this subsection, formulas that hold rigorously for the spatially averaged
distribution function and its associated entropy are derived without using
linear approximations based on small perturbation amplitudes. The distribution
function $f(x,v,t)$ is a periodic function of $x$ and the period length is
given by $L\equiv 2\pi/k$. We use

$\langle\cdots\rangle\equiv\frac{1}{L}\int_{-L/2}^{L/2}\cdots\;dx,$ (87)

to denote the average with respect to $x$. The Vlasov equation is averaged in
$x$ to yield

$\frac{\partial}{\partial t}\langle f\rangle(v,t)+\frac{\partial}{\partial
v}\left[\langle f\rangle(v,t)\,a(v,t)\right]=0,$ (88)

where $\langle f\rangle(v,t)\equiv\langle f(x,v,t)\rangle$ is the electron
distribution function averaged in $x$, and $a(v,t)$ represents the average
acceleration of electrons defined by

$\langle f\rangle(v,t)\,a(v,t)\equiv-\frac{e}{m}\langle
E(x,t)\,f(x,v,t)\rangle.$ (89)

From Eq. (88), we have

$\displaystyle\frac{\partial}{\partial t}\left[\langle
f\rangle(v,t)\frac{1}{2}mv^{2}\right]+\frac{\partial}{\partial v}\left[\langle
f\rangle(v,t)\frac{1}{2}mv^{2}a(v,t)\right]$ $\displaystyle=\langle
f\rangle(v,t)\,m\,v\,a(v,t),$ (90)

and

$\frac{d}{dt}\int_{-\infty}^{+\infty}dv\;\langle
f\rangle(v,t)\frac{1}{2}mv^{2}=\int_{-\infty}^{+\infty}dv\;\langle
f\rangle(v,t)\,m\,v\,a(v,t),$ (91)

where $m\,v\,a(v,t)$ represents the $x$-averaged rate of change of the kinetic
energy of an electron with velocity $v$ at time $t$ caused by the electric
field. The total energy conservation of the Vlasov-Poisson system is written
as

$\frac{d}{dt}\left[\int_{-\infty}^{+\infty}dv\;\langle
f\rangle(v,t)\frac{1}{2}mv^{2}+\frac{1}{8\pi}\langle E^{2}\rangle\right]=0.$
(92)

The Gibbs entropy $S[f]$ per unit length in the $x$-direction is defined as a
functional of the distribution function $f(x,v,t)$ by

$S[f]\equiv-\left\langle\int_{-\infty}^{+\infty}dv\;f(x,v,t)\log
f(x,v,t)\right\rangle,$ (93)

which is found to be an invariant,

$\frac{d}{dt}S[f]=0.$ (94)

We now use the $x$-averaged distribution function $\langle f\rangle(v,t)$ to
define the entropy density $S[\langle f\rangle]$ by

$S[\langle f\rangle]\equiv-\int_{-\infty}^{+\infty}dv\;\langle
f\rangle(v,t)\log\langle f\rangle(v,t),$ (95)

which is not an invariant but a function of time $t$. From the viewpoint of
information theory, Information  the entropy density $S[\langle f\rangle]$ is
given from the average of

$S(v,t)\equiv-\log\langle f\rangle(v,t)=-\log[\langle
f\rangle(v,t)(dv/n_{0})]+\log(dv/n_{0}),$ (96)

which represents the Shannon information content (or self-entropy)
$-\log[\langle f\rangle(v,t)(dv/n_{0})]=-\log[P_{V}(v_{j})]$ [see Eq. (199)]
plus an additional constant given by $\log(dv/n_{0})$ for the probability
$\langle f\rangle(v,t)(dv/n_{0})$ of finding the electron velocity in the
interval $[v-\frac{1}{2}dv,v+\frac{1}{2}dv]$ where $dv$ is regarded as an
infinitesimal positive constant. Supplementary explanations on information
entropies in the Vlasov-Poisson system are presented in Appendix D. We now
consider $S(v,t)$ as the information entropy (except an additional constant)
of the electron with the velocity $v$ at time $t$. From Eq. (88), we obtain

$\displaystyle\left(\frac{\partial}{\partial t}+a(v,t)\frac{\partial}{\partial
v}\right)S(v,t)$ $\displaystyle=$
$\displaystyle-\left(\frac{\partial}{\partial
t}+a(v,t)\frac{\partial}{\partial v}\right)\log f(v,t)$ (97) $\displaystyle=$
$\displaystyle\frac{\partial a(v,t)}{\partial v}.$

Here, we define $u(v_{0},t)$ which satisfies the differential equation,

$\frac{\partial}{\partial t}u(v_{0},t)=a(u(v_{0},t),t)$ (98)

with the initial condition

$u(v_{0},t=0)=v_{0}.$ (99)

Then, $u(v_{0},t)$ represents the velocity of the electron which has the
initial velocity $v_{0}$ and the history of acceleration $a(v,t^{\prime})$
$(0\leq t^{\prime}\leq t)$. The interval
$[v_{0}-\frac{1}{2}dv_{0},v_{0}+\frac{1}{2}dv_{0}]$ in the $v$-space evolves
to the interval
$[u(v_{0},t)-\frac{1}{2}du(v_{0},t),u(v_{0},t)+\frac{1}{2}du(v_{0},t)]$ at
time $t$ when the velocities of the electrons in the interval change are given
by $u(v^{\prime},t^{\prime})$ ($v_{0}-\frac{1}{2}dv_{0}\leq v^{\prime}\leq
v_{0}+\frac{1}{2}dv_{0}$, $0\leq t^{\prime}\leq t)$ in the $v$-space. We note
that the number (or probability) of the electrons found in the interval
$[u(v_{0},t)-\frac{1}{2}du(v_{0},t),u(v_{0},t)+\frac{1}{2}du(v_{0},t)]$ is
invariant in time,

$\langle f\rangle(u(v_{0},t),t),t)\;du(v_{0},t)=\langle
f\rangle(v_{0},0)\;dv_{0}.$ (100)

Using Eqs. (88) and (97), we find

$\displaystyle\frac{\partial}{\partial t}S(u(v_{0},t),t)$ $\displaystyle=$
$\displaystyle-\biggl{[}\biggl{(}\frac{\partial}{\partial
t}+a(v,t)\frac{\partial}{\partial v}\biggr{)}\log\langle
f\rangle(v,t)\biggr{]}_{v=u(v_{0},t)}$ (101) $\displaystyle=$
$\displaystyle\biggl{[}\frac{\partial a(v,t)}{\partial
v}\biggr{]}_{v=u(v_{0},t)}$

which implies that $\partial a(v,t)/\partial v$ represents the rate of change
in the information entropy $S(v,t)\equiv-\log\langle f\rangle(v,t)$ of the
electron with the velocity $v$ at time $t$ along the trajectory $u(v_{0},t)$
in the $v$-space. Then, the increase of $S(v,t)$ along the trajectory during
the time interval $[0,t]$ is given by

$\displaystyle\Delta S(u(v_{0},t),t)\equiv S(u(v_{0},t),t)-S(v_{0},0)$
$\displaystyle=-\log\biggl{[}\frac{\langle f\rangle(u(v_{0},t),t)}{\langle
f\rangle(v_{0},0)}\biggr{]}=\int_{0}^{t}dt^{\prime}\biggl{[}\frac{\partial
a(v,t^{\prime})}{\partial v}\biggr{]}_{v=u(v_{0},t^{\prime})},$ (102)

from which we obtain

$\displaystyle\langle f\rangle(u(v_{0},t),t)=\langle
f\rangle(v_{0},0)\;\exp\bigl{[}-\Delta S(u(v_{0},t),t)\bigr{]}$
$\displaystyle=\langle
f\rangle(v_{0},0)\;\exp\biggl{(}-\int_{0}^{t}dt^{\prime}\biggl{[}\frac{\partial
a(v,t^{\prime})}{\partial v}\biggr{]}_{v=u(v_{0},t^{\prime})}\biggr{)}.$ (103)

We also find that the rate of change in the entropy of $S[\langle f\rangle]$
is given by

$\displaystyle\frac{d}{dt}S[\langle f\rangle]$ $\displaystyle=$
$\displaystyle\frac{d}{dt}\int_{-\infty}^{+\infty}dv\;\langle
f\rangle(v,t)\;S(v,t)$ (104) $\displaystyle=$
$\displaystyle\int_{-\infty}^{+\infty}dv\;\langle
f\rangle(v,t)\;\frac{\partial a(v,t)}{\partial v}.$

Another formula is obtained as

$\frac{\partial}{\partial t}\log\langle
f\rangle(u(v_{0},t),0)=\biggl{[}a(v,t)\frac{\partial}{\partial v}\log\langle
f\rangle(v,0)\biggr{]}_{v=u(v_{0},t)},$ (105)

which is integrated in time $t$ to derive

$\log\biggl{[}\frac{\langle f\rangle(u(v_{0},t),0)}{\langle
f\rangle(v_{0},0)}\biggr{]}=\int_{0}^{t}dt^{\prime}\biggl{[}a(v,t^{\prime})\frac{\partial}{\partial
v}\log\langle f\rangle(v,0)\biggr{]}_{v=u(v_{0},t^{\prime})}$ (106)

and

$\displaystyle\langle f\rangle(u(v_{0},t),0)$ $\displaystyle=\langle
f\rangle(v_{0},0)\exp\biggl{(}\int_{0}^{t}dt^{\prime}\biggl{[}a(v,t^{\prime})\frac{\partial}{\partial
v}\log\langle f\rangle(v,0)\biggr{]}_{v=u(v_{0},t^{\prime})}\biggr{)}.\hskip
22.76219pt$ (107)

Then, using Eqs. (IV.1) and (IV.1), we obtain

$\log\biggl{[}\frac{\langle f\rangle(u(v_{0},t),t)}{\langle
f\rangle(v_{0},0)}\biggr{]}=\Omega_{t}(v_{0}),$ (108)

and

$\langle f\rangle(u(v_{0},t),t)=\langle
f\rangle(u(v_{0},t),0)\;\exp\Omega_{t}(v_{0}),$ (109)

where the function $\Omega_{t}(v_{0})$ is defined by

$\displaystyle\Omega_{t}(v_{0})$ $\displaystyle\equiv$
$\displaystyle\int_{0}^{t}dt^{\prime}\;\Omega(u(v_{0},t^{\prime}),t^{\prime}),$
(110)

with

$\Omega(v,t)\equiv-a(v,t)\frac{\partial}{\partial v}\log\langle
f\rangle(v,0)-\frac{\partial a(v,t)}{\partial v}.$ (111)

Now, we consider the case in which the initial $x$-averaged distribution
function $\langle f\rangle(v,0)=f_{0}(v)$ is given by the Maxwellian
equilibrium distribution function in Eq. (25). Then, we have

$\frac{\partial}{\partial v}\log f_{0}(v)=-\frac{m\;v}{T},$ (112)

which is substituted into Eq. (IV.1) to obtain

$f_{0}(u(v_{0},t))=f_{0}(v_{0})\;\exp\biggl{(}-\frac{\Delta{\cal
E}(u(v_{0},t),t)}{T}\biggr{)},$ (113)

where

$\Delta{\cal E}(u(v_{0},t),t)\equiv
m\int_{0}^{t}dt^{\prime}\;u(v_{0},t^{\prime})\;a(u(v_{0},t^{\prime}),t^{\prime})$
(114)

represents the change in the kinetic energy of the electron with the initial
velocity $v_{0}$ caused by the acceleration due to the electric field during
the time interval $[0,t]$. Then, we can rewrite Eqs. (111), (110), and (109)
as

$\Omega(v,t)=\frac{m}{T}v\;a(v,t)-\frac{\partial a(v,t)}{\partial v},$ (115)
$\Omega_{t}(v_{0})=\frac{\Delta{\cal E}(u(v_{0},t),t)}{T}-\Delta
S(u(v_{0},t),t),$ (116)

and

$\displaystyle\langle
f\rangle(u(v_{0},t),t)=f_{0}(u(v_{0},t))\;\exp\Omega_{t}(v_{0})$
$\displaystyle=f_{0}(u(v_{0},t))\exp\biggl{[}\frac{\Delta{\cal
E}(u(v_{0},t),t)}{T}-\Delta S(u(v_{0},t),t)\biggr{]},\hskip 22.76219pt$ (117)

respectively.

Here, we should note that $\Omega_{t}(v_{0})$ defined in Eq. (116) corresponds
to the dissipation function employed by Evans and Searles to present their
fluctuation theorem. FT ; Evans  The fluctuation theorem by Evans and Searles
is based on the time-reversible Liouville equation and gives the formula for
the ratio of the probabilities that the dissipation function takes the
positive and negative values with the same absolute value. The Vlasov equation
differs from the Liouville equation treated by Evans and Searles in that it
contains the electron’s acceleration term determined from the distribution
function through Poisson’s equation. Therefore, the fluctuation theorem cannot
be directly applied to the function $\Omega_{t}(v_{0})$ in our case. However,
as explained in Appendix D, we can show the non-negativity of the expected
value of $\Omega_{t}(v_{0})$, which leads to the inequality in the form of the
second law of thermodynamics.

###  IV.2 Analysis to second order in perturbation amplitude

In this subsection, we investigate the time evolution of the $x$-averaged
distribution function considered in Sec. IV.A up to the second order in the
perturbation amplitude $\alpha$ based on the results of the linear analysis in
Sec. II. Recall that the initial condition for the distribution function is
given by

$\displaystyle f(x,v,t=0)$ $\displaystyle=$ $\displaystyle f_{0}(v)+{\rm
Re}[f_{1}(k,v,t=0)\exp(i\;k\;x)]$ (118) $\displaystyle=$ $\displaystyle
f_{0}(v)[1+\alpha\cos(k\,x)],$

where $f_{0}(v)$ is the Maxwellian in Eq. (25) and $f_{1}(k,v,t=0)=\alpha
f_{0}(v)$ is used. We have already derived the linear solution of the
perturbed distribution function $f_{1}(x,v,t)={\rm
Re}[f_{1}(k,v,t)\exp(i\;k\;x)]$ for $t>0$, from which the electric field
$E(x,t)={\rm Re}[E(k,t)\exp(i\;k\;x)]$ is obtained. Hereafter, we represent
the order of magnitude of the small perturbation amplitude by $\alpha$. Then,
the linear solutions $f_{1}(x,v,t)={\rm Re}[f_{1}(k,v,t)\exp(i\;k\;x)]={\cal
O}(\alpha)$ and $E(x,t)={\rm Re}[E(k,t)\exp(i\;k\;x)]={\cal O}(\alpha)$ are
used to derive

$\displaystyle\langle f\rangle(v,t)a(v,t)$ $\displaystyle=$
$\displaystyle-\frac{e}{m}\langle E(x,t)f(x,v,t)\rangle$ (119)
$\displaystyle=$ $\displaystyle-\frac{e}{2m}{\rm Re}[E^{*}(k,t)f_{1}(k,v,t)],$

where higher order terms than ${\cal O}(\alpha^{2})$ are neglected. Now, we
can write

$\langle f\rangle(v,t)=f_{0}(v)+f_{2}(v,t),$ (120)

where we keep terms only up to ${\cal O}(\alpha^{2})$ and $f_{2}(v,t)$
represents the ${\cal O}(\alpha^{2})$ part given by

$\displaystyle f_{2}(v,t)=\int_{0}^{t}dt^{\prime}\frac{\partial}{\partial
v}\left[\frac{e}{m}\langle E(x,t^{\prime})f(x,v,t^{\prime})\rangle\right]$
$\displaystyle=\frac{e}{2m}\int_{0}^{t}dt^{\prime}\;{\rm
Re}\left[E^{*}(k,t^{\prime})\frac{\partial f_{1}(k,v,t^{\prime})}{\partial
v}\right].$ (121)

We can also write

$\langle f\rangle(v,t)a(v,t)=f_{0}(v)a(v,t)+{\cal O}(\alpha^{4}),$ (122)

and neglect the ${\cal O}(\alpha^{4})$ part to obtain

$a(v,t)=-\frac{e}{2m}{\rm
Re}\left[E^{*}(k,t)\frac{f_{1}(k,v,t)}{f_{0}(v)}\right],$ (123)

and

$\displaystyle f_{2}(v,t)=-\int_{0}^{t}dt^{\prime}\frac{\partial}{\partial
v}[f_{0}(v)a(v,t)]$
$\displaystyle=f_{0}(v)\int_{0}^{t}dt^{\prime}\left[\frac{mv}{T}a(v,t^{\prime})-\frac{\partial
a(v,t^{\prime})}{\partial v}\right]$
$\displaystyle=f_{0}(v)\left[\frac{\Delta{\cal E}(v,t)}{T}-\Delta
S(v,t)\right],$ (124)

where

$\displaystyle\Delta{\cal
E}(v,t)\equiv\int_{0}^{t}dt^{\prime}\;m\,v\,a(v,t^{\prime})$
$\displaystyle=-\frac{e}{2}v\int_{0}^{t}dt^{\prime}\;{\rm
Re}\left[E^{*}(k,t^{\prime})\frac{f_{1}(k,v,t^{\prime})}{f_{0}(v)}\right],$
(125)

and

$\displaystyle\Delta S(v,t)\equiv
S(v,t)-S(v,0)=\int_{0}^{t}dt^{\prime}\;\frac{\partial
a(v,t^{\prime})}{\partial v}$
$\displaystyle=-\frac{e}{2m}\int_{0}^{t}dt^{\prime}\;{\rm
Re}\left[E^{*}(k,t^{\prime})\frac{\partial}{\partial
v}\left(\frac{f_{1}(k,v,t^{\prime})}{f_{0}(v)}\right)\right].$ (126)

Since $a(v,t)={\cal O}(\alpha^{2})$, the variation in the electron velocity
during the time interval $[0,t]$ is of ${\cal O}(\alpha^{2})$. Neglecting
${\cal O}(\alpha^{4})$ effects, we can regard $\Delta{\cal E}(v,t)$ as the
change in the kinetic energy of the electron with the initial velocity $v$
caused by the acceleration due to the electric field during the time interval
$[0,t]$. Recall that $\partial a(v,t)/\partial v$ represents the rate of
change in the information content $S(v,t)\equiv-\log\langle f\rangle(v,t)$
associated with the distribution of electrons in the velocity space. Then,
$\Delta S(v,t)$ represents the change in the information content (or the
information entropy) of the electron with the initial velocity $v$) during the
time interval $[0,t]$.

From Eq. (IV.2), we find

$\int_{-\infty}^{+\infty}dv\;f_{2}(v,t)=\int_{-\infty}^{+\infty}dv\;f_{0}(v)\left[\frac{\Delta{\cal
E}(v,t)}{T}-\Delta S(v,t)\right]=0$ (127)

and

$\displaystyle\frac{\Delta{\cal
E}(t)}{T}\equiv\frac{1}{T}\int_{-\infty}^{+\infty}dv\;f_{0}(v)\;\Delta{\cal
E}(v,t)$ $\displaystyle=\int_{-\infty}^{+\infty}dv\;f_{0}(v)\;\Delta
S(v,t)\equiv\Delta S[\langle f\rangle](t)$ (128)

which implies that, to the second order in $\alpha$, the change in the
information (or Gibbs) entropy $\Delta S[\langle f\rangle](t)\equiv S[\langle
f\rangle](t)-S[\langle f\rangle](0)$ is equal to the inverse temperature $1/T$
multiplied by the energy transfer $\Delta{\cal E}(t)$ from the electric field
energy to the kinetic energy during the time interval $[0,t]$.

###  IV.3 Calculation of $\Delta{\cal E}(v,t)$, $\Delta S(v,t)$, and
$f_{2}(v,t)$

Neglecting ${\cal O}(\alpha^{4})$ effects, the change $\Delta{\cal E}(v,t)$ in
the kinetic energy and the change $\Delta S(v,t)$ in the information content
(or the information entropy) of the electron with the initial velocity $v$
during the time interval $[0,t]$ are given by Eqs. (IV.2) and (IV.2),
respectively. The electric field $E(k,t)$ and the perturbed distribution
function $f_{1}(k,v,t)$ appearing in Eqs. (IV.2) and (IV.2) are evaluated
using Eqs. (37) and (38), respectively, in which the series summations
involving an infinite number of the complex-valued eigenfrequencies
$\\{\omega_{\mu}\\}$ converge more quickly for larger time $t$. On the other
hand, the Taylor expansions of $E(k,t)$ and $f_{1}(k,v,t)$ about
$\tau=kv_{T}t=0$ are given by Eqs. (78) and (85), respectively, which converge
more quickly for smaller time $t$. Therefore, to evaluate $\Delta{\cal
E}(v,t)$ and $\Delta S(v,t)$, the expressions of $E(k,t)$ and $f_{1}(k,v,t)$
in Eqs. (78) and (85) and those in Eqs. (37) and (38) are separately used for
smaller and larger time, respectively, and truncation to finite terms is made
in these expressions so that the time integrals can be performed not
numerically but analytically because the time dependence of the integrands are
given in the form of the summation of exponential functions multiplied by
polynomials of time. Under the conditions of examples shown later, good
convergence is verified when including terms up to ${\cal O}(\tau^{13})$ in
the integrands in Eqs. (IV.2) and (IV.2) for $\tau\equiv kv_{T}t\leq 1$ and
using two pairs of the complex frequencies for Eqs. (37) and (38) to obtain
$E(k,t)$ and $f_{1}(k,v,t)$ for $\tau\equiv kv_{T}t\geq 1$. Once that
$\Delta{\cal E}(v,t)$ and $\Delta s_{V}(v,t)$ are obtained following the
procedures described above, we can use Eq. (IV.2) to calculate $f_{2}(v,t)$ by

$f_{2}(v,t)=f_{0}(v)\;\Omega_{t}(v),$ (129)

where

$\Omega_{t}(v)\equiv\frac{\Delta{\cal E}(v,t)}{T}-\Delta S(v,t).$ (130)

Figure 8 shows contours of $\Delta{\cal E}(v,t)/T$ and $\Delta S(v,t)$ on the
$(v,t)$ plane, obtained from the analytical solution using the aforementioned
procedure for $k\lambda_{D}=1/2$. In most of the $(v,t)$ plane, we find
$\Delta{\cal E}(v,t)>0$ which indicates the increase of the electron kinetic
energy. Here, the resonance velocity is given by $v_{res}\equiv{\rm
Re}(\omega_{\mbox{\o}})/k=2.831\;v_{t}$. As time progresses, $\Delta{\cal
E}(v,t)$ becomes more concentrated around the resonance velocities $v=\pm
v_{res}$ while $\Delta S(v,t)$ clearly shows the positive maximum (negative
minimum) at a slightly smaller (larger) absolute value $|v|$ than $v_{res}$.

Figure 8:  Contours of $\Delta{\cal E}(v,t)/(\alpha^{2}T)$ (top) and $\Delta
S(v,t)/\alpha^{2}$ (bottom) on the $(v,t)$ plane for $k\lambda_{D}=1/2$.

The top and bottom panels of Fig. 9 respectively show the distributions of
$f_{2}(v,t)/(\alpha^{2}n_{0}v_{t}^{-3})$ and $f_{2}(v,t)/(\alpha^{2}f_{0}(v))$
in the $(v,t)$ plane, calculated from the difference between $\Delta{\cal
E}(v,t)/T$ and $\Delta S(v,t)$ using Eqs. (129) and (130). As time progresses,
the distribution of $f_{2}(v,t)$ spreads from around $v=0$ to $v_{res}$ while
such spreading is not clearly seen in the distribution of
$f_{2}(v,t)/f_{0}(v)$.

Figure 9:  Contours of $f_{2}(v,t)/(\alpha^{2}n_{0}v_{t}^{-3})$ (top) and
$f_{2}(v,t)/(\alpha^{2}f_{0}(v))$ (bottom) on the $(v,t)$ plane for
$k\lambda_{D}=1/2$.

Figure 10 shows $f_{2}(v,t)/(\alpha^{2}n_{0}v_{t}^{-3})$ and
$f_{2}(v,t)/(\alpha^{2}f_{0}(v,t))$ as functions of $v$, obtained from the
analytical solution for $\omega_{p}t=0.1,0.5,1,5$, and 10. We see that
$f_{2}(v,t)$ oscillate along the $v$ direction and the number of oscillations
increases with increasing time. Positive peaks and negative troughs of
$f_{2}(v,t)/f_{0}(v)$ can be seen around the resonant velocities $v=\pm
v_{res}$. The red dots represent the results of CD simulations for
$\alpha=0.1$, which agree well with the analytical solution. However, for
$\omega_{p}t\geq 5$, a significant discrepancy between the CD simulation
results and the analytical solution appears near $v=0$. This discrepancy is
attributed to the fact that the distribution function treated in the CD
simulations is flat between a finite number of contour lines This results in
an underestimation of the amplitude of the perturbed distribution function
$f_{1}(k,v,t)$ driven by $\partial f_{0}(v)/\partial v$ which is set to zero
in the neighborhood of $v=0$ for the CD simulations. Then, the absolute value
of $f_{2}(v,t)$ around $v=0$ is also underestimated. Indeed, it is confirmed
that increasing the number of contour lines and narrowing the intervals
between them in CD simulations reduces the difference between the CD
simulation results and the analytical solution near $v=0$.

Profiles of $f_{2}(v,t)/(\alpha^{2}n_{0}v_{t}^{-3})$ and
$f_{2}(v,t)/(\alpha^{2}f_{0}(v))$ obtained by the analytical formulas for
$\omega_{p}t=10$, 20, 50, and 100 are shown by blue curves in the top and
bottom panels of Fig. 11, respectively. The red curves in Fig. 11 represent
the profiles in the long-time limit. As $t\rightarrow+\infty$, $f_{2}(v,t)$
and $f_{2}(v,t)/f_{0}(v)$ converge to the structures which have the positive
maximums (negative minimums) at $|v|$ slightly larger (smaller) than
$v_{res}=2.831\;v_{t}$. This is consistent with the well-known picture of the
increase and decrease in particles’ number around the resonant velocity in the
Landau damping process.

•

Figure 10:  Profiles of $f_{2}(v,t)/(\alpha^{2}n_{0}v_{t}^{-3})$ (top) and
$f_{2}(v,t)/(\alpha^{2}f_{0}(v))$ (bottom) for $\omega_{p}t=0.1$, 0.5, 1, 5,
and 10. Results obtained by the analytical formulas and the CD simulation for
$\alpha=0.1$ are shown by blue curves and red dots.  Figure 11:  Profiles of
$f_{2}(v,t)/(\alpha^{2}n_{0}v_{t}^{-3})$ (top) and
$f_{2}(v,t)/(\alpha^{2}f_{0}(v))$ (bottom) obtained by the analytical formulas
for $\omega_{p}t=10$, 20, 50, and 100 are shown by blue curves. The red curves
represent the profile in the limit of $t\rightarrow+\infty$.

For an arbitrary function $A(v)$, the change in its average value during
$[0,t]$ is evaluated using Eqs.(129) and (130) as

$\Delta\overline{A}(t)=\frac{1}{n_{0}}\int_{-\infty}^{+\infty}dv\;f_{2}(v,t)A(v)=\frac{1}{n_{0}}\int_{-\infty}^{+\infty}dv\;f_{0}(v)\;\Omega_{t}(v)A(v).$
(131)

By setting $A(v)=mv^{2}/2$, it can be shown that the above equation is equal
to the increase in kinetic energy $\Delta{\cal E}(t)$ over the time interval
$[0,t]$ given by Eq. (IV.2). Furthermore, by calculating Eq. (131) for
$A(v)=v^{2}$ and $v^{4}$, the kurtosis of the velocity as a random variable at
time $t$ is given by

$\displaystyle K(t)$ $\displaystyle=$
$\displaystyle\frac{\overline{v^{4}}(t)}{[\overline{v^{2}}(t)]^{2}}=\frac{\overline{v^{4}}(0)+\Delta\overline{v^{4}}(t)}{[\overline{v^{2}}(0)+\Delta\overline{v^{2}}(t)]^{2}}$
(132) $\displaystyle=$
$\displaystyle\frac{3+\Delta\overline{v^{4}}(t)/[\overline{v^{2}}(0)]^{2}}{[1+\Delta\overline{v^{2}}(t)/\overline{v^{2}}(0)]^{2}}=3+\Delta
K(t),$

where $\Delta K(t)$ is given up to the second order in $\alpha$ by

$\Delta
K(t)=\frac{\Delta\overline{v^{4}}(t)}{[\overline{v^{2}}(0)]^{2}}-6\frac{\Delta\overline{v^{2}}(t)}{\overline{v^{2}}(0)}=\frac{m^{2}\Delta\overline{v^{4}}(t)}{T_{2}}-6\frac{m\;\Delta\overline{v^{2}}(t)}{T}.$
(133)

Both $\Delta{\cal E}(t)$ and $\Delta K(t)$ are proportional to $\alpha^{2}$.
Figure 12 shows time evolutions of $\Delta{\cal E}(t)/(\alpha^{2}T)$ and
$\Delta K(t)/\alpha^{2}$ for $k\lambda_{D}=1/2$. It is seen that $\Delta{\cal
E}(t)/(\alpha^{2}T)$ converges to $\Delta{\cal
E}(\infty)/(\alpha^{2}T)=1/(4k^{2}\lambda_{D}^{2})=1$ as
$t\rightarrow+\infty$. The electric field is a standing wave of the form $\sin
kx$, and its amplitude becomes zero twice during one period of plasma
oscillation $2\pi/{\rm Re}(\omega_{\mbox{\o}})=4.438/\omega_{p}$. Considering
that the sum of the kinetic energy of electrons and the energy of the electric
field remains constant, it can be observed from the figure that $\Delta{\cal
E}(t)/(\alpha^{2}T)$ oscillates, reaching its maximum value, $\Delta{\cal
E}(\infty)/(\alpha^{2}T)=1/(4k^{2}\lambda_{D}^{2})=1$, twice during the period
$2\pi/{\rm Re}(\omega_{\mbox{\o}})$, and converges as $t\rightarrow+\infty$.
As seen from Figs. 10 and 11, $\Delta K(t)$ is negative at early times like
$\omega_{p}t<2$, which reflects the fact that $f_{2}(v,t)$ is more localized
near $v=0$ for small $t$ as seen in Figs. 9 and 10. As $t$ increases, peaks of
$f_{2}(v)$ near the resonant velocities $v=\pm v_{res}$ become more prominent,
and $\Delta K(t)$ takes positive values and increases, which indicates a
relative increase in the number of electrons with larger $|v|$ compared to a
Gaussian distribution. We find that $\Delta K(t)$ converges to $\Delta
K(\infty)=23.92\alpha^{2}$ as $t\rightarrow+\infty$.

Figure 12:  Time evolutions of $\Delta{\cal E}(t)/(\alpha^{2}T)$ (top) and
$\Delta K(t)/\alpha^{2}$ (bottom).

##  V Time evolution of information entropies

In this section, we consider the time evolution of information entropies in
the Vlasov-Poisson system. For that purpose, we use the solution of the
initial value problem in Sec. II that gives the perturbed density and the
electric field of ${\cal O}(\alpha)$ as

$n_{1}(x,t)=n_{1}(t)\cos(k\;x),\hskip 14.22636ptE(x,t)=E_{1}(t)\sin(k\;x),$
(134)

where

$E_{1}(t)=-\frac{4\pi e}{k}n_{1}(t).$ (135)

Here, $n_{1}(t)$ and $E_{1}(t)$ satisfy the initial conditions,

$n_{1}(t=0)=\alpha\;n_{0},\hskip 14.22636ptE_{1}(t=0)=-\alpha\frac{4\pi
en_{0}}{k},$ (136)

and they approach zero as $t\rightarrow+\infty$,

$n_{1}(x,\infty)=0,\hskip 14.22636ptE(x,\infty)=0.$ (137)

The electric field energy density at time $t$ is given by

$\displaystyle\frac{1}{8\pi}\langle
E(x,t)^{2}\rangle=\frac{1}{8\pi}\left(\frac{4\pi
e}{k}\right)^{2}\langle[n_{1}(x,t)]^{2}\rangle$
$\displaystyle=\frac{1}{16\pi}[E_{1}(t)]^{2}=\frac{1}{16\pi}\left(\frac{4\pi
e}{k}\right)^{2}[n_{1}(t)]^{2}$
$\displaystyle=\frac{n_{0}T}{2(k\lambda_{D})^{2}}\frac{\langle
n_{1}(x,t)^{2}\rangle}{n_{0}^{2}}=\frac{n_{0}T}{4(k\lambda_{D})^{2}}\frac{[n_{1}(t)]^{2}}{n_{0}^{2}}.$
(138)

Also, from the energy conservation law, we obtain

$\displaystyle n_{0}\Delta{\cal E}(t)+\frac{1}{8\pi}\langle
E(x,t)^{2}\rangle=n_{0}\Delta{\cal E}(\infty)$
$\displaystyle=\frac{1}{8\pi}\langle[E(x,0)]^{2}\rangle=\alpha^{2}\frac{n_{0}T}{4(k\lambda_{D})^{2}},$
(139)

where $\Delta{\cal E}(t)$ represents the change in the kinetic energy of the
electron during the time interval $[0,t]$ given by Eq. (IV.2).

We regard the electron’s position and velocity as random variables which are
represented by $X$ and $V$, respectively, as explained in Appendix D. The
joint probability density function of $X$ and $V$ are denoted by $p(x,v,t)$
[see Eq. (181)], which is integrated with respect to $x$ and $y$ to give the
marginal probability density functions $p_{X}(x,t)$ and $p_{V}(v,t)$,
respectively [see Eqs. (182) and (183)]. The entropy $S_{p}(X)\equiv S[p_{X}]$
is derived from the electron’s position distribution function $p_{X}(x,t)$ as

$\displaystyle S_{p}(X)\equiv
S[p_{X}]\equiv-\int_{-L/2}^{+L/2}dx\;p_{X}(x,t)\log p_{X}(x,t)$
$\displaystyle=-\int_{-L/2}^{+L/2}\frac{dx}{L}\;\left(\frac{n(x,t)}{n_{0}}\right)\log\left(\frac{n(x,t)}{n_{0}}\right)+\log(L)$
$\displaystyle\simeq\log(L)-\frac{1}{2}\int_{-L/2}^{+L/2}\frac{dx}{L}\left(\frac{n_{1}(x,t)}{n_{0}}\right)^{2}.$
(140)

In the last line of Eq. (V), terms of higher orders than $\alpha^{2}$ are
neglected. The increase in $S_{p}(X)\equiv S[p_{X}]$ during the time interval
from $0$ to $t$ is given by

$\displaystyle\Delta S[p_{X}](t)=S[p_{X}](t)-S[p_{X}](0)$
$\displaystyle\simeq\frac{1}{2}\int_{-L/2}^{+L/2}\frac{dx}{L}\left[\left(\frac{n_{1}(x,0)}{n_{0}}\right)^{2}-\left(\frac{n_{1}(x,t)}{n_{0}}\right)^{2}\right]$
$\displaystyle=\frac{\langle[n_{1}(x,0)]^{2}-[n_{1}(x,t)]^{2}\rangle}{2n_{0}^{2}}.$
(141)

As $t\rightarrow\infty$, it approaches to

$\Delta
S[p_{X}](\infty)=\frac{\langle[n_{1}(x,0)]^{2}\rangle}{2n_{0}^{2}}=\frac{\alpha^{2}}{4}.$
(142)

The velocity distribution function $p_{V}(v,t)$ is used to express the entropy
$S_{p}(V)\equiv S[p_{V}]$ as

$\displaystyle S_{p}(V)\equiv
S[p_{V}]\equiv-\int_{-\infty}^{+\infty}dv\;p_{V}(v,t)\log p_{X}(v,t)$
$\displaystyle=-\frac{1}{n_{0}}\int_{-\infty}^{+\infty}dv\;\langle
f\rangle(v,t)\log\langle f\rangle(v,t)+\log(n_{0}).$ (143)

The increase in $S_{p}(V)\equiv S[p_{V}]$ during the time interval from $0$ to
$t$ is represented by

$\displaystyle\Delta S[p_{V}](t)=S[p_{V}](t)-S[p_{V}](0)$
$\displaystyle\simeq\frac{1}{n_{0}}\int_{-\infty}^{+\infty}dv\;f_{2}(v,t)\frac{mv^{2}}{2T}=\frac{\Delta{\cal
E}(t)}{T},$ (144)

where terms of higher orders than $\alpha^{2}$ are neglected again. Using Eq.
Eqs. (V) and (V), $\Delta S[p_{V}](t)$ can be rewritten as

$\displaystyle\Delta S[p_{V}](t)$ $\displaystyle=$
$\displaystyle\frac{\Delta{\cal E}(t)}{T}=\frac{1}{8\pi
n_{0}T}\left\langle[E(x,0)]^{2}-[E(x,t)]^{2}\right\rangle$ (145)
$\displaystyle=$
$\displaystyle\frac{1}{2(k\lambda_{D})^{2}}\frac{\left\langle[n_{1}(x,0)]^{2}-[n_{1}(x,t)]^{2}\right\rangle}{n_{0}^{2}}.$

Its value in the limit of $t\rightarrow\infty$ is given by

$\Delta S[p_{V}](\infty)=\frac{\Delta{\cal
E}(\infty)}{T}=\frac{\alpha^{2}}{4k^{2}\lambda_{D}^{2}}.$ (146)

The mutual information of the random variables $X$ and $V$ is defined by

$I(X,V)=S_{p}(X)+S_{p}(V)-S_{p}(X,V).$ (147)

Under the initial condition given by Eq. (118), $X$ and $V$ are statistically
independent at $t=0$, and accordingly

$[I(X,V)]_{t=0}=0.$ (148)

When $f(x,v,t)=(n_{0}L)p(x,v,t)$ satisfies the Vlasov equation in Eq. (4), the
entropy defined by

$S_{P}(X,V)\equiv
S[p]\equiv-\int_{-L/2}^{+L/2}dx\int_{-\infty}^{+\infty}dv\;p(x,v)\log p(x,v)$
(149)

is known to be an invariant, and the increase in $I(X,V)$ during the time
interval from $0$ to $t$ is written as

$\displaystyle\Delta I(X,V)=I(X,V)=\Delta S_{p}(X)+\Delta S_{p}(V)$
$\displaystyle=\frac{\langle[n_{1}(x,0)]^{2}-[n_{1}(x,t)]^{2}\rangle}{2n_{0}^{2}}+\frac{\Delta{\cal
E}(t)}{T}$
$\displaystyle=\frac{\langle[n_{1}(x,0)]^{2}-[n_{1}(x,t)]^{2}}{2n_{0}^{2}}\left(1+\frac{1}{k^{2}\lambda_{D}^{2}}\right)$
$\displaystyle=\frac{\Delta{\cal E}(t)}{T}(1+k^{2}\lambda_{D}^{2}).$ (150)

In the limit of $t\rightarrow\infty$, we obtain

$[\Delta
I(X,V)]_{t=\infty}=\frac{\alpha^{2}}{4}\left(1+\frac{1}{k^{2}\lambda_{D}^{2}}\right).$
(151)

The time evolutions of $\Delta S_{p}(X)\equiv\Delta
S[p_{X}]=(k\lambda_{D})^{2}\Delta{\cal E}(t)/T$, $\Delta S_{p}(V)\equiv\Delta
S[p_{V}]=\Delta{\cal E}(t)/T$, and $\Delta I(X,V)=\Delta S_{p}(X)+\Delta
S_{p}(V)$ are shown in Figure 13, where contributions of higher orders than
${\cal O}(\alpha^{2})$ are neglected. As functions of time $t$, $\Delta
S_{p}(X)\equiv\Delta S[p_{X}]$, $\Delta S_{p}(V)\equiv\Delta S[p_{V}]$ and
$\Delta I(X,V)$ have the same form but different magnitudes.

Figure 13:  Time evolutions of $\Delta S_{p}(X)=(k\lambda_{D})^{2}\Delta{\cal
E}(t)/(\alpha^{2}T)$, $\Delta S_{p}(V)/\alpha^{2}=\Delta{\cal
E}(t)/(\alpha^{2}T)$, and $\Delta I(X,V)/\alpha^{2}=(\Delta S_{p}(X)+\Delta
S_{p}(V))/\alpha^{2}$.

Comparing the velocity distribution functions $p_{V}(v,t)=\langle
f\rangle(v,t)/n_{0}$ and $p_{V}(v,0)=\langle f\rangle(v,0)/n_{0}$, and using
Eqs. (211) and (D), we obtain the relative entropy $S(p_{V},t||p_{V},0)$ as

$\displaystyle S(p_{V},t||p_{V},0)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{+\infty}dv\;p_{V}(v,t)\log\left[\frac{p_{V}(v,t)}{p_{V}(v,0)}\right]$
(152) $\displaystyle\simeq$
$\displaystyle\frac{1}{2n_{0}}\int_{-\infty}^{+\infty}dv\;\frac{[f_{2}(v,t)]^{2}}{f_{0}(v)}$
$\displaystyle=$
$\displaystyle\frac{1}{2n_{0}}\int_{-\infty}^{+\infty}dv\;f_{0}(v)[\Omega_{t}(v)]^{2},$

which takes a small value of ${\cal O}(\alpha^{4})$. Figure 14 shows the time
evolution of $S(p_{V},t||p_{V},0)$.

Figure 14:  Time evolution of the relative entropy $S(p_{V},t||p_{V},0)$
normalized by $\alpha^{4}$.

As $t$ increases, $S(p_{V},t||p_{V},0)$ increases while showing modulation and
approaches to the limit value that corresponds to the limit function
$\lim_{t\rightarrow+\infty}f_{2}(v,t)$ in Fig. 11. As described in Appendix D,
the non-negativity of $S(p_{V},t||p_{V},0)$ implies that, even in the
collisionless Vlasov-Poisson system, the inequality in the form of the second
law of thermodynamics holds in the relation between the heat transfer from the
Maxwellian velocity distribution to the electric field and the conditional
entropy of the electron position variable for a given velocity distribution.

##  VI Collisional relaxation to thermal equilibrium

Coulomb collisions are neglected in the Vlasov equation (1). However, even if
the collision frequency is very small but finite, Coulomb collisions
eventually relax the distribution function to the Maxwellian

$f_{M\infty}=n_{M\infty}\left(\frac{m}{2\pi
T_{M\infty}}\right)^{1/2}\exp\left[-\frac{m}{2T_{M\infty}}(v-u_{M\infty})^{2}\right]$
(153)

at $t=+\infty$ while conserving total particles’ number, momentum, and energy.
The Maxwellian distribution function $f_{M\infty}(v)$ is the equilibrium
solution of the Boltzmann equation that is given by including the collision
term in the Vlasov equation. When the collision operator acts on the
Maxwellian, the collision term vanishes. Then, substituting $f_{M\infty}(v)$
into the left-hand side into the Vlasov equation, it also vanishes. It can be
proven from the above-mentioned fact that $n_{M\infty}$, $T_{M\infty}$, and
$u_{M\infty}$ are all independent of $t$ and $x$. Under the initial condition
in Eq. (118), we can use conservation laws for particles’ number, energy, and
momentum to derive

$n_{M\infty}=n_{0},\hskip 14.22636ptu_{M\infty}=0,\hskip
14.22636pt\frac{n_{0}T_{M\infty}}{2}=\frac{n_{0}T_{0}}{2}+\frac{\langle[E(x,0)]^{2}\rangle}{8\pi}.$
(154)

Thus, the distribution function in the thermal equilibrium state reached by
the collisional relaxation is given by

$f_{M\infty}=n_{0}\left(\frac{m}{2\pi
T_{M\infty}}\right)^{1/2}\exp\left(-\frac{mv^{2}}{2T_{M\infty}}\right).$ (155)

Here, ions are treated as uniform background positive charge with infinite
mass, so no energy exchange between electrons and ions due to collisions is
considered to occur. From Eq. (154), we also have

$\frac{\Delta
T_{M\infty}}{T}=\frac{T_{M\infty}-T}{T}=\frac{\langle[E(x,0)]^{2}\rangle}{4\pi
n_{0}T}=\frac{\alpha^{2}}{2(k\lambda)^{2}}=2\frac{\Delta{\cal E}(\infty)}{T}.$
(156)

From the equilibrium probability distribution function
$p_{M}(x,v)=f_{M\infty}(v)/(n_{0}L)$ and its marginal probability distribution
functions $p_{MX}(x)=1/L$ and $p_{MV}(v)=f_{M\infty}(v)/n_{0}$, we can define
the entropies $S[p_{M}]$, $S[p_{MX}]$, and $S[p_{MV}]$ in the thermal
equilibrium where

$p_{M}=p_{MX}\cdot p_{MV}$ (157)

holds so that the random variables $X$ and $V$ are statistically independent,
and accordingly, the mutual information of $X$ and $V$ vanishes,

$I_{M}(X,V)\equiv S[p_{M}]-S[p_{MX}]-S[p_{MV}]=0.$ (158)

The deviation of the entropy $S[p_{MX}]$ from the initial value
$S[p_{X}](t=0)$ is given by

$\displaystyle\Delta S[p_{MX}]\equiv S[p_{MX}]-S[p_{X}](t=0)$
$\displaystyle=-\left\langle[1+\alpha\cos(kx)]\log[1+\alpha\cos(kx)]\right\rangle\simeq\frac{\alpha^{2}}{4}.$
(159)

We see from Eqs. (142) and (VI) that $\Delta S[p_{MX}]$ and $\Delta
S[p_{X}](\infty)$ agree with each other up to ${\cal O}(\alpha^{2})$. Next,
the deviation of $S[p_{MV}]$ from $S[p_{V}](t=0)$

$\displaystyle\Delta S[p_{MV}]\equiv S[p_{MV}]-S[p_{V}](t=0)$
$\displaystyle=\frac{1}{2}(\log T_{M\infty}-\log T)\simeq\frac{\Delta
T_{M\infty}}{2T}=\frac{\Delta{\cal E}(\infty)}{T}.$ (160)

It is also found that $\Delta S[p_{MV}]=\Delta S[p_{V}](\infty)$ up to ${\cal
O}(\alpha^{2})$. The deviation of the entropy $S[p_{M}]$ from its initial
value $S[p](t=0)$ evaluated as

$\displaystyle\Delta S[p_{M}]$ $\displaystyle\equiv$ $\displaystyle
S[p_{M}]-S[p](t=0)=\Delta S[p_{MX}]+\Delta S[p_{MV}]$ (161)
$\displaystyle\simeq$ $\displaystyle\frac{\alpha^{2}}{4}+\frac{\Delta{\cal
E}(\infty)}{T}=\frac{\alpha^{2}}{4}\left(1+\frac{1}{k^{2}\lambda_{D}^{2}}\right)$
$\displaystyle=$ $\displaystyle\frac{\Delta{\cal
E}(\infty)}{T}\left(1+k^{2}\lambda_{D}^{2}\right).$

We note that $S[p]$ is invariant in the collisionless process although it
increases by $\Delta S[p_{M}]$ when the system reaches the thermal equilibrium
state due to collisions.

In Fig. 15 , the magnitudes of the entropies $S_{P}(X)$, $S_{P}(V)$,
$S_{P}(X,V)$, and the mutual information $I(X,V)$ are represented by the area
inside the corresponding contours for $t=0$, $t\rightarrow+\infty$ in the
collisionless process, and $t=+\infty$ in the collisional process. Here, the
entropies $S_{P}(X)\equiv S[P_{X}]$, $S_{P}(V)\equiv S[P_{X}]$, and
$S_{P}(X,V)\equiv S[P]$ takes non-negative values and they are related to
$S_{p}(X)\equiv S[p_{X}]$, $S_{p}(V)\equiv S[p_{V}]$, and $S_{p}(X,V)\equiv
S[p]$ by the relations shown in Eqs. (194) and (D). The entropy $S_{P}(X,V)$
does not change in the collisionless process although the Landau damping
increases $S_{P}(X)$ and $S_{P}(V)$ by $\Delta S_{P}(X)=\alpha^{2}/4$ and
$\Delta S_{P}(V)=\alpha^{2}/(4k^{2}\lambda_{D}^{2})$, respectively, as shown
in Eqs. (142) and (146), and the mutual information content $I(X,V)$ increases
by $\Delta I(X,V)=\Delta S_{P}(X)+\Delta S_{P}(V)$. Let us compare the limit
state at $t\rightarrow+\infty$ in the collisionless process and the thermal
equilibrium state reached by collisions. In the latter state, the values of
$S_{P}(X)$ and $S_{P}(V)$ remain the same as in the former up to the ${\cal
O}(\alpha^{2})$ accuracy. However, in the thermal equilibrium, the mutual
information quantity $I(X,V)$ vanishes and the entropy $S_{P}(X,V)$ of the
whole system increases by the amount that $I(X,V)$ decreases.

Figure 15:  The entropies $S_{P}(X)$, $S_{P}(V)$, $S_{P}(X,V)$ and the mutual
information $I(X,V)$ for $t=0$ (left), $t\rightarrow+\infty$ in the
collisionless process (middle), and $t=+\infty$ in the collisional process
(right). Their magnitudes are represented by the areas inside the
corresponding contours.

##  VII Conclusions and Discussion

In this paper, the one-dimensional Vlasov-Poisson system describing a plasma
consisting of electrons and uniformly distributed ions with infinite mass is
considered. Using analytical solutions and contour dynamics simulations, we
elucidate how the information entropies determined from the distribution
functions of the electron position and velocity variables evolve in the Landau
damping process. Under the initial condition given by the Maxwellian velocity
distribution with the perturbed density distribution in the form of the cosine
function, linear and quasilinear analytical solutions describing the time
evolutions of the electric field and the distribution function are obtained
and shown to be in good agreement with results from numerical simulations
based on contour dynamics.

A novel approximate integral formula including the effect of an infinite
number of complex eigenfrequencies to correctly evaluate the electric field is
presented. In addition, the linear analytical solutions for the electric field
and the distribution function near the initial time are expressed as series
expansions in time and velocity variables. The quasilinear analytical solution
describing the time evolution of the spatially averaged velocity distribution
function is obtained, and its validity is confirmed by the contour dynamics
simulation results. These analytical expressions of the linear and quasilinear
solutions are useful for verification of the accuracy of simulations of the
Vlasov-Poisson system using methods other than the contour dynamics as well.
Using the quasilinear analytical solution, it becomes possible to accurately
determine the time evolutions of the electron kinetic energy and the
background velocity distribution function associated with the Landau damping.
Furthermore, the time evolutions of the information entropies of the electron
position and velocity variables, and the mutual information are determined
with an accuracy of the order of the squared perturbation amplitude
$\alpha^{2}$. It is well known that, in a collisionless process, the
information entropy determined from the joint probability density distribution
function of position and velocity variables (or the phase-space distribution
function) is one of the Casimir invariants. On the other hand, the decrease in
the squared mean of spatial density fluctuations increases the information
entropy of the position variable, and the ratio of the increase in the
electron kinetic energy to the temperature equals the increase in the
information entropy of the velocity variable to the order of $\alpha^{2}$. The
sum of these increases in the information entropies of the position and
velocity variables yields the mutual information that is initially zero.

The relative entropy obtained by comparing the velocity distribution at time
$t$ with the initial distribution is a positive quantity of order of
$\alpha^{4}$. This leads to the fact that, even in the collisionless process,
the inequality in the form of the second law of thermodynamics holds in the
relation between the heat transfer from the Maxwellian velocity distribution
to the electric field and the conditional entropy of the electron position
variable for a given velocity distribution.

When Coulomb collisions are taken into account, they relax the distribution
function at $t\rightarrow+\infty$ in the collisionless process further to the
thermal equilibrium state. In this relaxation, the mutual information of the
position and velocity variables decreases to zero, although the information
entropies of the position and velocity variables do not change to the order of
$\alpha^{2}$. Then, the entropy determined from the phase-space distribution
increases by the amount of the decrease in the mutual information. It
indicates the validity of Boltzmann’s H-theorem. Future extensions of the
present work include studies on the position dependence of the phase-space
distribution function of order $\alpha^{2}$, which is not included in the
quasilinear solution, and the analysis of the information entropies and the
mutual information of the position and velocity variables to the order of
$\alpha^{4}$.

###### Acknowledgements.

This work is supported in part by the JSPS Grants-in-Aid for Scientific
Research (Grant Nos. 19H01879 and 24K07000) and in part by the NINS program of
Promoting Research by Networking among Institutions (Grant No. 01422301).
Simulations in this work were performed on “Plasma Simulator” (NEC SX-Aurora
TSUBASA) of NIFS with the support and under the auspices of the NIFS
Collaboration Research program (Grant Nos. NIFS23KIPT009 and NIFS24KISM007).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts of interest to disclose.

### Author Contributions

K. Maekaku: Data curation (lead); Investigation (lead); Visualization (lead);
Software (lead); Methodology (equal); Writing–original draft (equal);
Writing–review & editing (lead). H. Sugama: Conceptualization (lead); Formal
analysis (lead); Investigation (supporting); Methodology (equal);
Writing–original draft (equal); Writing-review & editing (supporting);
Supervision (lead). T.-H. Watanabe: Methodology (equal); Writing–review &
editing (supporting).

## DATA AVAILABILITY

The data that support the findings of this study are available from the
corresponding authors upon reasonable request.

##  Appendix A Plasma dispersion function

The plasma dispersion function is defined by Stix ; Miyamoto

$Z(\zeta)=\frac{1}{\sqrt{\pi}}\int_{-\infty}^{+\infty}dz\frac{e^{-z^{2}}}{z-\zeta},$
(162)

in the case of ${\rm Im}\zeta>0$ . Analytic continuation needs to be done to
define $Z(\zeta)$ in the case of ${\rm Im}\zeta\leq 0$. For ${\rm Im}\zeta=0$,
we have

$Z(\zeta)=\frac{1}{\sqrt{\pi}}P\int_{-\infty}^{+\infty}dz\frac{e^{-z^{2}}}{z-\zeta}+i\sqrt{\pi}e^{-\zeta^{2}}.$
(163)

The plasma dispersion function is also written as

$Z(\zeta)=i\sqrt{\pi}e^{-\zeta^{2}}\left[1+{\rm erf}(i\zeta)\right],$ (164)

where the error function ${\rm erf}$ is defined by

${\rm erf}\;z=\frac{2}{\sqrt{\pi}}\int_{0}^{z}e^{-s^{2}}ds.$ (165)

The plasma dispersion function satisfies

$[Z(-\zeta^{*})]^{*}=-Z(\zeta),$ (166)

where ∗ represents the complex conjugate. The derivative of $Z(\zeta)$ with
respective to $\zeta$ is given by

$Z^{\prime}(\zeta)\equiv\frac{dZ(\zeta)}{d\zeta}=-2[1+\zeta Z(\zeta)].$ (167)

The series expansion of $Z(\zeta)$ about $\zeta=0$ is given by

$Z(\zeta)=i\sqrt{\pi}e^{-\zeta^{2}}-\sum_{n=1}^{\infty}\frac{(-1)^{n-1}2\sqrt{\pi}}{\Gamma(n+1/2)}\zeta^{2n-1},$
(168)

and the asymptotic expansion of $Z(\zeta)$ for $|\zeta|\gg 1$ is written as

$Z(\zeta)=-\sum_{n=1}^{N}\frac{\Gamma(n-1/2)}{\sqrt{\pi}}\frac{1}{\zeta^{2n-1}}+\frac{\Gamma(N+1/2)}{\sqrt{\pi}}e^{-\zeta^{2}}\int_{i\infty}^{\zeta}\frac{e^{s^{2}}}{s^{2N}}ds.$
(169)

##  Appendix B Time reversibility of Vlasov-Poisson System

In the Vlasov-Poisson system, the time reversal map $T$ from the two-
dimensional phase space to itself is defined by

$T(x,v)=(x,-v),$ (170)

which satisfies

$T^{2}\equiv T\circ T=I,$ (171)

where $I$ represents the identity map. We can easily confirmed that, when
$f(x,v,t)$ is a solution of the Vlasov equation, $f(T(x,v),-t)=f(x,-v,-t)$ is
a solution as well. It should be noted that the property of the time reversal
map mentioned above holds whether the linear or nonlinear case is considered.
When the initial condition in Eq. (118) is employed, the initial distribution
function $f(x,v,t=0)$ satisfies the condition $f(x,-v,t=0)=f(x,v,t=0)$. Then,
for the solution $f(x,v,t)$ under this initial condition, $f(x,-v,-t)$ becomes
the solution satisfying the same initial condition. Thus, from the uniqueness
of the solution under the same initial condition, we find that
$f(x,-v,-t)=f(x,v,t)$ holds for any time $t$ and that $E(x,-t)=E(x,t)$ results
from Poisson’s equation.

##  Appendix C Contour-dynamics simulation for Vlasov-Poisson system

The Contour Dynamics (CD) method was proposed for solving inviscid and
incompressible motions in fluid mechanics in the two-dimensional space.
Zabusky . Here, we briefly describe the application of the CD method to the
solution of the one-dimensional Vlasov-Poisson system with the periodic
boundary condition [see Ref. Sato  for details].

Normalizing the time $t$, spatial variable $x$, velocity $v$, distribution
function $f(x,v,t)$, and electrostatic potential $\phi(x,t)$, appropriately,
the Vlasov-Poisson equations can be written as

$\frac{\partial f(x,v,t)}{\partial t}+v\frac{\partial f(x,v,t)}{\partial
x}+a(x,t)\frac{\partial f(x,v,t)}{\partial v}=0,$ (172)

and

$-\frac{\partial\phi(x,t)}{\partial
x^{2}}=1-\int_{-\infty}^{+\infty}f(x,v,t)dv\equiv F(x,t),$ (173)

respectively. In Eq. (172), the acceleration field is given by
$a(x,t)\equiv\partial\phi(x,t)/\partial x$.

In the CD method, $N_{max}$ contours $C_{m}(t)$ $(m=1,\cdots,N_{max})$ are
considered in the $(x,v)$-plane. Each point $(x,v)$ on $C_{m}(t)$ moves
according to the equations of motion

$\frac{dx}{dt}=v,\hskip 14.22636pt\frac{dv}{dt}=a(x,t).$ (174)

The solution $f(x,v,t)$ of the Vlasov equation is expressed as a piece-wise
constant distribution function,

$f_{pw}(x,v,t)=\sum_{m=1}^{N_{max}}\Delta f_{m}\;I[(x,v)\in S_{m}(t)]$ (175)

where $I[P]=1$ if $P$ is true and $I[P]=0$ otherwise. Here, $S_{m}(t)$
$(m=1,\cdots,N_{max})$ is the internal region of the $N_{max}$ contours
$C_{m}(t)$ in the phase $(x,v)$ space at time $t$, and $\Delta f_{m}$ denotes
the jump of the distribution function across the contour $C_{m}(t)$. The
electrostatic potential $\phi(x,t)$ can be calculated from

$\phi(x,t)=\int_{-L/2}^{L/2}G(x;{\xi})F({\xi},t)d{\xi}+\text{const}.$ (176)

Here, the Green function $G$ is given by

$G(x;\xi)=\frac{1}{2L}\left(|x-{\xi}|-\frac{L}{2}\right)^{2},$ (177)

which is the solution of

$\frac{\partial G(x;\xi)}{\partial x^{2}}=\frac{1}{L}-\delta(x-{\xi})$ (178)

with $L$ being the period length in the $x$ direction.

The acceleration $a(x,t)$ of each particle is obtained using the CD
representation as follows,

$a(x)=\sum_{m}\Delta f_{m}\oint_{C_{m}}G(x;{\xi})dv.$ (179)

By substituting this into Eq. (174) and integrating it in time, the time
evolutions of the contours $C_{m}(t)$ are obtained, and the distribution
function at each time $t$ is determined from Eq. (175). In practice, when
performing numerical simulations, the contours $C_{m}(t)$ are divided into a
finite number of nodes $(x_{i},v_{i})$. Between the nodes, the contour
integral in Eq. (179) is approximated by the integral along line segments, and
the numerical solution of the equations of motion for the finite number of
nodes $(x_{i},v_{i})$ is obtained.

In CD simulations, without using linear approximations, the solution of the
nonlinear Vlasov equation is obtained in the form of the piece-wise constant
distribution function as in Eq. (175). We use the piece-wise constant
distribution function $f_{Mpw}(v)$ corresponding to the Maxwellian equilibrium
velocity distribution, for which the outermost contours are placed at $v=\pm
5v_{t}$, and $f_{Mpw}(v)=0$ for $|v|>5v_{t}$. For all CD simulations in this
study, we follow the same procedure as in Ref. Sato  to give $f_{pw}(x,v,t=0)$
that corresponds to the initial condition in Eq. (118). The number of contours
used is 200, with each contour represented by 100 nodes. At time $t$, the
value of the distribution function $f(x,v,t)$ at any point $(x,v)$ in the
phase space is obtained by the linear Hermite interpolation from the values of
the distribution function on the contour to which the nodes $(x_{i},v_{i})$ in
the vicinity of $(x,v)$ belong, and is used for comparison with theoretical
predictions.

##  Appendix D Information entropies and mutual information in Vlasov-Poisson
system

We here regard the electron’s position and velocity as random variables
denoted by $X$ and $V$. The joint probability density function of $(X,V)$ is
represented by $p(x,v)$ which satisfies the normalization condition,

$\int_{-L/2}^{+L/2}dx\int_{-\infty}^{+\infty}dv\;p(x,v)=1.$ (180)

It is related to the distribution function in Eq. (4) by

$p(x,v)=\frac{1}{n_{0}L}f(x,v),$ (181)

where $n_{0}$ is the average number density of electrons. Here, the functions
$p$ and $f$ depends on time $t$ although $t$ is omitted from the arguments of
these functions for simplicity. The marginal probability distribution
functions for $X$ and $V$ are defined by

$p_{X}(x)=\int_{-\infty}^{+\infty}dv\;p(x,v)=\frac{1}{L}\frac{n(x)}{n_{0}}$
(182)

and

$p_{V}(v)=\int_{-L/2}^{+L/2}dx\;p(x,v)=\frac{1}{n_{0}}\langle f\rangle(v),$
(183)

respectively. Here, the electron density is denoted by $n(x)$.

The information (or Shannon) entropy Information  is originally defined for
probability distributions of discrete variables. Rigorously speaking, for
probability distributions of continuous variables, it should be called the
differential entropy. When dividing the interval $[-L/2,L/2)$ of the variable
$x$ is divided into $N_{x}$ intervals of equal width, the central value
$x_{i}$ in each interval is given by

$x_{i}=i\Delta x,\hskip 14.22636pt\Delta x=\frac{L}{N_{x}},\hskip
14.22636pti=-\frac{(N_{x}-1)}{2},\cdots,\frac{(N_{x}-1)}{2}.$ (184)

In the same way, when dividing the interval $[-v_{\rm max},v_{\rm max}]$ of
the variable $v$ into $N_{v}$ intervals of equal width, the central value
$v_{j}$ in each interval is given by

$v_{j}=j\Delta v,\hskip 14.22636pt\Delta v=\frac{2v_{\rm max}}{N_{v}},\hskip
14.22636pti=-\frac{(N_{v}-1)}{2},\cdots,\frac{(N_{v}-1)}{2}.$ (185)

Setting

$z_{ij}=(x_{i},v_{j}),\hskip 14.22636ptN_{z}=N_{x}N_{v},$ (186)

and letting $v_{\rm max}$ be sufficiently large, we consider the $(x,v)$ space
as a collection of $N_{z}$ cells, each of which is centered at
$z_{ij}=(x_{i},v_{j})$. Then, we approximate continuous variables $(x,v)$ by
discrete ones $z_{ij}=(x_{i},v_{j})$, and relate $p(x,v)$ to the probability
distribution function $P(x_{i},v_{j})$ of the discrete variables
$(x_{i},v_{j})$ by

$p(x,v)\Delta x\Delta v=P(x_{i},v_{j}),$ (187)

which satisfies the normalization condition,

$\int_{-L/2}^{+L/2}dx\int_{-\infty}^{+\infty}dv\;p(x,v)=\sum_{i,j}P(x_{i},v_{j})=1.$
(188)

Using $P(x_{i},v_{j})$, the information entropy is defined by

$S[P]\equiv\sum_{i,j}P(x_{i},v_{j})S(x_{i},v_{j})\equiv-\sum_{i,j}P(x_{i},v_{j})\log
P(x_{i},v_{j}),$ (189)

where

$S(x_{i},v_{j})\equiv-\log P(x_{i},v_{j})$ (190)

represents the self-entropy which is also called the self-information or
information content. Similarly, we employ $p(x,v)$ to define the differential
entropy by

$\displaystyle S[p]$ $\displaystyle\equiv$
$\displaystyle\int_{-L/2}^{+L/2}dx\int_{-\infty}^{+\infty}dv\;p(x,v)s(x,v)$
(191) $\displaystyle\equiv$
$\displaystyle-\int_{-L/2}^{+L/2}dx\int_{-\infty}^{+\infty}dv\;p(x,v)\log
p(x,v)$

where

$s(x,v)\equiv-\log p(x,v).$ (192)

We can easily derive the following relations,

$S(x_{i},v_{j})\equiv s(x,v)-\log(\Delta x\Delta v),$ (193)

and

$S[P]\equiv S[p]-\log(\Delta x\Delta v).$ (194)

Here, we should carefully note that $S(x_{i},v_{j})$ and $S[P]$ are non-
negative while $s(x,v)$ and $s[p]$ can take both positive and negative values
but satisfy the conditions,

$s(x,v)\geq\log(\Delta x\Delta v),\hskip 14.22636ptS[p]\geq\log(\Delta x\Delta
v).$ (195)

From the marginal probability distribution functions $p_{X}(x)$ and
$p_{V}(v)$, we can define the marginal distribution functions $P_{X}(x_{i})$
and $P_{V}(v_{j})$ of the discrete variables $x_{i}$ and $v_{j}$ by

$P_{X}(x_{i})=\sum_{j}P(x_{i},v_{j}),\hskip
14.22636ptP_{v}(v_{j})=\sum_{i}P(x_{i},v_{j}),$ (196)

which are related to $p_{X}(x)$ and $p_{V}(v)$ by

$P_{X}(x_{i})=p_{X}(x)\Delta x,\hskip 14.22636ptP_{V}(v_{j})=p_{V}(v)\Delta
v.$ (197)

Then, $P_{X}(x_{i})$ and $P_{V}(v_{j})$ are used to define the information
entropies $S[P_{X}]$ and $S[P_{V}]$ as

$\displaystyle S[P_{X}]$ $\displaystyle\equiv$
$\displaystyle\sum_{i}P_{X}(x_{i})S_{X}(x_{i})\equiv-\sum_{i}P_{X}(x_{i})\log
P_{X}(x_{i})$ $\displaystyle S[P_{V}]$ $\displaystyle\equiv$
$\displaystyle\sum_{i}P_{V}(v_{j})S_{V}(v_{j})\equiv-\sum_{i}P_{V}(v_{j})\log
P_{V}(v_{j}),\hskip 14.22636pt$ (198)

where the self-entropies $S_{X}(x_{i})$ and $S_{V}(v_{j})$ are defined by

$S_{X}(x_{i})\equiv-\log P_{X}(x_{i}),\hskip 14.22636ptS_{V}(v_{j})\equiv-\log
P_{V}(v_{j}).$ (199)

Similarly, we use $p_{X}(x)$ and $p_{V}(v)$ to the entropies $S[p_{X}]$ and
$S[p_{V}]$ by

$\displaystyle S[p_{X}]$ $\displaystyle\equiv$
$\displaystyle\int_{-L/2}^{+L/2}dx\;p_{X}(x)s_{X}(x)\equiv-\int_{-L/2}^{+L/2}dx\;p_{X}(x)\log
p_{X}(x)$ $\displaystyle S[p_{V}]$ $\displaystyle\equiv$
$\displaystyle\int_{-\infty}^{+\infty}dv\;p_{V}(v)s_{V}(v)\equiv-\int_{-\infty}^{+\infty}dv\;p_{V}(v)\log
p_{V}(v),$

where

$s_{X}(x)\equiv-\log p_{X}(x),\hskip 14.22636pts_{V}(v)\equiv-\log p_{V}(v).$
(201)

We can also confirm the following relations,

$\displaystyle S_{X}(x_{i})$ $\displaystyle\equiv$ $\displaystyle
s_{X}(x)-\log(\Delta x)\geq 0$ $\displaystyle S_{V}(v_{j})$
$\displaystyle\equiv$ $\displaystyle s_{V}(v)-\log(\Delta v)\geq 0$ (202)

and

$\displaystyle S[P_{X}]$ $\displaystyle\equiv$ $\displaystyle
S[p_{X}]-\log(\Delta x)\geq 0$ $\displaystyle S[P_{V}]$ $\displaystyle\equiv$
$\displaystyle S[p_{V}]-\log(\Delta v)\geq 0.$ (203)

The conditional entropies $S_{P}(V|X)$ and $S_{P}(X|V)$ are defined from the
probability distribution functions $P(x_{i},v_{j})$, $P_{X}(x_{i})$, and
$P(v_{j})$ by

$\displaystyle
S_{P}(V|X)\equiv-\sum_{i,j}P(x_{i},v_{j})\log\left(\frac{P(x_{i},v_{j})}{P_{X}(x_{i})}\right)$
$\displaystyle=S[P]-S[P_{X}]=S_{P}(X,V)-S_{P}(X)\geq 0,$ $\displaystyle
S_{P}(X|V)\equiv-\sum_{i,j}P(x_{i},v_{j})\log\left(\frac{P(x_{i},v_{j})}{P(v_{j})}\right)$
$\displaystyle=S[P]-S[P_{V}]=S_{P}(X,V)-S_{P}(V)\geq 0.$ (204)

Similarly, the conditional entropies $S_{p}(V|X)$ and $S_{p}(X|V)$ are defined
from $p(x,v)$, $p_{X}(x)$, and $p_{V}(v)$ by

$\displaystyle
S_{p}(V|X)\equiv-\int_{-L/2}^{L/2}dx\int_{\textendash\infty}^{+\infty}dv\;p(x,v)\log\left(\frac{p(x,v)}{p_{X}(x)}\right)$
$\displaystyle=S[p]-S[p_{X}]=S_{p}(X,V)-S_{p}(X)\geq\log(\Delta v),$
$\displaystyle
S_{p}(X|V)\equiv-\int_{-L/2}^{L/2}dx\int_{\textendash\infty}^{+\infty}dv\;p(x,v)\log\left(\frac{p(x,v)}{p_{V}(v)}\right)$
$\displaystyle=S[p]-S[p_{V}]=S_{p}(X,V)-S_{p}(V)\geq\log(\Delta x).$ (205)

To represent the mutual dependence of the random variables $X$ and $V$, the
mutual information $I(X,V)$ is defined by

$\displaystyle
I(X,V)\equiv-\sum_{i,j}P(x_{i},v_{j})\log\left(\frac{P(x_{i},v_{j})}{P_{X}(x_{i})P_{V}(v_{j})}\right)$
$\displaystyle=S_{P}(X)+S_{P}(V)-S_{P}(X,V)$
$\displaystyle=S_{P}(X)-S_{P}(X|V)=S_{P}(V)-S_{P}(V|X)$
$\displaystyle=-\int_{-L/2}^{L/2}dx\int_{\textendash\infty}^{+\infty}dv\;p(x,v)\log\left(\frac{P(x,v)}{P_{X}(x)P_{V}(v)}\right)$
$\displaystyle=S_{p}(X)+S_{p}(Y)-S_{p}(X,Y)$
$\displaystyle=S_{p}(X)-S_{p}(X|V)=S_{p}(V)-S_{p}(V|X)\geq 0.$ (206)

The mutual information takes the same non-negative value whether it is defined
from the distribution functions of the discrete variables or from those of the
continuous variables.

Hereafter, we express the time dependence of the distribution function by
writing $p(x,v,t)=f(x,v,t)/(n_{0}L)$,
$p_{X}(x,t)=\int_{-\infty}^{+\infty}dv\;p(x,v,t)=n(x,t)/(n_{0}L)$, and
$p_{V}(v,t)=\int_{-L/2}^{+L/2}dx\;p(x,v,t)=\langle f\rangle(v,t)/n_{0}$. It is
well known that, for the solution $f(x,v,t)$ of the Vlasov equation in Eq.
(4), the entropy $S[p]$ defined from $p(x,v,t)=f(x,v,t)/(n_{0}L)$ by

$S[p]\equiv
S_{p}(X,V)\equiv-\int_{-L/2}^{+L/2}dx\int_{-\infty}^{+\infty}dv\;p(x,v,t)\log
p(x,v,t)$ (207)

is an invariant,

$\frac{d}{dt}S[p]\equiv\ \frac{d}{dt}S_{p}(X,V)=0.$ (208)

On the other hand, the entropies defined from $p_{X}(x,t)$ and $p_{V}(v,t)$ by

$S[p_{X}]\equiv S_{p}(X)\equiv-\int_{-L/2}^{+L/2}dx\;p_{X}(x,t)\log
p_{X}(x,t)$ (209)

and

$S[p_{V}]\equiv S_{p}(V)\equiv-\int_{-\infty}^{+\infty}dv\;p_{V}(v,t)\log
p_{V}(v,t)$ (210)

depend on time $t$ in general.

The relative entropy (or Kullback-Leibler divergence) of the distribution
$p_{V}(v,t)$ at time $t$ relative to the initial distribution $p_{V}(t,0)$ is
defined as

$S(p_{V},t||p_{V},0)\equiv\int_{-\infty}^{+\infty}dv\;p_{V}(v,t)\log\left(\frac{p_{V}(v,t)}{p_{V}(v,0)}\right),$
(211)

which takes only non-negative values and vanishes if and only if
$p_{V}(v,t)=p_{V}(v,0)$. The non-negativity is derived from the inequality
$\log(x^{-1})=-\log x\geq 1-x$ for $x>1$. The relative entropy
$S(p_{V},t||p_{V},0)$ is the average difference between the information
quantities, $s_{V}(v,0)-s_{V}(v,t)=-\log p_{V}(v,0)-(-\log
p_{V}(v,t))=\log[p_{V}(v,t)/p_{V}(v,0)]$, and represents the amount of
information lost when approximating the correct velocity distribution
$p_{V}(v,t)$ at time $t$ by the initial distribution function $p_{V}(v,0)$.
From Eq. (108), we have

$\Omega_{t}(v_{0})=\log\left[\frac{p_{V}(u(v_{0},t),t)}{p_{V}(u(v_{0},t),0)}\right],$
(212)

and

$\displaystyle
S(p_{V},t||p_{V},0)=\int_{-\infty}^{+\infty}dv_{0}\;p_{V}(v_{0},0)\log\biggl{[}\frac{p_{V}(u(v_{0},t),t)}{p_{V}(u(v_{0},t),0)}\biggr{]}$
$\displaystyle=\int_{-\infty}^{+\infty}dv_{0}\;p_{V}(v_{0},0)\;\Omega_{t}(v_{0}).$
(213)

When the initial distribution function is the Maxwellian, using Eqs. (IV.1)
and (D) leads to

$S(p_{V},t||p_{V},0)=\frac{\Delta{\cal E}(t)}{n_{0}T}-\Delta S[P_{X}](t)\geq
0.$ (214)

Here, $\Delta{\cal E}(t)$ is the increase in the kinetic energy of electrons
per unit volume and equals the decrease in the electric field energy density,

$\displaystyle\Delta{\cal E}(t)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{+\infty}dv\;f_{0}(v_{0})\Delta{\cal
E}(u(v_{0},t),t)$ (215) $\displaystyle=$ $\displaystyle-\frac{\langle
E^{2}\rangle(t)-\langle E^{2}\rangle(0)}{8\pi}$ $\displaystyle=$
$\displaystyle-n_{0}\Delta Q(X|V).$

We now call two systems described by $p_{X}(x)$ and $p_{V}(v)$ as the
$X$-system and the $V$-system, respectively. Then, we can regard $\Delta
Q(X|V)=-\Delta{\cal E}/n_{0}$ in Eq. (215) as the energy transfer per electron
from the $V$-system to the $X$-system. Recalling that $\Delta S[P_{X}](t)$
represents the increase in the entropy $S[P_{X}]$ during the time interval
from 0 to $t$ and that $S_{P}(X,V)$ is time-independent, we can use Eq. (D) to
get

$\Delta S[P_{V}]=\Delta S_{P}(V)=-\Delta S_{P}(X|V).$ (216)

From Eqs. (214), (215), and (216), we find

$S(p_{V},t||p_{V},0)=\Delta S_{P}(X|V)-\frac{\Delta Q(X|V)}{T}\geq 0.$ (217)

Since the initial distribution $p_{V}(v,0)$ is given by the Maxwellian, we
regard the $V$-system here as the thermal reservoir with the temperature $T$.
Then, Eq. (217) implies that the heat transfer $\Delta Q(X|V)$ from the
thermal reservoir to the $X$-system and the change in the conditional entropy
of the $X$-system in contact with the thermal reservoir satisfy the inequality
in the form of the second law of thermodynamics.

When the values of $p_{V}(v,t)$ and $p_{V}(v,0)$ are close to each other, we
obtain

$\displaystyle
S(p_{V},t||p_{V},0)=-\int_{-\infty}^{+\infty}dv\;p_{V}(v,t)\log\left[\frac{p_{V}(v,0)}{p_{V}(v,t)}\right]$
$\displaystyle\simeq\int_{-\infty}^{+\infty}dv\;p_{V}(v,t)\left[-\left(\frac{p_{V}(v,0)}{p_{V}(v,t)}-1\right)+\frac{1}{2}\left(\frac{p_{V}(v,0)}{p_{V}(v,t)}-1\right)^{2}\right]$
$\displaystyle=\frac{1}{2}\int_{-\infty}^{+\infty}dv\;p_{V}(v,t)\left(\frac{p_{V}(v,0)}{p_{V}(v,t)}-1\right)^{2}$
$\displaystyle=\frac{1}{2}\int_{-\infty}^{+\infty}dv\;\frac{[p_{V}(v,t)-p_{V}(v,0)]^{2}}{p_{V}(v,t)}$
$\displaystyle\simeq\frac{1}{2}\int_{-\infty}^{+\infty}dv\;\frac{[p_{V}(v,t)-p_{V}(v,0)]^{2}}{p_{V}(v,0)}.$
(218)

Using Eq. (D) and the ordering parameter $\alpha$ for the perturbation
amplitude, $S(p_{V},t||p_{V},0)={\cal O}(\alpha^{4})$ is derived. We also find
from Eqs. (214) and (216) that $\Delta{\cal E}(t)/(n_{0}T)=\Delta S[P_{X}](t)$
and $\Delta S_{P}(X|V)=\Delta Q(X|V)/T$ hold up to ${\cal O}(\alpha^{2})$.

## References

  * (1) L. D. Landau, J. Exp. Theor. Phys. 16, 574 (1946); J. Phys. USSR 10, 26 (1946).
  * (2) K. M. Case, Ann. Phys. 7, 349 (1959).
  * (3) N. G. Van Kampen and B. U. Felderhof, Theoretical Methods in Plasma Physics (North-Holland, Amsterdam, 1967), Chap. 12.
  * (4) D. R. Nicholson, Introduction to Plasma Theory (John Wiley & Sons, New York, 1983), Chap. 6.
  * (5) R. D. Hazeltine and F. L. Waelbroeck, The Framework of Plasma Physics (Perseus books, Massacheusetts, 1998), Chap. 6.
  * (6) I. Prigogine and T. Petrosky, J. Plasma Phys. 59, 611 (1998).
  * (7) A. Zocco, J. Plasma Phys. 81, 905810402 (2015).
  * (8) F. F. Chen, Introduction to Plasma Physics and Controlled Fusion 3rd Ed. (Springer, 2016), Chap. 7.
  * (9) G. W. Hammett and F. W. Perkins, Phys. Rev. Lett. 64, 3019 (1990).
  * (10) C. Villani, Phys. Plasmas 21, 030901 (2014).
  * (11) F. Zonca, L. Chen, and R. A. Santoro, Plasma Phys. Control. Fusion 38, 2011 (1996).
  * (12) H. Sugama and T.-H. Watanabe, J. Plasma Phys. 72, 825 (2006).
  * (13) H. Sugama and T.-H. Watanabe, Phys. Plasmas 13, 012501 (2006).
  * (14) A. I. Smolyakov, M. F. Bashir, A. G. Elfimov, M. Yagi, and N. Miyato, Plasma Phys. Rep. 42, 407 (2016).
  * (15) T. M. Cover and J. A. Thomas, Elements of Information Theory, 2nd Ed. (Wiley, 2006), Chap. 2.
  * (16) P. J. Morrison, Z. Naturforsch. 42a, 1115 (1987).
  * (17) K. Maekaku and Z. Yoshida, Phys. Plasmas 29, 082303 (2022).
  * (18) W. E. Drummond and D. Pines, Suppl. Nucl. Fusion Part 3, 1049 (1962).
  * (19) A. A. Vedenov, E. P. Velikhov, and R. Z. Sagdeev, Suppl. Nucl. Fusion Part 2, 465 (1962).
  * (20) N. J. Zabusky, M. H. Hughes, and K.V. Roberts, J. Comput. Phys. 30, 96 (1979).
  * (21) H. Sato, T.-H. Watanabe, and S. Maeyama, J. Comput. Phys. Ann. Phys. 445, 110626 (2021).
  * (22) R. D. Hazeltine and J. D. Meiss, Plasma Confinement (Addison-Wesley, Redwood City, CA, 1992), Chap. 5.
  * (23) D. J. Evans and D. J. Searles, Adv. Phys. 51, 1529 (2002).
  * (24) D. J. Evans, D. J. Searles, and S. R. Williams, Fundamentals of Classical Statistical Thermodynamics (Wiley-VCH, Weinheim, 2016).
  * (25) T. H. Stix, Waves in Plasmas (American Institute of Physics, New York, 1992), p. 202.
  * (26) K. Miyamoto, Plasma Physics and Controlled Fusion, 2nd Ed. (Springer Berlin Heidelberg, 2016), p. 212.