Thermodynamic Explanation of Landau Damping by Reduction to Hydrodynamics

Landau damping is the tendency of solutions to the Vlasov equation towards spatially homogeneous distribution functions. The distribution functions, however, approach the spatially homogeneous manifold only weakly, and Boltzmann entropy is not changed by the Vlasov equation. On the other hand, density and kinetic energy density, which are integrals of the distribution function, approach spatially homogeneous states strongly, which is accompanied by growth of the hydrodynamic entropy. Such a behavior can be seen when the Vlasov equation is reduced to the evolution equations for density and kinetic energy density by means of the Ehrenfest reduction.


Introduction
Dynamics of self-interacting gas is well described by the Vlasov equation, where m is the mass of one particle and f (t, r, p) is the one-particle distribution function on phase space and F is the force exerted on the gas from outside or by the gas itself. The force can be derived from a Hamiltonian (or energy of the system) as follows. For simplicity of notation, the dependence of the distribution function on time or on position and momentum will not always be written explicitly. Energy of the system with a self-consistent long-range interaction is E = dr dp p 2 2m f + 1 2 dr dp dr dp f (r, p)V(|r − r |) f (r , p ), (2) where V(|r − r |) is an interaction potential, e.g., electrostatic or gravitational potential. This is the mean-field approximation of the electric energy of the plasma, which is sufficient for discussing the Landau damping below. Derivative of energy with respect to the distribution functions is then δE δ f (r, p) = p 2 2m + dr dp V(|r − r |) f (r , p ).
F(t, r) = − ∂ ∂r δE δ f = − dr dp ∂V(|r − r |) ∂r f (r , p ). (4) Note that this dependence of the force on the distribution function makes the Vlasov equation nonlinear. All terms in the Vlasov Equation (1) have now been specified.
In plasma physics, Landau damping is usually shown for the Vlasov equation coupled with the Poisson equation for electrostatic potential. In three dimensions, the fundamental solution of Poisson equations (Green's function) is proportional to the Coulomb potential, which means that one can write down the solution to the Poisson equation and plug it back into the Vlasov equation. This way, the mean-field approximation is obtained. However, the situation is more complex for periodic domains as discussed in [1], but the mean-field approximation can still be used. Note also that the mean-field Vlasov equation describes also gravitational interaction and was used by Jeans to study gravitational collapse, which occurs on distances larger than the Jeans length (see [2]).
Landau damping is a feature of solutions of the Vlasov equation-solutions approach spatially homogeneous distributions. This was first derived by L. D. Landau [3,4] when considering the linearized Vlasov equation. Villani and Mouhot then proved it without the linearization, i.e., in the fully nonlinear setting (see [1,5], where conditions on the initial perturbation and decay rates can be found). More precisely, Landau damping can be seen as a development of fast oscillation in the p-dependence of f (t, r, p) and a simultaneous decay of the dependence on r. The distribution function approaches a spatially homogeneous, i.e., r-independent, distribution function only weakly (i.e., not pointwise but in the mean) due to the fast oscillations in the p-dependence. However, the approach of integral over the p-coordinate (particle density ρ(r)) is already strong, i.e., f (r, p) w →f (p) and ρ(r) = dpm f (r, p) →ρ = const.
Constant m is the mass of one particle. These two properties of solutions to the Vlasov equation are referred to as the Landau damping.
A physical explanation of Landau damping was given for instance in [6], where particles "surfing" on the electromagnetic waves either give energy to the wave or take energy from it (damping). Landau damping can be also studied by means of statistical mechanics via observing coarse-grained entropies [7,8]. It can be also observed in continuum simulations of the Vlasov equation and electromagnetic field, e.g., [9]. The emergence of dissipative phenomena from reversible equations is also observed when proposing closure relations in the Grad hierarchy, e.g., [10][11][12][13].
Villani [5] summarizes the results of the analytic computations in the linear (the second term on the R.H.S. of Equation (1) is linearized). Landau damping as follows. Under certain assumptions (the data are analytic and the Penrose stability condition holds), solutions to the linearized Vlasov equation have the following properties: (i) all nonzero spatial Fourier modes of ρ 1 (t, r) = h(t, r, p)dp (h stands for the perturbation of a spatially homogeneous f 0 distribution function), converge to 0, hence ρ 1 converges to a constant in space; (ii) the force also converges to 0; (iii) h(t, r, p)ϕ(r, p)drdp converges to h i (r , p)ϕ(r, p)drdr dp for any smooth test function ϕ. Here, h i (r, p) = h(0, r, p) stands for the initial value of the perturbation, and the extra integration comes from the convergence to the mean.
The nonlinear Landau damping enjoys the same long-time properties if the initial perturbation is a perturbation of a stable equilibrium and some other assumptions are satisfied-particularly that the interaction potential be no more singular than Newtonian (Fourier transposeŴ(k) = O(1/|k| 2 ) for k ∈ Z d ). The main two observations are again that (i) f (t, r, p) → f +∞ (p) weakly, which implies convergence of smooth observables, that is f (t, r, p)ϕ(r, p)drdp → f +∞ (p)ϕ(r , p)drdr dp = V f +∞ (p)ϕ(r , p)dr dp), and F(t, r) → 0 strongly with time.
Landau damping is also interesting from the thermodynamic point of view. Although it clearly described decay of the distribution function to homogeneous equilibria, the Boltzmann entropy (see e.g., [14]), is kept constant by the Vlasov equation, as can be verified by direct calculation. Actually, any functional of the form dr dpη( f ), η being a function from real numbers to real numbers, is conserved for any energy. Such functionals are called Casimirs of the Poisson bracket. This seems to be in contradiction with the usual thermodynamic interpretation of relaxation to some equilibrium, which is usually associated with growth of entropy. This seemingly paradoxical situation will be resolved in this paper by showing that the hydrodynamic entropy, which is more macroscopic than the Boltzmann entropy, grows while Boltzmann entropy is conserved, which seems to be in agreement with the the answer of Villani [15] when we asked him.
In their rigorous analysis of solutions to the Vlasov kinetic equation, Villani et al. have observed that the emergence of irregularities of solutions plays an analogical role in the kinetic theory as the emergence of ergodicity of particle trajectories in the particle dynamics (connected to the KAM theorem). This observation supports our analysis of solutions of the Vlasov kinetic equation that uses the methods of statistical mechanics. All methods of statistical mechanics (including in particular the Zwanzig-Mori approach [16,17]) are devised to extract important features and ignore the unimportant details. The important features emerge as patterns in the collection of trajectories (i.e., solutions) of the time evolution equations. The longer the time period for which we know the trajectories, the larger is the chance that we recognize the pattern. For this reason, we turn in this paper to the Ehrenfest method that follows the trajectories beyond the infinitesimal period shown in the vector field.
In this paper, we follow [18] (In [18], we investigate Landau damping from a more detailed level of description. In particular, extended kinetic theory with induced dissipation is shown to reproduce Landau damping. In the present work, however, we analyze Landau damping from a less detailed (hydrodynamic) level of description) and investigate the Vlasov Equation (1) by applying methods of statistical mechanics. This investigation offers us a possibility to see both the statistical methods and the Landau damping in a new perspective. In [18], we have first extended the Vlasov equation to include explicitly the micro-turbulence and then we have shown that its decay brings about the Landau damping. In this paper, we apply the methods of non-equilibrium statistical mechanics, namely the Maximum Entropy Principle (MaxEnt) and the Ehrenfest reduction method (see more in Section 2) directly to the Vlasov equation. We reduce it to two equations governing the time evolution of two scalar fields, ρ(r) = dpm f (r, p) and (r) = dp p 2 2m f (r, p). From an analysis of their solutions, we are then able to see the Landau damping as the approach of the field ρ(r) to a spatially homogeneous field (in Section 3).

