*Article* **Landau Damping of Langmuir Waves: An Alternative Derivation**

## **Andreas Shalchi**

Department of Physics and Astronomy, University of Manitoba, Winnipeg, MB R3T 2N2, Canada; Andreas.Shalchi@umanitoba.ca

**Abstract:** In this paper, a discussion of the Landau damping of Langmuir waves is presented together with a simple derivation which does not require the application of methods of complex analysis. A general dispersion relation is derived systematically which corresponds to a nonlinear equation. The latter equation is solved numerically but asymptotic limits are also discussed.

**Keywords:** Langmuir waves; Landau damping; particle wave interactions

## **1. Introduction**

In the current paper, the damping of plasma waves, also known as Langmuir waves [1], is explored. Plasma waves are related to electron oscillations and, thus, they are sometimes called plasma oscillations. The standard description of those waves is based on the momentum equation (Euler equation) for the electrons, their continuity equation, and Gauss' law of electrodynamics. Furthermore, it is assumed that the electrons behave like an ideal gas. In order to find the dispersion relation of Langmuir waves, the aforementioned fluid equations are linearized leading to the Bohm–Gross dispersion relation (see [2–4]),

$$
\omega^2 = \omega\_p^2 + 3v\_e^2 k^2. \tag{1}
$$

Therein, the wave number *k*, the plasma frequency

$$
\omega\_p = \sqrt{\frac{n\_0 e^2}{m\_e \epsilon\_0}}\tag{2}
$$

and the electron thermal speed

$$w\_c = \sqrt{\frac{k\_B T}{m\_c}}\tag{3}$$

are used in terms of the electron particle density (particles per volume), *n*0, the elementary charge, *e*, the electron mass, *me*, the permittivity of free space, 0, Boltzmann's constant, *kB*, and the absolute temperature, *T* (in Kelvin). Note, the electron thermal speed comes straight from the equipartition theorem which states that each degree of freedom has the average kinetic energy *kBT*/2. Therefore,

$$
\frac{1}{2}m\_{\mathfrak{c}}v\_{\mathfrak{c}}^2 = \frac{1}{2}k\_B T\_{\mathfrak{c}} \tag{4}
$$

is a consequence of having only one degree of freedom due to the one-dimensional motion of the electrons. This relation can be straight rearranged to obtain Equation (3).

A useful parameter in the description of plasma waves is the Debye length corresponding to the distance over which significant charge separation occurs; it is given by

**Citation:** Shalchi, A. Landau Damping of Langmuir Waves: An Alternative Derivation. *Physics* **2021**, *3*, 940–954. https://doi.org/10.3390/ physics3040059

Received: 31 August 2021 Accepted: 8 October 2021 Published: 29 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

$$
\lambda\_D = \sqrt{\frac{k\_B T \epsilon\_0}{n\_0 \epsilon^2}}.\tag{5}
$$

Then, the electron thermal speed (3) reads:

$$
v\_{\mathfrak{e}} = \omega\_p \lambda\_{D\_{\mathfrak{e}}} \tag{6}$$

and the plasma wave dispersion relation (1) becomes:

$$
\omega^2 = \omega\_p^2 \left( 1 + 3\lambda\_D^2 k^2 \right). \tag{7}
$$

By considering linearized fluid and Maxwell equations, one can show that Langmuir waves are longitudinal in nature.

In this paper, Landau damping, describing the damping of plasma waves due to the interaction between the wave and the electrons in the plasma, is revisited. This is a special case of wave–particle interactions in a plasma. According to [5], the physical mechanism of Landau damping can be understood as follows: particles having velocities slightly less than the phase velocity of the wave are accelerated by the wave's electric field to move with the wave phase velocity. Therefore, the group of particles moving slightly slower than the phase velocity gain energy from the wave. In a collisionless plasma characterized by a Maxwellian distribution (see Section 3), the number of slower particles is greater than the number of faster particles. Therefore, energy gained from the waves by slower particles is more than the energy given to the waves by faster particles, thus leading to net damping of the waves. An example plot is given in Figure 1.

**Figure 1.** Maxwellian distribution. The arrow points to the location of the phase speed *ω*/|*k*| of the Langmuir wave. In the considered case, the wave speed is larger than the electron thermal speed. Particles to the left of the arrow are slower, and particles to the right of the arrow are faster than the wave. Here we have used the wave number *k*, the frequency *ω*, the electron speed *v*, and the electron thermal speed *ve*.

The first theoretical description of Landau damping was presented in [6] and was based on the Vlasov equation [7]. Although Landau damping was originally derived for plasma waves, it can also be considered in the context of magnetohydrodynamic waves [8].

The derivation of traditional Landau damping, as provided in [6], is based on Laplace transforms and complex contour integrations. In [9], it was shown that Landau damping can also be described by using a Fourier transform and normal mode expansion. The latter modes are usually called van Kampen modes. In [10], it was demonstrated that Landau's and van Kampen's treatments of plasma oscillations are equivalent. Alternative and simpler derivations of Landau damping can be found in textbooks, see, e.g., [11]. However, all the aforementioned descriptions are based on contour integrals and other tools of complex analysis.

In order to avoid using tools of complex analysis, an alternative derivation of Landau damping is presented in the current article. This approach can be particularly useful for teaching this topic to undergraduate and graduate students in introductory plasma physics and magnetohydrodynamics courses.

The reminder of this paper is organized as follows. Section 2 lists the basic equations used to describe Landau damping. This includes the Vlasov equation, Gauss' law, and linearizing those equations. Section 3 contains the derivation of a nonlinear equation for the dispersion relation, *ω* = *ω*(*k*). Numerical and analytical solutions of this nonlinear relation are discussed in Sections 4 and 5, respectively. Finally, Section 6 discusses and summarizes the findings.

