Adaptive spectral solution method for the Landau and Lenard-Balescu equations

TL;DR

This paper introduces an adaptive spectral method using Laguerre polynomials for solving the Landau and Lenard-Balescu equations, improving efficiency and applicability to complex plasma collision problems.

Contribution

The authors develop an adaptive Laguerre polynomial spectral method that reduces memory usage and extends the solvable range of plasma kinetic equations.

Findings

Successfully solves relaxation problems previously infeasible with spectral methods.

Demonstrates applicability to quantum Lenard-Balescu equation for hydrogen plasma.

Addresses memory limitations through adaptive re-projection of basis functions.

Abstract

We present an adaptive spectral method for solving the Landau/Fokker-Planck equation for electron-ion systems. The heart of the algorithm is an expansion in Laguerre polynomials, which has several advantages, including automatic conservation of both energy and particles without the need for any special discretization or time-stepping schemes. One drawback is the $O(N^3)$ memory requirement, where $N$ is the number of polynomials used. This can impose an inconvenient limit in cases of practical interest, such as when two particle species have widely separated temperatures. The algorithm we describe here addresses this problem by periodically re-projecting the solution onto a judicious choice of new basis functions that are still Laguerre polynomials but have arguments adapted to the current physical conditions. This results in a reduction in the number of polynomials needed, at the expense of increased solution time. Because the equations are solved with little difficulty, this added time is not of much concern compared to the savings in memory. To demonstrate the algorithm, we solve several relaxation problems that could not be computed with the spectral method without re-projection. Another major advantage of this method is that it can be used for collision operators more complicated than that of the Landau equation, and we demonstrate this here by using it to solve the non-degenerate quantum Lenard-Balescu equation for a hydrogen plasma.

Tables

Table 1. Table 1: The rate of change of the ion temperature for fully ionized hydrogen ( n e = n i subscript 𝑛 𝑒 subscript 𝑛 𝑖 n_{e}=n_{i} ) with T e = 3.0 × 10 7 ​ K subscript 𝑇 𝑒 3.0 superscript 10 7 𝐾 T_{e}=3.0\times 10^{7}K and T i = 1.0 × 10 5 ​ K subscript 𝑇 𝑖 1.0 superscript 10 5 𝐾 T_{i}=1.0\times 10^{5}K , calculated from Eq. ( 80 ) and compared with the results in Ref. [ 10 ] . At these conditions, coupled modes are important, and one must capture accurately the details of the dielectric function in the integration. Clearly, our method is doing this. In the fourth column, we show the calculation done when the B n ​ n ′ ​ n ′′ subscript 𝐵 𝑛 superscript 𝑛 ′ superscript 𝑛 ′′ B_{nn^{\prime}n^{\prime\prime}} of Eq. ( 78 ) are set to zero. Although some difference can be seen, it is minimal.

$n_i$	Eq. (80)	Ref. [10]	Eq. (80), $B_{n{n}^{\prime }{n}^{\prime \prime }}=0$	Landau
$1\times {10}^{23}$	0.0297	0.0296	0.0292	0.0200
$1\times {10}^{24}$	0.225	0.225	0.221	0.128
$1\times {10}^{25}$	1.54	1.54	1.49	0.568

Equations

$$\frac{\partial f_e}{\partial t}$$

$$\frac{\partial f_i}{\partial t}$$

$$C_{ei}(f_e,f_i)=-\frac{1}{4{\pi }^2{\hslash }^2}\int d^3{\mathbf{v}}^{\prime }\int d^3\mathbf{k}\frac{|{\varphi }_{ei}(k){|}^2}{{\left|ϵ\left(k,\mathbf{k}\cdot \mathbf{v}+\frac{\hslash k^2}{2m_e}\right)\right|}^2}$$

$$\times$$

$$\mu \equiv \frac{m_e{m}_i}{m_e+m_i},$$

$${\varphi }_{ie}(k)=\frac{4\pi Z_i{e}^2}{k^2},$$

$${ϵ}_{ie}(X,Y)=1+\frac{{\eta }_{ie}^2}{X^2}\left[\frac{\gamma }{Z}{w}_e\left({\alpha }_e^{(ie)}Y\right)+w_i\left({\alpha }_i^{(ie)}Y\right)\right],$$

$$f_j(v,t)=f_j^{\mathrm{eq}}(v)\sum _{n=0}^N{A}_n^j(t)L_n^{\left(\frac{1}{2}\right)}\left(\frac{{\beta }_j{m}_j{v}^2}{2}\right)$$

