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.


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 [17].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 [4,12].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.: [10,17,18,21,28].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 [17].
The model system of equations is highly nonlinear.Moreover, the Poisson equation and the drift-diffusion equations are strongly coupled for the computation of the electric field and the charge densities.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 [30].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 accurately.Here, we only explore linear diffusion tensors.Nonlinear diffusion tensors can be addressed by considering a generalisation of the SG scheme [3,6], but are not in the scope of this paper.Recent studies propose the use of high order finite difference gradient reconstruction methods, together with a magnetic field aligned mesh [29,34] in order to resolve the anisotropic diffusion.
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 [14] 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 on generic meshes.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.Also, the combination of the HMM with the SG scheme has been widely explored, and has several desirable properties for drift-diffusion equations, see e.g.[7,8].As a remark, the HMM method can also be viewed as a gradient reconstruction method or a gradient scheme [13], and can thus be used in combination with the methods proposed in [29,34].
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 driftdiffusion 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.

Mathematical model
Suppose that an electric field E is applied between two parallel electrodes and that an electron flux is forced through a uniform medium, where an oblique magnetic field B is applied.The electric field E = −∇V in a spatial domain Ω over a time interval (0, T ) is governed by the Poisson equation where ε is the permittivity of vacuum, V is the potential, and ρ is the space charge density given by Here, q is the elementary charge, n i and n e 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 driftdiffusion equations: where Γ i and Γ e are the drift-diffusion flux densities for the ions and electrons, respectively, defined by Γ e := µ e En e − D e ∇n e .Equations (1b) and (1c) are coupled to the Poisson equation (1a) due to the electric field E 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 p = i, e, referring to a particle (either an ion or an electron).Boltzmann constant q elementary charge T g temperature of the background gas We now discuss the parameters used in the model.Firstly, the ionisation rate k and the mobility coefficients µ p are computed via local field approximations, and are hence functions of the reduced electric field E /n g .In this paper, we consider helium gas, whose ion mobility µ i (Figure 2, left) is obtained from [15].The ionisation rate k (Figure 1) and electron mobility µ e (Figure 2, right) were computed in [19] using cross-sectional data for helium gas in [11,25] for low energies and [20] for high energies.Here and throughout the paper, • refers to the L 2norm.Secondly, we discuss the diffusion coefficients.These are governed by the Stokes-Einstein relation where, q i = q and q e = −q.For ions, we assume that T i = T g and so we have We assume that electrons follow a Maxwellian distribution so that k B T e = 2¯ /3 and thus Thirdly, the local electron density n e,loc is computed by the relation [19] n e,loc = M e Γ e µ e M e E .
Finally, the magnetic tensor M e is given by with Hall parameters β i = µ e B i , i = 1, 2, 3 and magnetic field B. We now note that C(β) is a cross product matrix such that for any vector At this stage, it is also important to note that the magnetic tensor M e is positive semidefinite.That is, v T M e v ≥ 0 for any vector v in R 3 .We also note that which shows that β is an eigenvector of M e corresponding to the eigenvalue 1.
Essentially, this means that M e v preserves the component of the vector v that is parallel to the magnetic field.
For dimension d = 2, the magnetic field can be specified at an angle θ with respect to the x-axis by setting where B 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 For two-dimensional models, we neglect the z-component, and hence we see that the components orthogonal to the magnetic field are scaled by a factor 1/( β 2 +1).These properties agree with the classical definition [22] of a magnetic tensor.
2.1.Boundary conditions.We now describe the boundary conditions for the model.At the electrodes (corresponding to either anodes or cathodes), Dirichlet boundary conditions are imposed for the Poisson equation (1a), V = V a anode, Correspondingly, flux boundary conditions are imposed for the drift-diffusion equations (1b) and (1c) at the electrodes, i.e., where v th,p is the thermal velocity, f + = max(f, 0) denotes the positive part of the function f , and n is the outward unit normal vector at the boundary.The first two terms in the right hand side of equations (6b) and (6c) 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 (6c) is the secondary electron emission, which describes the production of electrons as a result of ion fluxes bombarding the wall [5], 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, 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 n e (x, 0), n i (x, 0) at time t = 0, we are then interested in calculating the densities n e (x, T ), n i (x, T ) at time t = T .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.

