Numerical consideration of quasilinear electron cloud dynamics in plasma

TL;DR

This paper evaluates numerical schemes for modeling hot electron cloud dynamics in plasma, demonstrating a scheme that accurately captures asymptotic behavior and is suitable for solving initial value problems.

Contribution

It introduces and validates a finite difference scheme that correctly models quasilinear electron cloud dynamics in plasma, ensuring proper asymptotic behavior.

Findings

The proposed scheme matches analytical solutions in the hydrodynamic limit.

It effectively models initial value problems for various electron distributions.

The scheme demonstrates correct asymptotic behavior for plasma electron dynamics.

Abstract

The dynamics of a hot electron cloud in the solar corona-like plasma based on the numerical solution of kinetic equations of weak turbulence theory is considered. Different finite difference schemes are examined to fit the exact analytical solutions of quasilinear equations in hydrodynamic limit (gas-dynamic solution). It is shown that the scheme suggested demonstrates correct asymptotic behavior and can be employed to solve initial value problems for an arbitrary initial electron distribution function.

Equations

$$W/nT\ll (k{\lambda }_D{)}^2,$$

$$\frac{\partial f}{\partial t}+v\frac{\partial f}{\partial x}=\frac{e^2}{{ϵ}_0{m}^2}\frac{\partial }{\partial v}\frac{W}{v}\frac{\partial f}{\partial v},$$

$$\frac{\partial W}{\partial t}=\frac{\pi {\omega }_{pe}}{n}{v}^{2}W\frac{\partial f}{\partial v},{\omega }_{pe}=kv$$

$$F(v,x,t=0)=g_0(v)\text{mbox}exp(-x^2/d^2),$$

$$W(v,x,t=0)=10^{-8}m{n}^{\prime }{v}_0^3/{\omega }_{pe},$$

$$\frac{\partial F}{\partial t}+\alpha V\frac{\partial F}{\partial X}=\frac{1}{\tau }\frac{\partial }{\partial V}D\frac{\partial F}{\partial V},$$

$$\frac{\partial D}{\partial t}=\frac{1}{\tau }{V}^{2}D\frac{\partial F}{\partial V},$$

$$F(V,X,t=0)=G_0(V)\text{mbox}exp(-X^2),\text{mbox}where{\int }_0^{\infty }{G}_0(V)dV=1,$$

$$W(V,X,t=0)=10^{-8},$$

$$U_j=\frac{1}{\Delta y}{\int }_{y_{j-1/2}}^{y_{j+1/2}}U(y)dy,$$

$${\nabla }_j^{+}U\equiv \frac{U_{j+1}-U_j}{y_{j+1}-y_j},{\nabla }_j^{-}U\equiv \frac{U_j-U_{j-1}}{y_j-y_{j-1}}.$$

$$\frac{\partial F}{\partial t}=\frac{1}{\tau }\frac{\partial }{\partial V}D\frac{\partial F}{\partial V},$$

$$\frac{\partial D}{\partial t}=\frac{1}{\tau }{V}^{2}D\frac{\partial F}{\partial V},$$

$$\frac{\partial }{\partial t}\left[F-\frac{1}{V^2}\frac{\partial D}{\partial V}\right]=0,$$

