# Combining the hybrid mimetic mixed method with the Scharfetter-Gummel scheme for magnetised transport in plasmas

In this paper, we propose a numerical scheme for fluid models of magnetised plasmas. One important feature of the numerical scheme is that it should be able to handle the anisotropy induced by the magnetic field. In order to do so, we propose the use of the hybrid mimetic mixed (HMM) scheme for diffusion. This is combined with a hybridised variant of the Scharfetter-Gummel (SG) scheme for advection. The proposed hybrid scheme can be implemented very efficiently via static condensation. Numerical tests are then performed to show the applicability of the combined HMM-SG scheme, even for highly anisotropic magnetic fields.

## Authors

• 5 publications
• 3 publications
• 1 publication
• 1 publication
• 1 publication
03/23/2020

### A Generalised Complete Flux scheme for anisotropic advection-diffusion equations

In this paper, we consider separating the discretisation of the diffusiv...
04/09/2021

### Flux correction for nonconservative convection-diffusion equation

Our goal is to develop a flux limiter of the Flux-Corrected Transport me...
02/25/2022

### A Hybrid High-Order scheme for the stationary, incompressible magnetohydrodynamics equations

We propose and analyse a hybrid high-order (HHO) scheme for stationary i...
11/07/2021

### Mixed methods for large-strain poroelasticity/chemotaxis models simulating the formation of myocardial oedema

In this paper we propose a novel coupled poroelasticity-diffusion model ...
06/26/2019

### Lattice Boltzmann method for simulation of diffusion magnetic resonance imaging physics in heterogeneous tissue models

We report the first implementation of the Lattice Boltzmann method (LBM)...
09/01/2020

### A hybrid WENO method with modified ghost fluid method for compressible two-medium flow problems

In this paper, we develop a simplified hybrid weighted essentially non-o...
10/17/2019

### Visualizing the world's largest turbulence simulation

In this exploratory submission we present the visualization of the large...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

Magnetostatic fields play an important role in plasma physics. One useful property of a magnetic field is that it confines the plasma, limiting charged particle loss to the walls; hence, allowing magnetron discharges to operate at relatively low voltages. Some examples of these can be seen in Hall-effect thrusters (HETs) and electron-cyclotron-resonance (ECR) discharges [13]. These physical processes are expensive to set-up, and so we turn to mathematical models and perform numerical simulations in order to get an overview and general understanding of how particles behave in these situations.

Early numerical studies for plasmas involve particle-based models, which involve Monte Carlo simulations [3, 8]. However, the inclusion of a magnetic field complicates the physics and hence the numerical modeling of magnetised discharges. In this paper, we consider fluid models. These consist of drift-diffusion equations for the particles (electrons and ions), coupled to the Poisson equation for the electric field, see e.g.: [6, 13, 14, 17, 23]. This is an approximation for the real physics which is computationally more efficient and avoids some of the numerical complications and high computational costs associated to using Monte Carlo type methods on purely particle-based models [13].

The model system of equations is highly nonlinear due to the coupling between the Poisson equation and the drift-diffusion equations. In order to resolve this, we decouple the equations by using the transversal method of lines, where a semi-discrete system is obtained after discretising with respect to the time variable. Following this, we look at the spatial discretisation. The Poisson equation is elliptic, and can thus be solved by classical techniques, such as the central difference method. In the absence of a magnetic field, the drift-diffusion equations can be solved efficiently and precisely by the Scharfetter-Gummel (SG) scheme [25]. On the other hand, the presence of a magnetic field introduces anisotropy in the drift-diffusion equations; hence requiring a numerical scheme which captures the anisotropic diffusion precisely. Recent studies propose the use of high order finite difference gradient reconstruction methods, together with a magnetic field aligned mesh [24, 29] in order to resolve this.

The aim of this paper is to propose a numerical scheme which is computationally efficient, and is able to handle the anisotropy induced by the magnetic field. In order to do so, we propose, as in [2], to separate the discretisation of the diffusive and advective fluxes. We then use the hybrid-mimetic-mixed (HMM) method [10]

for discretising the diffusive flux, and a modification of the SG scheme for the advective flux. The main interest of considering the HMM method is its ability to handle the anisotropic diffusion tensor. Moreover, the method can be fully hybridised, leading to a purely local stencil for solving the system of equations, which allows a very efficient implementation of the scheme. As a remark, the HMM can also be viewed as a gradient reconstruction method or a gradient scheme