$$f_j^{\mathrm{eq}}(v)\equiv n_j{\left(\frac{m_j{\beta }_j}{2\pi }\right)}^{3/2}\exp \left(-\frac{m_j{\beta }_j{v}^2}{2}\right),$$

$$L_m^{\left(\frac{1}{2}\right)}\left(\frac{m_j{\beta }_j{v}^2}{2}\right)$$

$$\frac{d{A}_n^e}{dt}$$

$$\frac{d{A}_n^i}{dt}$$

$$A_0^j(t)=1$$

$$\frac{d{A}_0^e}{dt}=\frac{d{A}_0^i}{dt}=0$$

$$A_1^e=-\frac{n_i{\beta }_e}{n_e{\beta }_i}{A}_1^i,$$

$$\frac{d{A}_1^e}{dt}=-\frac{n_i{\beta }_e}{n_e{\beta }_i}\frac{d{A}_1^i}{dt}$$

$$C_{1n^{\prime }{n}^{\prime \prime }}^{(ee)}$$

$$C_{1n^{\prime }{n}^{\prime \prime }}^{(ii)}$$

$$C_{1n^{\prime }{n}^{\prime \prime }}^{(ei)}$$

$$A_n^{\prime }=n!\sum _{j=0}^n(-1{)}^{n+j}{A}_j\frac{{\left(\frac{{\beta }^{\prime }}{\beta }\right)}^j{\left(\frac{{\beta }^{\prime }}{\beta }-1\right)}^{n-j}}{j!(n-j)!}.$$

$$A_0^{\prime }=A_0=1$$

$${\beta }^{\prime }=\frac{\beta }{1-A_1(\beta )}.$$

$${\Delta }_{\sigma }\equiv \left|\frac{{\beta }_{\sigma }^{\prime }}{{\beta }_{\sigma }}-1\right|$$

$$C_{lk}^{(\sigma {\sigma }^{\prime })n}=-\frac{n_e\sqrt{2\pi }{b}_2^{(\sigma {\sigma }^{\prime })}}{2\sqrt{b_1^{(\sigma {\sigma }^{\prime })}}}({\beta }_{\sigma }{m}_{\sigma }{)}^{1/2}({\beta }_{{\sigma }^{\prime }}{m}_{{\sigma }^{\prime }}{)}^{1/2}{e}^4{Z}_{\sigma }^2{Z}_{{\sigma }^{\prime }}^2\frac{n!}{\Gamma (n+3/2)}$$

$$\times$$

$$q_{lk}^{(ie)n}(X,Y)\equiv \sum _{j=0}^{\min (l,n)}{L}_{n-j}^{\left(-\frac{1}{2}\right)}\left({\overline{Y}}_{-}^2\right)$$

$$\times$$

$$Y^2$$

$$X^2$$

$${\overline{Y}}_{\pm }$$

$$Y_{\pm }$$

$$b_1^{(\sigma {\sigma }^{\prime })}$$

$$b_2^{(\sigma {\sigma }^{\prime })}$$

$${\overline{P}}_{ie}$$

$${\overline{Q}}_{ie}$$

$$P_{ie}$$

$$Q_{ie}$$

$$\gamma \equiv \frac{{\beta }_e}{{\beta }_i}$$

Full text

††journal: Journal of Computational Physics

1 Introduction

The Landau equation, or its equivalent formulation in terms of the Fokker-Planck equation [1], is a primary tool in kinetic theory [2]. In plasma physics, for example, it describes a plasma at conditions for which the Coulomb logarithm dominates the collision integrals, i.e., at weak coupling. Here, the kinetic energy of the particles dominates the potential and small-angle grazing collisions tend to be the most important scattering events. The numerical solution of this system presents many challenges; in particular, great care must be taken when using velocity discretization to ensure conservation of particles and energy [3, 4]. Although discretization methods are standard in the field, a spectral method, consisting of an expansion in Laguerre polynomials, was recently used [5] to solve the quantum Lenard-Balescu equation and several variants, including the Landau equation. The motivation for using this technique is that it neatly handles the integration over the dielectric function, but another of its nice properties is that it identically conserves particles and energy without requiring any special method for time stepping. The main drawback is that the method has expensive memory requirements; if one expands the distribution in $N$ polynomials, $N^{3}$ coefficients must be computed from the equation. If there are $s$ particle species, this increases to $(s^{2}-s+1)N^{3}$, so one can quickly find oneself with resolution problems.

