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 that 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.

## I. INTRODUCTION

Landau damping is one of the most intriguing physical processes involving collective interactions between waves and particles in collisionless plasmas.^{1} 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.^{2–15}

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.^{2–5} 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. Incidentally, it is shown in Ref. 7 that the plasma dispersion function can be expressed in the form of an infinite continued fraction. From this fact, we can also understand that an infinite number of complex eigenfrequencies exist as zeros of the dielectric function. It is also shown in Ref. 8 that a very high number of poles are required for the correct calculation of density and pressure at an early stage in the Vlasov–Poisson system.

When treating the position *x* and velocity *v* of electrons as random variables, denoted as *X* and *V*, respectively, the information entropy^{16}^{,} $Sp(X,V)\u2261\u2212\u222b\u222bp(x,v,t)\u2009log\u2009p(x,v,t)\u2009dx\u2009dv$ 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,^{17,18} and it does not depend on time *t*. On the other hand, the entropies $Sp(X)\u2261\u2212\u222bpX(x,t)\u2009log\u2009pX(x,t)\u2009dx$ and $Sp(V)\u2261\u2212\u222bpV(v,t)\u2009log\u2009pV(v,t)\u2009dv$ determined from the marginal probability density functions $pX(x,t)=\u222bp(x,v)\u2009dv$ and $pV(v,t)=\u222bp(x,v)\u2009dx$ 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.^{19,20} Using these linear and quasilinear analytical solutions, we clarify the time evolution of the information entropies $Sp(X),\u2009Sp(V)$, and the mutual information $I(X,V)\u2261Sp(X)+Sp(V)\u2212Sp(X,V)$ associated with the Landau damping. We here note that, in Ref. 21, the energy conversion in phase space density moments in the Vlasov–Maxwell system is investigated by separating the kinetic entropy [which corresponds to $Sp(X,V)$] to the two parts [which correspond to $Sp(X)$ and $Sp(X,V)\u2212Sp(X)$]. Also, in Ref. 22, the role of the entropy production in momentum transfer in the Vlasov–Maxwell system is discussed based on the information theory.

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.^{23,24} We also obtain the velocity distribution *p _{V}* in the limit of $t\u2192+\u221e$ 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 $Sp(X,V)$, and thus validate Boltzmann's H-theorem.

^{25}

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\u2192+\u221e$, 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 THE VLASOV–POISSON SYSTEM

### A. The Vlasov–Poisson system

*e*and

*m*are the electron charge and mass, respectively, and

**E**(

**x**,

*t*) represents the electric field. The distribution function $f(x,v,t)$ of electrons is defined such that $f(x,v,t)d3xd3v$ represents the number of electrons in the phase space volume element $d3x\u2009d3v$ around $(x,v)$ at time

*t*. We ignore the motion of ions that have the uniform density

*n*

_{0}and infinite mass. Then, the electric field $E(x,t)$ is determined by Poisson's equation

*y*and

*z*, and that $Ey=Ez=0$. Integration of $f(x,v,t)$ with respect to

*v*and

_{y}*v*is done to define

_{z}*x*,

*v*) phase space as

*v*=

*v*and $E(x,t)=Ex(x,t)$ are used. Poisson's equation in Eq. (2) is rewritten as

_{x}*x*,

*v*) phase space instead of the six-dimensional $(x,v)$ space. The number of electrons in the area element $dx\u2009dv$ about the point (

*x*,

*v*) in the phase space at time

*t*is given by $f(x,v,t)dx\u2009dv$.

### B. Linearization

*E*(

*x*,

*t*) are assumed to be given by

*x*-direction is given by

*k*> 0. Substituting Eq. (6) into Eq. (4) and using Eq. (8), we obtain the linearized Vlasov equation

### C. Laplace transform

*E*(

*k*,

*t*) are denoted by $f1(k,v,\omega )$ and $E(k,\omega )$, respectively. Instead of the variable

*p*used in the conventional Laplace transform, we put $p=\u2212i\omega $ and employ the complex frequency

*ω*to perform the analysis in the manner similar to the Fourier transform as seen in Ref. 4. The inverse Laplace transform gives

