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
(1a) | |||
where is the permittivity of vacuum, is the potential, and is the space charge density given by | |||
Here, is the elementary charge, and 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: |
(1b) | |||
(1c) |
where and are the drift-diffusion flux densities for the ions and electrons, respectively, defined by
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).
ionisation rate | |
mobility coefficient | |
diffusion coefficient | |
magnetic tensor | |
local electron density |
density of the background gas | |
---|---|
eV | mean electron energy |
permittivity of vacuum | |
Boltzmann constant | |
elementary charge | |
temperature of the background gas |
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
where, and . For ions, we assume that and so we have
(2a) | |||
We assume that electrons follow a Maxwellian distribution so that and thus | |||
(2b) |
Thirdly, the the local electron density is computed by the relation [15]
(3) |
Finally, the magnetic tensor is given by
(4a) | ||||
with Hall parameters , and magnetic field . We now note that is a cross product matrix such that for any vector , . Taking note that | ||||
we then see that is an eigenvector of corresponding to the eigenvalue 1. Essentially, this means that 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
(4b) | ||||
where is the strength of the magnetic field. Considering now a vector
orthogonal to the magnetic field, we see that
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),
(5a) | ||||
Then, flux boundary conditions are imposed for the drift-diffusion equations (1b) and (1c), i.e., | ||||
(5b) | ||||
(5c) | ||||
where is the thermal velocity, denotes the positive part of the function, and 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 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:
(5d) | ||||
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) |
(6b) |
(6c) |
Here,
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
(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
(8) |
We then note that the time derivative can be computed from the drift-diffusion equations, and can be written as
(9) |
Using now (8) for the charge density in the Poisson equation and evaluating in (9) at time , we then seek, for , such that
(10) | ||||
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
or equivalently,
Denoting the correction flux at time level by
(11) |
the above equation can be rewritten as
(12) | ||||
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 ,
respectively. By applying the divergence theorem and using , these can be rewritten as
(13a) | ||||
(13b) | ||||
(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
(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
(15a) | |||
(15b) | |||
(15c) | |||
(15d) | |||
(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 ,
(16a) | |||
(16b) | |||
(16c) | |||
where is the average ionisation rate in cell at time , and is the value of computed at cell . |
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.
(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 ,
(18a) | ||||
(18b) | ||||
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:
(19) | ||||
where is a discrete gradient that is defined as follows: For , we set
(20) |
(21a) | |||
We see in (20) that the discrete gradient consists of two components. The first component, in (21a), is a linearly exact reconstruction of the gradient. The second component, , is a stabilisation term, which is piecewise defined such that for each convex hull of and (see Figure 3, right), , with | |||
(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
Moreover, we obtain the following simplified expression for the stabilisation term.
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
with
We then compute
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
Here, we see that is an approximation of
in the sense that , are approximated by and , respectively. ∎
4.4. 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 [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
(22) |
the Scharfetter-Gummel flux across the eastern edge of cell shared with cell is given by
(23) |
where is the Bernoulli function with
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
(24) |
with
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
The advective fluxes are then given by
(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