Here, we aim to address this difficulty by proposing an adaptive method in which the solution is periodically re-projected onto a new set of basis functions that are better suited to representing the current physical conditions. In this way, we greatly reduce the number of coefficients required, at the expense of increasing the time it takes to solve the actual differential equations. We test the algorithm on electron-ion temperature equilibration, as this problem encapsulates the challenge to the spectral method, namely representing distributions whose temperatures vary widely between species and undergo large variations in time. With the original method, one is limited in how far one can separate the temperatures of the two species; increasing the temperature ratio means needing more polynomials to represent the distribution. Here, we show that this restriction is removed with the adaptive method and we can simulate an arbitrarily large temperature separation. We show also that the algorithm is effective in like-mass cases, where the particles do not remain Maxwellian over the course of the equilibration.

In addition to the Landau equation, our method can be adapted to solve a range of other kinetic equations, and we demonstrate this here by solving the non-degenerate quantum Lenard-Balescu equation. This equation has rarely been attempted numerically [6, 5] because it requires an integration in frequency and wavenumber over a dielectric function. This dielectric function adds another degree of non-linearity to the system because it contains the distribution function, although this is not the main source of difficulty. The real problem is that under some circumstances, the integrand contains very sharp peaks that are difficult to capture by discretizating either in velocity or over wavenumber and frequency. The spectral method, along with a fairly innocuous approximation in the dielectric function, allows us to circumvent this problem and compute these integrals to good accuracy and quickly enough to be used a numerical solution of the equation. We demonstrate the accuracy of the method by using it to compute the electron-ion equilibration rate in the coupled-mode regime, where missing the sharp peak in the integrand can significantly affect the answer. Our results perfectly match those obtained with sophisticated techniques [7, 8] that were developed to handle this problem. Of course, we can do more than just compute the rate, which is just a single number, but we can compute entire time evolutions under the Lenard-Balescu equation. We give a few examples in section 4.2. We find, however, that there is not much difference between calculating the distributions with the full Lenard-Balescu collision operator versus using Landau with a well-chosen Coulomb logarithm, even in cases where the temperature separation is large. This is consistent with previous numerical comparisons of the two equations [5, 6].

2 Equation

To get at the Landau equation, we start with the non-degenerate quantum Lenard-Balescu equation. To use the spectral method, we must calculate the ordinary differential equations satisfied by the expansion coefficients, and this is somehow easier with $\hbar$ non-zero. The Landau equation then appears in the limit $\hbar\rightarrow 0$ and as $\epsilon$, the dielectric function (to be defined below), goes to $1$. The quantum version of the Lenard-Balescu equation features a completely convergent collision operator; a dielectric function provides convergence at large scales, $k\rightarrow 0$, and quantum diffraction ensures it at small scales, $k\rightarrow\infty$. The kinetic equations for the electrons and ions are given by

[TABLE]

where

[TABLE]

and $\hbar$ is Planck’s constant, $\mu$ the reduced mass,

[TABLE]

and $m_{e}$ and $m_{i}$ are the electron and ion masses. We use the Coulomb potential,

[TABLE]

where $e$ is the electron charge and $Z_{i}$ is the charge state of species $i$. Note that (4) defines all the collision integrals in (1) and (2); to get $C_{ee}$, for example, one just sets $i=e$. In the dielectric function we take the limit $\hbar\rightarrow 0$, which has almost no effect on the solution [9, 10] and results in the form

[TABLE]

where the functions $w(x)$ have been derived elsewhere [5] and are given in section 3.6. As already mentioned, we get the Landau equation by taking $\hbar\rightarrow 0$ everywhere in the equation and setting $\epsilon=1$. This renders the $k$-integral divergent at both of its limits; these must then be cut off, resulting in a Coulomb logarithm. We now discuss our solution method and its implementation.

3 Method

3.1 Expansion

We use the expansion in Laguerre polynomials used in Ref. [5] to solve the quantum Lenard-Balescu equation. For the species $j$,

[TABLE]

where $f_{j}^{\mathrm{eq}}(v)$ is the Maxwell distribution at the temperature $\beta_{j}$,

[TABLE]

and $L_{n}^{(1/2)}(x)$ is the $n^{\mathrm{th}}$ Laguerre polynomial of order $1/2$. We stress that the temperature $\beta_{j}$ is just an arbitrary expansion parameter and need not be related to any real property of the system. This parameter plays a central role in the adaptive method.