*L*and

_{f}*L*need to pass above all poles of $f1(k,v,\omega )$ and $E(k,\omega )$ on the complex

_{E}*ω*-plane, respectively. Now, the Vlasov equation and Poisson's equation in Eqs. (9) and (11) are represented in terms of $f1(k,v,\omega )$ and $E(k,\omega )$ by

*E*(

*k*,

*t*) and $f1(k,v,t)$ for

*t*> 0 as

*ω*is represented by $\u2202\omega \u2261\u2202/\u2202\omega $.

### D. Conditions of equilibrium and initial perturbation

*T*represents the equilibrium temperature and

*λ*is defined by

_{D}*ω*is given by

*α*is a small constant. Then, from Poisson's equation, we have

*E*(

*k*,

*t*) and $f1(k,v,t)$ given in Eqs. (23) and (24), respectively.

### E. Solution of initial value problem

*t*> 0 as

Figure 2 shows *E*(*k*, *t*) calculated as a function of time using Eq. (37) for $k\lambda D=1/2$. The summation over *μ* in Eq. (37) is done by including a finite number of pairs of the complex eigenfrequencies $(\omega \mu ,\u2212\omega \mu *)$ from $(\omega \xf8,\u2212\omega \xf8*)$ 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 *μ* 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}$.

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

## III. APPROXIMATE REPRESENTATION OF LINEAR SOLUTIONS IN EARLY TIME

### A. Approximate expression for complex eigenfrequencies with large absolute values

*n*is an integer. Then, we define $\delta \mu $ by

*n*needs to be large, $n\u226b1$, to satisfy the condition $|\zeta \mu |\u226b1$. Recall that when $\omega \mu $ is a complex-valued eigenfrequency, $\u2212\omega \mu *$ is so, too. Then, the approximate values of complex-valued eigenfrequency with large absolute values are denoted by

*ω*and $\u2212\omega n*$ which are given by

_{n}*n*is used instead of

*μ*to specify different approximate eigenfrequencies. The approximation is good for $n\u226b1$.

In Fig. 4, the exact complex eigenfrequencies and the approximate eigenfrequencies for $k\lambda D=1/2$ are plotted as red and dots, respectively, in the region where $Re(\omega )>0$. Here, Eqs. (53) and (55) with $n=0,1,2,\u2026$ 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.

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

*C*is a positive constant and $\tau \u2261kvTt$. This form of the infinite series is included in Eq. (37) for

*E*(

*k*,

*t*). When $n\u226b1$ and $|\zeta n|\u226bC$, we use Eqs. (53) and (55) to make the following approximations:

*τ*= 0, using Eq. (56) and $|\delta n|\u226a1$ gives

*N*,

*E*(

*k*,

*t*), diverges to infinity as $t\u2192+0$. On the other hand, when

*v*is fixed, the series expansion for representing $f1(k,v,t)$ uniformly converges. It is because the extra factor $1/(\omega n\u2212kv)$ included in the series expansion accelerates the convergence as $n\u2192\u221e$. Now, we use

*E*(

*k*,

*t*) as

*τ*= 0, we obtain

*ζ*-plane where the radius of

*C*is given by $R\u226b1$ and the orientation of

*C*is taken clockwise by varying the argument

*θ*of

*ζ*from

*π*to $\u2212\pi $. Comparing Eqs. (67) and (68) shows that changing the order of operations $lim\tau \u2192+0$ and $\u2211\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 \u2192+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 $14(3\u2212\kappa 2)$. We use Eq. (68) to confirm that the approximate solution $EN(k,t)$ given in Eq. (65) can correctly approach to the correct initial value as $t\u2192+0$,

*N*= 100. There, the approximate solution $EN(k,t)$ for

*N*= 100 and the CD simulation results for $\alpha =0.01$ are plotted by 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 $limN\u2192\u221eE<N(k,0)/E(k,0)=(3\u2212\kappa 2)/4=11/16=0.6875$ for $\kappa \u2261k\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),\u2009E>N(k,t)$, and $EN(k,t)=E<N(k,t)+E>N(k,t)$ for