Ehrenfest Reduction
A reduction of a dynamical system (DS1) to another dynamical system (DS2) is a pattern recognition in the phase portrait corresponding to DS1. The phase portrait corresponding to DS1 is a collection of the trajectories that arise in DS1 (i.e., solutions to the governing equations of DS1). The recognized pattern in the DS1 phase portrait is the phase portrait corresponding to DS2. The first step in the reduction is thus an information about the phase portrait. How shall we obtain such information? We recall two examples.
First, it is the Gibbs reduction of the Liouville equation to the equilibrium thermodynamics. The phase portrait is assumed to be ergodic and the recognized pattern is the Gibbs equilibrium distribution function obtained by maximizing the Gibbs entropy with constraints (MaxEnt principle). In the Gibbs analysis, the constraints are the state variables of the equilibrium thermodynamics (i.e., the total mass and the total energy). The Gibbs entropy is a measure of disorder. Since the constraints remain constant during the time evolution, the reduced time evolution is no time evolution.
The second example is the Ehrenfest reduction of DS1 to DS2. The constraints are now fields that serve as state variables in DS2. These two fields do not remain constant in the DS1 time evolution. The pattern is again obtained by the MaxEnt principle and its time evolution by following the DS1 time evolution in a small time interval. Below, we describe the Ehrenfest reduction in detail (in the rest of this section) and then we apply it on the Vlasov Equation (1) in Section 3.