The equations for the $A_{m}(t)$ are found by multiplying both sides of the equation by

[TABLE]

and integrating over $d^{3}{\bf v}$. The resulting coupled ordinary differential equations are then of the form

[TABLE]

where, for the Landau equation, the $C^{(jk)}_{nn^{\prime}n^{\prime\prime}}$ are constants proportional to $\log\Lambda$. These will be given in section 3.5 and those of the quantum Lenard-Balescu equation in section 3.6.

3.2 Conservation

Particles will be conserved exactly provided

[TABLE]

is satisfied throughout the calculation (this simple relation is the reason for the choice of $1/2$ for the order of the Laguerre polynomials). Equations (11) and (12) have the property that

[TABLE]

identically, so that (13) is trivially satisfied if it holds initially. As for conservation of energy, its expression in terms of the $A_{n}$ is

[TABLE]

which can easily be derived by computing the kinetic energy integrals for each species and demanding their sum be constant. The equations (11) and (12) have the property that

[TABLE]

which arises from the identities

[TABLE]

so that the method guarantees conservation of energy.

3.3 Adaptive temperature

The parameters $\beta_{j}$ are (effective) inverse temperatures and can be chosen arbitrarily; it is this freedom that we exploit in the adaptive method. Clearly, if the distribution we wish to resolve is in fact Maxwellian at the temperature $\beta$, then choosing $\beta_{j}=\beta$ means that the expansion is simply $A_{n}^{j}=\delta_{n0}$ and we need only one coefficient to represent the distribution. Choosing a different $\beta_{j}$ would require more coefficients, with the number depending on how far $\beta$ is from $\beta_{j}$. Although this is a somewhat trivial example, it is generally the case that some choices of $\beta_{j}$ are better than others from the standpoint of minimizing $N$. Moreover, a good choice at one point in the solution may be very poor at a later time. For example, if we consider electron-ion temperature equilibration, the species are roughly Maxwellian over their entire evolution with only their temperatures changing in time. For species $j$, if we set $\beta_{j}$ to its initial temperature then we can easily resolve the initial condition, but as time goes on the actual temperature of the distribution becomes very different from this $\beta_{j}$. If the initial temperature separation of the two species was large enough, we will run out of resolution with $N$ coefficients. The obvious remedy to this problem is to recalculate, at various intervals, the expansion coefficients using the current temperature, with the result that most will return to zero. Although this only sounds useful when the distributions remain approximately Maxwellian, the latter is often the case in systems of particles with large mass disparities such as electrons and ions. Furthermore, using an effective temperature in (8) that corresponds to the kinetic energy of the distribution is usually the most efficient choice even when particles are not Maxwellian [5].

If we are given a set of coefficients $A_{i}$, for any species, corresponding to the temperature parameter $\beta$, and we wish to find the coefficients $A_{i}^{\prime}$ corresponding to $\beta^{\prime}$, the formula is

[TABLE]

The derivation is given in A. Note that this transformation has no effect on particle or energy conservation as

[TABLE]

and if (15) is satisfied, then conservation of kinetic energy and (20) guarantee that (15) is also satisfied in the primed variables.

The next issue to consider is when to actually perform these reprojections. The most obvious thing to do is to find the $\beta^{\prime}$ associated with the current kinetic energy and check its ratio with the $\beta$ used in the expansion, $\beta^{\prime}/\beta$. When its deviation from unity is larger than some predefined tolerance, we trigger a reprojection. Computing this $\beta^{\prime}$ is a simple matter, and we derive the formula in B,

[TABLE]

For each species $\sigma$, we define the quantity

[TABLE]

and if $\Delta_{\sigma}>c$ we reproject. The smaller we make $c$, the more frequent will be the reprojections, which will generally lower the number of polynomials needed. Of course, there are other alternatives. One is to look at the size of the highest-order coefficients and to reproject when they appear to be getting too large. Another is simply to reproject at regular intervals. We will give examples of all three schemes in action below.

3.4 Equation coefficients

Deriving the coefficients $C_{nn^{\prime}n^{\prime\prime}}$ is a straightforward but tedious task. We do not outline the steps here, but refer the reader to Ref. [5], where the calculation was done in detail for the one-component, single-temperature case. The generalization to multiple species with different temperatures is fairly obvious and requires nothing conceptually different from what was presented there. Although the collision operators defined in (4) are 6-fold integrals, and to get the coefficients we must integrate these further over $d^{3}{\bf v}$, the final 9-fold integrals can be reduced by various identities to

