Genetic Exponentially Fitted Method for Solving Multi-dimensional Drift-diffusion Equations

A general approach was proposed in this article to develop high-order exponentially fitted basis functions for finite element approximations of multi-dimensional drift-diffusion equations for modeling biomolecular electrodiffusion processes. Such methods are highly desirable for achieving numerical stability and efficiency. We found that by utilizing the one-one correspondence between continuous piecewise polynomial space of degree $k+1$ and the divergence-free vector space of degree $k$, one can construct high-order 2-D exponentially fitted basis functions that are strictly interpolative at a selected node set but are discontinuous on edges in general, spanning nonconforming finite element spaces. First order convergence was proved for the methods constructed from divergence-free Raviart-Thomas space $RT_0^0$ at two different node sets

1. Introduction. We will propose in this article a general approach to construct exponentially fitted methods for numerically solving the drift-diffusion equation where u is the density of charged particles, D is the diffusion coefficient, β is a problem-dependent constant assumed positive in this study, φ is the given electrostatic potential field. The source function f models the generation or recombination of the particles, and is assumed here to be independent of u. The drift-diffusion equation is widely applied in semiconductor device modeling, where u(x) is the density of charge carriers [18], and in biomolecular simulations, where u(x) can be the density of ions, ligands, lipids, or other diffusive molecules [21,22,36]. Generalizations of Eq.(1.1) to model quantum effects, finite particle sizes, particle-particle correlations and interactions have been made recently [15,23,36], in many cases by replacing the electrostatic potential with an effective potential that models these non-electric interactions. The drift-diffusion equations arising in biomolecular simulations are in general multidimensional due to the intrinsic 3-D structure of macromolecules. For drift diffusion in bulk solution or through ion channels 3-D drift-diffusion equations are usually adopted so the charge density can be solved at sufficient temporal and spatial accuracy. For lateral transport of charged particles on surfaces, such as lateral motion of lipids, lipid clusters, or proteins on membrane surfaces, 2-D drift-diffusion will be very useful for the reduction of dimensionality [19]. Surface drift-diffusion equation rather than 2-D drift-diffusion equation has to be used for lateral diffusion on smoothly curved surfaces [3,2,36], for which the gradient and divergence operators in (1.1) have to be replaced by surface gradient and surface divergence operators, respectively. This surface drift-diffusion equation can be solved using surface finite element methods, which share many properties with 2-D finite element methods. It is known that numerical solutions of drift-diffusion equations may suffer instability in the advection-dominated regimes that usually appear near the strongly charged molecular surfaces such as ion channels, deoxyribonucleic acid (DNA), or membrane surfaces with proteins bound or in the vicinity. Numerical solutions of drift-diffusion equation in these biomolecular systems necessitates stable and efficient numerical methods. Exponentially-fitted methods, known originally as the Scharfetter-Gummel method [32], is a class of finite element methods for solving Eq.(1.1) through an exponential approximation to the density function u(x) on individual mesh elements. This is made possible by the introduction of the Slotboom variable ρ = ue βφ , (1.2) with which Eq.(1.1) can be symmetrized to be The density function u(x) can then be exponentially fitted through locally constant or linear approximations of the density current J = De −βφ ∇ρ. With exponentially fitted methods the gradient of the electric potential field can be incorporated to the basis functions, and this makes the methods intrinsically stable against strong drift when the electric field is strong. This desirable feature has motivated many studies of the exponentially fitted methods [1,32,10,16,30,11,26,5]. Exponential fitting for 1-D drift-diffusion equation is obtained by exactly solving a two-point boundary value problem locally in a divergence-free space of Lagrange basis functions [32]. The resultant nodal basis functions for u(x) involve Bernoulli functions. In multi-dimensional space, however, it is generally difficult to find a divergence-free space where the basis functions are interpolative at nodes and continuous on edges, for the local problem does not admit an analytical solution on triangles, quadrilaterals, or tetrahedra. To circumvent this difficulty a number of averaging schemes, such as exact average, volume harmonic average, and edge harmonic average, were introduced to the exponentially transformed diffusion coefficient in the symmetrized equation (1.3). 2-D exponentially fitted methods of this type were constructed in [10]. Such a construction of multi-dimensional exponential fitting methods is further analyzed in [16], where it is found that the classical Scharfetter-Gummel exponentially fitted scheme can be exactly reproduced by the harmonic averaging of the piecewise linear Galerkin method. Linear convergence in the exponentially weighted norm is proved as well. While the convergence and stability of these multi-dimensional methods in numerical computations were not reported in all of these studies, it was shown that they can well resolve the sharp gradient of the charge density without incurring any unphysical oscillations. Averaging techniques, such as exact average and volume harmonic average, may suffer from overflow or lead to over-smeared solutions. Some of the best results are obtained using edge harmonic average, which utilizes approximated exponential fitting along edges of triangles.
There are continual efforts to obtain genetic exponential fitting in multi-dimensional spaces rather than through averaging. These include approximate spline fitting [24,25,30,5,4], nonconforming finite element [29], finite volume methods [20], and exponentially fitted difference methods [17]. The approximate spline methods strongly resemble the 1-D exponential fitting as the locally divergence free current field solutions are sought. Such solutions, however, are not unique, and thus approximations will be in place to achieve uniqueness. Nevertheless, locally constant approximation to the density current in 2-D violates the nodal interpolation constraints [28]. This major drawback is overcome by using locally linear approximation to the density current, and nodal interpolative basis functions for density functions are obtained [31].
These basis functions are discontinuous on edges in general so the finite element space is nonconforming. Extension of this approach to 3-D equations seems feasible. The work in [5,4] are among the few attempts to establishing 3-D exponential fitting methods, where 1-D exponential fitting is sought along individual edges connecting to the vertices, and the nodal basis functions are defined to be the linear combinations of local 1-D exponentially fitted basis functions. The resulting finite element space is conforming. Despite these efforts there are very few applications of exponentially fitted methods to realistic multi-dimensional drift-diffusion equations. The analysis and application in [26], for example, is on 1-D drift-diffusion equation with quantumcorrected potential. Sometimes Streamline-Upwind/Petrov-Galerkin (SUPG) methods that were developed for advection-diffusion equations have to be adopted [12] for stabilization. Moreover, it seems unclear whether the current methodologies can be extended for constructing high order exponentially fitting methods.
We are motivated to propose in this article a general approach for constructing exponential fitting methods in multi-dimensional spaces. We start with the divergencefree vector space that is obtained by taking the curl of the nodal basis of k th order Raviart-Thomas conforming space P k+1 , to define a prototype of the basis functions for the Slotboom variable ρ. The final form of the exponentially fitted basis functions for ρ(x) and u(x) will be computed by enforcing the nodal interpolative constraints. Our approach is different from that in [31] as the latter seeks a divergence-free approximation of the density current and restrictive assumptions are introduced to get analytical representations of the basis functions. We do not enforce the density current to be divergence-free, and this relaxation enables our approach to reproduce the standard Lagrange spaces for solving the Poisson equation at the limit of vanishing electric potential in the drift-diffusion equation. With our approach 2-D exponentially fitted methods of arbitrary high order can be readily constructed by starting with high-order Raviart-Thomas conforming spaces. We will show that our method does not entail some of the crucial pitfalls indicated in previous studies. First-order 3-D exponentially fitted method can be established similarly, but the extension of our approach to high-order methods for 3-D problems is difficult due to the lack of one-one correspondence between P k+1 and the divergence-free vector space in 3-D.
The rest of the paper is organized as follows. In the next section we review some of the most important exponential fitting methods and summarize their major features. In Section 2 we define the divergence-free spaces and their basis functions, and describe the procedure to compute exponentially fitted charge density u(x) using these basis functions. Here we shall show that our construction is general, can reproduce 2-D exponentially fitted method previously developed, and degenerates to the standard conforming Lagrange finite element spaces. In Section 3 first order convergence is proved for constructions based on RT 0 0 . We make conclusion remarks and outline the future work in Section 4.