[9], and can thus be used in combination with the methods proposed in [24, 29].

The outline of this paper is as follows. We start by presenting the mathematical model for magnetised transport in Section 2. This model consists of a set of drift-diffusion equations, coupled to the Poisson equation. We then describe in Section 3 the choice of time discretisation, and how to decouple the system. Following this, we employ a finite volume type spatial discretisation in Section 4, using the HMM for the diffusive fluxes, and a modified SG for the advective fluxes. Details of the implementation are then given in Section 5. Numerical tests are provided in Section 6 to show the effectiveness of the proposed method. Finally, in Section 7 we provide a summary of the results and some recommendations for future research.

## 2. Mathematical model

Suppose that an electric field is applied between two parallel electrodes and that an electron flux is forced through a uniform medium, where an oblique magnetic field is applied. The electric field in a spatial domain over a time interval is governed by the Poisson equation

 ∇⋅(εE)=ρonΩ×(0,T), (1a) where ε is the permittivity of vacuum, V is the potential, and ρ is the space charge density given by ρ=q(ni−ne). Here, q is the elementary charge, ni and ne are the ion and electron densities, respectively. Under the assumption that ions are not magnetised, the model which describes the evolution of the ion and electron densities is then given by the drift-diffusion equations:
 ∂ni∂t+∇⋅Γi−kne,locng=0onΩ×(0,T), (1b) ∂ne∂t+∇⋅MeΓe−kne,locng=0onΩ×(0,T), (1c)

where and are the drift-diffusion flux densities for the ions and electrons, respectively, defined by

 Γi :=μiEni−Di∇ni, Γe :=μeEne−De∇ne.

Equations (1b) and (1c) are coupled to the Poisson equation (1a) due to the electric field and the space charge density . The parameters and constants used in the model are enumerated in Tables 1 and 2. For quantities which are common to both ions and electrons, we use the subscript , referring to a particle (either an ion or an electron).

We now discuss the parameters used in the model. Firstly, the ionisation rate and the mobility coefficients are computed via local field approximations, and are hence functions of the reduced electric field . In this paper, we consider helium gas, whose ion mobility (Figure 2, left) is obtained from [11]. The ionisation rate (Figure 1) and electron mobility (Figure 2, right) were computed in [15] using cross-sectional data for helium gas in [7, 20] for low energies and [16] for high energies. Here and throughout the paper, refers to the -norm.

Secondly, we discuss the diffusion coefficients. These are governed by the Stokes-Einstein relation

 Dp=kBTpμpqp,

where, and . For ions, we assume that and so we have

 Di=kBTgμiq. (2a) We assume that electrons follow a Maxwellian distribution so that kBTe=2/3¯ϵ and thus De=−¯ϵ2μe3q. (2b)

Thirdly, the the local electron density is computed by the relation [15]

 ne,loc=∥MeΓe∥μe∥MeE∥. (3)

Finally, the magnetic tensor is given by

 Me =1∥β∥2+1(I+ββT−C(β)), (4a) C(β) =⎡⎢⎣0−β3β2β30−β1−β2β10⎤⎥⎦, with Hall parameters βi=μeBi, i=1,2,3 and magnetic field B. We now note that C(β)is a cross product matrix such that for any vector v∈R3, C(β)v=β×v. Taking note that Meβ =1∥β∥2+1(β+β(βTβ)−β×β) =β, we then see that βis an eigenvector of Mecorresponding to the eigenvalue 1. Essentially, this means that Me preserves the component of the vector that is parallel to the magnetic field.

For dimension , the magnetic field can be specified at an angle with respect to the -axis by setting

 B =[B1,B2,0]T (4b) =∥B∥[cosθ,sinθ,0]T,

where is the strength of the magnetic field. Considering now a vector

 β⊥=μe∥B∥[−sinθ,cosθ,0]T

orthogonal to the magnetic field, we see that

 Meβ⊥=1∥β∥2+1(β⊥−∥β∥2ez).

For two-dimensional models, we neglect the -component, and hence we see that the components orthogonal to the magnetic field are scaled by a factor . These properties agree with the classical definition [18] of a magnetic tensor.