Time discretisation
We form a partition of the time interval (0, T ) by taking 0 = t (0) < t (1) < . . .t (M ) = T .A semi-implicit Euler method will be considered for the time integration.This will be obtained by starting off with a fully implicit Euler method, for which some of the implicit terms will be made explicit in order to have a linear scheme.Even though the Euler method is only first-order accurate in time, it is sufficient since we are only interested in the steady-state solution.Using a uniform time step of ∆t = t (m+1) −t (m) for m = 0, . . .M −1, and given the densities n Here, p is the drift-diffusion flux density of particle p (ions or electrons) at time t (m+1) .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 ∆t (m) = t (m+1) − t (m) .We note that due to the electric field, the Poisson equation (7a) is coupled to the continuity equations (7b) and (7c).This results in a nonlinear scheme, which is computationally very expensive.To avoid this, we consider instead a semi-implicit scheme which is obtained by decoupling the system of equations (7a)-(7c).
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 m = 0, . . .M − 1, The advantage of using ( 8) 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 O(∆t).Hence, a very small time step, at most equal to the dielectric relaxation time [32, 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 t (m+1) , we may write We then note that the time derivative ∂ρ ∂t can be computed from the drift-diffusion equations, and can be written as Using now (9) for the charge density in the Poisson equation and evaluating ∂ρ ∂t in (10) at time t (m+1) , we then seek, for m = 0, . .
In this case, the quantities E, M e , µ i , µ e , D i , D e , n i , and n e at time t (m+1) 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 [32,33].That is, we keep E implicit and take first-order approximations to M e , µ i , µ e , D i , D e , n i , n e , which leads to or equivalently, Denoting the correction flux at time level m + 1 by the above equation can be rewritten as We now look at the coercivity properties of the modified Poisson equation which results from the second-order approximation of the charge density.Due to the correction terms, the matrix P = εI associated to the Poisson equation is now modified and we have e is positive semidefinite, and thus we have that v T Pv ≥ εv T v for any vector v, which ensures the coercivity of the matrix P associated to the modified Poisson equation.

Spatial discretisation
For the spatial discretisation, we start by forming Ω h , a partition of the domain Ω into N K control volumes.Each of the control volumes K has edges σ ∈ E K such that max σ∈E K ,K∈Ω h |σ| = h, and we denote by N σ the total number of edges contained in this partition.We then denote by X D the space which contains the collection of unknowns (one on each cell, and one on each edge).For this work, we focus on finite volume type schemes.

Finite volume discretisation.
To write the equations (7a),(7b), and (7c) in their finite volume forms, we take their integrals over all control volumes K ∈ Ω h .This results in, for each control volume K, e,loc dA = 0, e,loc dA = 0, respectively.By applying the divergence theorem and using E (m+1) = −∇V (m+1) , these can be rewritten as e,loc dA = 0, (14b) e,loc dA = 0. (14c) Here, n K,σ is the outward unit normal vector of K at σ. Due to the definition (3) of n e,loc , the term n is nonlinear in n e .Here, we choose a first-order approximation and take the local electron density at time m so that 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 − σ Λ∇c • n K,σ ds for some diffusion tensor Λ and scalar c; on the other hand, the advective fluxes take the form σ cV • n K,σ ds, where V is a velocity field.In equations (14a)-(14c), three diffusive fluxes and two advective fluxes are present.We then denote the discrete diffusive and advective fluxes for each edge σ in a cell K ∈ Ω h by Replacing the fluxes in (14a)-(14c) with their discrete counterparts (16a)-(16e) at the appropriate time level and under the assumption that n e,K respectively in a cell K, we then obtain our numerical scheme.For the Poisson equation, we start by considering a first-order approximation of the charge density, given in (8).That is, for m = 0, . . ., M , given the densities n |K| n Starting with m ≥ 1, we then use the second-order approximation for the charge density in solving the Poisson equation.That is, instead of using (17a), we use a semi-implicit discretisation based on (13), given by where is a discretisation of the semi-implicit correction flux Γ c defined in (12).
We note here that we cannot immediately use (17d) at the first time step due to the yet unknown mobility coefficients, correction and diffusive fluxes at time t = 0.One natural option would be to use (17a) for solving the Poisson equation at the first time step when m = 0.Alternatively, we may opt to use a predictorcorrector method, e.g., starting with (17a), we get an initial estimate of the electric field, which will then be used for an initial estimate of the diffusion and mobility coefficients.These coefficients can then be used for computing the diffusive and correction fluxes, which in turn allows us to use (17d) for solving the Poisson equation at time t = 0.However, we notice in our numerical tests that this does not have a significant impact on the quality of the numerical solutions.Hence, for the simulations, we start with (17a).4.2.Boundary conditions and conservation of fluxes.We now see that there are N K equations (corresponding to the number of cells) involved in each of (17a)-(17d).However, we have N K + N σ unknowns for each of these equations: N K corresponding to the number of cells, and N σ corresponding to the number of edges.This is standard for hybrid finite volume methods, and the equations for the N σ 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 K and L, we should have that the total (diffusive and advective) flux F K,σ going from K to L via σ is equal to the negative total flux F L,σ going from L to K via σ, i.e.
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 (6b) and (6c).This is done by setting, for the corresponding boundary edges σ ∈ E K , ).
(19b) Here, n e,σ and n i,σ are the values (to be determined) of n e and n i 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.[14,16].Denoting by v ∈ X D a collection of discrete values, one on each cell, and one on each face (see Figure 3, left), we then define the diffusive fluxes by the relation: where ∇ D is a discrete gradient that is defined as follows: For K ∈ Ω h , we set ∀w ∈ X D , ∀x ∈ K , We see in (21a) that the discrete gradient consists of two components.The first component, ∇ K w in (21b), is a linearly exact reconstruction of the gradient.The second component, S K w(x), is a stabilisation term, which is piecewise defined such that for each convex hull (T K,σ ) σ∈E K of σ and x K (see Figure 3, right), S K w(x) = S K,σ w, with We now illustrate the HMM method for Cartesian meshes.We note that even though the HMM is only illustrated for Cartesian meshes, it is also applicable to generic types of meshes, which will be demonstrated in the tests in Section 6.2.Since we work on square cells, the expressions in (21b) and (21c) can be simplified.Firstly, we have for every edge σ ∈ E K , |σ| = h.By adopting the compass notation and denoting by σ N , σ S , σ E , σ W the north, south, east and west edges of cell K, respectively, we find that Moreover, we obtain the following simplified expression for the stabilisation term.
Finally, owing to (20), the diffusive flux F D K,σ along the edge σ of cell K can be uniquely determined by choosing v ∈ X D such that v K = 0, v σ = −1 and v σ = 0 for σ = σ.(20) and denoting by Λ K the value of the diffusion tensor at cell K, we then have We then compute Using definition (21a) of the discrete gradient and denoting the values of c at the cell K and its edges by c K , c σ N , c σ S , c σ E , c σ W , respectively, we then compute the expression corresponding to ∇ D c.Using the appropriate values of ∇ D c and ∇ D v for each of the regions T K,σ , we obtain Here, we see that in the sense that c x , c y are approximated by 2 h (c σ E − c K ) and 1 h (c σ N − c σ S ), respectively.