General Formulation
Let us first recall the Ehrenfest reduction, which was developed in [19,20]. The starting point is to expand solutions to the DS1 time evolution equation (denoted by state variables x), in Taylor series in time (i.e., close to the initial condition x(t + τ)| τ=0 = x(t)), By the symbol x, we denote the state variables of the DS1 dynamics. Substituting Equation (7) into the expansion leads to where •, • stands for the usual L 2 -product, i.e., integration over the respective space or phase-space domain.
A less detailed level of description, i.e., DS2, has state variables where π, • is the projection operator, derivatives of which δπ a δx j will be denoted simply as π a i . For linear constant-in-time projection operators, the derivatives constitute a constant matrix, and the angular brackets have the same meaning, i.e., scalar product, as above.
Exact evolution of variables y can be obtained bẏ but, in order to evaluate it, it would be necessary to solve Equation (7). Since the aim is to obtain evolution equation for y only in terms of y, a reduction has to be carried out. In other words, some information has to be forgotten. The approximation can be sought by following two routes: (i) projecting the more detailed evolution to the less detailed state variables; (ii) expanding the sought evolution equation on the coarser level. Comparison of these two approaches yields an assessment of Taylor series expansion of the sought evolution equation on the coarser level (to arbitrary order).
In the latter route, we search for the unknown operator φ that governs the evolution on the coarser level DS2ẏ = φ(y).
Expanding its solution around an arbitrary initial condition y(t) yields To identify the asymptotic expansion of the unknown operator φ, we write Now, we compare this coarser DS2 level expansion to the DS1 time evolution. In order to achieve that, however, the DS1 time evolution, has to be made dependent on y. We choose the mapping from y to x to be the MaxEnt mappingx(y), i.e., to find x such that the detailed entropy S(x) is maximal subject to the constraints represented by y = π(x). Then, by substituting (13) into (12)) with projection of Equation (9), we obtain π, x(t + τ) = π, x(t) + τ π, J(x(y)) + τ 2 2 π, δJ(x(t)) δx where the quasi-equilibriumx(y) initial conditions were chosen so that only y appears in the expansion. The comparison leads to and which completes the estimate of the sought evolution equation on the coarser level via (13) Note that the τ 0 term has to be strictly satisfied due to the requirement that both the more and less detailed evolution coincide at τ = 0 on the quasi-equilibrium. The zeroth correction R (0) k can be also naturally understood as the least biased way to express x in terms of y is the MaxEnt [21] mapping from y to x,x(y), which is exactly what the first approximation is. The first correction R (1) k is typically non-zero as the detailed evolution (7) carries x from the quasi-equilibriumx(y) to values that are not in the image of the MaxEnt mapping, i.e., out of the quasi-equilibrium (or Legendre) submanifold, and corrections have to be introduced.

Hamiltonian Version of Ehrenfest Reduction
The purpose of this section is to reformulate the Ehrenfest reduction in the case when the detailed evolution DS1 is Hamiltonian.

Hamiltonian Structure of the Vlasov Equation
Vlasov Equation (1) is a Hamiltonian evolution, since it is generated by the Boltzmann-Poisson bracket which is of course antisymmetric and fulfills both the Leibniz rule and Jacobi identity. The Hamiltonian structure of plasma evolution equations was first found by Phil Morrison and John Greene in [22]. See also [23] or [24]. By rewriting the bracket as dr dp the terms rhs are the right-hand side of the evolution equation for f , which is the Vlasov Equation (1). This Poisson bracket can be derived from the Liouville Poisson bracket by projection or it can be seen as the Lie-Poisson bracket corresponding to symplectic transformations on the one-particle cotangent bundle (see, e.g., [25]). Any real-valued function of the distribution function σ : R → R then generates a Casimir of the bracket, This means in particular thatṠ and Boltzmann entropy is thus indeed conserved by the Vlasov equation. The Boltzmann entropy is a plausible entropy (a Casimir of the Boltzmann Poisson bracket), but it is not the only one possible. Nevertheless, in the case of ideal gases, it follows by MaxEnt from the Liouville entropy (see e.g., [14]), which can be seen as the phase-space analogue of the Shannon entropy −k B ∑ i p i ln p i . Shannon entropy is uniquely determined by Shannon's axioms defining uncertainty [26]. In this sense, the Boltzmann entropy S (B) is unique. Moreover, energy is conserved automatically due to the antisymmetry of the Poisson bracket. Each Poisson bracket can be seen as constructed by means of a Poisson bivector L, Antisymmetry of the Poisson bracket is expressed by antisymmetry of the Poisson bivector, and Jacobi identity for the bracket is also inherited from a corresponding formula for the bivector (see e.g., [27] p. 332). Casimirs of the Poisson bracket are functionals that are not evolving regardless of the choice of energy, i.e.,Ṡ = {S, E} = 0 ∀E, (24) which means that This relation is satisfied for all Casimirs. In particular, entropy is always taken as a Casimir of the Poisson bracket in order to be conserved by the Hamiltonian evolution, which means that the above relation holds for S being entropy. More generally, this allows a clear separation of reversible and irreversible part of evolution dynamics.