[TABLE]

where

[TABLE]

$\Delta\beta_{ie}\equiv\beta_{e}-\beta_{i}$, and the various dimensionless variables are

[TABLE]

where

[TABLE]

For the dielectric function we use the classical approximation; that is, we neglect quantum diffraction [5]. The main benefit of this is that the dielectric function takes the comparatively simple form (7). That is, the variables $X$ and $Y$ separate, facilitating analytical evaluation of the $X$-integral without losing any important physical content. The constant $\gamma$ is the temperature ratio

[TABLE]

and the $\eta$ are given by

[TABLE]

with

[TABLE]

The complex function, $w(Y)$, is given by [5],

[TABLE]

and the constants $\alpha$ are

[TABLE]

3.5 Landau equation

So far, we have kept $\hbar$ in the calculation. Now, we can set it to zero by making $X=0$ except in $dX/X$ where the $\hbar$ will cancel. If we also set $\epsilon(X,Y)=1$, the result is the Landau equation. Only the terms of $q(X,Y)$ that are even in $Y$ survive the integration in (25) and so we write

[TABLE]

Now, we separate this function into a piece proportional to $X^{2}$ and one containing all higher-order terms in $X$,

[TABLE]

where $G(Y)$ is a finite polynomial in $Y$ and $R(X,Y)$ is $\mathrm{O}(X^{4})$. Clearly, the terms $X^{4}$ and higher vanish in the integrand in (27) when $X$ is set to zero, but the $X^{2}$ term gives us

[TABLE]

i.e., the Coulomb logarithm. Thus we have

[TABLE]

where

[TABLE]

The polynomials $G(Y)$ are given by

[TABLE]

where

[TABLE]

The bulk of the work in the adaptive spectral method is in computing the coefficients $S_{nn^{\prime}n^{\prime\prime}}$; the resulting ordinary differential equations are quicky solved with a standard Runge-Kutta method. To determine the coefficients, we use Gaussian quadrature to compute their defining integral, equation (53). To do this, we first observe that all $G(Y)$ are even, as can easily be deduced by glancing at (56). This suggests the substitution $x=Y^{2}$ and we then have

[TABLE]

Clearly, the quadrature scheme should be based on associated Laguerre polynomials of order $1/2$ as the weight function is $x^{1/2}e^{-x}$. To determine the number of points, $N_{p}$, we consider that Gaussian quadrature is exact for integrands consisting of polynomials of degree up to $2N_{p}-1$ that multiply the weight. From (56) it can easily be inferred that the maximum degree of $G_{nn^{\prime}n^{\prime\prime}}(\sqrt{x})$ is around $3N$ (i.e., when $n=n^{\prime}=n^{\prime\prime}=N$), although in practice it is usually slightly less. We therefore have

[TABLE]

Once $N_{p}$ has been chosen, we need only the weights, $W_{j}$, and quadrature points, $x_{j}$. The latter are given by the roots of the $N_{p}^{\mathrm{th}}$ Laguerre polynomial of order $1/2$,

[TABLE]

and the weights are

[TABLE]

The integral determining $S_{nn^{\prime}n^{\prime\prime}}$ is then given by

[TABLE]

provided $N_{p}$ is large enough, otherwise (62) is an approximation. The final ingredient needed is an efficient way to determine $G_{nn^{\prime}n^{\prime\prime}}(\sqrt{x})$ at the quadrature points. For this, we need to determine $L^{\left(\frac{1}{2}\right)}_{n}(\xi x_{j})$ and $L^{\left(-\frac{1}{2}\right)}_{n}(\xi x_{j})$ for all $n$ up to $N$, all $j$ up to $N_{p}$, and a variety of $\xi$. The best way to do this is to use the formulas

[TABLE]

and the recurrence relation

[TABLE]

for each $\xi x_{j}$. The formula (65) is numerically stable so there is no problem going to large $n$. Equation (62) must be evaluated every time there is a re-projection because of the temperature dependences in various parts of (56).

We now have what we need to compute the coefficients of the ordinary differential equations for the Landau system. To solve these equations, we use the fourth-order Runge-Kutta method as implemented in the Boost C++ library [11]. Numerical results are given in Sec. 4.1.