$$F_{\infty}(V)=\left\{\begin{array}[]{ll}1,&\mbox{$V<1$}\\ 0,&\mbox{$V>1$}\end{array}\right.\;\;\;\;\;D_{\infty}(V)=\left\{\begin{array}[]{ll}V^{3},&\mbox{$V<1$}\\ 0,&\mbox{$V>1$}\end{array}\right.$$

$$F_{\infty}(V)=\left\{\begin{array}[]{ll}1,&\mbox{$V<1$}\\ 0,&\mbox{$V>1$}\end{array}\right.\;\;\;\;\;D_{\infty}(V)=\left\{\begin{array}[]{ll}V^{3}(1-V),&\mbox{$V<1$}\\ 0,&\mbox{$V>1$}\end{array}\right.$$

$${\int }_{V_{j-1/2}}^{V_{j+1/2}}\left[F(V,t_{k+1})-F(V,t_k)\right]dV={\int }_{t_k}^{t_{k+1}}\left[w(V_{j+1/2},t)-w(V_{j-1/2},t)\right]dt,$$

$$w(V,t)=\frac{D}{\tau }\frac{\partial F}{\partial V},$$

$$w_{j+1/2}=a_j{\nabla }_j^{+}{F}_j,$$

$$\frac{1}{\Delta t}{\int }_{t_k}^{t_{k+1}}w(V_{j+1/2},t)dt=\sigma w_j^{k+1}+(1-\sigma )w_j^k,$$

$$F_j^{k+1}-F_j^k={\nabla }_j^{-}{a}_j{\nabla }_j^{+}\left(\sigma F_j^{k+1}+(1-\sigma )F_j^k\right),$$

$$a_j={\left(\frac{\tau }{\Delta V}{\int }_{V_j}^{V_{j+1}}\frac{dV}{D(V)}\right)}^{-1}={\left[\tau {\int }_0^1\frac{ds}{D(V+\Delta Vs)}\right]}^{-1},$$

$$F_j^{k+1}-F_j^k=\frac{\Delta t}{\tau }{\nabla }_j^{-}{D}_j^k{\nabla }_j^{+}{F}_j^{k+1}$$

$$D_j^{k+1}-D_j^k=\frac{\Delta t}{\tau }{V}_j^2{D}_j^k{\nabla }_j^{+}{F}_j^{k+1}$$

$$F_j^{k+1}-F_j^k=\frac{\Delta t}{\tau }{\nabla }_j^{-}{D}_j^k{\nabla }_j^{+}{F}_j^k$$

$$D_j^{k+1}-D_j^k=\frac{\Delta t}{\tau }{V}_j^2{D}_j^k{\nabla }_j^{+}{F}_j^k$$

$$F_j^{k+1}-F_j^k=\frac{\Delta t}{2\tau }{\nabla }_j^{-}(D_{j+1}^k+D_j^k){\nabla }_j^{+}{F}_j^k$$

$$D_j^{k+1}-D_j^k=\frac{\Delta t}{\tau }{V}_j^2{D}_j^k{\nabla }_j^{-}{F}_j^k$$

$$\Delta t\le \text{mbox}min\left[\tau /(2V_j^2{\nabla }_j^{+}{F}_j^{k+1}),\tau \Delta V^2/2a_j(V_j)\right]\text{mbox}for1\le j\le M,k>0$$

$$\frac{\partial F}{\partial t}+\gamma \frac{\partial F}{\partial X}=0,\gamma =\text{mbox}const>0$$

Full text

1 Introduction

Accelerated particle beams occur in the wide range of astrophysical situations as solar flares, cosmic rays, radio jets, magnetospheres of pulsars, planetary atmospheres, etc [1, 2]. The bright signatures of the electron beams in a plasma are the solar type III bursts [3, 4]. In accordance with current understanding of these bursts an electron beam propagating along open magnetic field lines from the Sun toward the Earth generate Langmuir waves, which are partly transformed into observable radio emission via nonlinear plasma processes [3, 5]. The typical density of electron beams is low and Coulomb collisions have no influence on beam dynamics. The main process of beam interaction with the surrounding plasma is resonant Cherenkov’s generation and absorption of plasma waves [5]. The Langmuir waves excited by the electron beam flatter the electron distribution function [6]. Thus, for the characteristic time of electron-wave interaction (quasilinear time) $\tau\approx n^{\prime}/\omega_{pe}n$ (where $n^{\prime}$, $n$ are the beam an plasma density, and $\omega_{pe}=\sqrt{e^{2}n/\epsilon_{0}m}$ is the electron plasma frequency) plateau is formed at the electron distribution function [6]. Propagation of electrons disturbs local equilibrium and in the next spatial point generation of waves repeats. Many authors considered the problem of electron beam propagation analytically [7, 8, 9, 10, 11] as well as numerically [12, 13, 14, 15, 16, 17]. However, the results obtained are far from quantitative agreement [5]. This is mainly connected with the fact that the system of kinetic equations describing the problem is nonlinear with stiff relaxation terms. The problem can be significantly simplified if the smallness of quasilinear time can be taken into account as it was suggested by Ryutov and Sagdeev [7]. Thus, implying that plateau is established at the electron distribution function and high level of plasma waves is generated at every spatial point one can turn from kinetic to gas-dynamic description. The gas-dynamic system of equation describing electron cloud dynamics in hydrodynamic limit has been recently derived by Mel’nik [11]. The solution obtained in [11, 18] for the initial distribution function $\partial g_{0}(v)/\partial v>0$ is a compact object propagating in a plasma with conservation of the particle number, energy, and wave energy. However, the solutions obtained were not supported by numerical calculations. It is argued in [19] that it is even impossible to obtain the stable solution for initially unstable electron beam. On the contrary, gas-dynamic description [11, 18] demonstrates that initially unstable electron beam can lead to interesting solutions.

Numerical consideration of quasilinear equations with initial distribution function $\partial g_{0}(v)/\partial v<0$ has been conducted several times [10, 12, 13, 14, 15, 18].

In the given paper the numerical solution of kinetic equation of quasilinear theory is discussed. The dynamics of the electron beam is considered at the scale much larger that the size of electron beam. We show that the transport term should be approximated by the higher than first order finite difference operator. The different presentation of collision terms are examined. The results of numerical solution are presented for typical parameters of electron beams and solar plasma are presented.

2 The statement of the problem

Let us consider the propagation of the electron beam cloud when the energy density of excited Langmuir waves is much less than that of surrounding plasma

[TABLE]

where $W$ is the energy density of Langmuir waves, $T$ is the temperature of surrounding plasma, $k$ is the wave number, and $\lambda_{D}=\sqrt{k_{B}T\epsilon_{0}/ne^{2}}$ is electron Debye length. Our analysis is limited by one-dimensional kinetic equations following [7, 8]. One-dimensional beam propagation is supported by numerical solution of 3D equations [22]. In the application to solar burst III type those electrons propagate along magnetic field which energy $\mu_{0}H^{2}/2\gg nmv_{0}^{2}/2$ [3], that ensures one-dimensional character of electron propagation.

In the case of type III bursts, as it was shown by Vedenov, Velikhov and Sagdeev [20] and by Drummond and Pines [21], one can use equations of quasilinear theory [7]

[TABLE]

[TABLE]

where $f(v,x,t)$ is the electron distribution function, $W(v,x,t)$ is the spectral energy density of Langmuir waves. $W(v,x,t)$ plays the same role for waves as the electron distribution function does for particles. The system (2,3) describes the resonant interaction $\omega_{pe}=kv$ of electrons and Langmuir waves, i.e. electron with the velocity $v$ can emit or absorb a Langmuir wave with the phase velocity $v_{ph}=v$. The group velocity of Langmuir waves is small as $v_{g}\approx v_{Te}^{2}/v\ll v$ and therefore the corresponding term in the left side of the equation (3) is omitted [7].

The clouds of fast electrons are formed in the spatially limited regions of solar corona where acceleration takes place. Therefore, spatially bounded beam is taken for consideration. The initial electron distribution function is

[TABLE]

where $d$ is the characteristic size of the electron cloud and $g_{0}(v)$ is the initial distribution of electrons in the velocity space and $\int\limits_{0}^{\infty}g_{0}(v)dv=n^{\prime}$. It is also implied that initially the spectral energy density of Langmuir waves is of the thermal level and homogeneously distributed in space

[TABLE]

where $v_{0}$ is some characteristic velocity of the electron cloud. The system of kinetic equations (2,3) is nonlinear with two characteristic time scales. The first is the quasilinear time $n/n^{\prime}\omega_{pe}$ that is determined by the interaction of particles and waves. The second scale is the time length of the electron cloud $d/v_{0}>>n/n^{\prime}\omega_{pe}$. And we are interested in dynamics of the electron cloud at the time scale $t>>d/v_{0}$.

3 Numerical method

The problem we are confronted with is an initial value problem. There are a variety of techniques available for the numerical solution of such partial differential equations. We use finite differencing (see, for example, [23, 24, 25]). For our further consideration we rewrite partial differential equations together with initial conditions in the following form

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where we use normalized velocity $V=v/v_{0}$, distance $X=x/d$, quasilinear time $\tau=n/n^{\prime}\omega_{pe}\pi$, electron distribution function $F(V,X,t)=f(v,x,t)v_{0}/n^{\prime}$, and $D(V,X,t)=W(v,x,t)\omega_{pe}/vmn^{\prime}v_{0}^{3}$, $\alpha=v_{0}/d$. All terms in equations (6,7) are presented in 1/s units.

For numerical solution of the equation (6,7) we will introduce the grid points $V_{j}$, for $j=1,...,M$, and $X_{i}$ for $i=1,...,N$ with uniform mesh width $\Delta V=V_{j+1}-V_{j}$, $\Delta X=X_{i+1}-X_{i}$. The discrete time level $t_{k}$ is also uniformly spaced with the time step $\Delta t$. Given any function $U$ we denote its nodal value by $U_{j+1/2}=U(y_{j+1/2})$ and its cell average values by

[TABLE]

where $y$ is an variable ($X$ or $V$). We also introduce operators

[TABLE]

The use of operator splitting [26] allows us to solve equations (6, 7) by determining the finite difference operators for each individual term of the equation. We now describe the finite difference schemes used to advance each term.

3.1 Quasilinear relaxation

It case the electron beam is homogeneously distributed in space transport term in (6) should be omitted. Thus we obtain the system of equations that describes quasilinear relaxation of an electron beam in velocity space

[TABLE]

[TABLE]

The equations(12,13) can be solved analytically [6]. Thus substituting (13) into (12) one obtains quasilinear integral

[TABLE]

that allows us to obtain final distribution of particles and waves via initial conditions. The initially unstable electron distribution function $\partial G_{0}(V)/\partial V>0$ leads to generation of plasma waves and flattering of the electron distribution function. Quasilinear relaxation continues until $\partial F(V)/\partial V=0$, and plateau is formed at the electron distribution function. The characteristic time of beam-plasma interaction is $\tau$ (for this time a half of initial electron beam energy is transferred into waves).

For initially monoenergetic electron beam $G_{0}(V)=\delta(V-1)$ using (14) [6] one finds the following steady state solution (the solution of the system at $t\rightarrow\infty$)

[TABLE]

Similar to (15) for the initial distribution function $G_{0}(V)=2V$, for $V<1$ which has been considered in [18] one obtains the following steady state solution

[TABLE]

The solution of quasilinear equations (12, 13) shows that the maximum of the spectral energy density depends on the initial electron distribution function. The maximum of $D(V)$, $D=1$ is reached at $V=1$ when the initial electron distribution function is an monoenergetic beam.

To construct an conservative finite difference scheme we follow [24]. For equation (12) one writes the equation of balance in the cell $V_{j-1/2}\leq V\leq V_{j+1/2}$, and $t_{k}\leq t\leq t_{k+1}$

[TABLE]

where

[TABLE]

is the particle flux in velocity space. Integrating each term in (17) and using that

[TABLE]

[TABLE]

one obtains the following general finite-difference equation

[TABLE]

where $\sigma$ is a number $0\leq\sigma\leq 1$,and $a_{j}(D)$ is a functional

[TABLE]

which can be approximated in a number of different ways [24].

Let us consider a few interesting cases. For $\sigma=1$ and $a_{j}=D_{j}/\tau$ we rederive the fully implicit scheme used by Grognrad [14]

[TABLE]

[TABLE]

Equation (23) is unconditionally stable whereas second one (24) is stable only when $\Delta t\leq\tau/(2V_{j}^{2}\nabla_{j}^{+}F_{j}^{k+1})$. Using (23) and (24) together one obtains unbounded growth of the Langmuir waves [14]. To solve the problem spontaneous terms have been added to the right hand side of equations (23) and (24) [14].

The scheme (23,24) is not acceptable when spontaneous terms are omitted [14] (as in our case). Indeed, in absence of spontaneous terms we know that the relaxation is non-linear diffusion in velocity space which decreases $|\partial F/\partial V|$ in a finite domain but increases $|\partial F/\partial V|$ to very large values at the limits of this domain.

For $\sigma=0$ one obtains fully explicit schemes. If the functional is approximated as $a_{j}=D_{j}/\tau$ we obtain the scheme is used by Takakura [13](hereafter scheme I)

[TABLE]

[TABLE]

and more accurate scheme (scheme II), when $a_{j}=(D_{j+1}+D_{j})/2\tau$

[TABLE]

[TABLE]

Both schemes are conditionally stable. The criteria of stability is

[TABLE]

Using that $D$ is always less than $1$, and $\nabla_{j}^{+}F_{j}^{k+1}\leq 1/\Delta V^{2}$ for an arbitrary initial distribution function timestep should be as small as $\Delta t\leq\tau\Delta V^{2}/2$.

To test the schemes I and II we use typical beam-plasma parameters that ensures smallness of $\tau$. The results of numerical tests are presented in fig. 1-2. Test run shows that both schemes correctly approximate the process of quasilinear relaxation. The scheme II gives better approximation for both initial electron distribution functions presented in this section. In fig. 2 we obtain the best coincidence between numerical results and the analytical solution (16). In case of initially monoenergetic electron beam (actually $G_{0}(V)=2\mbox{exp}(-(V-1)^{2}/\Delta V_{0}^{2})/\sqrt{\pi}\Delta V_{0}$ for $V\leq 1$, $\Delta V_{0}\ll 1$) we have the correct asymptotic behavior. Decreasing initial dispersion in the beam $\Delta V_{0}$ we obtain the spectral energy density approaching to the analytical solution (13) (see 3). Scheme I gives higher plateau but more accurate drop of distribution function at $V=1$, whereas scheme II better approximates plateau but has more smooth border at $V=1$. From fig. 1, 2 we can see that scheme II lead to appearance of ”accelerated” particles ($F$ is different from zero at $V=1+\Delta V$). Consequently, if we take into account transport of particles the maximum velocity of the plateau will unphysically grow at beam propagation. Thus, scheme I may be applied for short-time beam dynamics at the time scale $t\leq d/v_{0}$ whereas scheme II is more suitable for long time dynamics $t\gg d/v_{0}$.

3.2 Transport of particles

Propagation term plays an important role for initial distribution function $\partial G_{0}(V)/\partial V<0$. The equation to consider is

[TABLE]

Since dynamic calculations at $t>>\tau$ consume much computer time and therefore finite difference schemes for (6,7) are usually taken as simple as possible. First order upwind representation of transport operator (30) is taken in the majority of cases [10, 12, 13, 14, 15]

[TABLE]

However, this first order scheme seems to be not enough to caught the correct asymptotic behavior of the system at $t\gg d/v_{0}$. It is shown [27, 28, 29] that monotonic transport is the best finite difference method for equation (30). Using this method one finds that

[TABLE]

where

[TABLE]

where $F_{i}^{k}$ is the value of $F$ at position $X_{i}$ and the time $t=k\Delta t$.

In fact we are interested in the evolution of the initial distribution. Therefore, a good test of the numerical solution is the propagation of spatially finite distribution. We take an initial distribution $F(X,t=0)=\mbox{exp}(-X^{2})$. In fig. 2 the numerical solutions are compared to the exact analytical solution of (30) (which is just a Gaussian moving with velocity $\beta$) for upwind (31) and monotonic transport (32) methods. For all methods, the numerical solution conserves particle number, so that in general, the height of the numerical solution is a good measure of accuracy.

The fact that upwind scheme has ”worser” approximation can drastically influence the electron beam dynamics. Although, the quasilinear terms are large enough in comparison with the propagation term, they are of the same order for some points of $(X,V)$ plain [30]. In these areas of $(X,V)$ plain the error of the first order upwind operator can be comparable to the quasilinear diffusion in the velocity space

[TABLE]

when $D\approx 0$ or $\partial F/\partial V\leq 0$.

On our opinion this can lead to a wrong asymptotic behavior of numerical solution. Thus, when the initial electron distribution function is stable as to generation of plasma waves $\partial G_{0}(V)/\partial V<0$ the transport term of kinetic equations plays the main role. The electron propagation causes the change of electron distribution function and finally the unstable electron distribution $\partial F/\partial V>0$ appears. The rate of the quasilinear relaxation strongly depends on beam density at this point. Therefore, the underestimated electron density will change the initial point of relaxation [30].

4 Numerical results

In this section we consider the long time dynamics of electron cloud for typical parameters of the beam and plasma and compare to the analytical solution found in hydrodynamic limit (gas-dynamic solution [18]).

For the time of quasilinear relaxation $\tau$ plateau is formed at the electron distribution function and the high level of plasma waves is generated at every spatial point [20]. The smallness of quasilinear relaxation time can be used to turn from kinetic to gas-dynamic description of the problem that was suggested by Ryutov and Sagdeev [7] and done by Mel’nik [11]. Plateau at the electron distribution function and high level of Langmuir turbulence are implied to be formed. The method is similar to ordinary hydrodynamics, where we integrate kinetic equations implying that Maxwell’s distribution is assumed at every point. Following [11, 18] we can obtain analytical solution when $\tau\ll t$ for $G_{0}(V)=2V,\;\;\;v\leq 1$:

[TABLE]

[TABLE]

where

[TABLE]

(see [31] for details.)

Electrons (35) accompanied by Langmuir waves (36) propagate in a plasma as a beam-plasma structure with the constant velocity $\gamma/2$.

The electron distribution function $F(V,X,t)$ and the spectral energy density of Langmuir waves $D(V,X,t)V$ are presented in fig. 4 at the time moment $t=10.0$s. As it is implied [11, 18] plateau is established in the wide area of velocities from $1$ down to $V\approx 55$ at every spatial point. We notice that the plateau height exponentially increases on the forward front (which is presented by the points $X>55$ in fig. 3) and exponentially decreases at the back front ($X<55$ ). In the point ( $X\approx 55$ ) the plateau height reaches maximum value. The symmetry of the form of the initial electron beam conserves, but now it is a form of distribution of electron stream (whereas electrons move with various velocities). In fig. 4 these Langmuir waves are presented at the moment $t=10$s. Like electrons, Langmuir waves are concentrated near $X\approx 55$. The spectral energy density $DV$ reaches its maximum value at $V\approx 0.8$ that coincides with the theoretical value (36).

The observable difference between the profile of beam-plasma structure and the theoretical solution is explained by the fact that despite the smallness of quasilinear time it is a finite value $\tau>0$. Some electrons on the tails of the structure do not take part in quasilinear relaxation and therefore propagate freely away from the structure. Decreasing the quasilinear time we can reach better agreement with the analytical solution (fig. 5).

5 Conclusions

The propagation of an electron cloud has complex nonlinear properties. In order to describe electron propagation in plasma correctly one needs attention in choose of the method of numerical solution. The system of partial differential equations has stiff terms and finite difference scheme should correctly describe different scales of the system. The different schemes have different optimal approximations for various time scales. The analysis of the system shows that first order upwind transport of the particles may lead to a wrong asymptotic regime. To ensure sufficient accuracy monotonic scheme has been suggested for numerical consideration. The monotonic scheme is found to be accurate enough for the problem considered.

To design schemes with the correct behavior one should use a discrete analogous of the asymptotic limit of continuous system. Therefore, the results in asymptotic regimes are compared with exact analytical solution. Test calculation demonstrate optimistic agreement for different parameters of plasma and a beam. However, the difference between numerical solution and gas-dynamic solution is observed. The main assumption of the gas-dynamic approach is that plateau is form at every spatial point. This is not true for regions of $(X,V)$ plain where the fast electron density is low. Quasilinear relaxation for these particles is not a fast process and electron propagate in a plasma almost freely. The velocities of these electron differ from the speed of beam plasma structure and the electrons move away from the structure. This leads to the loss of particles and as a result the structure becomes lower and wider. However, the mentioned difference becomes smaller if we consider the systems with smaller $\tau$.