Formal Solution of Hamiltonian Evolution
The first step in Ehrenfest reduction was to expand the evolution of detailed variables x as a Taylor series in time. Assume that evolution of x is purely Hamiltonian, which means thaṫ Solution to the Hamiltonian evolution equations can be expanded in series as (see e.g., [27] pp. 334) When expanding the Poisson brackets by means of the Poisson bivector and when taking into account only linear functionals A, i.e., A(x) = A i x i , this last equation becomes which is the analogue of Equation (9). Taylor expansion (28) can be interpreted as a new evolution equation for x i where the first term on the right-hand side generates reversible evolution while the second irreversible.
Such an evolution equation can be regarded as a self-regularized evolution of x i corresponding to time step τ. A similar idea has already been brought up in symplectic integration [28]. The original Hamiltonian evolution, Equation (26), conserved both energy (due to antisymmetry) and entropy (due to degeneracy of the Poisson bracket). None of these properties is kept in the self-regularized evolution Equation (29). Indeed, evolution of entropy becomes The inequality follows from concavity of S(x) (and thus negative semi-definiteness of d 2 S). Entropy is thus produced by the self-regularized equations. Similarly, evolution of energy becomes due to convexity of energy (d 2 E positive semi-definite). Energy is thus dissipated by the selfregularized evolution.

Projection of the Poisson Bracket
Assume now that when functionals dependent only on the less detailed variables y are plugged into the Poisson bracket, the resulting expression depends only on variables y and it is referred to as the less-detailed Poisson bracket ↓ {•, •}, i.e., That is the case for example when y are the hydrodynamic fields of density, momentum density and entropy density and x the one-particle distribution functions. Poisson brackets can often be derived from more detailed Poisson brackets by such projections (see, e.g., [25]).
The projections are naturally made in the so-called energetic representation (see [29]), where, for example, in the case of hydrodynamic field entropy density is among the state variables instead of energy density. On the other hand, in the entropic representation, the field of energy density is among the state variables. In that case, it is natural to employ the MaxEnt principle to find the least biased estimate of the detailed variable x based on the knowledge of the less detailed variables y. In the case of the energetic representation, the MaxEnt principle is replaced by the principle of minimal energy: find x such that energy E(x) is minimal subject to the knowledge of field y that include the less detailed entropy density. This can be expressed by Legendre transformation which gives the dependencex(y * ) and a new functional the conjugate lower energy, which is a functional of conjugate less detailed (lower) variables y * . The energy on the lower level is then obtained by a subsequent Legendre transformation and These equations imply, in particular, that From Equation (33), it then follows that which will be useful later.
Assuming that we have already projected the Poisson bracket to the lower level of description, evolution of the less detailed state variables is given bẏ The Taylor expansion in time on the less-detailed level of description then yields A a y a (t + τ) = A a y a (t) This last equation is the analogue of Equation (12), and can be regarded as evolution equation where the right-hand side consists of a reversible (first) and an irreversible (second) term. This is again a self-regularization of evolution (39).