3.6 Quantum Lenard-Balescu equation

To solve the QLB equation, we need a way to compute efficiently the integral over the dielectric function in Eq. (25). The strategy is to exploit the separation between the variables $X$ and $Y$ in the Eq. (7) and to do the $X$-integral analytically. In fact, we will be interested only in the order-unity term and the one proportional to $\ln\eta$. At weak coupling, these are by far the dominant terms in the collision operator.

Unlike in the Landau equation, we do not take the limit $X\rightarrow 0$. As demonstrated in Ref. [5], the terms $\mathrm{O}(X^{4})$ and higher in the expansion (50) do not require the dielectric function for their integrals to converge. This means that we may drop the dielectric function from those integrals as it will contribute terms higher order in $\eta$ than are of interest to us here. However, the integral over $R(X,Y)$ does contribute an order-unity term, which we call $B_{nn^{\prime}n^{\prime\prime}}^{(\sigma\sigma^{\prime})}$ and is defined by

[TABLE]

where $R^{(\sigma\sigma^{\prime})}_{nn^{\prime}n^{\prime\prime}}(X,Y)$ can be calculated from (50). To evaluate this, we use Gauss-Laguerre quadrature with a small modification from what was done in the previous section. Making the substitutions $y=X^{2}$ and $x=Y^{2}$,

[TABLE]

Because the lowest order term in $R(\sqrt{y},\sqrt{x})$ is $y^{2}$, the $1/y^{2}$ in Eq. (67) will cancel, leaving the weight $e^{-y}$. The $1/x$ will be similarly cancelled and thus we see that the best choice here is to use the order [math] Laguerre polynomials for the $y$-integral and order $1/2$ for $x$,

[TABLE]

where $x_{i}$ are given by Eq. (61), $y_{i}$ are the zeros of the order zero (or unassociated) Laguerre polynomials,

[TABLE]

and the weights $W_{j}^{(0)}$ are

[TABLE]

The double integration makes $B_{nn^{\prime}n^{\prime\prime}}^{(\sigma\sigma^{\prime})}$ comparatively expensive to calculate. However, as we will show in Sec. 4.2, these coefficients tend to have an insignificant impact on the calculation, at least for the relaxation problems we consider.

The piece of the integral that does require the dielectric function for convergence we write as

[TABLE]

The expansion of such integrals in $\eta$ was given in Appendix C of Ref. [5] and is easily adapted to the multi-species case,

[TABLE]

where

[TABLE]

and

[TABLE]

$\gamma_{E}\approx 0.5772156649...$ is the Euler constant, and $w_{R}^{(\sigma\sigma^{\prime})}(Y)$ and $w_{I}^{(\sigma\sigma^{\prime})}(Y)$ are given in (7). The coefficients $S_{nn^{\prime}n^{\prime\prime}}^{(\sigma\sigma^{\prime})}$ are the same as those defined in Sec. 3.5.

We know of no simple closed form for the integrals $F^{(\sigma\sigma^{\prime})}_{nn^{\prime}n^{\prime\prime}}$, but despite appearances they are fairly simple to evaluate numerically. We can once again use the Gauss-Laguerre scheme used to compute $S_{nn^{\prime}n^{\prime\prime}}^{(\sigma\sigma^{\prime})}$, so that

[TABLE]

with the weights given in (61). Evaluating $F^{(\sigma\sigma^{\prime})}(\sqrt{x_{i}})$ is fairly simple; although it contains hypergeometric functions, these need only be evaluated during a reprojection and their values may be stored for subsequent time steps. The GNU Scientific Library [12] contains functions that quickly and accurately perform these evaluations. The Laguerre polynomials that appear in $F^{(\sigma\sigma^{\prime})}(\sqrt{x_{i}})$ can be evaluated using the recurrence relation given in Sec. 3.5.

Putting the various pieces together,

[TABLE]

where

[TABLE]

4 Numerical results

4.1 Landau Equation

4.1.1 Electron-proton system

To begin, we demonstrate the resolution problems that can arise when we have too few polynomials and do not use reprojection. We consider a hydrogen system at density $n_{e}=n_{p}=10^{25}\ \mathrm{cm}^{-3}$ and with initial electron and proton temperatures $T_{e}=10\ \mathrm{keV}$ and $T_{p}=10\ \mathrm{eV}$. This is a highly unrealistic temperature separation, but it serves to illustrate the numerical method. In Fig. 1, we plot the proton distribution at the initial time and at $t=1\times 10^{-14}\ \mathrm{s}$ using six polynomials ($n_{\mathrm{max}}=5$). We can see that as the solution progresses, and the effective temperature becomes more widely separated from the expansion temperature, we become unable to represent the distribution with only six polynomials.