Genetic
Multi-dimensional Exponential Fitting. We first introduce notations that will be used in the discussion below. Let Ω ∈ R 2 be a bounded domain with Lipschitz continuous boundary. We introduce for s > 0 the standard Sobolev space H s (Ω) to the functions in Ω, and use the standard inner product (·, ·) s , the norm · s , and the semi-norm | · | s . The inner product and the norm in L 2 (Ω) are denoted by · and (·, ·). The subspace of L 2 (Ω) that consists of functions with zero mean value is denoted by L 2 0 (Ω). With these we define the space H(div; Ω): where n is the outer normal direction on the boundary ∂Ω, and V k is the space of vector-valued polynomials of degree k on the element K. To fulfill the condition of the vanishing normal component on the boundary ∂Ω one usually chooses V k to be the classical Raviart-Thomas element of order k (RT k ) [27] or the Brezzi-Douglas-Marini element of order k (BDM k ) [8]. For other H(div; Ω) spaces constructed more recently the readers may refer to [35,33]. We will work on the divergence-free subspace D h of V h , which is defined to be [35,9] that if V k is chosen to be RT k then In case that V h = RT k , this subspace is also denoted RT 0 K . According to the Helmholtz decomposition, any divergence-free vector v can be given as the curl of a potential function ψ: where the subscripts denote the partial derivatives. Eq.(2.1) suggests that we can construct a divergence-free subspace by taking the curl of some appropriate space. This is actually possible, as stated by the following well-known result concerning the RT k and BDM k triangular elements [9,34]: There exists a one-to-one map curl: Similar correspondence between H(curl) spaces and the divergence free subspaces D h in 3-D is also indicated [6]. Two sample spaces of D h = RT 0 k on the reference triangle are given here: 3) We are now in a position to construct the exponentially fitted methods for solving the following boundary value problem for Eq.
in Ω, where Γ D and Γ N are the Dirichlet and Neumann subsets of the boundary ∂Ω, respectively. We denote by u h and J h (u h ) = D(∇u h + βu h ∇φ) the finite element approximations of u and J(u). Let {v i } be the basis of D h on a triangular element K, then the approximate density current function J h (u h ) in K can be given as a linear combination of {v i }: It remains to find the basis functions ρ j of the certain finite element space on K in which ρ h can be approximated.