Comparing the Solutions on Different Levels
The next step in the Ehrenfest reduction is to alter the right-hand side of the less-detailed evolution equation so that the solution is closer to the solution on the detailed level of description.
More precisely, Taylor expansion (28) gives an approximation of x i (t + τ) provided x i (t) is known. Similarly, Taylor expansion (40) gives an approximation of y a (t + τ) provided y a (t) is known. The value y a (t + τ) can be, however, approached by two ways: (i) Projection from x(t) to y(t) and subsequent evolution by Taylor expansion (40) or (ii) evolution by expansion (28) and subsequent projection from x(t + τ) to y(t + τ). The second route is of course more precise, but one has to solve the detailed evolution equations for x. Therefore, the first route is more suitable, but it needs a further correction to make it closer to the second route. An additional term must thus be added to expansion (40), which makes the two routes equivalent up to the given order of τ. This procedure is summarized in Figure 1. Step 2: This solution, x(t + τ), is then projected to the less detailed level to obtain π(x(t + τ)).
Step 1 : Alternative route is to first project x(t) to y(t).
Step 2 : The less detailed evolution equation (generated by the projection of the Poisson bracket) then takes y(t) and gives y(t + τ).
Step 2 : We have thus π(x(t + τ)) and y(t + τ), which should ideally be equal, but they are typically not. The value π(x(t + τ)) is of course more precise because it is constructed from the detailed evolution equations. To make the value y(t + τ) more precise, the less detailed evolution equations are altered by adding the difference between the self-regularized detailed and less-detailed equations. Such equations for y are then the reduced Ehrenfest equations. Let us thus compare the self-regularized evolutions (28) and (40). The first terms of the expansions are already related by the projection because (using Equation (38)) The second terms of the expansions are no longer related by such a projection. Therefore, one has to add the difference between these two second terms. The regularized evolution equation for y then becomesẏ which is the less-detailed evolution equation obtained by the Ehrenfest reduction, analogue of Equation (18) The canonical Poisson bivector is a constant matrix.
The Poisson bracket for hydrodynamic density and momentum density (ρ, u) is not a canonical one. However, it can be casted into the canonical form by transformation from (ρ, u) into the so-called Clebsch variables (see [30] and [31], or even (ρ, u, s) as in [32]). Constant Poisson bivectors can be thus often met.
If the bivector is constant, expansions of the solutions (28) and (40) obtain a particularly simple form, which gives evolution equationsẋ Matrices M and ↓ M are clearly symmetric. Moreover, it follows from convexity of energy, the second differential of which is thus a positive definite matrix, that which is the analogue of Equation (18).

Canonical Hamiltonian System
For instance, the canonical Poisson bivector (44) is constant. The canonical Poisson bracket, which provides natural kinematics on cotangent bundles, is formed from the Poisson bivector.
Let the variables on a cotangent bundle be denoted by (r, p), which can be interpreted as position and momentum of a particle. Evolution of these variables is then given by with H(r, p) being energy of the particle in a static external field, e.g., The self-regularized evolution of (r, p) is then, according to Equation (47), For the particular choice of Hamiltonian (52), the evolution equations become This is the self-regularized evolution of the particle. Assuming that the external field V(r) is convex and that it has a minimum, the Hamiltonian (or energy) is also convex. Consequently, the dissipative matrix M is positive-semidefinite, and energy is dissipated in course of the evolution aṡ and the evolution stops at the position and momentum satisfying V r = 0 and which is the equilibrium position of the particle within the field as well as the position of lowest total energy. Let us now consider a less detailed level of description where only the position r plays the role of state variable. The mapping from (r, p) to (r) is the projection between the two levels of description. The Poisson bivector on the less detailed description is a 1 × 1 matrix, i.e., a number, and it is given by The less detailed Poisson bivector is identically zero. The less detailed evolution is thus zero as well as the self-regularized less detailed evolution.
However, the less detailed evolution derived by the Ehrenfest reduction, Equation (50), is not zero and readsṙ This reduced evolution tends again to the minimum of the external field. Lower-level energy ↓ E(r) = H(r,p(r)) = V(r) is being damped by the reduced evolution equation until the state with minimal energy is reached. In order to restore energy conservation, another state variable (entropy or energy) should be taken into account so that mechanical energy can be transformed into thermal (see e.g., [14]).

From Vlasov to Mechanical Equilibrium
Now, we are coming to the main point of the paper-the reduction of the Vlasov equation to mechanical equilibrium. Mechanical equilibrium is the level of description where only density and kinetic energy density play the role of state variables. Since those two fields can still evolve in time, mechanical equilibrium is still a generally non-equilibrium description. The reason for the name is that momentum density is already in a quasi-equilibrium state given by the acting mechanical forces.
The reduced equations will naturally exhibit the tendency to spatially homogeneous equilibrium, which can be interpreted as a manifestation of Landau damping. However, let us first define the projection from the Vlasov level of description to the level of mechanical equilibrium.