#### **2. Basic Equations**

In order to describe the statistical behavior of a physical system not in a state of equilibrium the Boltzmann equation is used. Let us assume that the plasma electrons are described via the distribution function *f*(*x*,*v*, *t*) which depends on position, *x*, and velocity, *v*, of the particles, as well as the time, *t*. Furthermore, the considered particles experience an external force, *F*, and collisions. With the help of Liouville's theorem [12], one finds:

$$\frac{d}{dt}f\left(\vec{x},\vec{v},t\right) = S\left(\vec{x},\vec{v},t\right),\tag{8}$$

where *S x*,*v*, *t* describes the collisions. For the total time-derivative of the distribution function,

$$\frac{d}{dt}f\left(\vec{x},\vec{v},t\right) = \frac{\partial f}{\partial t} + \sum\_{n} \dot{\mathfrak{x}}\_{n} \frac{\partial f}{\partial \mathfrak{x}\_{n}} + \sum\_{n} \psi\_{n} \frac{\partial f}{\partial \upsilon\_{n}},\tag{9}$$

can be used. Here, the dot defines time derivative.

Equations (8) and (9) can be combined to find the Boltzmann equation. The collisionless Boltzmann equation, on the other hand, is given by

$$\frac{\partial f}{\partial t} + \sum\_{n} \dot{\mathbf{x}}\_{n} \frac{\partial f}{\partial \mathbf{x}\_{n}} + \sum\_{n} \dot{\upsilon}\_{n} \frac{\partial f}{\partial \upsilon\_{n}} = 0. \tag{10}$$

In what follows, only this collisionless case is considered. Furthermore, *v*˙*<sup>n</sup>* is replaced by using Newton's second law, *Fn* = *mv*˙*n*, to find:

$$\frac{\partial f}{\partial t} + \sum\_{n} v\_{n} \frac{\partial f}{\partial x\_{n}} + \frac{1}{m} \sum\_{n} F\_{n} \frac{\partial f}{\partial v\_{n}} = 0,\tag{11}$$

with *m* being the particle mass.

For the case considered in this paper, where *Fn* describes Coulomb interactions, Equation (11) is usually called the Vlasov equation [7]. The Vlasov equation, discussed here, is commonly used to investigate the interaction between particles and fields. In the case of Landau damping, this equation describes the interaction between Langmuir waves and the electrons of the plasma. However, the Vlasov equation is also used as a starting point to describe the interaction between magnetic turbulence and energetic particles such as cosmic rays (see, e.g., [13–15]).

For the investigations, presented here, it is sufficient to consider the one-dimensional case for which the Vlasov equation reads:

$$\frac{\partial f}{\partial t} + v \frac{\partial f}{\partial x} + \frac{1}{m} F \frac{\partial f}{\partial v} = 0. \tag{12}$$

In the case considered, the force is just the electric force and, thus,

$$F = qE = q\delta E,\tag{13}$$

can be used, where *q* is the electric charge of the particle. The electric field, *δE*, therein is the wave field of the plasma wave. Using this in Equation (12) yields:

$$\frac{\partial f}{\partial t} + v \frac{\partial f}{\partial x} + \frac{q}{m} \delta E \frac{\partial f}{\partial v} = 0. \tag{14}$$

In order to solve the Vlasov equation, a perturbation approach of the form,

$$f(\mathbf{x}, v, t) = f\_0(v) + \delta f(\mathbf{x}, v, t), \tag{15}$$

is employed.

Here, *f*<sup>0</sup> is the unperturbed distribution often assumed to be Maxwellian (see Section 3 below). The quantity *δ f* describes the fluctuations corresponding to the deviations from the unperturbed distributions. These fluctuations are a consequence of the interactions with the electric field of the wave and are assumed to be small so that the perturbation approach can be justified. Using Equation (15) in Equation (14) yields in first-order:

$$\frac{\partial \delta f}{\partial t} + v \frac{\partial \delta f}{\partial x} + \frac{q}{m} \delta E \frac{\partial f\_0}{\partial v} = 0. \tag{16}$$

Equation (16) is called the linearized Vlasov equation. As a second equation, Gauss' law,

$$
\vec{\nabla} \cdot \vec{E} = \frac{1}{\epsilon\_0} \rho\_{\varepsilon \prime} \tag{17}
$$

is employed.

Note that SI units are used here rather than Gaussian or cgs units. In the onedimensional case considered, Equation (17) becomes:

$$\frac{\partial}{\partial \mathbf{x}} \delta E = \frac{1}{\epsilon\_0} \delta \rho\_\varepsilon = \frac{q}{\epsilon\_0} \int\_{-\infty}^{+\infty} dv \,\delta f(\mathbf{x}, v, t),\tag{18}$$

where the electric charge density *δρ<sup>e</sup>* is replaced by the velocity integral over the fluctuations *δ f* . In order to solve Vlasov and Gauss equations, a Fourier ansatz,

$$\begin{array}{llll}\delta f\left(\mathbf{x},\boldsymbol{\upsilon},t\right) & \approx & \delta f\left(\mathbf{k},\boldsymbol{\upsilon},\omega\right)e^{i\left(k\mathbf{x}-\omega t\right)},\\\delta E\left(\mathbf{x},t\right) & \approx & \delta E\left(\mathbf{k},\omega\right)e^{i\left(k\mathbf{x}-\omega t\right)},\end{array} \tag{19}$$

is used.

Using Equation (19) in Equations (16) and (18) yields:

$$-i\omega \delta f + ivk \delta f + \frac{q}{m} \delta E \frac{\partial f\_0}{\partial v} = 0,\tag{20}$$

and

$$ik\delta E = \frac{q}{\epsilon\_0} \int\_{-\infty}^{+\infty} dv \,\delta f(k, v, \omega). \tag{21}$$