Now we perform the same calculation using a scheme in which a reprojection is triggered when $|A^{j}_{4}|>10^{-2}$ for either species $j$. In Fig. 2 we have plotted the proton distribution at the initial time and at every second reprojection. The electron distribution is shown in Fig. 3. Clearly, we have solved the resolution issue as the calculation now proceeds all the way to equilibrium without any problems.

4.1.2 Electron-electron system

In the disparate-mass system, the distributions are basically Maxwellians throughout their whole evolution. The adaptive scheme is clearly ideal for such problems, but we demonstrate here that a wide mass separation is not essential for the method to be useful. We consider electrons at a density $n_{e}=2\times 10^{23}\ \mathrm{cm}^{-3}$ where half the particles start in a Maxwell distribution with temperature $T_{1}=10^{7}K$ and the other half have temperature $T_{2}=2.5\times 10^{8}K$, a separation by a factor of 25. In this case, the system will certainly not maintain Maxwellian profiles over the course of the equilibration and we must thus ensure that we have enough polynomials to capture the entire evolution. We take $n_{\mathrm{max}}=20$ for this calculation, which is not enough to complete the calculation without reprojections. We use a scheme of reprojecting at regular intervals, which, as a general rule, should probably be used for cases in which the species do not remain roughly Maxwellian. In fact, for this case we reproject on every time step. This is not really necessary, but running this way minimizes the number of polynomials needed. In practice, one should find a balance between the number of polynomials used and the number of time steps between reprojections. The result is shown in Figure 4, where various times are shown along with the ultimate equilibrium. Clearly, our algorithm handles this case perfectly well, providing some encouraging evidence of its wider applicability.

4.2 Lenard-Balescu equation

Relaxation problems in which the electrons are much hotter than the ions can exhibit the coupled mode effect in which an ion acoustic mode is excited, leading to a lowering of the energy transfer rate between electrons and ions. This effect has been discussed in many places, from both theoretical [13, 7, 8, 14] and computational [15, 10] perspectives. Near the frequency of an ion acoustic mode, the dielectric function contributes a sharp peak to the integrand in Eq. (25) which can be a challenge to capture accurately [8, 15] in rate calculations, let alone quickly enough to be used in a kinetic solver. In Ref. [10], it was shown that if one performs the integration over the wavenumber ($X$ in our case) first, the remaining integral over the frequency ($Y$) can be accurately handled with Gaussian quadrature. This is exactly the approach we use here. At $t=0$, when the distributions are Maxwellian and the expansions use the true values of $\beta$, the equilibration rate, i.e., the derivative of the ion temperature, is given by

[TABLE]

This rate has been thoroughly studied with the quantum Lenard-Balescu equation, with theories including the effects of quantum degeneracy [7, 15, 10, 16] and strong coupling [16]. In Table 1, we compare the rate calculated with Eq. (80) with the computations done in Ref. [10] in which both a careful numerical integration over the double integral and an analytic method similar to ours were used. Although that work included the effects arising from Fermi statistics, the conditions we are considering here are sufficiently non-degenerate that this should have no impact on the answer. Indeed, as can be seen in the table, our computation exactly agrees with that of Ref. [10].

We show in the fourth column of Table 1 the effect of setting $B_{nn^{\prime}n^{\prime\prime}}=0$; the change is negligible, and is similar to what one finds by neglecting Eq. (64) of Ref. [10]. That the $B$ coefficients do not change the answer significantly is rather convenient; their computation is the most expensive part of the reprojections. However, we cannot rule out the possibility that they may be more important for some problems.