*N*= 1. Even in this case of

*N*= 1, the difference between $EN(k,t)$ and the CD simulation result near

*t*= 0 is as small as about 5%.

As described in Appendix D, 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 D. Accordingly, $f1(k,\u2212v,\u2212t)=f1(k,v,t)$ and $E(k,\u2212t)=E(k,t)$ hold for any *t* and *v*. Therefore, *E*(*k*, *t*) and $f1(k,v,t)$ for *t* < 0 are immediately given from Eqs. (37) and (38) with making use of the time reversal.

### C. Approximate solution in early time

*ζ*-plane. Here, the radius of

*C*is given by $R\u226b1$ and the orientation of

*C*is taken clockwise by varying the argument

*θ*of

*ζ*from

*π*to $\u2212\pi $. Supplementary explanations about the derivation of Eq. (76) are given after Eq. (B10) in Appendix B. Equation (37) is rewritten here for $\tau \u2261kvTt>0$ as

*E*(

*k*,

*t*) about $\tau =kvTt=0$ is written as

*ζ*yields

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 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 12$. 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 pt<2$ (which corresponds to $\tau \u2261kvTt\u22612(k\lambda D)\omega pt<1$). Figure 7 shows the plot of $f1(x,v,t)/f0(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 pt\u22641$.

## IV. ANALYTICAL SOLUTION BASED ON QUASILINEAR THEORY AND ITS VERIFICATION BY NUMERICAL SIMULATION

### A. Spatially averaged distribution function

*x*and the period length is given by $L\u22612\pi /k$. We use

*x*. The Vlasov equation is averaged in

*x*to yield

*x*, and

*a*(

*v*,

*t*) represents the average acceleration of electrons defined by

*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

*x*-direction is defined as a functional of the distribution function $f(x,v,t)$ by

*x*-averaged distribution function $\u27e8f\u27e9(v,t)$ to define the entropy density $S[\u27e8f\u27e9]$ by

*t*. From the viewpoint of information theory,

^{16}the entropy density $S[\u27e8f\u27e9]$ is given from the average of

*dv*is regarded as an infinitesimal positive constant. Supplementary explanations on information entropies in the Vlasov–Poisson system are presented in Appendix E. 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

*v*

_{0}and the history of acceleration $a(v,t\u2032)(0\u2264t\u2032\u2264t)$. The interval $[v0\u221212dv0,v0+12dv0]$ in the

*v*-space evolves to the interval $[u(v0,t)\u221212du(v0,t),u(v0,t)+12du(v0,t)]$ at time

*t*when the velocities of the electrons in the interval are given by $u(v\u2032,t\u2032)$ ( $v0\u221212dv0\u2264v\u2032\u2264v0+12dv0$, $0\u2264t\u2032\u2264t)$. We note that the number (or probability) of the electrons found in the interval $[u(v0,t)\u221212du(v0,t),u(v0,t)+12du(v0,t)]$ is invariant in time,

*v*at time

*t*along the trajectory $u(v0,t)$ in the

*v*-space. Then, the increase in

*S*(

*v*,

*t*) along the trajectory during the time interval $[0,t]$ is given by

*t*to derive

*x*-averaged distribution function $\u27e8f\u27e9(v,0)=f0(v)$ is given by the Maxwellian equilibrium distribution function in Eq. (25). Then, we have

*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

Here, we should note that $\Omega t(v0)$ defined in Eq. (116) corresponds to the dissipation function employed by Evans and Searles to present their fluctuation theorem.^{26,27} 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(v0)$ in our case. However, as explained in Appendix E, we can show the non-negativity of the expected value of $\Omega t(v0)$, which leads to the inequality in the form of the second law of thermodynamics.

### B. Analysis to second order in perturbation amplitude

*x*-averaged distribution function considered in Sec. IV A up to the second order in the perturbation amplitude

*α*based on the results of the linear analysis in Sec. II. Recall that the initial condition for the distribution function is given by

*t*> 0, from which the electric field $E(x,t)=Re[E(k,t)\u2009exp(i\u2009k\u2009x)]$ is obtained. Hereafter, we represent the order of magnitude of the small perturbation amplitude by

*α*. Then, the linear solutions $f1(x,v,t)=Re[f1(k,v,t)\u2009exp(i\u2009k\u2009x)]=O(\alpha )$ and $E(x,t)=Re[E(k,t)exp(i\u2009k\u2009x)]=O(\alpha )$ are used to derive

*v*caused by the acceleration due to the electric field during the time interval $[0,t]$. Recall that $\u2202a(v,t)/\u2202v$ represents the rate of change in the information content $S(v,t)\u2261\u2212log\u2009\u27e8f\u27e9(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]$.

*α*, the change in the information (or Gibbs) entropy $\Delta S[\u27e8f\u27e9](t)\u2261S[\u27e8f\u27e9](t)\u2212S[\u27e8f\u27e9](0)$ is equal to the inverse temperature $1/T$ multiplied by the energy transfer $\Delta E(t)$ from the electric field energy to the kinetic energy during the time interval $[0,t]$.

### C. Calculation of $\Delta E(v,t),\u2009\Delta S(v,t)$, and $f2(v,t)$

*v*during the time interval $[0,t]$ are given by Eqs. (125) and (126), respectively. The electric field

*E*(

*k*,

*t*) and the perturbed distribution function $f1(k,v,t)$ appearing in Eqs. (125) and (126) 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 $f1(k,v,t)$ about $\tau =kvTt=0$ are given by Eqs. (78) and (85), respectively, which converge more quickly for smaller time

*t*. Therefore, to evaluate $\Delta E(v,t)$ and $\Delta S(v,t)$, the expressions of

*E*(

*k*,

*t*) and $f1(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 $O(\tau 13)$ in the integrands in Eqs. (125) and (126) for $\tau \u2261kvTt\u22641$ and using two pairs of the complex frequencies for Eqs. (37) and (38) to obtain

*E*(

*k*,

*t*) and $f1(k,v,t)$ for $\tau \u2261kvTt\u22651$. Once that $\Delta E(v,t)$ and $\Delta sV(v,t)$ are obtained following the procedures described above, we can use Eq. (124) to calculate $f2(v,t)$ by

Figure 8 shows contours of $\Delta 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 E(v,t)>0$ which indicates the increase in the electron kinetic energy. Here, the resonance velocity is given by $vres\u2261Re(\omega \xf8)/k=2.831\u2009vt$. As time progresses, $\Delta E(v,t)$ becomes more concentrated around the resonance velocities $v=\xb1vres$ while $\Delta S(v,t)$ clearly shows the positive maximum (negative minimum) at a slightly smaller (larger) absolute value $|v|$ than *v _{res}*.

The top and bottom panels of Fig. 9 respectively show the distributions of $f2(v,t)/(\alpha 2n0vt\u22123)$ and $f2(v,t)/(\alpha 2f0(v))$ in the (*v*, *t*) plane, calculated from the difference between $\Delta E(v,t)/T$ and $\Delta S(v,t)$ using Eqs. (129) and (130). As time progresses, the distribution of $f2(v,t)$ spreads from around *v* = 0 to *v _{res}* while such spreading is not clearly seen in the distribution of $f2(v,t)/f0(v)$. Figure 10 shows $f2(v,t)/(\alpha 2n0vt\u22123)$ and $f2(v,t)/(\alpha 2f0(v,t))$ as functions of

*v*, obtained from the analytical solution for $\omega pt=0.1,0.5,1,5$, and 10. We see that $f2(v,t)$ oscillate along the

*v*direction and the number of oscillations increases with increasing time. Positive peaks and negative troughs of $f2(v,t)/f0(v)$ can be seen around the resonant velocities $v=\xb1vres$. The red dots represent the results of CD simulations for $\alpha =0.1$, which agree well with the analytical solution. However, for $\omega pt\u22655$, 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 $f1(k,v,t)$ driven by $\u2202f0(v)/\u2202v$ which is set to zero in the neighborhood of

*v*= 0 for the CD simulations. Then, the absolute value of $f2(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 $f2(v,t)/(\alpha 2n0vt\u22123)$ and $f2(v,t)/(\alpha 2f0(v))$ obtained by the analytical formulas for $\omega pt=10$, 20, 50, and 100 are shown by 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\u2192+\u221e,\u2009f2(v,t)$ and $f2(v,t)/f0(v)$ converge to the structures which have the positive maxima (negative minima) at $|v|$ slightly larger (smaller) than $vres=2.831\u2009vt$. 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.

*A*(

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

*v*

^{4}, the kurtosis of the velocity as a random variable at time

*t*is given by

*α*by

*v*= 0 for small

*t*as seen in Figs. 9 and 10. As

*t*increases, peaks of $f2(v)$ near the resonant velocities $v=\xb1vres$ 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(\u221e)=23.92\alpha 2$ as $t\u2192+\u221e$.

## V. TIME EVOLUTION OF INFORMATION ENTROPIES

*t*is given by

*X*and

*V*, respectively, as explained in Appendix E. The joint probability density function of

*X*and

*V*are denoted by $p(x,v,t)$ [see Eq. (E2)], which is integrated with respect to

*x*and

*y*to give the marginal probability density functions $pX(x,t)$ and $pV(v,t)$, respectively [see Eqs. (E3) and (E4)]. The entropy $Sp(X)\u2261S[pX]$ is derived from the electron's position distribution function $pX(x,t)$ as

*t*is given by

*t*is represented by

*X*and

*V*is defined by

*X*and

*V*are statistically independent at

*t*= 0, and accordingly

*I*(

*X*,

*V*) during the time interval from 0 to

*t*is written as

*t*, $\Delta Sp(X)\u2261\Delta S[pX],\u2009\Delta Sp(V)\u2261\Delta S[pV]$ and $\Delta I(X,V)$ have the same form but different magnitudes.

*t*increases, $S(pV,t||pV,0)$ increases while showing modulation and approaches to the limit value that corresponds to the limit function $limt\u2192+\u221ef2(v,t)$ in Fig. 11. As described in Appendix E, the non-negativity of $S(pV,t||pV,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

*t*= +

*∞*while conserving total particles' number, momentum, and energy. The Maxwellian distribution function $fM\u221e(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 $fM\u221e(v)$ into the left-hand side into the Vlasov equation, it also vanishes. It can be proven from the above-mentioned fact that $nM\u221e,\u2009TM\u221e$, and $uM\u221e$ 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

*X*and

*V*are statistically independent, and accordingly, the mutual information of

*X*and

*V*vanishes,

In Fig. 15, the magnitudes of the entropies $SP(X),\u2009SP(V),SP(X,V)$, and the mutual information *I*(*X*, *V*) are represented by the area inside the corresponding contours for *t* = 0, $t\u2192+\u221e$ in the collisionless process, and *t* = + *∞* in the collisional process. Here, the entropies $SP(X)\u2261S[PX],\u2009SP(V)\u2261S[PX]$, and $SP(X,V)\u2261S[P]$ takes non-negative values and they are related to $Sp(X)\u2261S[pX],Sp(V)\u2261S[pV]$, and $Sp(X,V)\u2261S[p]$ by the relations shown in Eqs. (E15) and (E24). The entropy $SP(X,V)$ does not change in the collisionless process although the Landau damping increases $SP(X)$ and $SP(V)$ by $\Delta SP(X)=\alpha 2/4$ and $\Delta SP(V)=\alpha 2/(4k2\lambda D2)$, respectively, as shown in Eqs. (142) and (146), and the mutual information content *I*(*X*, *V*) increases by $\Delta I(X,V)=\Delta SP(X)+\Delta SP(V)$. Let us compare the limit state at $t\u2192+\u221e$ in the collisionless process and the thermal equilibrium state reached by collisions. In the latter state, the values of $SP(X)$ and $SP(V)$ remain the same as in the former up to the $O(\alpha 2)$ accuracy. However, in the thermal equilibrium, the mutual information quantity *I*(*X*, *V*) vanishes and the entropy $SP(X,V)$ of the whole system increases by the amount that *I*(*X*, *V*) decreases.

## 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\u2192+\u221e$ 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$.

## ACKNOWLEDGMENTS

This work was 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

^{28,29}

*ζ*is given by

*ζ*= 0 is given by

### APPENDIX B: CONTOURS FOR INTEGRALS IN THE COMPLEX PLANE

This appendix presents supplementary explanations about contours used for the integrals in the complex plane, which appear in Secs. II and III. Figure 16 shows the contours explained in this appendix.

*ζ*-plane from $Re(\zeta )=\u2212\u221e$ to $Re(\zeta )=+\u221e$. We let

*L*denote the real axis as the contour for integration. We also consider

_{E}*τ*to take a positive fixed value. Then, to apply the residue theorem for deriving the third line of Eq. (37), we need to deform the integration contour from

*L*to the closed one by adding to

_{E}*L*the semicircle $C\u2212:\zeta =Rei\theta $ in the lower half-plane $Im(\zeta )<0$ where the radius is given by $R\u226b1$ and the orientation of $C\u2212$ is taken clockwise by varying the argument

_{E}*θ*of

*ζ*from 0 to $\u2212\pi $. We also choose $C\u2212$ such that $A(\zeta )$ has no poles on $C\u2212$. The deformation of the contour to $C\u2212$ mentioned above is justified because the part of the contour integral along the semicircle $C\u2212$ in the lower half-plane $Im(\zeta )<0$ vanishes in the limit of $R\u2192+\u221e$ as shown below.

*θ*of

*ζ*from

*π*to 0. Then, $C=C++C\u2212$ represents the circle with the radius $R\u226b1$ and the clockwise orientation. For $|\zeta |\u2261R\u226b1$, we have

*R*limit. Therefore, a closed integration contour including $C\u2212$ can also be used to apply the residue theorem in deriving Eq. (38).

*L*with $C+$ in Eq. (B8) and obtain

_{E}*R*limit. Therefore, when we substitute Eq. (B10) into Eq. (B9), it is still valid even though the contour

*C*contains the arc $CII$ where Eq. (B10) does not hold. Thus, we obtain the second line of Eq. (76) where the residue theorem is also used to derive the third line. Taking the limit of $\tau \u2192+0$, Eqs. (B9) and (76) give

*τ*is assumed to be positive. We finally consider the case of

*τ*= 0 in which $e\u2212i\zeta \tau =1$ and Eq. (B7) does not hold. Again, using Eq. (B3) in the large

*R*limit, we have

*C*is taken clockwise. Then, Eq. (B12) is used to derive Eq. (68) where the summation of the residues about the poles of $A(\zeta )$ are also shown.

### APPENDIX C: CONTOUR-DYNAMICS SIMULATION FOR THE VLASOV–POISSON SYSTEM

The Contour Dynamics (CD) method was proposed for solving inviscid and incompressible motions in fluid mechanics in the two-dimensional space.^{23} 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. 24 for details).

*t*, spatial variable

*x*, velocity

*v*, distribution function $f(x,v,t)$, and electrostatic potential $\varphi (x,t)$, appropriately, the Vlasov–Poisson equations can be written as

*N*contours $Cm(t)(m=1,\u2026,Nmax)$ are considered in the (

_{max}*x*,

*v*)-plane. Each point (

*x*,

*v*) on $Cm(t)$ moves according to the equations of motion

*P*is true and $I[P]=0$ otherwise. Here, $Sm(t)(m=1,\u2026,Nmax)$ is the internal region of the

*N*contours $Cm(t)$ in the phase (

_{max}*x*,

*v*) space at time

*t*, and $\Delta fm$ denotes the jump of the distribution function across the contour $Cm(t)$. The electrostatic potential $\varphi (x,t)$ can be calculated from

*G*is given by

*L*being the period length in the

*x*direction.

*a*(

*x*,

*t*) of each particle is obtained using the CD representation as follows:

*t*is determined from Eq. (C4). In practice, when performing numerical simulations, the contours $Cm(t)$ are divided into a finite number of nodes (

*x*,

_{i}*v*). Between the nodes, the contour integral in Eq. (C8) is approximated by the integral along line segments, and the numerical solution of the equations of motion for the finite number of nodes (

_{i}*x*,

_{i}*v*) is obtained.

_{i}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. (C4). We use the piece-wise constant distribution function $fMpw(v)$ corresponding to the Maxwellian equilibrium velocity distribution, for which the outermost contours are placed at $v=\xb15vt$, and $fMpw(v)=0$ for $|v|>5vt$. For all CD simulations in this study, we follow the same procedure as in Ref. 24 to give $fpw(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*) in the vicinity of (

_{i}*x*,

*v*) belong, and is used for comparison with theoretical predictions.

### APPENDIX D: TIME REVERSIBILITY OF THE VLASOV–POISSON SYSTEM

*T*from the two-dimensional phase space to itself is defined by

*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),\u2212t)=f(x,\u2212v,\u2212t)$ 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,\u2212v,t=0)=f(x,v,t=0)$. Then, for the solution $f(x,v,t)$ under this initial condition, $f(x,\u2212v,\u2212t)$ 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,\u2212v,\u2212t)=f(x,v,t)$ holds for any time