The first equation can be rearranged to read:

$$
\delta f = -i \frac{q}{m} \frac{\delta E}{\omega - vk} \frac{\partial f\_0}{\partial v}. \tag{22}
$$

Using this in Equation (21) yields:

$$ik\delta E = -i\frac{q^2}{\epsilon\_0 m} \int\_{-\infty}^{+\infty} dv \, \frac{\delta E}{\omega - vk} \frac{\partial f\_0}{\partial v}.\tag{23}$$

As the electric field does not depend on the particle velocity *v*, one can cancel *δE* in the latter equation. Therefore, the following dispersion relation,

$$1 = -\frac{q^2}{\epsilon\_0 m k} \int\_{-\infty}^{+\infty} dv \, \frac{1}{\omega - vk} \frac{\partial f\_0}{\partial v},\tag{24}$$

is obtained.

The function *f*<sup>0</sup> therein corresponds to the unperturbed particle density per volume. In the following, the replacement *f*<sup>0</sup> → *f*0*n*<sup>0</sup> is made and the plasma frequency, defined via Equation (2), is used. This leads Equation (24) to be rewritten as:

$$1 + \frac{\omega\_p^2}{k^2} \int\_{-\infty}^{+\infty} dv \, \frac{1}{\omega/k - v} \frac{\partial f\_0}{\partial v} = 0. \tag{25}$$

Note, *f*<sup>0</sup> is now the number of particles per velocity meaning that due to normalization, one needs:

$$\int\_{-\infty}^{+\infty} dv \, f\_0(v) = 1,\tag{26}$$

which fixes constants in the distribution function *f*0.

#### **3. The Dispersion Relation**

Let us now discuss the obtained dispersion relation as given by Equation (25). In order to evaluate this further, the Maxwellian distribution,

$$f\_{\mathbb{D}}(v) = \frac{1}{\sqrt{2\pi}} \frac{1}{v\_{\varepsilon}} e^{-v^2/(2v\_{\varepsilon}^2)}\,\,\,\,\,\tag{27}$$

is employed.

The latter function is visualized via Figure 1. Note, this is a Maxwellian distribution in one dimension and, therefore, it is in coincidence with a Gaussian distribution. The used form is correctly normalized meaning that it satisfies Equation (26). Furthermore, the electron thermal speed, given by Equation (3), is used. From Equation (27), it follows that:

$$\frac{\partial f\_0}{\partial v} = -\frac{1}{\sqrt{2\pi t}} \frac{v}{v\_\varepsilon^3} e^{-v^2/(2v\_\varepsilon^2)}.\tag{28}$$

Using this, Equation (25) reads:

$$1 = \frac{\omega\_p^2}{\sqrt{2\pi}kv\_\epsilon^3} \int\_{-\infty}^{+\infty} dv \, \frac{\upsilon}{\omega - \upsilon k} e^{-\upsilon^2/(2v\_\epsilon^2)}.\tag{29}$$

Note that it looks like there is a singularity at *ω* = *vk*. However, this is not so as long as *ω* is a complex number with non-vanishing imaginary part.

The aim here is to derive the dispersion relation *ω* = *ω*(*k*). To this end, let us rewrite the integral occurring in Equation (29) as:

$$I(\omega, k) = \int\_{-\infty}^{+\infty} dv \, \frac{\upsilon}{\omega - \upsilon k} e^{-v^2/(2v\_\varepsilon^2)}.\tag{30}$$

To continue,

$$\begin{array}{rcl} \int\_{-\infty}^{+\infty} dv \, e^{-v^2/(2v\_\varepsilon^2)} & = & \int\_{-\infty}^{+\infty} dv \, \frac{\omega - v\mathbf{k}}{\omega - v\mathbf{k}} e^{-v^2/(2v\_\varepsilon^2)} \\ & = & \omega \int\_{-\infty}^{+\infty} dv \, \frac{1}{\omega - v\mathbf{k}} e^{-v^2/(2v\_\varepsilon^2)} \\ & - & k \int\_{-\infty}^{+\infty} dv \, \frac{v}{\omega - v\mathbf{k}} e^{-v^2/(2v\_\varepsilon^2)} \end{array} \tag{31}$$

is considered.

The integral on the left-hand-side of Equation (31) is a usual Gaussian integral which can be solved via

$$\int\_{-\infty}^{+\infty} d\mathbf{x} \, e^{-c\mathbf{x}^2} = \sqrt{\frac{\pi}{c}} \qquad \text{if} \qquad \text{Re}(c) > 0. \tag{32}$$

The second integral on the right-hand-side of Equation (31) is the desired integral *I*(*ω*, *k*) as defined via Equation (30).

Thus, Equation (31) can be rewritten as:

$$
\sqrt{2\pi}v\_{\varepsilon} = \omega \int\_{-\infty}^{+\infty} dv \, \frac{1}{\omega - vk} e^{-v^2/(2v\_{\varepsilon}^2)} - kI\{\omega, k\}.\tag{33}
$$

Alternatively, this can be written as;

$$I(\omega, k) = \frac{\omega}{k} \int\_{-\infty}^{+\infty} dv \, \frac{1}{\omega - vk} e^{-v^2/(2v\_\varepsilon^2)} - \frac{\sqrt{2\pi}v\_\varepsilon}{k}. \tag{34}$$

Therewith the dispersion relation (29) becomes:

$$1 = \frac{\omega\_p^2}{\sqrt{2\pi}kv\_\varepsilon^3} \left[ \frac{\omega}{k} \int\_{-\infty}^{+\infty} dv \, \frac{1}{\omega - vk} e^{-v^2/(2v\_\varepsilon^2)} - \frac{\sqrt{2\pi}v\_\varepsilon}{k} \right].$$