Projection
Since the particle density approaches the equilibrium in a strong sense, we can expect that if entropy were a function of the density, it could approach its maximum. Therefore, it is sensible to regard the Vlasov equation from a less detailed (or lower) level of description-for example, the level of mechanical equilibrium where state variables are particle density: kinetic energy density: ε(r) = dp The level of description where ρ and ε play the role of state variables is referred to as the level of mechanical equilibrium (motivated by [33]) because momentum is not present among the state variables.
Let us justify the reduction to the level of mechanical equilibrium from the perspective of rigorous mathematical results summarized in Section 1. Consider a particular choice of the test function that is a tensor product ϕ(r, p) = p ⊗φ(r), whereφ is smooth. Then, we have that h(t, r, p)ϕ(r, p)drdp = drφ(r) dp ph(t, r, p) → drφ(r) dr dp ph i (r , p), where the extra integration with respect to r gives the average towards which the perturbation converges. As this relation holds for all smoothφ and particularly for smooth cutoff functions (convolution of a characteristic function with a mollifier), we can observe that lim t→+∞ dv vh(t, x, v) = dydv vh i (y, v), which is a constant in space. Therefore, not only the macroscopic density, but also the macroscopic momentum, tends, as t → ∞, to a spatially homogeneous state. Hence, it is natural to consider the mechanical equilibrium as a reasonable level of description for long time behavior of the system. We shall explore this Landau damping via pertinent Ehrenfest reduction that is able to reveal irreversible behavior on the reduced level of description. Total energy (2) can be expressed as Mass and kinetic energy density are thus obtained from the distribution function by projection (59) and (60).
On the other hand, knowledge of fields ρ and ε is insufficient for exact reconstruction of the distribution function. What is the estimate of the distribution function when only those two fields are accessible? The answer is the maximum entropy principle (MaxEnt), which reveals that the least biased estimate of the distribution function is the solution to equation where fields ρ * and ε * play the role of Lagrange multipliers ensuring constraints (59) and (60), but can be also given geometrical interpretation [14]. Solving this equation leads tõ where the h −3 prefactor comes ensures the correct units of the distribution function. Using constraints (59) and (60), the Lagrange multipliers can be expressed as functions of ρ and ε, which finally leads tõ which is the MaxEnt estimate of the distribution function subject to the knowledge of field ρ and ε.
Boltzmann entropy evaluated at the MaxEnt distribution function then becomes the entropy on the level of mechanical equilibrium, which is in fact the local equilibrium version of the Sackur-Tetrode equation [29].
Having constructed the projection from f onto the level of mechanical equilibrium, a mapping from mechanical equilibrium back into the space of distribution functions and the implied mechanical-equilibrium entropy, we have specified the static relations between the two levels of description (Vlasov and mechanical equilibrium). Let us now have a look at dynamical relations between the levels, i.e., how evolution on one level corresponds with evolution on the another level.

Construction of the Reduced Evolution
Let us now apply the Ehrenfest reduction to the passage from the Vlasov level of description onto the level of mechanical equilibrium. State variables are identified as and y = (ρ(r), ε(r)).
The more detailed evolution equation is the Vlasov Equation (1). The sought reduced evolution equations are evolution equations for ρ and ε. The zeroth approximation is zero, which readily follows from the observation that wheref (ρ, ε; p) is the quasi-equilibrium state calculated in (65). The zeroth approximation of the macroscopic density then follows from the corresponding projection of the last expression for J(x(y)), which is The first term vanishes as it is an odd function in all p k and the second term can be simply integrated to zero as F(t, r) is independent on momentum p. Analogously, one can show that the projection yielding kinetic energy density ε vanishes as well on J(x(y)) and hence the zeroth approximation R (0) ε is zero. The first approximation that corresponds to irreversible evolution can again be calculated from the relations (17) identified above. As the zeroth approximation is zero, i.e., R (0) = π, J(x(y)) = 0, only the first terms contribute to the first order correction. Let us start by calculating the variation of where the first integral on the last but one line vanishes unless k = j due to symmetry off (an even function in p i ). Its nonzero value is 2εm/3 as can be seen from the projections. Note that ∆ stands for the Laplacian.
In exactly the same manner, one can proceed in obtaining the first (irreversible) approximation of the macroscopic evolution equation for kinetic energy density. Finally, the reduced evolution equations become ∂ρ ∂t ∂ε ∂t where the force in mechanical equilibrium F (ME) is defined as Equations (73) and (74) are the evolution equations approximating the Vlasov equation on the level of mechanical equilibrium, and they have been derived by the Ehrenfest reduction. The equations are clearly irreversible in the sense of time-reversal transformation [34].
In summary, the reduced evolution on the level of mechanical equilibrium is irreversible although the detailed evolution (Vlasov Equation (1)) was completely reversible, and entropy on the level of mechanical equilibrium is expected to grow in time although the Boltzmann entropy is conserved by the Vlasov equation. Let us further analyze the evolution equations in the next section.