*t*and that $E(x,\u2212t)=E(x,t)$ results from Poisson's equation.

### APPENDIX E: INFORMATION ENTROPIES AND MUTUAL INFORMATION IN THE VLASOV–POISSON SYSTEM

*X*and

*V*. The joint probability density function of (

*X*,

*V*) is represented by

*p*(

*x*,

*v*) which satisfies the normalization condition,

*n*

_{0}is the average number density of electrons. Here, the functions

*p*and

*f*depend 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

*n*(

*x*).

^{16}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 $[\u2212L/2,L/2)$ of the variable

*x*is divided into

*N*intervals of equal width, the central value

_{x}*x*in each interval is given by

_{i}*v*into

*N*intervals of equal width, the central value

_{v}*v*in each interval is given by

_{j}*x*,

*v*) space as a collection of

*N*cells, each of which is centered at $zij=(xi,vj)$. Then, we approximate continuous variables (

_{z}*x*,

*v*) by discrete ones $zij=(xi,vj)$, and relate

*p*(

*x*,

*v*) to the probability distribution function $P(xi,vj)$ of the discrete variables (

*x*,

_{i}*v*) by

_{j}*p*(

*x*,

*v*) to define the differential entropy by

*s*(

*x*,

*v*) and $s[p]$ can take both positive and negative values but satisfy the conditions,