In the remaining integral, the substitution,

$$\alpha = \frac{v}{\sqrt{2}v\_{\varepsilon}}\tag{35}$$

is employed to derive

$$1 = \frac{\omega\_p^2}{\sqrt{2\pi}k\upsilon\_\varepsilon^3} \left[ \frac{\omega}{k} \sqrt{2}\upsilon\_\varepsilon \int\_{-\infty}^{+\infty} d\mathbf{x} \frac{e^{-\mathbf{x}^2}}{\omega - \sqrt{2}\upsilon\_\varepsilon k\mathbf{x}} - \frac{\sqrt{2\pi}\upsilon\_\varepsilon}{k} \right]. \tag{36}$$

To continue, the parameter,

$$\alpha := \frac{\omega}{\sqrt{2}v\_{\varepsilon}k},\tag{37}$$

is defined.

Note that the latter parameter is a complex number in the general case. Therewith, one can write Equation (36) as:

$$1 = \frac{\omega\_p^3}{\sqrt{2\pi}k^3v\_\varepsilon^3} \left[ \frac{\omega}{\omega\_p} \int\_{-\infty}^{+\infty} dx \, \frac{1}{a-x} e^{-x^2} - \frac{\sqrt{2\pi}v\_\varepsilon k}{\omega\_p} \right]. \tag{38}$$

In the result obtained, Equation (6) can be now used to finally arrive at

$$1 + \frac{1}{\sqrt{2\pi} (\lambda\_D k)^3} \left[ \sqrt{2\pi} (\lambda\_D k) - \frac{\omega}{\omega\_p} J(\mathfrak{a}) \right] = 0,\tag{39}$$

where the integral,

$$J(\mathbf{a}) = \int\_{-\infty}^{+\infty} d\mathbf{x} \, \frac{1}{\alpha - \mathbf{x}} e^{-\mathbf{x}^2},\tag{40}$$

is used.

Let us solve this integral. First, the integral is generalized via

$$J(\mathfrak{a}, \mathfrak{z}) = e^{-\mathfrak{a}^2} \int\_{-\infty}^{+\infty} d\mathfrak{x} \, \frac{1}{\mathfrak{a} - \mathfrak{x}} e^{\mathfrak{z}(\mathfrak{a}^2 - \mathfrak{x}^2)}. \tag{41}$$

Note, we are looking for *J α*, *β* = 1 ≡ *J α* . One can derive that

$$\begin{array}{rcl} \frac{d}{d\boldsymbol{\beta}} J(\boldsymbol{a}, \boldsymbol{\beta}) &=& \boldsymbol{e} \quad e^{-\boldsymbol{a}^2} \int\_{-\infty}^{+\infty} \boldsymbol{d}x \, \frac{\boldsymbol{a}^2 - \boldsymbol{x}^2}{\boldsymbol{a} - \boldsymbol{x}} e^{\boldsymbol{\beta}(\boldsymbol{a}^2 - \boldsymbol{x}^2)} \\\\ &=& \boldsymbol{e} \quad e^{-\boldsymbol{a}^2} \int\_{-\infty}^{+\infty} \boldsymbol{d}x \, (\boldsymbol{a} + \boldsymbol{x}) e^{\boldsymbol{\beta}(\boldsymbol{a}^2 - \boldsymbol{x}^2)} \\\\ &=& \boldsymbol{a} e^{-\boldsymbol{a}^2} \int\_{-\infty}^{+\infty} \boldsymbol{d}x \, e^{\boldsymbol{\beta}(\boldsymbol{a}^2 - \boldsymbol{x}^2)} .\end{array} \tag{42}$$

Here, the symmetry of the integral is used and, at the end, one would only need to solve a Gaussian integral. With the help of Equation (32), one finally derives:

$$\frac{d}{d\beta}J(\alpha,\beta) = a\sqrt{\pi}e^{-\alpha^2}\beta^{-1/2}e^{\beta\alpha^2}.\tag{43}$$

Note, here and above, *β* is assumed to be a positive real number. To continue, the result is integrated over *β*:

$$J(a,\beta=1) = J(a,\beta=0) + a\sqrt{\pi}e^{-a^2} \int\_0^1 d\beta \,\beta^{-1/2} e^{\beta a^2}.\tag{44}$$

In the integral therein the substitution *β* = *ξ*<sup>2</sup> is used:

$$\frac{d\beta}{\sqrt{\beta}} = 2d\overline{\xi}.\tag{45}$$

Therewith:

$$J(a,\beta=1) = J(a,\beta=0) + 2a\sqrt{\pi}e^{-a^2} \int\_0^1 d\xi \, e^{a^2\xi^2}.\tag{46}$$

In Appendix A, it is demonstrated that

$$J(a, \beta = 0) = -i\pi e^{-a^2}.\tag{47}$$

Furthermore, *J α*, *β* = 1 = *J α* and, thus, Equation (46) reads:

$$J(a) = 2a\sqrt{\pi}e^{-a^2} \int\_0^1 d\xi \, e^{a^2\xi^2} - i\pi e^{-a^2}.\tag{48}$$

Then, using the substitution *t* = *αξ*:

$$J(\mathfrak{a}) = 2\sqrt{\pi}e^{-\mathfrak{a}^2} \int\_0^{\mathfrak{a}} dt \, e^{t^2} - i\tau e^{-\mathfrak{a}^2}. \tag{49}$$

The remaining integral corresponds to an imaginary error function (see, e.g., [16]):

$$\text{Erfi}(\alpha) = \frac{2}{\sqrt{\pi t}} \int\_0^{\alpha} dt \, e^{t^2}. \tag{50}$$

Using this in Equation (49) finally yields:

$$J(a) = \pi e^{-a^2} \left[ \text{Erfi}\left(a\right) - i \right]. \tag{51}$$

After using Equation (51) in Equation (39), one finds for the dispersion relation:

$$1 + \frac{1}{\sqrt{2\pi} \left(\lambda\_D k\right)^3} \left\{ \sqrt{2\pi} \left(\lambda\_D k\right) - \pi \frac{\omega}{\omega\_p} \left[ \text{Erfi}\left(a\right) - i \right] e^{-a^2} \right\} = 0,\tag{52}$$

where the parameter *α* is given by Equation (37). Alternatively, the result obtained can be written as:

$$\begin{split} 1 + \frac{1}{\left(\lambda\_{\rm D} k\right)^{2}} & \quad - \quad \sqrt{\frac{\pi}{2}} \frac{\omega}{\omega\_{p} \left(\lambda\_{\rm D} k\right)^{3}} \text{Erfi}\left(\alpha\right) e^{-\alpha^{2}} \\ & \quad + \quad i \sqrt{\frac{\pi}{2}} \frac{\omega}{\omega\_{p} \left(\lambda\_{\rm D} k\right)^{3}} e^{-\alpha^{2}} = 0. \end{split} \tag{53}$$

Furthermore, the parameter *α* can be rewritten as:

$$
\kappa = \frac{1}{\sqrt{2}} \frac{\omega}{\omega\_p} \frac{1}{\lambda\_D k}. \tag{54}
$$

Equation (53) is a nonlinear equation for the dispersion relation *ω* = *ω*(*k*). It depends only on two parameters: *ω<sup>p</sup>* and *λD*. Below, numerical and analytical solutions of Equation (53) are considered.

#### **4. Numerical Solution for the General Case**

Instead of dealing with Equation (53), let us return and use Equation (39) as starting point for numerical investigations. Equation (48) can be written as:

$$J(\alpha) = 2\alpha \sqrt{\pi} e^{-\alpha^2} K(\alpha) - i \pi e^{-\alpha^2},\tag{55}$$

where

$$\mathcal{K}(\mathfrak{a}) = \int\_0^1 d\xi^\pi e^{\mathfrak{a}^2 \xi^2} \tag{56}$$

is used.

Substituting these equations into Equation (39) gives:

$$\begin{split} \mathbf{1} &\quad + \quad \frac{1}{\sqrt{2\pi} \left(\lambda\_D k\right)^3} \Big[ \sqrt{2\pi} \left(\lambda\_D k\right) \\ &\quad - \quad \frac{\omega}{\omega\_p} 2a\sqrt{\pi}e^{-a^2} K(a) + i\pi \epsilon \frac{\omega}{\omega\_p} e^{-a^2} \Big] = 0. \end{split} \tag{57}$$

For numerical investigations, it is useful to define the dimensionless quantities,

*x* := *λDk* and *y* := *ω*/*ωp*, (58)

so that the parameter *α*, defined via Equation (54), becomes:

$$
\alpha = \frac{y}{\sqrt{2}\chi}.\tag{59}
$$

Therewith, Equation (57) can be written as:

$$1 + \frac{1}{x^2} - \frac{y^2}{x^4} e^{-y^2/(2x^2)} K\left(\frac{y}{\sqrt{2}x}\right) + i\sqrt{\frac{\pi}{2}} \frac{y}{x^3} e^{-y^2/(2x^2)} = 0. \tag{60}$$

The solution of the latter equation is the complex quantity *y* as a function of the real variable *x*. It is straightforward to solve Equation (60) by using a standard Newton-solver in combination with MATLAB; the details are given in Appendix B. The numerical solutions are shown in Figure 2, where the analytical solution obtained for asymptotic limits is also given. This analytical solution is derived in Section 5 below.

**Figure 2.** Visualized is the numerical solution of Equation (60). Shown are the real part of *y* = *ω*/*ω<sup>p</sup>* (top two lines) as well as its imaginary part (bottom two lines) versus the variable *x* = *λDk*. Compared are the numerical solution (solid lines) with the analytical limits (dotted lines) given by Equations (82) and (91), respectively. *ω<sup>p</sup>* is the plasma frequency and *λ<sup>D</sup>* is the Debye length.

#### **5. Analytical Solutions in Asymptotic Limits**

In Section 3, the nonlinear Equation (53) was derived and solved numerically. In order to obtain pure analytical solutions, we consider the case of

$$|\mathfrak{a}| \equiv \frac{|\omega|}{\sqrt{2}v\_{\mathfrak{e}}k} \gg 1,\tag{61}$$