Advective fluxes. 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 [30,31].Here, we briefly recall the Scharfetter-Gummel flux, and for simplicity of exposition, we consider an isotropic diffusion tensor Λ = εI and velocity field V = [v 1 , 0] T .Defining the Péclet number the Scharfetter-Gummel flux across the eastern edge σ E of cell K shared with cell E is given by where B is the Bernoulli function with .
We now note that the definition of the Scharfetter-Gummel flux ( 23) is based on the Péclet number P .Hence, for a good definition of advective fluxes F A K,σ to be used in (17b) and (17c), we need a proper definition of a grid-based Péclet number for anisotropic advection-diffusion problems.In particular, for each K ∈ Ω h , σ ∈ E K , we define as in [9, Section 3.2.2] a local grid-oriented Péclet number with Here, the quantities V K,σ and λ K,σ measure the strength of advection and diffusion across the outward normal n K,σ , respectively; hence allowing the Péclet number P K,σ to capture the relative strength of advection over diffusion along n K,σ 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,7,8], where instead of the Bernoulli function, we extract the advective component of the Scharfetter-Gummel fluxes by defining The advective fluxes are then given by As a remark, the accuracy and convergence of numerical solutions from the combination of the HMM with the Scharfetter-Gummel scheme for advection-diffusion equations have already been studied in [2].Remark 4.2 (Charge conservation).Looking at (17d), we see that due to (10), the integral of the term ∂ρ ∂t over a cell can actually be expressed in terms of the integral of the charge density fluxes over the cell's faces.This suggests that our usage of a finite volume discretisation in space ensures that a positive contribution of a given cell is balanced by a negative contribution from its neighboring cell, i.e. the charge is always conserved for internal faces.We note however that before reaching steady state, charge is not necessarily conserved at the boundary, due to the fact that Dirichlet boundary conditions (6a) are imposed for the potential.This means that a net charge is present when integrating the charge density over the control volume at the boundary.

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 m = 0 and initial ion and electron density profiles of n (0) i and n (0) e , respectively, the idea is to sequentially solve, the potential V (m+1) , followed by the ion density n .This process is repeated until we reach a time T = t (m+1) such that the solution is already at steady state.Numerically, the solution is said to be at steady state whenever (26) We now elaborate the details of the computations involved in the first few steps of the scheme.We start by computing the potential V (1) via solving (17a) with diffusive fluxes (16a) defined as in (20).Under the assumption that the initial ion and electron densities are equal, we then solve at time m = 0, for each K ∈ Ω h , the balance of fluxes − 2ε(V (1)  σ W − 2V (1) This is solved together with the conservation of fluxes for interior edges σ between cells K and L, given by Finally, to complete the set of equations for the potential V (1) , we impose, for the appropriate edges, Dirichlet boundary conditions, and Neumann boundary conditions, given by K . (27d) The linear system (27a)-(27d) then fully defines V (1) .As a remark, we note that due to the absence of anisotropy in the Poisson equation, the HMM method on Cartesian meshes boils down to a standard central difference scheme.We also note that using only the consistent part (21b) of the discrete gradient will yield an electric field E (1) that does not satisfy the homogeneous Neumann boundary condition (6d).Hence, we reconstruct the HMM discrete gradient as in (21a) in order to obtain the value of the electric field E (1) on each cell K ∈ Ω h .At first glance, it may seem that the hybrid scheme is very expensive to implement, since it involves solving a system of N K + N σ equations in N K + N σ unknowns.However, by looking at the first N K equations (27a) corresponding to the cells K ∈ Ω h , we see that they are fully local in the sense that the value of V on a cell K only depends upon the values of V on its edges.Hence, static condensation can be applied for an efficient implementation of the scheme, i.e., in terms of the matrix system Ax = b that needs to be solved, we actually realise that upon partitioning and writing where A 11 and A 12 contains the first N K equations (27a) for the unknowns V K , we see that A 11 is actually diagonal, and hence, the main cost of the scheme only involves solving a system of N σ equations for N σ unknowns V σ .
Using the value of E (m+1) , the ionisation rates k and the mobility coefficients µ e and µ i are computed by the lookup tables as in Figures 1 and 2. The diffusion coefficients D i and D e are then calculated according to the definitions (2a) and (2b), respectively.Following this, we compute the local electron density on each cell K ∈ Ω h as in (15).These quantities allow us to compute the diffusive and advective fluxes (16b)-(16e) by using the formulations ( 20) and (25).We then substitute the appropriate expressions corresponding to (16b) and (16d) into (17b) to find n (m+1) i . The system of equations (17b) used to compute n i are given by, for each e,loc,K n g , where, denoting by E x and E y the x-and y-components of the electric field, respectively, the Péclet numbers are given by P Here, we also used the fact that P K,σ W = −P K,σ E , P K,σ S = −P K,σ N .Following this, the conservation of flux conditions are obtained by substituting, for interior edges, the appropriate expressions for (16b) and (16d) in (18).
For the boundary edges corresponding to electrodes, the flux boundary conditions (19a) are applied.Finally, the homogeneous Neumann boundary conditions require that We then use n (m+1) i to compute the secondary electron emission Γ (m+1) i •n for (19b).Lastly, a similar process may be used to assemble the system of equations needed to compute n (m+1) e , i.e., by assembling the N K equations (17c), and imposing the appropriate flux conservation conditions (18) and boundary conditions.The main difference between the equations for n e and n i comes from the presence of the magnetic tensor M e in (17c).Due to anisotropy, the expression for the diffusive flux of n e is different from the one obtained for n i .Moreover, this requires us to use the Péclet number as defined in (24) for the advective fluxes of n e .
For m ≥ 1, we use (17d) instead of (17a) to solve the Poisson equation for the potential V and electric field E, which will be used for computing the ion and electron densities in the next time steps.

Numerical tests
In this section, we perform numerical tests using the scheme (17) described above for the model system of equations (1a), (1b), and (1c).Here, the domain involves 2 parallel plates that are 0.02m apart, which is given by Ω = (0, 0.02) × (0, 0.02).We start with a simple test case, for which there is an anode on the left wall, and a cathode on the right wall.In particular, we set For all test cases, we use the following parameters in Table 3.The thermal velocities of the ions and electrons presented below were calculated by the formula v th,p = 8k B T p πm p with m i = 6.6462 × 10 −27 kg and m e = 9.1093 × 10 −31 kg.
Table 3. Parameters for the test.
temperature of the background gas T g 300K density of the background gas n g 1.9 × 10 23 m −3 thermal velocity, electron v th,e 7.7277 × 10 5 m/s thermal velocity, ion v th,i 1.2598 × 10 3 m/s emission coefficient γ i 0.2 initial density, electron n e (x, 0) 5 × 10 13 m −3 initial density, ion n i (x, 0) 5 × 10 13 m −3 6.1.1D test.We start with a 1D test case, obtained by setting M e = I, and consider Ω = (0, 0.02), partitioned uniformly.This will be referred to as test case 1.We compare in Figure 4 the steady-state density profiles between our proposed scheme and the classical Scharfetter-Gummel scheme [30] on a fine mesh with 512 cells, whilst taking a small time step of ∆t = 2e-9s.Here, we set the tolerance value to indicate that the densities have reached steady state in (26) to be tol = 0.01.Typically, this is achieved at time t = 1e-5s.In Figure 4, we see that the numerical solutions from both methods give a similar shape for the density profiles.To be specific, the ion densities peak on the interval (0, 0.002), and a quasi-neutral state is reached for x > 0.002.However, upon having a closer look at Figure 4, right, we see that the numerical solution obtained from the Scharfetter-Gummel scheme exhibits a steeper dip in the ion density before reaching quasi-neutrality.Moreover, the ion and electron densities in this region are slightly larger than that of scheme (17).In particular, we have that the maximum ion and electron densities in this region are 3.6 × 10 15 m −3 for the Scharfetter-Gummel scheme, which is approximately 10% larger than 3.2 × 10 15 m −3 , obtained from scheme (17).We now compare the solution profiles on a coarser mesh, now consisting of 64 cells in Figure 5.In Figure 5, we observe that the numerical solutions on the coarse mesh still exhibit the same behavior as the solutions on the fine mesh.Namely, the ion densities still peak in the interval (0, 0.002) and a quasi-neutral state is still observed for x > 0.002.However, we see for both schemes an increase in the ion and electron densities compared to the solutions obtained on the fine mesh.This is not unexpected, since the solutions on a coarse mesh are in most cases less accurate than the ones on a fine mesh.We note however that the solution peaks only rise by around 12.5%, which is an acceptable trade-off, since the computational costs are reduced by a factor of 8.
Due to the use of a second-order approximation of the charge density in the Poisson equation by solving (17d), we may take a larger time step of ∆t = 4e-8s for both methods, and still obtain results that are essentially the same as those presented in Figures 4 and 5.
6.2.2D tests.We now proceed to a 2D test.For these tests, we will present our results on Cartesian meshes (Figure 6, left).In order to illustrate that our scheme can also be applied to more generic types of meshes, we also ran tests on perturbed Cartesian meshes (Figure 6, right), constructed following the guidelines provided in [23].That is, starting with a uniform Cartesian mesh (Figure 6, left), if the maximum diameter of the cells are given by h, then the internal nodes (x, y) are perturbed randomly by taking x := x + 0.4β x h, ŷ := y + 0.4β y h, where β x , β y are random values from a uniform distribution on (−0.5, 0.5).We then consider a two-dimensional magnetic field.That is, we take B as in (4b), with B(x) = 0.2 min 1, 2x 0.02 min 1, 2y 0.02 .
Physically, this means that the magnetic field grows stronger as we get closer to the upper and right parts of the domain, with a maximum strength of 0.2.We then perform numerical tests with magnetic field angles of θ = 45 • , and θ = 90 • in (4b), respectively on a uniform and perturbed grid consisting of 64 × 64 cells, with a time step of ∆t = 4e-8s.We will refer to this as test case 2. We will only present one figure for each of the tests, as the results obtained on the perturbed mesh are essentially the same as those obtained on the uniform mesh.In both Figures 7 and 9, we see that the electron and ion densities are sparse on the upper right region of the domain.This is due to the fact that the electrons did not have sufficient time to go to the opposite electrode.Physically, this is expected, because the magnetic field creates a barrier with an angle of 45 • and 90 • , respectively, in the upper right half of the domain.This barrier separates the electrons that have a smaller loss term at the upper boundary than at the electrode for a 45 • and 90 • angle.We now discuss the behaviour of the density profiles for the 45 • angle in more detail.In particular, we discuss separately the effects of the magnetic field at the upper right, central, and left regions, respectively (see Figure 7, right).Firstly, the loss at the upper right corner of the domain is equal to the wall loss.Going towards the upper center, this wall loss at the upper boundary decreases because the electrons are not fully magnetised.Additionally, the influx from other parts of the domain is high at this location.Finally, at the upper left,  the electron flux density from other parts of the domain is lower due to the smaller magnetisation and therefore the density decreases.Figures 8 and 10 show the effect of the magnetic field angle on the fluxes Γ e .In particular, we see in Figure 8 that after magnetisation, the fluxes are now oriented in an angle of 45 • at the upper right region of the domain.On the other hand, Figure 10 shows that by applying a 90 • magnetic field, the fluxes which were initially traveling from left to right horizontally has been decreased significantly.As discussed in (5), this is expected since the flux is orthogonal to the magnetic field.The discussion above and the results in Figures 7-10 show that the scheme (17) can capture the shape of the solution profile and the physics properly on both uniform and distorted meshes.Of course, as seen in the 1D test, if we want a more accurate depiction of the peaks of the solution profile, we need to work on a finer mesh.We also present in Figure 11 the potential V in both tests.
In order to further test the limits of the scheme, we consider a test case for which the boundary is not fully a cathode.To do so, we slightly alter the boundary conditions of the Poisson equation.In particular, instead of purely having a cathode  Note that since the upper part of the right wall is no longer a cathode, the corresponding boundary conditions for the ion and electron continuity equations at this region should also be changed.That is, we set, for p = i, e, ∇n p • n = 0 on x = 0.02, y ∈ 1 4 × 0.02, 0.02 .
We will refer to this as test case 3.As with test case 2, we perform our numerical simulations on a grid with 64 × 64 cells, with a time step of ∆t = 4e-8s.Note that upon changing the boundary conditions from test case 2, Figures 12 and 13 exhibit two-dimensional potential and density profiles, even with no magnetisation (setting M e = I).
Finally, we illustrate in Figures 14 and 15 that our numerical scheme is able to handle test case 3, even with the most extreme magnetic field with θ = 45 • .This is the most extreme test in the sense that a 45 • angle maximizes the term  β 1 β 2 in the magnetic tensor M e defined in (4a), which translates to an anisotropic diffusion tensor with very strong cross-diffusion.We now give a detailed discussion about the results in Figures 14 and 15.We observe in Figure 12, right and Figure 15, left that without magnetisation, most of the electrons travel from the left wall to the right wall, parallel to the x-axis at the lower part of the domain, and with approximately an angle of −45 • with respect to the x-axis at the upper part of the domain.Upon enabling a magnetic field of θ = 45 • , we see that the flux in −45 • direction is greatly decreased (see Figure 15, right), which is expected since the flux is perpendicular to the magnetic field (see Eq.( 5)).We also observe in Figure 14 that the maximum electron density is located near the center of the wall.This is due to the fact that above this region, there is a loss of electrons towards the upper boundary, whereas below this region, the electrons are attracted to the bottom right.

Conclusion
In this work, we propose the use of a combined hybrid mimetic mixed (HMM) method with the Scharfetter-Gummel (SG) scheme for magnetised transport problems encountered in plasma physics.The main reason behind the employment of the HMM method is due to its ability to handle anisotropic diffusion tensors on generic meshes, which is mainly encountered in the presence of magnetic fields.Another advantage of hybrid schemes is that the stencils involved for solving the linear system of equations are purely local, which allows the use of static condensation for a very efficient implementation.Another important factor we introduced is a grid-based Péclet number which captures the relative strength of advection over diffusion properly, which was used for the definition of the advective fluxes.
Numerical tests show that on a standard one-dimensional model, the proposed scheme is slightly better than the Scharfetter-Gummel scheme, both on fine and coarse meshes.It is also important to note that even though the peaks of the solution profiles are not captured properly on the coarse mesh, both schemes still capture the shape (location of density peaks and quasi-neutral regions) of the solution correctly.For two-dimensional models, the numerical results obtained from our proposed scheme agree with what is physically expected, even for distorted meshes and for extreme tests which involve highly anisotropic magnetic fields with strong cross-diffusion terms.On the other hand, it is not straightforward to extend the classical Scharfetter-Gummel scheme in order for it to handle the anisotropic diffusion that arises from the presence of the magnetic field.
Numerical tests showed the robustness and applicability of the HMM-SG method on both uniform and distorted quadrilateral meshes.Hence, an interesting and natural follow-up would be the application of the proposed method onto the magnetic field aligned meshes in [29,34].Another interesting aspect to consider would be the application of the proposed scheme to models in cylindrical coordinates, which is commonly encountered in applications.Another avenue that may be explored for future research would involve using the complete flux scheme [24,31] or hybrid high order schemes [26,27] for the spatial discretisation, which will allow us to resolve the peaks of the solution profile, even on coarse meshes.We also note that due to the use of a first-order Euler time integration scheme for the continuity equations, we might not have been able to fully take advantage of the second-order approximation of the charge densities.Hence it might also be interesting to investigate the use of second-order time integration schemes for the continuity equations, especially if it allows us to preserve the accuracy of the solution whilst taking much larger time steps.These are quite challenging aspects to explore since a proper study and choice of flux limiters would be needed to ensure that the numerical solutions do not exhibit nonphysical oscillations.

Figure 1 .Figure 2 .
Figure 1.Ionisation rate 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.

Figure 3 .
Figure 3. Notation for a Cartesian cell in dimension d = 2.

Example 4 . 1 (
Computing the diffusive flux along an edge).As an example, to compute the diffusive flux along the eastern edge σ

Figure 13 .
Figure 13.Steady state density profiles, test case 3, no magnetic field (Left: n i , Right: n e ).

Table 1 .
Parameters in the model.

Table 2 .
Constants in the model.