### 2.1. Boundary conditions

We now describe the boundary conditions for the model. In the electrodes (corresponding to either anodes or cathodes), Dirichlet boundary conditions are imposed for the Poisson equation (1a),

 V =Vaanode, (5a) V =Vccathode. Then, flux boundary conditions are imposed for the drift-diffusion equations (1b) and (1c), i.e., Γi⋅n =(μiE⋅n)+ni+14vth,ini, (5b) (MeΓe)⋅n =(μeMeE⋅n)+ne+14vth,ene−γiΓi⋅n, (5c) where vth,p is the thermal velocity, f+=max(f,0) denotes the positive part of the function, and n is the outward unit normal vector at the boundary. The first two terms in the right hand side of equations (5b) and (5c) describe the loss of electrons and ions, respectively, through the wall via the drift and thermal fluxes. On the other hand, the third term in (5c) is the secondary electron emission, which describes the production of electrons as a result of ion fluxes bombarding the wall [4], where γi is the average number of electrons emitted per incident ion.

Finally, homogeneous Neumann boundary conditions are imposed over the other boundaries of the domain. That is, for (1a), (1b), and (1c) we impose:

 ∇V⋅n =0, (5d) ∇ni⋅n =0, ∇ne⋅n =0.

For the drift-diffusion equations, this means that there is no diffusion of ions and electrons across these boundaries. For the Poisson equation, this means that the electric field lines are parallel to the boundary.

Given the electron and ion densities at time , we are then interested in calculating the densities at time . In this work, we consider a numerical scheme based on the transversal method of lines. That is, we start by presenting a semi-discrete version of the model system of equations with respect to the time variable. Following this, we then describe the discretisation in space.

## 3. Time discretisation

We form a partition of the time interval by taking . In order to avoid issues with the stability of the numerical solution, we use the implicit Euler time integration. Using a uniform time step of , for , and given the densities , we seek such that

 (6a)
 n(m+1)i−n(m)iΔt+∇⋅Γ(m+1)i−k(m+1)n(m+1)e,locng=0, (6b)
 n(m+1)e−n(m)eΔt+∇⋅M(m+1)eΓ(m+1)e−k(m+1)n(m+1)e,locng=0. (6c)

Here,

 Γ(m+1)p=μ(m+1)pE(m+1)n(m+1)p−D(m+1)p∇n(m+1)p

is the drift-diffusion flux density of particle (ions or electrons) at time . As a remark, we note that even though the scheme was presented with a uniform time stepping, it can be easily adapted into a scheme with an adaptive time stepping .

We note that due to the electric field, the Poisson equation (6a) is coupled to the continuity equations (6b) and (6c). We now describe a few ways for decoupling the system of equations (6a)-(6c).

### 3.1. First-order approximation of the charge density

The simplest way to decouple the Poisson equation from the continuity equations is to use a first-order approximation of the charge density. This leads to seeking, for , such that

 ∇⋅(εE(m+1))=q(n(m)i−n(m)e). (7)

The advantage of using (7) is that the charge density will then only depend on the electron and ion densities at the previous time step. On the other hand, this is a first-order approximation, which means that the leading order of error in time is . Hence, a very small time step, at most equal to the dielectric relaxation time [27, Section II.C.], is needed in order to ensure that the electric field does not change sign during one time step, and that the approximation of the electric field is good enough to be used for computing the electron and ion densities in the drift-diffusion equations.

### 3.2. Second-order approximation of the charge density

We now present a second-order approximation by performing a Taylor expansion on the charge density. Centering the Taylor expansion at time , we may write

 ρ(x,t(m+1))=ρ(x,t(m))+Δt∂ρ(x,t(m+1))∂t+O(Δt2). (8)

We then note that the time derivative can be computed from the drift-diffusion equations, and can be written as

 ∂ρ∂t=−q(∇⋅Γi−∇⋅MeΓe). (9)

Using now (8) for the charge density in the Poisson equation and evaluating in (9) at time , we then seek, for , such that

 ∇⋅(εE(m+1))= q(n(m)i−n(m)e)−qΔt∇⋅Γ(m+1)i (10) +qΔt∇⋅(M(m+1)eΓ(m+1)e).