Conservation of Total Energy
Evolution of the total energy on the level of mechanical equilibrium can be obtained by using the reduced evolution Equations (73) and (74) througḣ Total energy is thus conserved by the reduced evolution equations.

Dissipativity of Reduced Evolution
To assess the growth of the reduced entropy, we evaluate the following expressioṅ and obtaiṅ The first term on the right-hand side is always non-negative. Using the Young inequality |a · b| ≤ (a 2 + b 2 )/2, the remaining three terms on the first line can be estimated from below as which means that the entropy production is bounded from below bẏ Unlike the first two positive terms, the last term has no definite sign, which means that it has to be estimated using the positive terms. By Hölder inequality | dr f g| ≤ dr f 2 drg 2 and the Young inequality, the last term can be bounded from above as where c is an arbitrary constant. Setting c 2 /2 = 3/4, i.e., 1/2c 2 = 1/3, the last expression becomes equal to the two positive terms on the right-hand side of Equation (81), and we obtain eventually thaṫ S (ME) ≥ 0. The reduced entropy production is always non-negative, and the reduced evolution equations are thus compatible with the second law of thermodynamics.

Homogeneous Equilibrium Solution
A natural question now arises whether there exists a (spatially) homogeneous solution of the macroscopic evolution Equations (73) and (74) describing the system in mechanical equilibrium, i.e., a solution to Hence, one immediately arrives at a necessary condition for the existence of a homogeneous solution that states that F ME = 0, perfectly corresponding to the observation of Villani et al. that F ME strongly tends to zero with long times (as perturbations decay). With the knowledge of the particular form of the quasi-equilibriumf and with the assumption of homogeneous ρ, ε we may proceed further with to observe that a homogeneous solution exists if or Ω is periodic as a system (e.g., a torus) The choice of the homogeneous solution is not, however, arbitrary. It follows from the initial condition and the evolution equation accompanied by entropy functional guiding this process. To prove mathematically that the reduced evolution equations tend to this homogeneous solution is beyond the scope of this paper (if it can be done at all: strongly coupled nonlinear parabolic PDEs). However, the approach of ρ, precise evolution of which is given by the Vlasov equation, to a homogeneous equilibrium was proven by Villani and Mouhot as a manifestation of nonlinear Landau damping. The approach of ρ to a homogeneous equilibrium by means of the reduced evolution Equation (73) and (74) can be anticipated.
Is the homogeneous solution the state where a maximum of the mechanical equilibrium entropy is reached? The system of evolution equations is accompanied by entropy functional (66) indicating the evolution of mechanical equilibrium variables ρ and ε towards the equilibrium values. Maximization of the reduced entropy with respect to constraints can be written as These equations lead to Therefore, there is a relation between particle density and kinetic energy density (in the homogeneous equilibrium) via a constant α and consequently total energy becomes for homogeneous distribution of particles where ν = drV(|r|),ρ = M/|Ω|,ε = E/|Ω| and Ω is the spatial integration domain. The MaxEnt problem has thus a solution homogeneous in space, and the solution has been determined explicitly (for known total energy and mass). Is the solution unique?
To proceed further and comply with the definitions of energies above, we need to restrict ourselves to finite domains. Hence, we are forced further to assume that Ω is a finite domain, e.g., a torus, in line with the assumptions in the work of Mouhot and Villani [1]. This identifies constant α as function of E, M, Ω and of the potential V. Moreover, Equation (91) gives the Lagrange multiplier M * as a function of E, M and Ω as Equation (91) thus becomes which is a nonlinear integral equation for ρ(r). Equation (95) can be linearized aroundρ throughρ := ρ ρ − 1 1 in magnitude, which yields which is a linear Fredholm integral equation. By applying Fourier transform to the equation, we obtain thatρ whereρ(ξ) = drρ(r) exp(−2πir · ξ). The transformed equation means thatρ = 0 and thusρ has to be constant. From the definition ofρ, it follows that it must be zero everywhere, which means that the only solution to the linearized problem is the homogeneous in space solution.
We have shown using physical (not mathematical) tools that there is a (at least to some extent unique) homogeneous solution to the macroscopic problem (83) and (84), which is a MaxEnt equilibrium (maximizes the reduced entropy). Moreover, entropy production is non-negative until this homogeneous density distribution is reached. Therefore, this approach is providing thermodynamic arguments for linear Landau damping.