Finally, in Figures 5 to 7, we compare temperature evolutions computed with the QLB equation and the Landau equation. For the Coulomb logarithm, we choose a $\Lambda$ that varies with time to match the logarthmic term we calculated for the Lenard-Balescu equation. That is, we use equation (78) with $B_{nn^{\prime}n^{\prime\prime}}=F_{nn^{\prime}n^{\prime\prime}}=0$ and we recompute $\eta$ from (39) at every time step based on the current conditions. In Figure 5 we show the temperature evolution for hydrogen with $n_{i}=n_{e}=10^{25}\ \mathrm{cm}^{-3}$ and initial temperatures $T_{e}=10^{8}K$ and $T_{i}=10^{7}K$. The electrons are shown in black and the protons in green. Both the Landau and Lenard-Balescu solutions are plotted but there are completely indistinguishable from each other. Lowering the density to $n_{i}=n_{e}=10^{23}\ \mathrm{cm}^{-3}$, we show the temperature evolution for $T_{i}=10^{5}K$ and $T_{e}=10^{6}K$ in Figure 6, and $T_{e}=5\times 10^{6}K$ in Figure 7. In this case, there are visible differences between the Landau and quantum Lenard-Balescu equations but, even when the temperatures are separated by a factor of $50$, these differences are very small. It is possible that by going to even more unrealistically large temperature separations, like those listed in Table 1, one might see more pronounced differences between the solutions. However, it becomes difficult to integrate the ordinary differential equations for the Lenard-Balescu equation at these conditions. The reason for this is unclear. It could perhaps be that we are beginning to violate the time scale assumptions underlying the derivation of the equation [17] and solutions are becoming unphysical. Perhaps there is a more mundane explanation, but as these conditions are highly unrealistic physically we do not pursue the issue here.

5 Conclusion

We have presented an adaptive spectral method for solving kinetic equations, such as the Landau and quantum Lenard-Balescu equations. Our method identically conserves both energy and particles, while the adaptive scheme greatly mitigates the large memory requirements normally attendant on this spectral approach.

We demonstrated the potential of this algorithm by using it to solve the quantum Lenard-Balescu equation, a problem that would be very difficult with standard Fokker-Planck methods [4, 3]. We compared the solutions of electron-ion equilibration problems calculated with the Landau and Lenard-Balescu equations and find that, consistent with other studies into these equations, there is little advantage to solving the full Lenard-Balescu equation, at least for these problems.

This method can be generalized to distributions that are spatially dependent and non-isotropic in velocity, and we hope that it will be used for such problems in the future. It is possible that among the expanded range of problems that will then be accessible, some will provide a good reason to go beyond the Landau equation.

Acknowledgments

We thank Ömer F. Tekin and Hai Le for helpful discussions. Part of this research was performed while the authors were visiting the Institute for Pure and Applied Mathematics (IPAM), which is supported by the National Science Foundation. This work was performed under the auspices of the U.S. Department of Energy at the Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344.

Appendix A Reprojection formula

We wish to take a distribution $f(v)$ in the form:

[TABLE]

where $\{A_{n}(\beta)\}$ are already calculated, and express it in the form

[TABLE]

so that

[TABLE]

We begin evaluating this by multiplying on both sides by $L_{n}^{(1/2)}\left(m\beta^{\prime}v^{2}/2\right)$ and integrating both sides by $\int d^{3}v$. Our equation becomes

[TABLE]

Evaluating the right hand side gives us

[TABLE]

We simplify this by changing the bounds of the integral (using the fact that our function is even), and by using the substitution $x^{2}\equiv m\beta v^{2}/2$. Our equation becomes:

[TABLE]

We can now replace our Laguerre polynomials with Hermite polynomials using

[TABLE]

[TABLE]

.

This turns the left hand side of (88) into

[TABLE]

When $j>n$ our integral is equal to zero. Otherwise (92) becomes

[TABLE]

Now, using the identity

[TABLE]

we find (20).

Appendix B Calculating the energy

In order to do the reprojection, we need to know the value of $\beta^{\prime}$ at time $t$ in terms of the expansion parameter $\beta$. Assume that at time $t$ we know the coefficients ${A_{i}}$ for the Laguerre expansion of $f$ in terms of $\beta$. The kinetic energy ($\mathrm{KE}$) of a species is given by

[TABLE]

where $n_{i}$ is the number density of the species under consideration. Using the fact that $L_{0}^{\left(\frac{1}{2}\right)}(x)=1$ and $L_{1}^{\left(\frac{1}{2}\right)}(x)=\frac{3}{2}-x$, (96) becomes

[TABLE]

Noting that $A_{0}(\beta)=1$ and $\Gamma{(5/2)}=(3/2)\Gamma{(3/2)}=\frac{3}{4}\sqrt{\pi}$ this simplifies to

[TABLE]

The energy of a particle in a Maxwell distribution at $\beta^{\prime}$ is given by

[TABLE]

and setting Eq. (102) equal to Eq. (103) we get Eq. (22).