2.1.
Attempts to approximate density current exactly in divergencefree spaces. If each basis function ρ j has the representation for some set of real numbers {m j,i }, then (i) ∇ρ j and ∇ρ h will be exactly divergencefree, and (ii) the right-hand side of Eq.(2.8), as the gradient of a scalar, is curl free. Consequently, its integration along any curve l(P 1 , P 2 ) connecting two given points in K is path-independent, for l(P1,P2) For this reason one can choose the path to be the line connecting (x, y) and the starting point (x 0 , y 0 ), c.f. Figure 2.1. Then the basis function ρ j (x, y) can be computed from with u j (x 0 , y 0 ) = ρ j (x 0 , y 0 )e −βφ(x0,y0) . These functions {u j (x, y)} form a basis function of u h in the element K. On the standard reference triangle one can choose a two-section path (x 0 , y 0 ) → (x, y 0 ) → (x, y), each section parallel to the axes, giving rise to are the two components of the vector v i , and correspondingly To uniquely determine the total N k coefficients m ji and the constant u j (x 0 , y 0 ) in the expression (2.11) we will enforce the interpolative constraints at N k + 1 nodes. For example, for D h = RT 0 0 on the triangle in Figure (2.1) we can expand equation (2.13) with j = 2 at three vertices (1, 2, 3) to get the following (2.14) If we choose (x 1 , y 1 ) to be (x 0 , y 0 ) then we must have from the interpolative condition that These three equations form a linear system where F 1i and F 2i are defind as Therefore, with the similar expansions of the other two basis functions u 1 (x, y), u 3 (x, y), one can arrive at the general linear system below: where I is the identity matrix, and The following theorem proves that the matrix F is invertible for either set of the nodes, i.e., the vertices (1, 2, 3) or the edge centers (4, 5, 6).
Theorem 2.2. The matrix F is nonsingular for the basis {v i } of RT 0 0 . Proof. We notice that every component F ji represents a linear transformation of Consequently it is sufficient to prove that F j (v) = 0 for j = 1, 2 if and only if v = 0. We only need to prove that if F j (v) = 0 for j = 1, 2 then v = 0, as the reverse statement is trivial. We denote v by v = (a, b) because any v ∈ D h | K is a constant vector. If the vertices are chosen to be the node set for enforcing the interpolative constraints, we have m ji v i is enforced along with the interpolative constraints we will have an over-determined problem for m ji , suggesting that one can not find a set of interpolative basis functions for u h that constitute a constant density current. Since this curl-free condition is not satisfied, different basis functions ρ j , u j can be generated if different paths are chosen for the integrations. Figure  (2.2) shows the difference between the two basis functions that are obtained along different integral paths. The differences on the edges suggest the basis functions on neighboring triangles are discontinuous in general. The finite element space is hence nonconforming. In addition, this method can not be generalized to construct high order exponentially fitted method. For example, direct computation on D h = RT 0 1 in the reference triangle with φ = 0, along the integration path (0, 0) → (x, 0) → (x, y), gives rise to a singular matrix F : is constructed in [31], by assuming local linear approximations of the electric potential and the density current J h . The six parameters of the latter and one constant in u j are determined by applying seven conditions, namely one divergence-free condition, three linear conditions arising in the curl-free nature of ∇ρ h , and three interpolation constraints. The linear approximation of the density current provides more degrees of freedom than the constant approximation in RT 0 0 presented above. These additional degrees of freedom make it possible to enforce the curl-free condition along with the interpolative constraints. The component (∇φ×x)·e z of space allows the simulation of the diffusion flow orthogonal to the electric field lines. Interestingly, we can prove that the nonconforming space DF(K) is a special case of our construction. Consider the standard reference triangle and let D h | K = RT 0 0 then D h | K = span{(1, 0) T , (0, 1) T }. Let Path 1 be defined by (1, 0) → (x, 0) → (x, y) and Path 2 be defined by (1, 0) → (0, y) → (x, y). Note that by averaging the results for the two symmetric paths we will arrive at the final solution u j (x, y) and if we assume the same linear potential in the element φ| K = ax + by + c, then for Path 1 we have Assuming that the integration in Eq. (2.22) is path-independent, then the results from Path 1 and Path 2 are the same, and equal to their average e −β(ax+by) (2.27) Applying the approximation e x ≈ 1 + x we arrive at