Some Qualitative Insight into Macroscopic Evolution Equations: Linearization.
It was shown in the preceding section that homogeneous equilibrium is a state ("the" state when near to the equilibrium) for which the mechanical equilibrium entropy S (ME) attains its maximum. The approach of ρ and ε by means of evolution Equations (73) and (74) towards the homogeneous equilibrium can be thus expected.
We can gain some analytical insight into the governing macroscopic evolution equations from their analytical solution, which is available only under certain assumptions. We shall thus solve a linearized version of the equations around constant solutions.
The macroscopic evolution Equations (73) and (74) can posses a constant (both in time and space) solution¯ ,ρ only when F (ME) = 0 =ρ dr ∂V(|r −r |) ∂r j , which requires the system Ω to be periodic as discussed above. Then, actually arbitrary constant¯ ,ρ is a solution, but the thermodynamically admissible single one is chosen by the entropy functional as described above.
Denoting small perturbations around this constant solution with tildes, we may approximate the evolution equations to the leading (linear) order as follows: (F (ME) stands forF Fourier transform in space from (ρ(r),ε(r)) to (ρ(ξ),ε(ξ)) turns the evolution equations to Taking a radially symmetric potential of form V(|r|) = Kr −k , K being a constant, Fourier transform of the potential is d being the dimension. Assuming 0 < k < d, the Gamma functions are always positive, which means that the Fourier transform is real and positive for K > 0 (electrostatic case) while negative for K < 0 (gravitational case). In the electrostatic case, the first principal minor of the matrix M is positive and determinant of the matrix, det M = 80 27 is also positive. The matrix is then positive definite, and the linearized system of ordinary differential Equations (100) for each Fourier mode has a solution decaying exponentially fast to zero. Therefore, the potential (or electrostatic) energy also decays exponentially. The potential energy is the L 1 norm of the convolution of the density, potential and density. By Young-Hölder inequality, it is lower than the L 1 norm of the potential multiplied by a square of the L 2 norm of density. The latter is equal to the square of the L 2 norm of the Fourier image of density, which decays exponentially in time. Therefore, potential energy also decays exponentially in time. In the gravitational case, however, the definiteness of the matrix is lost for long enough Fourier modes as in the Jeans instability.
Due to linearity of Fourier transform, the decay of perturbations occurs iff F −1 e −τπ 2 ξ 2 λ ± t (t, r) decay with time, where we denoted λ ± the eigenvalues of matrix M. With known potential V, we may readily assess this condition. To illustrate this, we further assume that the dependence of λ ± onV(ξ) can be dropped, i.e.,ε ρ 2V (ξ), and hence λ ± ξ 2 ∝ −ξ 2 . Thence, ε, ρ are a linear combination of two fundamental solutions to the diffusion operator. Particularly, we have withλ ± = 4ε 9ρ (5 ± √ 10). Therefore, the initial delta pulse (at the origin) is flattening out with time as (κtτ) −3/2 , which can be used for identification of the essential parameter in the Ehrenfest reduction, the time-scale τ. Experiments measuring the Landau damping report characteristic times of these decays, hence assuming exponential decay. Note, however, that C 1 t −3/2 + C 2 is practically indistinguishable for exponential decay over finite time intervals and real data. Finally, note that one should convolve the initial condition with the above identified responses to delta pulses, but we argue that such a response is already providing a qualitative insight into long-time behavior (decay).
In summary, linearization of Equations (73) and (74) around homogeneous solutions leads to the conclusion that the equations approach the homogeneous solutions exponentially, which is compatible with maximization of entropy at the homogeneous solutions and Landau damping. Landau damping can thus be manifestly visible on the less detailed level of mechanical equilibrium.

Conclusions
Landau damping is the property of the distribution function, governed by reversible Vlasov equations, that approaches spatially homogeneous equilibria in a weak sense (i.e., in the mean). Particle density (integral of the distribution function) then approaches the homogeneous equilibrium strongly. The reversible Vlasov equation does not alter the Boltzmann entropy, which depends on the distribution function. Conservation of entropy and the approach to the homogeneous equilibrium (Landau damping) are thus seemingly in contrast.
Here, we have shown (by means of the Ehrenfest reduction) that the Vlasov equation can be approximated by evolution equations for particle density ρ and kinetic energy density ε, (73) and (74). These state variables are equipped with their own entropy, (66), which is given by maximization of the Boltzmann entropy. The evolution equations approach homogeneous equilibrium, where the reduced entropy is maximized. It is thus the reduced hydrodynamic entropy (66) that is maximized during the Landau damping, and Landau damping has been given an alternative thermodynamic interpretation. The paradox of thermodynamics of Landau damping is resolved this way.
Author Contributions: M.G. initiated the research. All authors developed the concepts and methods used. M.P. and V.K. carried out the calculations wrote most of the paper.