meaning that the phase speed of the wave, *ω*/*k*, is assumed to be much larger than the electron thermal speed, *ve*. Due to this assumption, one can employ the following asymptotic expansion (see, e.g., [16]),

$$\text{Erfi}\left(\alpha\right) \approx \frac{1}{\sqrt{\pi\alpha}}e^{\alpha^2} \left[1 + \left(2\alpha^2\right)^{-1} + 3\left(2\alpha^2\right)^{-2}\right].\tag{62}$$

Using the expansion (62) in Equation (53) yields:

$$\begin{split} 1 + \frac{1}{\left(\lambda\_D k\right)^2} &\quad - \quad \sqrt{\frac{\pi}{2}} \frac{\omega}{\omega\_p \left(\lambda\_D k\right)^3} \frac{1}{\sqrt{\pi a}} \left[1 + \frac{1}{2a^2} + \frac{3}{4a^4}\right] \\ &\quad + \quad i\sqrt{\frac{\pi}{2}} \frac{\omega}{\omega\_p \left(\lambda\_D k\right)^3} e^{-a^2} = 0. \end{split} \tag{63}$$

With the help of Equation (54) this equation can be rewritten as:

$$1 - \frac{1}{\left(\lambda\_D k\right)^2} \frac{1}{2a^2} \left[1 + \frac{3}{2a^2}\right] + i\sqrt{\pi} \frac{\alpha}{\left(\lambda\_D k\right)^2} e^{-a^2} = 0. \tag{64}$$

Using definition (58) of *x* in Equation (64) gives:

$$
\alpha^2 - \frac{1}{2\alpha^2} - \frac{3}{4\alpha^4} + i\sqrt{\pi}\alpha e^{-\alpha^2} = 0.\tag{65}
$$

To further evaluate Equation (65), let us write

$$
\omega = \omega\_{\mathbb{R}} + i\delta\omega\_{\prime} \tag{66}
$$

and, therefore,

$$
\mathfrak{a} = \mathfrak{a}\_{\mathbb{R}} + i\delta\mathfrak{a}.\tag{67}
$$

In the following, it is assumed that *δω* and *δα* are small. As *δω* corresponds to the imaginary part of *ω*, it describes the damping of the Langmuir wave. Thus, the assumption of small *δω* and *δα* corresponds to weak damping.

To continue, let us Taylor-expand and take into account only terms linear in *δα*:

$$\begin{array}{rcl} \frac{1}{\mathfrak{a}^2} & = & \frac{1}{\left(\mathfrak{a}\_R + i\delta\mathfrak{a}\right)^2} \approx \frac{1}{\mathfrak{a}\_R^2} - \frac{2i}{\mathfrak{a}\_R^3}\delta\mathfrak{a}\_\prime \\\frac{1}{\mathfrak{a}^4} & = & \frac{1}{\left(\mathfrak{a}\_R + i\delta\mathfrak{a}\right)^4} \approx \frac{1}{\mathfrak{a}\_R^4} - \frac{4i}{\mathfrak{a}\_R^5}\delta\mathfrak{a}\_\prime \\\ e^{-\mathfrak{a}^2} & = & e^{-\left(\mathfrak{a}\_R + i\delta\mathfrak{a}\right)^2} \approx \left(1 - 2i\mathfrak{a}\_R\delta\mathfrak{a}\right)e^{-\mathfrak{a}\_R^2}. \end{array} \tag{68}$$

Using the above equations in Equation (65) yields:

$$\begin{aligned} &\mathbf{x}^2 - \frac{1}{2a\_R^2} + \frac{i\delta\mathbf{q}}{a\_R^3} - \frac{3}{4a\_R^4} + \frac{3i\delta\mathbf{q}}{a\_R^3} \\ &+ -i\sqrt{\pi} \left(a\_R - 2a\_R^2 i\delta\mathbf{a} + i\delta\mathbf{a}\right) e^{-a\_R^2} = 0. \end{aligned} \tag{69}$$

Taking the real part of the latter equation gives:

$$\left(\alpha^2 - \frac{1}{2a\_R^2} - \frac{3}{4a\_R^4} + \left(2\sqrt{\pi}a\_R^2\delta a - \sqrt{\pi}\delta a\right)e^{-a\_R^2} = 0. \tag{70}$$

As soon as *α<sup>R</sup>* 1 is focused on (see, e.g., Equation (61)), the term with the exponential function can be neglected. Thus, in the lowest-order:

$$
\alpha^2 - \frac{1}{2\alpha\_R^2} - \frac{3}{4\alpha\_R^4} = 0.\tag{71}
$$

Due to Equation (37), one gets:

$$
\alpha\_R = \frac{\omega\_R}{\sqrt{2}v\_c k}.\tag{72}
$$

Furthermore, one finds

$$
\infty = \lambda\_D k = \frac{\upsilon\_c}{\omega\_p} k \tag{73}
$$

as a consequence of Equation (6).

Combining all the above findings, allows us to finally derive from Equation (71):

$$1 - \frac{\omega\_p^2}{\omega\_R^2} \left( 1 + 3 \frac{v\_\epsilon^2 k^2}{\omega\_R^2} \right) = 0. \tag{74}$$

This can be rewritten as:

$$
\omega\_{\mathbb{R}}^4 - \omega\_p^2 \omega\_{\mathbb{R}}^2 - 3\omega\_p^2 v\_c^2 k^2 = 0,\tag{75}
$$

which can be understood as a quadratic equation for *ω*<sup>2</sup> *<sup>R</sup>*. It has the solutions

$$
\omega\_R^2 = \frac{1}{2} \left[ \omega\_p^2 \pm \sqrt{\omega\_p^4 + 12v\_c^2 k^2 \omega\_p^2} \right]. \tag{76}
$$

Therein, taking the plus-sign and employing Equation (6), one finds:

$$
\omega\_R^2 = \frac{1}{2}\omega\_p^2 \left[1 + \sqrt{1 + 12\lambda\_D^2 k^2}\right]. \tag{77}
$$

So far it was assumed that the restriction, given by Equation (61), is valid. Furthermore, the damping effect was assumed to be small meaning that *δω* is small. This essentially means that *ω<sup>R</sup> δω*. Furthermore, the condition, given by Equation (61), turns into

$$
\omega\_R^2 \equiv \frac{\omega\_R^2}{2v\_\varepsilon^2 k^2} \gg 1. \tag{78}
$$

Using therein Equations (6) and (77) allows us to rewrite this condition:

$$1 + \sqrt{1 + 12\lambda\_D^2 k^2} \gg 4(\lambda\_D k)^2. \tag{79}$$

It is clear that this can only be valid if

$$
\lambda\_D k \ll 1.\tag{80}
$$

Using condition (80) in Equation (77) allows us to Taylor-expand so that one obtains:

$$
\omega\_R^2 = \omega\_p^2 \left( 1 + 3\lambda\_D^2 k^2 \right). \tag{81}
$$

One can find that the result obtained concides with the plasma wave dispersion relation of Equation (7). In dimensionless quantities (58):

$$y\_R = \sqrt{1 + 3x^2}.\tag{82}$$

Considering the imaginary part of Equation (69), one finds:

$$
\frac{\delta\alpha}{\alpha\_R^3} + \frac{3\delta\alpha}{\alpha\_R^5} + \sqrt{\pi}u\_R e^{-\alpha\_R^2} = 0.\tag{83}
$$

This equation can be rearranged:

$$
\delta \mathfrak{a} = -\sqrt{\pi} \frac{\mathfrak{a}\_{\mathbb{R}}^{\delta}}{3 + \mathfrak{a}\_{\mathbb{R}}^{2}} e^{-\mathfrak{a}\_{\mathbb{R}}^{2}}.\tag{84}
$$

As soon as *α<sup>R</sup>* 1:

$$
\delta\mathfrak{a} = -\sqrt{\pi} \mathfrak{a}\_{\mathbb{R}}^{4} e^{-\mathfrak{a}\_{\mathbb{R}}^{2}}.\tag{85}
$$

From Equation (54) it follows that

$$
\alpha\_R = \frac{1}{\sqrt{2}} \frac{\omega\_R}{\omega\_p} \frac{1}{\lambda\_D k} \qquad \text{and} \qquad \delta\alpha = \frac{1}{\sqrt{2}} \frac{\delta\omega}{\omega\_p} \frac{1}{\lambda\_D k}. \tag{86}
$$

Therefore, one finds:

$$
\delta\omega = -\sqrt{\frac{\pi}{8}} \frac{\omega\_R^4}{\omega\_p^3} \frac{1}{\left(\lambda\_D k\right)^3} e^{-a\_R^2}.\tag{87}
$$

Combining Equation (81) and condition (80), one gets *ω<sup>R</sup>* ≈ *ω<sup>p</sup>* in the lowest order. Therewith, the result obtained for *δω* reads:

$$
\delta\omega = -\sqrt{\frac{\pi}{8}} \frac{\omega\_p}{\left(\lambda\_D k\right)^3} e^{-a\_R^2}.\tag{88}
$$

Therein one can use Equations (81) and (86):

$$\begin{array}{rcl} \alpha\_R^2 &=& \frac{\omega\_R^2}{2\omega\_P^2} \frac{1}{\left(\lambda\_D k\right)^2} \\\\ &=& \frac{1}{2\left(\lambda\_D k\right)^2} \left(1 + 3\lambda\_D^2 k^2\right) \\ &=& \frac{1}{2\left(\lambda\_D k\right)^2} + \frac{3}{2} \end{array} \tag{89}$$

and finds:

$$
\delta\omega = -\sqrt{\frac{\pi}{8}} \frac{\omega\_p}{\left(\lambda\_D k\right)^3} e^{-\frac{1}{2\lambda\_D^2 k^2} - \frac{3}{2}}.\tag{90}
$$

In dimensionless quantities (58), one gets:

$$
\delta y = -\sqrt{\frac{\pi}{8}} \frac{1}{x^3} e^{-\frac{1}{2x^2} - \frac{3}{2}}.\tag{91}
$$

Equations (82) and (91) are visualized in Figure 2 together with the numerical solution discussed in Section 4. One can immediately see that the agreement between the numerical solution and the analytical one is much better for smaller values of *x* corresponding to small wave numbers. The assumption of small wave numbers is part of the analytical calculation presented above (see, e.g., Equation (80)) and, therefore, the deviation between the analytical and numerical solutions for larger wave numbers was predictable.

Note that as soon as the case *λDk* 1 is considered, one can further approximate Equation (90) by

$$
\delta\omega = -\sqrt{\frac{\pi}{8}} \frac{\omega\_p}{\left(\lambda\_D k\right)^3} e^{-\frac{1}{2\lambda\_D^2 k^2}},\tag{92}
$$

which is in perfect agreement with Equation (17) of [6]. A discussion about the validity and possible improvements of asymptotic formulas for the Landau damping rates of electrostatic waves can be found in [17].

Using the form (66) in the electric field given by Equation (19), yields:

$$
\delta E(\mathbf{x}, t) \propto \delta E(k, \omega) e^{i(kx - \omega t)} \propto e^{-i\omega\_R t + \delta \omega t}.\tag{93}
$$

As *δω* < 0 was derived, one finds damping of the electric field associated with the wave. This damping is the renowned Landau damping and the corresponding damping rate is given by Equations (90) or (92) depending on the desired accuracy.

#### **6. Discussion and Conclusions**

The derivations and discussions, presented in the current paper, are based on lecture notes I have developed for a course about magnetohydrodynamics. The aim was to obtain a mathematically simpler description of Landau damping compared to what Landau originally developed in his seminal paper [6]. This means, in particular, that tools of complex analysis, such as contour integrals, are avoided and one just solves standard integrals in combination with imaginary error functions.

Although using methods of complex analysis was avoided throughout the current paper, one still arrives at the correct Langmuir dispersion relation as given by Equation (81) as well as the correct Landau damping rate as given by Equation (90). These results were obtained analytically by employing Equations (61) and (80). The more general dispersion relation is given by Equation (57) corresponding to a nonlinear equation. However, it can be solved by using standard tools of computational physics. All findings are visualized in Figure 2.

Although the main intention behind the current paper was to present a simpler derivation of known results, Equation (53) is quite general since it does not require to consider asymptotic limits. Therefore, the method, proposed in this paper, could lead to an improved understanding of Landau damping in some other limits. Furthermore, it could be possible that the method developed here can be used to explore other effects such as the dispersion relation of electron-acoustic waves.

It needs to be emphasized that the exploration of Landau damping is still an active field of research pursued both in mathematics and in theoretical physics. A comprehensive overviews of Landau damping which includes some historical remarks, derivations, and mathematical details can be found in [18,19].

**Funding:** This research was funded by the Natural Sciences and Engineering Research Council (NSERC) of Canada.

**Data Availability Statement:** This study does not report any data.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **Appendix A. Some Mathematical Details**

Equation (41) can be written as

$$J(\alpha, \beta) = e^{(\beta - 1)\alpha^2} \int\_{-\infty}^{+\infty} dz \, \frac{1}{\alpha - z} e^{-\beta z^2} \,. \tag{A1}$$

Note that *α* is a complex number, given by Equation (37), and *β* is a positive real number. Equation (A1) can be rewritten as:

$$\begin{split} J(\boldsymbol{\alpha}, \boldsymbol{\beta}) &= \quad \epsilon^{(\beta - 1)\alpha^2} \int\_{-\infty}^{+\infty} dz \, \frac{\boldsymbol{\alpha} + z}{\boldsymbol{\alpha}^2 - z^2} e^{-\beta z^2} \\ &= \quad \epsilon^{(\beta - 1)\alpha^2} \int\_{-\infty}^{+\infty} dz \, \frac{\boldsymbol{\alpha}}{\boldsymbol{\alpha}^2 - z^2} e^{-\beta z^2} \end{split} \tag{A2}$$

where the symmetry properties of one of the integrals is used. In the limit *β* → 0,

$$J(a, \beta \to 0) = e^{-a^2} \int\_{-\infty}^{+\infty} dz \, \frac{a}{a^2 - z^2}. \tag{A3}$$

The remaining integral can be solved via substitution:

$$z = \mathfrak{a}\tanh\left(y\right),\tag{A4}$$

and then,

$$dz = \alpha \left[ 1 - \tanh^2 \left( y \right) \right] dy,\tag{A5}$$

in the integral. Furthermore, for the inverse hyperbolic tangent function, the relations,

$$
\tanh^{-1}\left(+\infty\right) = -\frac{1}{2}i\pi\,,\tag{A6}
$$

and

$$
\tanh^{-1}\left(-\infty\right) = \frac{1}{2}i\pi,\tag{A7}
$$

are employed. Therefore, one gets:

$$\int\_{-\infty}^{+\infty} dz \, \frac{a}{a^2 - z^2} = -i\pi \,\tag{A8}$$

and, thus:

$$J(\mathfrak{a}, \mathfrak{z} \to 0) = -i\pi e^{-\mathfrak{a}^2}.\tag{A9}$$

#### **Appendix B. The Newton-Solver**

The task is to solve the nonlinear relation given by Equation (60). For convenience, let us write this equation down again:

$$1 + \frac{1}{\mathbf{x}^2} - \frac{y^2}{\mathbf{x}^4} e^{-y^2/(2\mathbf{x}^2)} K(a) + i\sqrt{\frac{\pi}{2}} \frac{y}{\mathbf{x}^3} e^{-y^2/(2\mathbf{x}^2)} = 0. \tag{A10}$$

In the latter, the function,

$$K(\alpha) = \int\_0^1 dz \, e^{\alpha^2 z^2} \,\tag{A11}$$

is used, as well as:

$$\alpha = \lambda\_D k, \qquad y = \omega / \omega\_{p\prime} \qquad \text{and} \qquad \alpha = \frac{y}{\sqrt{2}\alpha}. \tag{A12}$$

Therewith, one can write:

$$\mathcal{K}(x,y) = \int\_0^1 dz \, e^{\frac{y^2}{2\tilde{w}^2}z^2}. \tag{A13}$$

The nonlinear relation above is solved via a standard Newton method; see, e.g., in [20]. The task is to find *y* for a given *x*; therefore, it is needed to compute *y* for each *x*. This *y* is obtained via the iteration method,

$$y\_{n+1} = y\_n - \frac{f(y\_n)}{f'(y\_n)'} \tag{A14}$$

where, in the case considered here,

$$f(y) = 1 + \frac{1}{\mathbf{x}^2} - \frac{y^2}{\mathbf{x}^4} e^{-y^2/(2\mathbf{x}^2)} K(\mathbf{x}, y) + i\sqrt{\frac{\pi}{2}} \frac{y}{\mathbf{x}^3} e^{-y^2/(2\mathbf{x}^2)}.\tag{A15}$$

Furthermore, one also needs to compute the derivative of the latter function with respect to *y*. One derives:

$$\begin{split} f^{\prime}(y) \equiv \frac{\partial f}{\partial y} &= -2\frac{y}{\lambda^{4}}e^{-y^{2}/(2\mathbf{x}^{2})}\mathcal{K}(\mathbf{x},y) + \frac{y^{3}}{\lambda^{4}}e^{-y^{2}/(2\mathbf{x}^{2})}\mathcal{K}(\mathbf{x},y) - \frac{y^{2}}{\lambda^{4}}e^{-y^{2}/(2\mathbf{x}^{2})}\mathcal{K}^{\prime}(\mathbf{x},y) \\ &+ \quad i\sqrt{\frac{\pi}{2}}\frac{1}{\lambda^{3}}e^{-y^{2}/(2\mathbf{x}^{2})} - i\sqrt{\frac{\pi}{2}}\frac{y^{2}}{\lambda^{5}}e^{-y^{2}/(2\mathbf{x}^{2})}\, \end{split} \tag{A16}$$

where

$$K'(x,y) = \int\_0^1 dz \, \frac{yz^2}{x^2} e^{\frac{y^2}{2x^2}z^2} \tag{A17}$$

is used.

Combining Equations (A14)–(A17) determines *y* for a given *x* numerically. Performing these calculations for a set of *x* values yields the dispersion relation *y*(*x*) which was looked for.

#### **References**