In this case, the quantities and at time are unknown and hence the Poisson equation is still coupled to the drift-diffusion equations. One way to deal with this is to consider a semi-implicit treatment of these quantities [27, 28]. That is, we keep implicit and take first-order approximations to , which leads to

 ∇⋅(εE(m+1))= q(n(m)i−n(m)e)−qΔt∇⋅(μ(m)iE(m+1)n(m)i−D(m)i∇n(m)i) +qΔt∇⋅M(m)e(μ(m)eE(m+1)n(m)e−D(m)e∇n(m)e),

or equivalently,

 ∇⋅((εI+qΔt =q(n(m)i−n(m)e)+qΔt∇⋅(D(m)i∇n(m)i) −qΔt∇⋅(M(m)eD(m)e∇n(m)e).

Denoting the correction flux at time level by

 Γ(m+1)c:=(μ(m)in(m)iI−μ(m)en(m)eM(m)e)E(m+1), (11)

the above equation can be rewritten as

 ∇⋅(εE(m+1))+qΔt∇⋅Γ(m+1)c =q(n(m)i−n(m)e)+qΔt∇⋅(D(m)i∇n(m)i) (12) −qΔt∇⋅(M(m)eD(m)e∇n(m)e).

## 4. Spatial discretisation

For the spatial discretisation, we start by forming , a partition of the domain into control volumes. Each of the control volumes has edges such that , and we denote by the total number of edges contained in this partition. We then denote by the collection of unknowns (one on each cell, and one on each edge). For this work, we focus on finite volume type schemes.

### 4.1. Finite volume discretisation

To write the equations (6a),(6b), and (6c) in their finite volume forms, we take their integrals over all control volumes . This results to, in each control volume ,

 ∫K∇⋅(εE(m+1))dA=q∫K(n(m+1)i−n(m+1)e)dA,
 1Δt∫K(n(m+1)i−n(m)i)dA+∫K∇⋅Γ(m+1)idA−ng∫Kk(m+1)n(m+1)e,locdA=0,
 1Δt∫K(n(m+1)e−n(m)e)dA+∫K∇⋅M(m+1)eΓ(m+1)edA−ng∫Kk(m+1)n(m+1)e,locdA=0,

respectively. By applying the divergence theorem and using , these can be rewritten as

 −∑σ∈EK∫σε∇V(m+1)⋅nK,σds=q∫K(n(m+1)i−n(m+1)e)dA, (13a) 1Δt∫K(n(m+1)i−n(m)i)dA +∑σ∈EK∫σΓ(m+1)i⋅nK,σds (13b) 1Δt∫K(n(m+1)e−n(m)e)dA +∑σ∈EK∫σM(m+1)eΓ(m+1)e⋅nK,σds (13c)

Here, is the outward unit normal vector of at . Due to the definition (3) of , the term is nonlinear in . Here, we choose a first-order approximation and take the local electron density at time so that

 n(m)e,loc:=∥M(m)eΓ(m)e)∥μ(m)e∥M(m)eE(m)∥. (14)

Key to the definition of a finite volume scheme is the choice of how to discretise the fluxes. In order to be able to treat the anisotropy, we propose to split the fluxes into a diffusive and an advective component. Diffusive fluxes are of the form for some diffusion tensor and scalar ; on the other hand, the advective fluxes take the form , where is a velocity field. In equations (13a)-(13c), three diffusive fluxes and two advective fluxes are present. We then denote the discrete diffusive and advective fluxes for each edge in a cell by

 FDVK,σ≈−∫σε∇V⋅nK,σds, (15a) FDiK,σ≈−∫σDi∇ni⋅nK,σds, (15b) FDeK,σ≈−∫σMeDe∇ne⋅nK,σds, (15c) FAiK,σ≈∫σniμiE⋅nK,σds, (15d) FAeK,σ≈∫σneμeMeE⋅nK,σds. (15e)

Replacing the fluxes in (13a)-(13c) with their discrete counterparts (15a)-(15e) at the appropriate time level and under the assumption that and are piecewise constant with values and respectively in a cell , we then obtain our numerical scheme. That is, for , given the densities , we seek such that for each ,

 ∑σ∈EKFDV,(m+1)K,σ=q|K|(n(m)i,K−n(m)e,K), (16a) |K|n(m+1)i,KΔt+∑σ∈EKFDi,(m+1)K,σ+∑σ∈EKFAi,(m+1)K,σ=|K|n(m)i,KΔt+ng|K|k(m+1)Kn(m)e,loc,K, (16b) |K|n(m+1)e,KΔt+∑σ∈EKFDe,(m+1)K,σ+∑σ∈EKFAe,(m+1)K,σ=|K|n(m)e,KΔt+ng|K|k(m+1)Kn(m)e,loc,K, (16c) where k(m+1)K is the average ionisation rate in cell K at time t(m+1), and n(m)e,loc,K is the value of n(m)e,loc computed at cell K.

Starting with , we then use the second-order approximation for the charge density in solving the Poisson equation. That is, instead of using (16a), we use a semi-implicit discretisation based on (12), given by

 ∑σ∈EKFDV,(m+1)K,σ +qΔt∑σ∈EKFCV,(m+1)K,σ (16d) =q|K|(n(m)i,K−n(m)e,K)+qΔt(∑σ∈EK(FDe,(m)K,σ−FDi,(m)K,σ)),

where

 FCV,(m+1)K,σ≈−∫σΓ(m+1)c⋅nK,σds

is a discretisation of the semi-implicit correction flux defined in (11).

We note here that we cannot immediately use (16d) at the first time step due to the yet unknown mobility coefficients, correction and diffusive fluxes at time . Hence, we choose to use (16a) for solving the Poisson equation at the first time step when .

### 4.2. Boundary conditions and conservation of fluxes

We now see that there are equations (corresponding to the number of cells) involved in each of (16a)-(16d). However, we have unknowns for each of these equations: corresponding to the number of cells, and corresponding to the number of edges. This is standard for hybrid finite volume methods, and the equations for the unknowns along the edges are obtained by first imposing conservation of fluxes for interior edges. That is, if is an edge shared by two cells and , we should have that the total (diffusive and advective) flux going from to via is equal to the negative total flux going from to via , i.e.

 (FDK,σ+FAK,σ)+(FDL,σ+FAL,σ)=0. (17)

Finally, we impose boundary conditions on the remaining edges, which are located along the boundary of the domain. The Dirichlet and Neumann boundary conditions are straightforward to impose, so we only describe here how to impose the flux boundary conditions (5b) and (5c). This is done by setting, for the corresponding boundary edges ,

 FDi,(m+1)K,σ+FAi,(m+1)K,σ=|σ|((μ(m+1)iE(m+1)⋅nK,σ)+n(m+1)i,σ+14vth,in(m+1)i,σ). (18a) FDe,(m+1)K,σ+FAe,(m+1)K,σ= (18b) −γi(FDi,(m+1)K,σ+FAi,(m+1)K,σ).

Here, and are the values (to be determined) of and at , respectively. To complete the definition of the numerical scheme, we are then left to describe the choice of how to compute the discrete diffusive and advective fluxes.

### 4.3. Diffusive fluxes

In this section, we discuss the discretisation of the diffusive fluxes . Here, we choose to use the hybrid mimetic mixed (HMM) method, due to its ability to handle anisotropic diffusion tensors on generic meshes, see e.g. [10, 12]. For this paper, we only present the HMM method for Cartesian meshes, since the tests and models are performed on square cells. Denoting by the collection of all discrete values, one on each cell, and one on each face (see Figure 3, left), we then define the diffusive fluxes by the relation:

 ∀K∈Ωh,∀v∈XD, (19) ∑σ∈EKFDK,σ(vK−vσ)=∫KΛK∇Dc(x)⋅∇Dv(x)dA,

where is a discrete gradient that is defined as follows: For , we set

 ∀w∈XD,∀x∈K, ∇Dw(x) =¯¯¯¯¯∇Kw+SKw(x), (20)
 ¯¯¯¯¯∇Kw=1|K|∑σ∈EK|σ|(wσ−wK)nK,σ. (21a) We see in (20) that the discrete gradient consists of two components. The first component, ¯¯¯¯¯∇Kw in (21a), is a linearly exact reconstruction of the gradient. The second component, SKw(x), is a stabilisation term, which is piecewise defined such that for each convex hull (TK,σ)σ∈EK of σ and xK (see Figure 3, right), SKw(x)=SK,σw, with SK,σw=√2dK,σ[wσ−wK−¯¯¯¯¯∇Kw⋅(¯¯¯xσ−xK)]nK,σ. (21b)

Since we work on square cells, the expressions (21a) and (21b) can be simplified. Firstly, we have for every edge , . By adopting the compass notation and denoting by , , , the north, south, east and west edges of cell , respectively, we find that

 ¯¯¯¯¯∇Kw=1h[wσE−wσWwσN−wσS].

Moreover, we obtain the following simplified expression for the stabilisation term.

 SK,σ=⎧⎪⎨⎪⎩√2h(wσE−2wK+wσW)nK,σforσ=σEorσ=σW,√2h(wσN−2wK+wσS)nK,σforσ=σNorσ=σS.

Finally, owing to (19), the diffusive flux along the edge of cell can be uniquely determined by choosing such that and for .

###### Example 4.1 (Computing the diffusive flux along an edge).

As an example, to compute the diffusive flux along the eastern edge of a cell , we choose such that , , . Substituting into (19) and denoting by the value of the diffusion tensor at cell , we then have

 FDK,σE=∫KΛK∇Dc⋅∇DvdA=∑σ∈EK∫TK,σΛK∇Dc⋅∇DvdA,

with

 ΛK=[λK,11λK,12λK,21λK,22].

We then compute

 ∇Dv=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩−1hexforσ=σNorσ=σS,(−1h−√2h)exforσ=σE,(−1h+√2h)exforσ=σW.

Using definition (20) of the discrete gradient and denoting the values of at the cell and its edges by , respectively, we then compute the expression corresponding to . Using the appropriate values of and for each of the regions , we obtain

 FDK,σE=−1(2λK,11(cσE−cK)+λK,12(cσN−cσS)).

Here, we see that is an approximation of

 −∫σEΛ∇c⋅nK,σEds=−∫σE(λK,11cx+λK,12cy)ds

in the sense that , are approximated by and , respectively. ∎

We now discuss the discretisation of the advective fluxes . The definition of the discrete advective fluxes is motivated by the Scharfetter-Gummel or exponential scheme [25, 26]. Here, we briefly recall the Scharfetter-Gummel flux, and for simplicity of exposition, we consider an isotropic diffusion tensor and velocity field . Defining the Péclet number

 P=hv1ε, (22)

the Scharfetter-Gummel flux across the eastern edge of cell shared with cell is given by

 FK,σ=−ε|σE|h(B(P)cE−B(−P)cK), (23)

where is the Bernoulli function with

 B(z)=zez−1.

We now note that the definition of the Scharfetter-Gummel flux (23) is based on the Péclet number . Hence, for a good definition of advective fluxes to be used in (16b) and (16c), we need a proper definition of a grid-based Péclet number for anisotropic advection-diffusion problems. In particular, for each , we define as in [5, Section 3.2.2] a local grid-oriented Péclet number

 PK,σ=hVK,σ2λK,σ, (24)

with

 VK,σ =1|σ|∫σV⋅nK,σds, λK,σ =nTK,σΛKnK,σ.

Here, the quantities and measure the strength of advection and diffusion across the outward normal , respectively; hence allowing the Péclet number to capture the relative strength of advection over diffusion along properly. As a remark, we note that for isotropic diffusion, the expression (24) of the Péclet number boils down to the classical definition (22) of the Péclet number.

Since the diffusive fluxes have already been defined in Section 4.3, we consider only the advective component of the Scharfetter-Gummel flux. This is done by following the ideas in [1, 2], where instead of the Bernoulli function, we extract the advective component of the Scharfetter-Gummel fluxes by defining

 Asg(z)=B(−z)−1.

The advective fluxes are then given by

 FAK,σ (25)

## 5. Summary and implementation of numerical scheme

In this section, we present a summary of the numerical scheme for the magnetised transport problem described by the model system of equations (1a), (1b), and (1c). Starting with and initial ion and electron density profiles of and , respectively, the idea is to sequentially solve, the potential , followed by the ion density and then the electron density . This process is repeated until we reach a time such that the solution is already at steady state. Numerically, the solution is said to be at steady state whenever

 ∥V