aa r X i v : . [ phy s i c s . p l a s m - ph ] F e b Quantum Fokker-Planck Modeling of Degenerate Electrons
Hai P. Le ∗ Lawrence Livermore National Laboratory, Livermore, California 94551, USA
Abstract
An implicit and conservative numerical scheme is proposed for the isotropic quantum Fokker-Planck equation describing the evolution of degenerate electrons subject to elastic collisions withother electrons and ions. The electron-ion and electron-electron collision operators are discretizedusing a discontinuous Galerkin method, and the electron energy distribution is updated by animplicit time integration method. The numerical scheme is designed to satisfy all conservationlaws exactly. Numerical tests and comparisons with other modeling approaches are shown todemonstrate the accuracy and conservation properties of the proposed method.
Keywords:
1. Introduction
Coulomb collisions play a central role in plasma physics as the main driving mechanism fortemperature equilibration and transport phenomena. They are typically modeled by a Fokker-Planck (FP) collision operator with the underlying assumptions that the particles are non-degenerate,weakly coupled and that small angle scattering dominates [1]. This approach works well for lowdensity, high temperature ‘classical’ plasmas. At higher densities and lower temperatures, strongCoulomb coupling and quantum e ff ects becomes significant and the classical treatment is no longerapplicable. Quantum kinetic theory for dense plasmas has developed over the last couple decades,but the numerical solution of these equations remains a challenge due to the complex mathematicalstructure involved [2, 3].In the present work, we focus on developing numerical solutions for an isotropic quantumFokker-Planck (qFP) model of degenerate electrons. The electron distribution function is assumedto be spatially homogeneous and isotropic, such that the qFP model can be formulated in termsof the electron energy distribution function. The qFP model is a direct extension of the classicalmodel, in which the e ff ects of quantum degeneracy are taken into account [4, 5, 6], includingthe Pauli exclusion principle, which forbids transitions into occupied states. The qFP equationcan be derived from the grazing collision limit of the quantum Boltzmann equation (also knownas the Boltzmann-Uehling-Uhlenbeck equation) [7]. The domain of applicability for the qFPequation was examined by Daligault [8]. It can be applied to electron collisions in high energy ∗ Corresponding author.
Email address: [email protected] (Hai P. Le)
Preprint submitted to Elsevier January 3, 2021 ensity physics experiments, where the electrons remain weakly coupled over a wide range oftemperature. In addition, since the formulation of the qFP equation is very similar to its classicalversion, existing numerical methods developed for the classical FP equation can be readily adapted[9, 10, 11, 12, 13].The classical FP model has been studied extensively in the context of plasma transport [14, 15]and kinetic simulations based on the FP model are routinely used to study nonequilibrium systems[16, 17, 18]. The qFP model receives much less attention and very few studies discuss numericalmethods for solving the qFP equation. Most of the studies dicussing about the qFP equation arein the context of plasma transport theory [19, 20, 21]. We mention here those studies that discussnumerical solutions to qFP and / or quantum kinetic equations in general. Hu et al. [22] developedan asymptotic preserving scheme for the qFP equation. In this work, the collision operator isdiscretized using the spectral method of Pareschi et al. [11], originally designed for the classicalFP equation. A similar scheme was also proposed for the quantum Boltzmann equation [23].Kitamura [24] studied ultra-fast thermalization dynamics of primary and secondary electrons inmetals by numerically solving the quantum Lenard-Balescu (LB) equation. A similar e ff ort atsolving the non-degenerate quantum LB equation was due to Scullard et al. [25, 26]. We alsonote that Monte Carlo collision algorithms which simulate Coulomb collisions can be extended toinclude quantum degeneracy e ff ects, and these provide an alternative to the qFP approach. [27, 28]The numerical scheme presented here is built upon previous work on nonequilibrium electronkinetics modeling [13, 29]. For simplicity, we consider the case of a uniform and spatially ho-mogeneous plasma. We further assume that the electron distribution is isotropic and formulatethe qFP operators for both electron-ion (ei) and electron-electron (ee) collisions. Degeneracy andstrong coupling e ff ects of ions are neglected, and the ion distribution is assumed to be Maxwellian.The resultant qFP equation is discretized and solved by a fully implicit Discontinuous Galerkin(DG) method. The discretization method is high-order, and can be applied to a non-uniform en-ergy grid. The fully implicit time stepping scheme allows us to take time steps larger than themean collision time. These features are attractive for problems with a large dynamical range ofconditions. Other physics can be incorporated. For example, the same numerical formulation wasapplied to the case of inelastic collisions, e.g., excitation and ionization, which are modeled byBoltzmann collision operators [13]. The extension to include quantum degeneracy for inelasticcollisions will be examined in future work. Transport e ff ects due to spatial inhomogeneity can beincorporated by employing the full spherical Harmonic decomposition of the velocity distributionfunction [1, 30, 31]. This was demonstrated for the case of classical FP equation [29], and can beextended to the qFP equation in a straightforward manner.The rest of the paper is organized as follows. In section 2, we summarize basic propertiesof degenerate electrons and the qFP model for ei and ee collisions. Section 3 describes the DGdiscretization the collision operators, its conservation properties and time integration method. Sev-eral numerical tests are presented in section 4 to validate the qFP model. Comparisons with MonteCarlo simulations and simplified model are also shown. Finally we make some concluding remarksin section 5. The formulation of the qFP equation is discussed in Appendix A.2 . Quantum Fokker-Planck Model Before introducing the qFP model, we review some basic definitions and properties of de-generate electron distributions. Consider a spatially homogeneous and isotropic electron gas ina fully ionized plasma free of external fields. The electron energy distribution function (EEDF)evolves due to collisions among the electrons and with surrounding ions. Electron thermalizationis primarily due to ee collisions, while ei collisions are responsible for temperature equilibrationbetween the ions and electrons. In thermal equilibrium, the mean occupation number of a state ofenergy ε , denoted as ˜ f ( ε ), follows a Fermi-Dirac distribution:˜ f ⋆ ( ε ) = + e ( ε − µ ) / T e , (1)where µ is the chemical potential and T e is the electron temperature. Here we use a superscript( ⋆ ) to indicate thermal equilibrium condition. The EEDF f ( ε ) is related to the mean occupationnumber ˜ f ( ε ) as follows: f ( ε ) = π mh ! / √ ε ˜ f ( ε ) . (2)The total electron density N e and energy E e are given by the zeroth and first order energy momentsof the distribution function: N e = Z ∞ f ( ε ) d ε, (3a) E e = Z ∞ ε f ( ε ) d ε. (3b)For a Fermi-Dirac distribution (1), we have: N e = √ π λ − I / ( µ/ T e ) , (4a) E e = N e T e I / ( µ/ T e ) I / ( µ/ T e ) , (4b)where λ ≡ h √ π mT e is the electron thermal deBroglie wavelength and I n ( y ) ≡ R ∞ x n + e x − y dx is theFermi-Dirac integral of order n . Unlike the classical Maxwell distribution, T e (and µ ) cannot betrivially expressed in terms of N e and E e , i.e., inverting Eq. (4) requires an iterative procedureor numerical fit [32]. Temperature and chemical potential are uniquely defined only when theelectrons are in thermal equilibrium. In the case of a nonequilibrium distribution with fixed N e and E e , Eq. (4) can be used to determine the corresponding equilibrium distribution.Two asymptotic limits can be associated with the Fermi-Dirac distribution: classical (high tem-perature) and strongly degenerate (low temperature). These can be characterized by the degener-acy parameter, defined as T e / E F where E F ≡ ~ m (3 π N e ) / is the Fermi energy. In the classicallimit ( T e / E F ≫ e µ/ T → N e λ and E e → N e T e and the distribution approaches the classical3axwellian. For a strongly degenerate system ( T e / E F ≪ µ → E F and E e → N e E F . At thezero temperature limit, all energy states up to E F are fully occupied, i.e,. ˜ f ( ε ) = ε ≤ E F and0 otherwise.In a nonequilibrium system, the time evolution of the EEDF subject to ei and ee collisionsobeys the following qFP equation (see Appendix A for the mathematical formulation): ∂ f ( ε, t ) ∂ t = ∂∂ε ( J ei + J ee ) . (5)Both ei and ee collision operators are written as a divergence of particle fluxes in energy space.The flux due to ei collisions is written as: J ei = γ ei " f √ ε (cid:16) − ˜ f (cid:17) + T i ∂∂ε f √ ε ! , (6a) γ ei = mM r m π N e Ze ln Λ ei , (6b)where Z is the ion charge state, M is the ion mass, and T i is the ion temperature. Eq. (6) assumesthat the ions are non-degenerate and thermalized, i.e., their energy distribution is a Maxwelliandistribution at T i . It also assumes that the ions’ drift velocities are negligible and their thermalvelocities are much lower than those of the electrons, i.e., T i M ≪ T e m . The particle flux due to eecollisions is: J ee = γ " L f √ ε (1 − ˜ f ) + K ∂∂ε f √ ε ! , (7a) γ = r m π e ln Λ ee , (7b) L = Z ε f d ε ′ , (7c) K = Z ε f (1 − ˜ f ) ε ′ d ε ′ + ε / Z ∞ ε f (1 − ˜ f ) ε ′− / d ε ′ , (7d)The two logarithmic terms ln Λ ei and ln Λ ee in Eqs. (6) and (7) arise from the integral of theCoulomb di ff erential cross section over all impact parameters. It is well-known that this integraldiverges at both the lower and upper limits, so cut-o ff s must be imposed. It is often assumed, in thecase of a dense plasma, that the lower cut-o ff is determined by the thermal deBroglie wavelength,and the upper cut-o ff by Coulomb screening. However, at high enough density, the interparticleseparation can be comparable or greater than the screening length, and therefore is also a suitablechoice for the upper cut-o ff . Di ff erent choices of these cut-o ff s result in di ff erent formula for ln Λ (see [6, 33] and the references therein). For the purpose of this work, we will assume that ln Λ ei andln Λ ee are constant and focus on the numerical discretization of the collision operators. The qFPEq. (5) must be solved subjected to the boundary condition that the particle fluxes J ei and J ee vanishat ε = ∞ . The blocking factors (1 − ˜ f ) in (6) and (7) are direct consequences of the Pauliexclusion principle. These blocking factors increase the complexity and non-linearity of the kinetic4quation, i.e., the ei flux is at most quadratic and the ee flux cubic in f . Similar to the classicalcase, the qFP collision operators (6) and (7) conserve density and energy. The flux divergenceform of the collision operators and imposed boundary conditions guarantee conservation of totalelectron density for both ei and ee collisions. As a result, most finite volume discretizations easilysatisfy this condition. It is more di ffi cult to achieve energy conservation since it involves a high-order moment quantity. For ei collisions, the statement of energy conservation means that the totalenergy (electron + ion) is constant in time:32 N i ∂ T i ∂ t + Z ∞ ε ∂ J ei ∂ε d ε = . (8)Hence, to impose energy conservation, one can update the ion energy (or temperature) accordingto the discrete form of the integral in Eq. (8). Energy conservation for the ee collision operator isexpressed as: Z ∞ ε ∂ J ee ∂ε d ε = . (9)It is less straightforward to satisfy this condition at the discrete level, even for the classical FPequation. We address this issue in the next section.
3. Numerical Method
In this section, we describe a DG method for solving the qFP equation (5). Let us discretizethe solution domain into N non-overlapping energy groups, ε ∈ (cid:2) ε i − / , ε i + / (cid:3) = Ω i , where i ± / f as an expansionin terms of an orthonormal set of basis functions U p ( ε ): ε ∈ Ω i : f i ( ε, t ) = p max X p = ˆ f i , p ( t ) U p ( ε ) . (10)It is convenient to introduce a mapping from Ω i to [ − ,
1] using the transformation x = ε − ε i ) ∆ i ,where ε i denotes the energy at the center of group i and ∆ i the width of the group. Hereafter tosimply the expressions, we shall occasionally interchange the variables x and ε . In this work wechoose the normalized Legendre polynomials for the basis functions, i.e., U p ( x ) = q p + P p ( x )where P p ( x ) is the regular Legendre polynomial. The orthogonality relationship for the basisfunctions is: Z − U p ( x ) U q ( x ) dx = δ pq . (11)From (10) and (11), we can show that any discrete p -th order moment of the distribution functioncan be evaluated exactly as linear combination of expansion coe ffi cients ˆ f q ( q ≤ p ). For example,5t is straightforward to show that number density and energy of a group are linear combinations ofonly the first two coe ffi cients: ¯ n i = Z Ω i f d ε = ∆ i √ f i , , (12a)¯ e i = Z Ω i f ε d ε = ∆ i ε i √ f i , + ∆ i √ f i , , (12b)where R Ω i denotes an integral over group i . Eq. (12) indicates that the time evolution of numberdensity and energy of the system can be correctly described if the expansion in Eq. (10) is car-ried to at least second order ( p max = p max = J . The di ff erence in the treatment of ei and ee collisions will be pointed out whereverappropriate. The time rate of change of the expansion coe ffi cients is obtained by multiplyingEq. (5) by the basis function and integrating over group i : ∆ i d ˆ f i , p dt = h ˆ J i + / U p (1) − ˆ J i − / U p ( − i − Z Ω i J dU p d ε d ε. (13)Here ˆ J is the numerical flux at an interface, which can be decomposed into a convective anddi ff usive parts: ˆ J i + / = ˆ J Ci + / + ˆ J Di + / . (14)The convective and di ff usive fluxes for ei collision term are: h ˆ J Ci + / i ei = γ ei ε − / i + / f Ci + / (cid:16) − ˜ f Ci + / (cid:17) , (15a) h ˆ J Di + / i ei = γ ei T i ∂∂ε f D √ ε ! i + / . (15b)The convective and di ff usive fluxes for ee collision term are: h ˆ J Ci + / i ee = γ L i + / ε − / i + / f Ci + / (cid:16) − ˜ f Ci + / (cid:17) , (16a) h ˆ J Di + / i ee = γ K i + / ∂∂ε f D √ ε ! i + / . (16b)6n Eqs. (15) and (16), we introduce average values of the distribution function at an interface forconvective and di ff usive fluxes, denoted by the superscript C and D respectively. For convectivefluxes, the average value at the interface is defined using the Chang-Cooper method [9]: f Ci + / = δ i + / f − i + / + (cid:0) − δ i + / (cid:1) f + i + / , (17)where f − i + / denotes the value of f evaluated at the right boundary of cell i and f + i + / the value of f evaluated at the left boundary of cell i +
1. The weighting coe ffi cient δ i + / is defined as: δ i + / = ω i + / − e ω i + / − , (18)where ω i + / is related to the ratio of the convective to the di ff usive coe ffi cient. Note that f C isdefined di ff erently for ei and ee collisions. Here we use a version of Chang-Cooper average forenergy grid similar to Buet and Le Thanh [35]. For ei collision, ω i + / = ε i + − ε i T i and for ee colli-sion, ω i + / = ( ε i + − ε i ) L i + / K i + / . We point out that the blocking factor is omitted from these definitions,because it directly depend on the solution at the interface. This simplification does not introducesignificant error when the distribution is weakly degenerate. For a degenerate distribution, wewould expect some errors in the energy range ε . µ + T e . However, a Chang-Copper type ofaverage in this range is not justifiable so any averaging of the form (17) where δ i + / ∈ [0 ,
1] is rea-sonable. We emphasize that this simplification is only applied to the calculation of the weightingfactors ω i + / and δ i + / .We note that a standard DG scheme does not guarantee positive distribution inside all groupswhen p max ≥
1. This is a known issue for the use of DG schemes to solve kinetic equations,where the particle distribution always remains a non-negative scalar. Application of a positivitypreserving limiter [34] can alleviate this problem but for the case of FP equation it will destroyconservation properties. The adaptation of the Chang-Cooper flux, motivated by its use in stan-dard FV schemes for FP equation, is a way to mitigate this problem, by taking into account theexponential fall o ff of the distribution to avoid overestimating the flux.The average value of f for the di ff usive flux is calculated according to the recovery-based DGscheme of van Leer et al. [36]. The procedure is briefly summarized here. For each interface i + /
2, we first recover a polynomial g ≡ g ( f i , f i + ) that is continuous across two adjacent cells i and i +
1, i.e., over Ω i ∪ Ω i + , from a L minimization: Z Ω i ( g − f i ) U p d ε = , (19a) Z Ω i + ( g − f i + ) U p d ε = . (19b)If f is represented by a polynomial of order p max , the recovered polynomial g is of order 2 p max + Z Ω i J dU p d ε d ε ≃ X n w n J ( x n ) dU p dx , (20)7here n denotes a quadrature point and w n its weight. The approximation in Eq. (20) is used for eicollisions. For ee collisions, it is more advantageous to rewrite Eq. (20) for p = p max +
1, while the order of accuracy of thedi ff usive term, using the recovery-based approach, is 2 p max + [36]. This exponential accuracy isonly obtained in one dimension which is the case here. For higher dimensions, the order of ac-curacy for recovery-based DG schemes is 2( p max + p max =
1, we can expect overall second orderaccuracy for the solution of the qFP equation. This will be checked numerically in Sec. 4.1.
It is important to ensure that the numerical method respect any conservation law associatedwith the collision operators. It is easy to show that the present numerical scheme is density con-serving. Summing Eq. (13) for p = dN e dt = N X i = d ¯ n i dt = ˆ J N + / − ˆ J / . (21)Since the imposed boundary conditions are that particle fluxes vanish at domain boundaries, thesum in Eq. (21) is exactly zero. This is true for both ei and ee collision operators. For ee collisions,we can show that total electron energy is conserved in the continuous limit. At the discrete level,the time evolution of the total energy can be written as: " dE e dt ee = N X i = d ¯ e i dt = ˆ J ee , N + / ε N + / − ˆ J ee , / ε / − X i Z Ω i J ee d ε. (22)Again the flux terms are zeros due to the imposed boundary conditions. To conserve total energy,it is essential to ensure that the discrete sum in Eq. (22) is exactly zero. Note that the integralin Eq. (22) shows up on the right hand side of Eq. (13) for p = J ee into the integral in Eq. (22) and integrating by part, we obtain: Z Ω i J ee d ε = γ h K ε − / f i i + / i − / − γ Z Ω i " ∂ K ∂ε − L (1 − ˜ f ) ε − / f d ε. (23)Let us examine the second term on the right hand side of Eq. (23). We first di ff erentiate Eq. (7d)with respect to ε to obtain an expression for ∂ K ∂ε . The derivative of the integrals can be evaluatedusing the first fundamental theorem of calculus. After some simplification, we arrive at: ∂ K ∂ε = ε / Z ∞ ε f (1 − ˜ f ) ε ′− / d ε ′ . (24)8o simplify the notation, let us define h ( ε ) ≡ − ˜ f ( ε ). Using (24) and (7c), the second term on theright hand side of Eq. (23) becomes: Z Ω i J ee d ε ∝ − γ Z Ω i " f Z ∞ ε f h ε ′− / d ε ′ − f h ε − / Z ε f d ε ′ d ε = − γ Z Ω i " f Z ∞ ε f h ε ′− / d ε ′ + ∂∂ε Z ∞ ε f h ε − / ! Z ε f d ε ′ ! d ε. (25)The second term on the right hand side of (25) can be integrated by part. This results in two terms,one of which cancels the first term on the right hand side of (25) exactly. Hence, we arrive at: Z Ω i J ee d ε ∝ − γ " Z ∞ ε ε ′ / f h d ε ′ ! Z ε f d ε ′ ! i + / i − / . (26)Putting this result back to Eq. (23), we obtain: Z Ω i J ee d ε = γ h K ε − / f i i + / i − / − γ " Z ∞ ε ε ′ / f (cid:16) − ˜ f (cid:17) d ε ′ ! Z ε f d ε ′ ! i + / i − / . (27)We can see that the rate of change in total energy expressed in Eq. (22) with the integral given byEq. (27) is not exactly zero. When summing over all groups, there is a single term remained in theexpression: N X i = d ¯ e i dt = − γ K N + / ε − / N + / f N + / (28)Eq. (28) gives the error in total energy at the discrete level. This is a direct consequence of thefinite truncation of the solution domain. In order to exactly conserve total energy, we propose thefollowing approximate expression for the flux integral: Z Ω i J ee d ε ≃ γ h K ε − / (cid:0) f − f N + / (cid:1)i i + / i − / − γ " Z ∞ ε ε ′ / f (cid:16) − ˜ f (cid:17) d ε ′ ! Z ε f d ε ′ ! i + / i − / , (29)where f N + / is the solution at the last interface. This additional term acts as a correction to theenergy flux to ensure conservation. Another acceptable approximation is: Z Ω i J ee d ε ≃ γ h K (cid:16) ε − / f − ε − / N + / f N + / (cid:17)i i + / i − / − γ " Z ∞ ε ε ′ / f (cid:16) − ˜ f (cid:17) d ε ′ ! Z ε f d ε ′ ! i + / i − / . (30)It is straightforward to verify that the sum of Eq. (29) or (30) over all energy groups is exactlyzero. Hence, total energy is exactly conserved. 9 .3. Time Integration The semi-discrete system (13) is integrated by a backward Euler method. It can be put into thefollowing form: ˆ f n + i , p − ˆ f ni , p ∆ t = R i , p ( ˆ f n + ) , (31)where R i , p contains the right hand side of Eq. (13) evaluated at the next time step. Note that thefactor ∆ i / R i , p . Linearizing Eq. (31) about thesolution at the current time step leads to: R i , p ( ˆ f n + ) ≃ R i , p ( ˆ f n ) + X j , q ∂ R i , p ∂ f j , q (cid:16) ˆ f n + j , q − ˆ f nj , q (cid:17) , (32)where P j , q denotes a sum over all group and polynomial indices. Since the convective and di ff usivecoe ffi cients for ee collisions involve integrals of the distribution function over all energy groups,the Jacobian ∂ R i , p ∂ f j , q is generally a dense matrix. Appendix B gives detailed analytical expressionsfor the these coe ffi cients and their derivatives with respect to the solution vector, which are usedto construct the Jacobian. We also describe in Appendix B an e ffi cient method to evaluate thesecoe ffi cients and derivatives during the computation.Substituting Eq. (32) into (31), we obtain a linear system to be solved at each time step: X j , q ∆ t − δ i j δ pq − ∂ R i , p ∂ f j , q ! (cid:16) ˆ f n + j , q − ˆ f nj , q (cid:17) = R i , p ( ˆ f n ) , (33)where δ i j and δ pq are Kronecker deltas. We note that inverting Eq. (33) constitutes the mostcomputationally intensive part of the calculation since this operation scales as O ( N ) where N dof = N ( p max + ff errors. The results shown in this paper are obtained from solving the linearized system (33).We note that the backward Euler method is first order accurate in time. Since the method isfully implicit, we can relax the stability constraint and only impose time step restriction based onaccuracy requirement. The time step size is estimated by limiting maximum rate of change ofdensity within an energy group using the following expression: ∆ t new = min ≤ i ≤ N ¯ n n + i >β N e α ∆ t | − ¯ n ni / ¯ n n + i | , (34)where ∆ t new is an estimated time step for the next cycle. The min operation is only performed overgroups that have non-negligible densities. Here we set α = . β = − . For all simulationsshown in this work, the time step is estimated using Eq. (34). The only time we use constant timestep size is at the beginning of Sec. 4.1 to check numerical convergence.10 . Simulation Results In this section, we present several numerical tests to validate the numerical method. We firstconsider the thermalization of a monoenergetic distribution to demonstrate the order of accuracyof the scheme. Here only the ee collision term is examined. Since both ee and ei collision termsare discretized the same way, we expect similar order of accuracy. The electron distribution in unitcm − -eV − is initialized as follows: f ( ε, t = = √ π N e δ h exp " − ( ε − ε h ) δ h (35)where N e = cm − , ε h = . δ h = .
12 eV. This problem is similar to the one fromEpperlein[10], but here the electrons are degenerate due to high density and low temperature. Atthis density, the Fermi energy E F = .
34 eV is comparable to the total energy E e = . τ − ee = e m π ~ + e − µ ⋆ / T ⋆ e ln Λ ee , (36)where µ ⋆ / T ⋆ e = . T ⋆ e = .
28 eV are the equilibrium values calculated from the initialdistribution. In the classical limit, we recover the standard expression for τ ee [37].We use an energy grid from 0 to 3 eV with constant energy spacings. Fig. 1 shows the typicalevolution of the occupation number ˜ f ∝ f / √ ε at di ff erent times as the distribution relaxes towarda Fermi-Dirac distribution. The classical Maxwellian distribution is also shown to highlight thedi ff erence in the equilibrium spectrum when degeneracy is taken into account. The numericalresults in Fig. 1 are obtained using 256 groups and a constant time step size of ∆ t = . τ ee . Tocheck the order of accuracy of the scheme, Figs. 2a and 2b show the L p ( p = , , ∞ ) norms of theerror, defined as the di ff erence between the numerical and a reference solution of the distribution.Fig. 2a shows the error as a function of the number of energy groups N using a constant ∆ t = . τ ee . The reference solution is calculated at t = . τ ee using N = N = t = . τ ee using ∆ t = − τ ee . The results from Figs. 2a and 2b clearlyindicate second order convergence in energy group size and first order temporal convergence.The next test case is a more challenging problem requiring a nonuniform energy grid for ef-ficiency. This problem is a simple model describing thermalization of excited electrons in metalsproduced by x-ray absorption [24]. The relaxation of these high energy electrons due to inter-action with the conduction electrons are modeled using the qFP equation. Here the conductionelectrons are treated as a uniform electron gas at zero temperature. The electron distribution inunit cm − -eV − is initialized as: f ( ε, t = = π mh ! / √ ε ˜ f ( ε ) , (37a)˜ f ( ε ) = ε ≤ E F E / F x h √ πδ h √ ε exp (cid:20) − ( ε − ε h ) δ h (cid:21) if ε > E F , (37b)11 F ̃ f t/τ ττ = 0t/τ ττ = 0.1t/τ ττ = 0.3t/τ ττ = . ̃t/τ ττ → ∞→ermi∞axwell Figure 1: Relaxation of a monoenergetic electron energy distribution functions at di ff erent times indicated by di ff erentcolors. The symbols represent the expected equilibrium Fermi-Dirac distribution while the dashed line is the classicalequilibrium Maxwellian distribution −3 −2 −4 −3 −2 −1 L , L , L ∞ ∞N −2 L L L ∞ (a) −3 −2 Δt/τ ee −3 −2 τ ,Δ τ ,Δ τ ∞ ∞Δtτ τ τ ∞ (b) Figure 2: The L , L and L ∞ errors of the distribution function for the problem shown in Fig. 1 as function of (a)number of groups N and (b) time step size ∆ t . E F = x h = . × − , ε h = δ h =
10 eV. The initial distribution consists of apopulation of cold dense electrons at zero temperature and a small fraction of high energy electronswith a total electron density of 8 . × cm − . The high energy electrons are initialized as aGaussian distribution centered at 1 keV. The ee mean collision time is defined using Eq. (36) with µ ⋆ / T ⋆ e = .
66 and T ⋆ e = .
36 eV. In this test, ei collision is turned o ff .In this simulation, we use an energy grid with nonuniform spacing and clustered points near theGaussian to ensure that the particle distributions at low and high energies are adequately resolved.The main strategy is to divide the physical domain into a low and high energy portion. For eachportion, we utilize the formula describe in Appendix C to generate a grid that has grid spacingssuccessively increased by a constant ratio. The energy grid for this problem is defined as: ε i + / = , for i = ε i − / + r i − ∆ ε , for i = , , · · · , M ε i − / + r i − ∆ ε , for i = M + , M + , · · · , N (38)We ran the simulation with successively refined grids (64, 128 and 256 energy groups). Table 1lists values of M , r and r for di ff erent values of N . N M r r
64 52 1 .
18 1 . .
05 1 . .
02 1 . Table 1: Grid generation parameters for simulations shown in Figs. 3 - 6
In Eq. (38), ∆ ε is calculated from Eq. (C.4) using N ε = M , ε min = − eV and ε max = ∆ ε is calculated using N ε = N − M , ε min =
900 eV and ε max = ff er from numerical noises, and the features in the high energytail of the EEDF are better resolved with the qFP approach. At steady-state, the EEDF’s fromboth methods approach the correct equilibrium limit. To demonstrate that the numerical resultsare grid-independent, Fig. 4 shows the distributions at di ff erent times using three di ff erent of grids(64, 128 and 256 groups). It is evident that the results obtained with 128 groups are already veryclose to the converged solution. Fig. 5 shows that the relative H function, calculated from thedistribution, is always decreasing and approaching its minimum at equilibrium. Here we define13 ε/E F −7 −6 −5 −4 −3 −2 −1 ̃ f t/τ ττ = 0t/τ ττ = 200t/τ ττ = 400t/τ ττ = 500t/τ ττ → ∞Fermĩ→2.36̃eV∞ (a) ε/E F −7 −6 −5 −4 −3 −2 −1 ̃ f (b) Figure 3: Thermalization of high energy electrons in a uniform electron gas at zero temperature. N e = . × cm − , E F =
7. Time evolution of the EEDF using (a) qFP equation and (b) Monte Carlo collision method. Di ff erentcolors indicate di ff erent times, and the symbols represent an exact Fermi-Dirac distribution at µ/ T e = .
66 and T e = .
36 eV. The qFP simulation uses 256 groups, and the Monte Carlo simulation uses 2 million particles. ε/E F −7 −6 −5 −4 −3 −2 −1 ̃ f Figure 4: Electron energy distribution functions at t /τ ee = ff erent energy grids (64,128 and 256 groups). Di ff erent colors indicates di ff erent times similarly to Fig. 3. Di ff erent linestyles indicate resultsfrom di ff erent grids. The solid and dashed lines are almost indistinguishable.
200 400 600t/τ ee −2 −1 τ Figure 5: Evolution of the relative H function. the relative H function similar to [22]: H = N − e Z ∞ f ln " ˜ f − ˜ f − ˜ f ⋆ ˜ f ⋆ d ε, (39)This result is consistent with H-theorem [22, 6]. Fig. 6 shows the normalized energy errors bothglobally (solid) and per time step (dashed). The errors are extremely small ( < − ), confirmingthat the numerical scheme is energy conserving. This also illustrates the point that the linearizedsystem (33) is a robust approximation to (31). The next two test cases model electron-ion equilibration process in a dense hydrogen plasma( Z =
1) using both ei and ee collision operators. For simplicity, we assume constant ion tem-perature, so we expect that at steady state the EEDF becomes a Fermi-Dirac distribution at thistemperature. The first test is a heating problem where the distribution is evolved from a degenerateto a classical distribution. The second test is the opposite where the EEDF goes from a classicalto a degenerate distribution. The energy grid are generated according to the formula: ε i + / = , for i = (cid:16) ε max ε min (cid:17) N ε i − / , for i = , , · · · , N (40)where N is the number of groups, ε min = − eV and ε max = eV. Both simulations use thesame energy grid consisting of 64 groups.In the limit where the ee collision time is much faster than ei equilibration time, the EEDFcan be approximated as a Fermi-Dirac distribution at all times (instantaneous thermalization) and15 −1 t/τ ττ −16 −15 −14 −13 e n e r gy e rr o r Figure 6: Total energy error (solid) and error per time step (dash). the ei energy transfer rate can be obtained from direct evaluation of the energy moment of the fullqFP collision integral [6]. In many situations, we can assume that the ion thermal velocity is muchsmaller than the electron thermal velocity ( mT i MT e ≪ dE e dt = − N e ( T e − T i ) τ Eei , (41a)( τ Eei ) − = mM Ze m π ~ + e − µ/ T e ln Λ ei . (41b)Eq. (41) was first derived by Brysk in 1974 [38]. In the classical limit, the relaxation time τ Eei approaches the Spitzer form ( τ Eei ∝ T − / e ). In the opposite limit, it becomes a constant.For the heating test, the initial electron distribution is a Fermi-Dirac distribution at N e = cm − and T e = E F = . µ/ T e = .
2. For the ions, we set N i = N e and T i =
60 eV. Since the electron distribution is dynamically evolved, the correspondingei collision time also changes as a function of time. For convenience, we define a reference timebased on the ei collision time evaluated at initial time τ Eei = τ Eei ( t = Λ ei = ln Λ ee . Figs. 7a and 7b show the evolution of T e , µ and the EEDF during the heatingprocess. It is evident that the temperature evolution from the qFP model (solid line) agrees verywell with Brysk’s model (symbols). This is expected since for Hydrogen, τ ee /τ Eei ∝ m / M ≪ µ is decreasing in time. The steady-state distribution (thickest line in Fig. 7b) matches theanalytical Fermi-Dirac distribution at 60 eV (symbols).Fot the cooling test, we initialize the electron distribution as a Fermi-Dirac distribution at N e = cm − and T e =
60 eV, which corresponds to µ/ T e = − .
9. For the ions, we set N i = N e τei0 T e / τ F −202468 μ / T e (a) F ̃ f t/τ Eτi0 ̃=̃0.01t/τ
Eτi0 ̃=̃0.1t/τ
Eτi0 ̃=̃0.5t/τ
Eτi0 ̃=̃1t/τ
Eτi0 ̃=̃22Fermĩ(60̃eV) (b)
Figure 7: Evolution of (a) T e and µ and (b) EEDF during heating test. The symbols in (a) are solutions obtainedfrom Brysk’s model (Eq. (41)), while the dashed-doted line indicates the equilibrium temperature ( T i =
60 eV). Thesymbols in (b) correspond to an exact Fermi-Dirac distribution at equilibirum temperature. τei0 T e / τ F −202468 μ / T e (a) F ̃ f t/τ Eτi0 ̃=̃0.1t/τ
Eτi0 ̃=̃0.5t/τ
Eτi0 ̃=̃0.8t/τ
Eτi0 ̃=̃1t/τ
Eτi0 ̃=̃2Fermĩ(5̃eV) (b)
Figure 8: Evolution of (a) T e and µ and (b) EEDF during cooling test. The symbols in (a) are solutions obtained fromBrysk’s model (Eq. (41)), while the dashed-doted line indicates the equilibrium temperature ( T i = .0 0.5 1.0 1.5 2.0 2.5 3.0ε/E F ̃ f Figure 9: Steady state electron energy distribution functions for the cooling problem obtained from four di ff erentenergy grids (16, 32, 64 and 128 groups). and T i = T e , µ and the EEDF during the cooling process. Similar to the heating test, the temperature evolution isin very good agreement with Brysk’s results, and the distribution approaches to the correct Fermi-Dirac equilibrium at 5 eV. The steady-state distribution shows that the states below the Fermienergy are almost fully occupied, which resembles a strongly degenerate distribution. Fig. 9 showsthe steady state distributions from four di ff erent grids (16, 32, 64 and 128 groups), and indicatesthat the results using 64 groups are grid-independent. Although not yet converged, the resultsobtained from 16 and 32 groups are very reasonable, given the large energy range we need tocover to resolve the initial distribution. This highlights the advantage of nonuniform griddingcapability for problems with large temperature separation.
5. Conlusions
We construct a conservative numerical scheme for the isotropic quantum Fokker-Planck equa-tion describing the evolution of degenerate electrons due to electron-ion and electron-electronelastic collisions. A discontinuous Galerkin method is employed to discretize the collision oper-ators and the electron energy distribution is time integrated using a backward Euler method. Thehigh-order discontinuous Galerkin discretization o ff ers additional degrees of freedom to describeenergy transfer and enforce energy conservation. The fully implicit time integration relaxes therestriction on the time step, making it feasible to include other physics occurring at the slowertime scales. The extension to include inelastic collisions is planned for future work. We anticipatethat this capability will be useful for modeling electron collisions in high energy density physicsexperiments. 18 cknowledgement I would like to thank the anonymous referees whose comments / suggestions helped improveand clarify this manuscript. I wish to thank H.A. Scott for many helpful discussions. This workwas performed under the auspices of the U.S. Department of Energy by Lawrence LivermoreNational Laboratory under Contract DE-AC52-07NA27344.This document was prepared as an account of work sponsored by an agency of the UnitedStates government. Neither the United States government nor Lawrence Livermore National Se-curity, LLC, nor any of their employees makes any warranty, expressed or implied, or assumes anylegal liability or responsibility for the accuracy, completeness, or usefulness of any information,apparatus, product, or process disclosed, or represents that its use would not infringe privatelyowned rights. Reference herein to any specific commercial product, process, or service by tradename, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorse-ment, recommendation, or favoring by the United States government or Lawrence Livermore Na-tional Security, LLC. The views and opinions of authors expressed herein do not necessarily stateor reflect those of the United States government or Lawrence Livermore National Security, LLC,and shall not be used for advertising or product endorsement purposes. Appendix A. Collision Operators
In this appendix, we formulate the ei and ee collision operators for an isotropic EEDF. Fol-lowing the work of Brown and Haines [5], the qFP equation for an isotropic electron velocitydistribution f v ( v ) due to collisions with another particle distribution F v ( V ) (either electron or ion)can be written as: 1 Y ∂ f v ∂ t = v ∂∂ v " mM f v (1 − ˜ f v ) I + v (cid:0) I − I ′ + J − − J ′− (cid:1) ∂ f v ∂ v , (A.1)where Y = π Z e m ln Λ . Here Z and M denote charge and mass of the colliding partner, and ln Λ iswritten generically for both ei and ee collisions. Also ˜ f v ≡ h m g f v is the mean occupation numberwhere g is the spin degeneracy. The set of integrals in Eq. (A.1) are defined as follows: I p = π v p Z v V p + F v ( V ) dV , I ′ p = π v p Z v V p + F v ( V ) ˜ F v ( V ) dV , (A.2a) J p = π v p Z ∞ v V p + F v ( V ) dV , J ′ p = π v p Z ∞ v V p + F v ( V ) ˜ F v ( V ) dV . (A.2b)The velocity distribution f v ( v ) is related to the EEDF f ε ( ε ) via the change of variable f ε ( ε ) d ε = π v f v ( v ) dv ( ε = mv / ∂ f ε ∂ t = m / Y √ ∂∂ε " mM f ε √ ε (1 − ˜ f ε ) I + ε (cid:0) I − I ′ + J − − J ′− (cid:1) ∂∂ε f ε √ ε ! (A.3)where ˜ f ε = ˜ f v . 19hroughout this work the ions are assumed to be non-degenerate and thermalized. If we makeanother assumption that the ion thermal velocity is much smaller than the electron velocity, thenthe set of integrals in (A.2) for ei collisions reduces to: I → N i , (A.4a) I → N i T i Mv , (A.4b) J − → , (A.4c) I ′ , J ′− → , (A.4d)where N i and T i are the ion density and temperature.Inserting these results into Eq. (A.3), the ei collision operator for the EEDF can be written as: ∂ f ε ∂ t ! ei = ∂∂ε γ ei " f ε √ ε (cid:16) − ˜ f ε (cid:17) + T i ∂∂ε f ε √ ε ! , (A.5)where γ ei = mM q m π N e Ze ln Λ ei . Here we also apply the charge neutrality condition, i.e., N e = ZN i .For ee collisions, we first rewrite the integrals (A.2) in term of f ε : I = π Z v f v v ′ dv ′ = Z ε f ε d ε ′ , (A.6) I − I ′ = π v Z v f v (cid:16) − ˜ f v (cid:17) v ′ dv ′ = ε − Z ε ε ′ f ε (cid:16) − ˜ f ε (cid:17) d ε ′ , (A.7) J − − J ′− = π v − Z ∞ v f v (cid:16) − ˜ f v (cid:17) v ′ dv ′ = ε / Z ∞ ε ε ′− / f ε (cid:16) − ˜ f ε (cid:17) d ε ′ . (A.8)Let us define the following set of integrals: L = I , (A.9a) K = ε (cid:0) I − I ′ + J − − J ′− (cid:1) . (A.9b)The ee collision operator can then be written as: ∂ f ε ∂ t ! ee = γ ∂∂ε " L f ε √ ε (1 − ˜ f ε ) + K ∂∂ε f ε √ ε ! , (A.10)where γ = q m π e ln Λ ee . For ease of notation, in the main text we drop the subscript ε and consis-tently refer to f as the electron energy distribution function. Appendix B. Fast Evaluation of Electron-electron Flux Coe ffi cients K and L In this appendix, we give discrete expressions of the coe ffi cients K i + / and L i + / appeared inthe expression of the ee collision fluxes (Eq. (16)), and describe an e ffi cient method to evaluate20hem during the computation. The fact that these coe ffi cients are non-local gives rise to the denseJacobian matrix in Eq. (33). Analytical expressions of their derivatives with respect to the solu-tion vector are also shown here since they are required for constructing the Jacobian matrix. Forconvenience, let us define the following group-wise moment of order k : M ( k ) i = Z Ω i f (1 − ˜ f ) ε ′ k d ε ′ . (B.1)The evaluation of Eq. (B.1) with the Pauli blocking factor is nontrivial. A practical way is toapproximate the integral using Gaussian quadrature: M ( k ) i ≃ ∆ i X n [ i ] p max X p = w n ˆ f i , p ε ki , n U p ( ε i , n ) − C p max X q = ˆ f i , q √ ε i , n U q ( ε i , n ) , (B.2)where C = π (cid:16) h m (cid:17) / . Here n [ i ] denotes the quadrature point n with weight w n and energy ε i , n lying within the interval [ ε i − / , ε i + / ]. We point out that Eq. (B.1) does not need to be reevaluatedat every computational cycle. Instead, for each group i we can define two sets of coe ffi cient S ( k ) i , p and T ( k ) i , pq ( p , q = , · · · , p max ): S ( k ) i , p = ∆ i X n [ i ] w n ε ki , n U p ( ε i , n ) , (B.3a) T ( k ) i , pq = C ∆ i X n [ i ] w n ε k − / i , n U p ( ε i , n ) U q ( ε i , n ) . (B.3b)Note that these coe ffi cients only need to be calculated once at the beginning. During the calcula-tion, we can utilize these coe ffi cients to update the moment: M ( k ) = p max X p = S ( k ) p ˆ f p − p max X q = T ( k ) pq ˆ f p ˆ f q . (B.4)Here we have omitted the group index for ease of notation. The derivatives of M ( k ) with respect tothe expansion coe ffi cient ˆ f p can be evaluated as follows: ∂ M ( k ) ∂ ˆ f p = S ( k ) p − p max X q = (1 + δ pq ) T ( k ) pq ˆ f q . (B.5)Using (B.1) and (12a), the coe ffi cients K i + / and L i + / defined in Eq. (7) can now be writtenas: K i + / = i X j = M (1) j + ε i + / ) / N X j = i + M (1 / j , (B.6a) L i + / = √ i X j = ∆ j ˆ f i , . (B.6b)21he derivatives of K i + / and L i + / with respect to a group coe ffi cient ˆ f j , p can be evaluated as: ∂ K i + / ∂ ˆ f j , p = ∂ M (1) j ∂ ˆ f j , p , for j ≤ i ε i + / ) / ∂ M (1 / j ∂ ˆ f j , p , otherwise (B.7a) ∂ L i + / ∂ ˆ f j , p = δ p √ ∆ j , for j ≤ i , otherwise (B.7b)We remark that since K i + / and L i + / depend on solution of all energy groups, the Jacobiancontaining the derivatives ∂ K i + / ∂ ˆ f j , p and ∂ L i + / ∂ ˆ f j , p is a dense matrix. Appendix C. Energy grid generation
We define an energy grid that has the grid spacings increased by a constant ratio r : ε i + / = ε i − / + r i − ∆ ε for i = , , · · · , N ε (C.1)where N ε is the number of energy groups. By recursively applying (C.1), we obtain an expressionrelating the first and last grid points: ε N + / = ε / + ∆ ε N ε − X i = r i (C.2)By rewriting P N ε − i = r i − = (1 − r N ε ) P ∞ i = r i = (1 − r N ε ) / (1 − r ), we obtain an expression for deter-mining ∆ ε : ∆ ε = ( ε N ε + / − ε / )(1 − r )1 − r N ε (C.3)Hence for a energy grid consisting of N ε groups spanning from ε min to ε max , we can estimate theinitial grid spacing ∆ ε by: ∆ ε = ( ε max − ε min )(1 − r )1 − r N ε (C.4)In practice, we choose ε min to be a small but non-zero quantity to retain resolution in the lowenergy portion of the spectrum. However, the physical domain should always start from zero, soone can either add an zero energy grid point or simply set the first value of the energy grid to zero. References [1] I. P. Shkarofsky, T. W. Johnston, M. P. Bachynski,
The Particle Kinetics of Plasmas (Addison-Wesley,, 1966).[2] D. Kremp, M. Schlanges, W.-D. Kraeft,
Quantum statistics of nonideal plasmas (Springer, Berlin, 2005).[3] M. Bonitz,
Quantum kinetic theory (Springer, Cham, 2. ed edition, 2016).[4] P. Danielewicz,
Physica A (1), 167 (1980).[5] S. R. Brown and M. G. Haines,
Journal of Plasma Physics (4), 577 (1997).[6] J. Daligault, Physics of Plasmas (3), 032706 (2016).[7] E. A. Uehling and G. E. Uhlenbeck, Physical Review (7), 552 (1933).
8] J. Daligault,
Phys. Rev. Lett. , 045002 (2017).[9] J. S. Chang and G. Cooper,
Journal of Computational Physics (1), 1 (1970).[10] E. M. Epperlein, Journal of Computational Physics (2), 291 (1994).[11] L. Pareschi, G. Russo, G. Toscani,
Journal of Computational Physics (1), 216 (2000).[12] W. Taitano, L. Chac´on, A. Simakov,
Journal of Computational Physics , 391 (2016).[13] H. P. Le and J.-L. Cambier,
Physics of Plasmas (12), 122105 (2017).[14] A. B. Langdon, Physical Review Letters (9), 575 (1980).[15] E. M. Epperlein and M. G. Haines, Phys. Fluids , 1029 (1986).[16] E. M. Epperlein, G. J. Rickard, A. R. Bell, Phys. Rev. Lett. (21), 2453 (1988).[17] O. Larroche, Physics of Fluids B: Plasma Physics (8), 2816 (1993).[18] A. L. Milder et al. , Phys. Rev. Lett. (2), 025001 (2020).[19] M. Lampe,
Physical Review (1), 306 (1968).[20] S. R. Brown and M. G. Haines,
Journal of Plasma Physics (2), 129 (1999).[21] J. Daligault, Physics of Plasmas (8), 082703 (2018).[22] J. Hu, S. Jin, B. Yan, Communications in Computational Physics (5), 1541 (2012).[23] F. Filbet, J. Hu, S. Jin, ESAIM: M2AN (2), 443 (2012).[24] H. Kitamura, Journal of Electron Spectroscopy and Related Phenomena , 45 (2019).[25] C. R. Scullard et al. , Physics of Plasmas (9), 092119 (2016).[26] C. R. Scullard et al. , Journal of Computational Physics , 109110 (2020).[27] A. Turrell, M. Sherlock, S. Rose,
Journal of Computational Physics , 13 (2013).[28] P. Borowik, J.-L. Thobel, L. Adamowicz,
Journal of Computational Physics , 397 (2017).[29] H. P. Le, M. Sherlock, H. A. Scott,
Phys. Rev. E , 013202 (2019).[30] A. R. Bell, A. P. L. Robinson, M. Sherlock, R. J. Kingham, W. Rozmus,
Plasma Physics and Controlled Fusion (3), R37 (2006).[31] A. G. R. Thomas et al. , Journal of Computational Physics , 1051 (2012).[32] H. M. Antia,
ApJS , 101 (1993).[33] D. O. Gericke, M. S. Murillo, M. Schlanges, Phys. Rev. E (3), 036418 (2002).[34] X. Zhang and C. W. Shu, Proc. R. Soc. A , 2752 (2011).[35] C. Buet and K.-C. Le Thanh, Technical report CEA-R6172 (2008).[36] B. van Leer, M. Lo, M. van Raalte, In
Miami, Florida(2007). American Institute of Aeronautics and Astronautics.[37] J. D. Huba,
NRL PLASMA FORMULARY (Naval Research Laboratory, Washington, DC, 2016).[38] H. Brysk,
Plasma Physics (10), 927 (1974).(10), 927 (1974).