*x*and

_{i}*v*by

_{j}*p*(

*x*,

*v*), $pX(x)$, and $pV(v)$ by

*X*and

*V*, the mutual information

*I*(

*X*,

*V*) is defined by

*t*in general,.

*t*relative to the initial distribution $pV(t,0)$ is defined as

*t*by the initial distribution function $pV(v,0)$. From Eq. (108), we have

*X*-system and the

*V*-system, respectively. Then, we can regard $\Delta Q(X|V)=\u2212\Delta E/n0$ in Eq. (E36) as the energy transfer per electron from the

*V*-system to the

*X*-system. Recalling that $\Delta S[PX](t)$ represents the increase in the entropy $S[PX]$ during the time interval from 0 to

*t*and that $SP(X,V)$ is time-independent, we can use Eq. (E26) to get

*V*-system here as the thermal reservoir with the temperature

*T*. Then, Eq. (E38) 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.

*α*for the perturbation amplitude, $S(pV,t||pV,0)=O(\alpha 4)$ is derived. We also find from Eqs. (E35) and (E37) that $\Delta E(t)/(n0T)=\Delta S[PX](t)$ and $\Delta SP(X|V)=\Delta Q(X|V)/T$ hold up to $O(\alpha 2)$.

## REFERENCES

*Theoretical Methods in Plasma Physics*

*Introduction to Plasma Theory*

*The Framework of Plasma Physics*

*The Vlasov Equation I. History and General Properties*

*Introduction to Plasma Physics and Controlled Fusion*

*Elements of Information Theory*

*Plasma Confinement*

*Fundamentals of Classical Statistical Thermodynamics*

*Waves in Plasmas*

*Plasma Physics and Controlled Fusion*