28)
for some reorganized constants α j , η j , γ j . This last formula (2.28) is of the exact same form as in [31]. However, the restrictive assumption a, b = 0 in the local linear representation of the electric potential is no longer necessary in our construction here because we do not compute the analytical integration of the function e βφ and therefore the analytical form of the potential φ is not needed. When u j of the representation in Eq.(2.27) is used to compute the current J h , the first term will produce a current component in the direction of the electric field; the second term does not contribute to J h , and the last two terms will provide additional fluxes that are not necessarily aligned with the electric field. Eq.(2.28) indicates that these additional fluxes are either not necessarily orthogonal to the electric field in general. Under very special circumstances, for instances a = b in Eq.(2.27), one can solve the same m j,1 , m j,2 , and this crosswind diffusion flux will disappear in Eq.(2.27) as it will be absorbed into the second term. It is worth noting the divergence-free condition was only approximately enforced to get a constant current field in [28].

General construction of exponentially fitted methods by giving up
divergence-free condition. We would note that it is not fully justified to enforce the divergence-free condition in constructing the exponentially fitted finite element space. Divergence-free condition is seldom enforced in solving the Laplace equation, which is the limit of the drift-diffusion equation (1.1) at the limit of a vanishing electric field and with f (x) = 0. Divergence-free condition shall not be used for solving drift-diffusion equations with inhomogeneous source function f (x), or unsteady driftdiffusion equations.
We will give up the divergence-free condition to seek an alternative construction. Our construction for the drift-diffusion equation must regenerate the standard finite element space P k+1 for solving the Poisson equation at the limit of vanishing electric field. These standard basis functions ψ j of P k+1 can be derived from the known basis of D h through the following correspondence 2.29) and the enforcement of the interpolative constraints. We are motivated to construct the basis ρ j for ρ h in a similar manner On the reference triangleK, we can evaluate ρ j at an arbitrary point (x, y) via different paths (2.32) The expansion coefficients m ji and the constants ρ j (x 0 , y 0 ) can be solved from the following linear system similar to (2.21) resulting from the enforcement of the interpolative constraints at selected node sets, with 34) We can then evaluate u j (x, y) = ρ j (x, y)e −βφ(x,y) . (2.35) Proper scaling can make u j (x, y) to be interpolative as well. It turns out in general that N k + 1 = dimP k+1 (K), meaning the dimension of the finite element space for u h | K corresponding to the divergence-free subspace D h | K = RT 0 k | K is equal to the dimension of the space S h | K = P k+1 (K) whose curl produces D h | K . The one-one correspondence between RT 0 k and P k+1 indicates that the matrix F corresponding to Eq.(2.29), i.e., F with must be invertible. We shall note that the finite element spaces constructed here are nonconforming and depend on the integral path, similar to the spaces obtained in subsection 2.1, because the solved parameters m ji may not fit Eq.(2.30). They do fit when φ = 0. Two examples of the node set for the finite element space u h | K = span{u j } constructed on D h | = RT 0 0 on the reference triangle are shown in Figure 2.3. For a weak potential φ = e −2 √ x 2 +y 2 the exponentially fitted basis functions are very similar to the basis of P 1 (K) at the same set of nodes, c.f. Figure 2.4 and 2.6. The nature of exponential fitting is more clearly illustrated when the electric field is strong, c.f. Figure 2.5 and 2.7. It can be seen that the finite element spaces spanned by ρ j , u j are invariant with respect to affine linear maps, regardless of the choice of node set. This feature is of particular importance to the estimation of interpolation errors through mappings to the reference triangles. The basis functions u j constructed on D h = RT 0 1 with φ = exp −4 √ x 2 +y 2 are shown in Figure 2.8, where the characters of interpolation at nodes and exponential fitting are also well depicted. The method developed here thus provides for the first time a feasible approach to systematically construct highorder exponentially fitted methods for solving 2-D drift-diffusion equation. Our construction of the 2-D exponentially fitted methods can be illustrated by the left diagram in Figure 2.9. Such construction does not apply to 3-D problems because the mapping P k+1 →Q h is not one-to-one. Consequently, even a computable basis ofV h has been found and exponentially fitted, it will be difficult in general to compute the basis for u h from these fitted basis functions, in contrast to 2-D cases.
and following the procedure identical to that in subsection 2.1 to finally get the following basis functions for u h |K: x x0 e βφ(s,y0,z0) v x i (s, y 0 , z 0 )ds + y y0 e βφ(x,t,z0) v y i (x, t, z 0 )dt+ z z0 e βφ(x,y,r) v z i (x, y, r)dr , (2.37) where The total four parameters, i.e., u j (x 0 , y 0 , z 0 ) and the three coefficients m j,i , will be determined by enforcing the interpolative constraints at four nodes, which can be the four nodes of the tetrahedral or the centers of its four faces. A linear system similar to Eq.(2.21) can be formed and its solvability can be established through a direct computation. The space span{u j } is nonconforming similar to the 2-D cases.
3. Convergence of Exponentially Fitted Methods based on RT 0 0 . In this section we shall give a convergence analysis of the Galerkin finite element approximation of the problem (2.5) using the nonconforming exponentially fitted basis functions (2.11) or (2.35). We consider the symmetrized form of the drift-diffusion equation, and adopt the approach in [31] Fig. 2.9. Illustration of the construction of 2-D exponentially fitted methods (left). In 3-D (right), the mapping between P k+1 andQ h is not one-to-one so constructing basis for u h from exponentially fittedV h for arbitrary k is difficult.
nonconforming finite element spaces and will be reported separately. Notice that the interpolative basis functions of the nonconforming finite element spaces for the Slotboom variable ρ h can be obtained from u j (x, y) through scaling. Let V = H 1 0 (Ω) the weak solution of the problem (2.5) can be uniquely solved from the following weak formulation where the coercive, symmetric bilinear form a(ρ, v) : V × V → R is given For the purpose of Galerkin finite element approximation of this weak formulation we define the finite element space By construction ρ j is continuous only at selected nodes but in general discontinuous on edges, hence V h is not a subspace of V in general. Following the standard treatments for nonconforming finite element methods (E.g., Chapter 10 on [7]) we introduce the space V + V h in which we can define discrete bilinear and linear form, as It can be seen that this norm becomes the semi-norm |v| 1,Ω for v ∈ V , which is itself a norm due to the Poincáre inequality. To check that v h = 0 if and only if v = 0 for v ∈ V h we notice that v must be a constant for |v| 1,K = 0. Since v ∈ V h is continuous at nodes by construction the constant must be the same for all elements, and hence shall be zero for v = 0 on ∂Ω. It follows immediately that a h (ρ h , v h ) is continuous and coercive on V + V h because The nonconforming Galerkin finite element approximation of problem (1.3) finally reads: There are two major components in the convergence analysis of the nonconforming exponentially fitted finite element method, one is the interpolation error in V h , and the other is the consistency error arising from the nonconformity of the method, according to the second Strang lemma. The analysis in [31] follows this standard approach. Estimation of the interpolation error is obtained through a decomposition where Π h ρ is the interpolation of ρ in V h , and Π 1 ρ is the P 1 conforming interpolation of ρ in T h . This decomposition is applicable to our constructions with interpolative constraints enforced at the same set of nodes as those for conforming Lagrange finite element spaces. If the edge nodes are chosen for the constraint enforcement we will have to adopt a decomposition alternative to (3.6) (Ω) and let {T h } be a family of regular triangulations ofΩ. Then there exist positive constants C 1 (φ), C 2 (φ) and C 3 (φ) that depend only on φ, such that ρ − ρ h h ≤ C 1 (φ)h|ρ| 2,Ω + C 2 (φ) ∇φ ∞,Ω |Π 1 ρ| 1,Ω + C 3 (φ)h ∇φ ∞,T h ρ 2,Ω . (3.9) Proof. In [31].
4. Summary and Future Work. We proposed a general approach to construct exponentially fitted methods for solving 2-D drift-diffusion equations. Our constructions are based on the divergence-free subspace of H(div; Ω), usually RT 0 k . The basis functions u j for the density function are computed from the basis of RT 0 k , and through enforcement of interpolative constraints at selected node sets, giving rise to nonconforming finite element spaces. We showed that our approach can reproduce some of the previous constructions of the exponentially fitted methods, and will recover the standard conforming or nonconforming finite element spaces when the electric potential is zero. Thanks to the one-one mapping between RT 0 k and P k+1 , our approach can be used to develop arbitrary high-order 2-D exponentially fitted methods from RT 0 k . Extension of our approach for high-order 3-D methods is difficult because there does not exist a similar one-one correspondence between k th order divergence-free vector space and P k+1 , but construction of first-order 3-D exponentially fitted method is possible.
General convergence theory for the high-order methods can be established following decomposition (3.6), the estimate of P k conforming interpolation errors, the estimate of Π h ρ − Π k ρ, and the consistency error of high-order nonconforming finite element methods. These results will be reported in a forthcoming paper. We will also conduct extensive numerical experiments of the methods, in particular the high order methods, and apply them to realistic biomolecular drift-diffusion problems.