Steady-state Simulation of Semiconductor Devices using Discontinuous Galerkin Methods

Design of modern nanostructured semiconductor devices often calls for simulation tools capable of modeling arbitrarily-shaped multiscale geometries. In this work, to this end, a discontinuous Galerkin (DG) method-based framework is developed to simulate steady-state response of semiconductor devices. The proposed framework solves a system of Poisson equation (in electric potential) and drift-diffusion equations (in charge densities), which are nonlinearly coupled via the drift current and the charge distribution. This system is decoupled and linearized using the Gummel method and the resulting equations are discretized using a local DG scheme. The proposed framework is used to simulate geometrically intricate semiconductor devices with realistic models of mobility and recombination rate. Its accuracy is demonstrated by comparing the results to those obtained by the finite volume and finite element methods implemented in a commercial software package.


I. INTRODUCTION
Simulation tools capable of numerically characterizing semiconductor devices play a vital role in device/system design frameworks used by the electronics industry as well as various related research fields [1]- [8].Indeed, The authors are with the Division of Computer, Electrical, and Mathematical Sciences and Engineering, King Abdullah University of Science and Technology (KAUST), Thuwal 23955-6900, Saudi Arabia (e-mails:{liang.chen,hakan.bagci}@kaust.edu.sa).
in the last several decades, numerous commercial and open source technology computer aided design (TCAD) tools, which implement various transport models ranging from semi-classical to quantum mechanical models, have been developed for this purpose [9].Despite the recent trend of device miniaturization that requires simulators to account for quantum transport effects, many devices with larger dimensions (at the scale of 1µm) and with more complex geometries are being designed and implemented for various applications.Examples of these nanostructured devices range from photodiodes and phototransistors to solar cells, light emitting diodes, and photoconductive antennas [10].Electric field-charge carrier interactions on these devices can still be accurately accounted for using semi-classical models, however, their numerical simulation in TCAD raises challenges due to the presence of multi-scale and intricate geometric features.
Among the semi-classical approaches developed for modeling charge carrier transport, the drift-diffusion (DD) model is among the most popular ones because of its simplicity while being capable of explaining many essential characteristics of semiconductor devices [1]- [3].
One well-known challenge in using the DD model is the exponential variation of carrier densities, which renders standard numerical schemes used for discretizing the model unstable unless an extremely fine mesh is used.This challenge traces back to the convection-dominated convection-diffusion equations, whose solutions show sharp boundary layers.Various stabilization techniques have been proposed and incorporated with different discretization schemes [11]- [21].The Scharfetter-Gummel (SG) method [11] has been one of the workhorses in semiconductor device modeling; it uses exponential functions to approximate the carrier densities so that the fine mesh requirement can be alleviated.The SG method has been first proposed for finite difference discretization, and then generalized to finite volume method (FVM) [12]- [17] and finite element method (FEM) [18]- [21].
Although significant effort has been put into the numerical solution of the convection-dominated convection-diffusion problem in the last three decades, especially in the applied mathematics community, a fully-satisfactory numerical scheme for general industrial problems is yet to be formulated and implemented, for example see [27], [28], [31]- [33] for surveys.
Recently, the discontinuous Galerkin (DG) method has attracted significant attention in several fields of computational science [34]- [38].DG can be thought of as a hybrid method that combines the advantages of FVM and FEM.It uses local high-order expansions to represent/approximate the unknowns to be solved for.
Each of these expansions is defined on a single mesh element and is "connected to other expansions defined on the neighboring elements via numerical flux.This approach equips DG with several advantages: The order of the local expansions can be changed individually, the mesh can be non-conformal (in addition to being unstructured), and the numerical flux can be designed to control the stability and accuracy characteristics of the DG scheme.More specifically, for semiconductor device simulations, the instability caused by the boundary layers can be alleviated without introducing much numerical diffusion.We should note here that for a given order of expansion p, DG requires a larger number of unknowns than FEM.However, the difference decreases as p gets larger, and for many problems, DG benefits from hand/or p-refinement schemes [36], [38] and easily compensate for the small increase in the computational cost.
Those properties render DG an attractive option for multi-scale simulations [29], [34]- [39], and indeed, time domain DG has been recently used for transient semiconductor simulations [40]- [42].However, in device TCAD, the non-equilibrium steady-state response of semiconductor devices is usually the most concerned case and it is computationally very costly to model in time domain because the simulation has to be executed for a very large number of time steps to reach the steady-state [2], [43].
The steady-state simulation calls for solution of a nonlinear system consisting of three coupled second-order elliptic partial differential equations (PDEs).The first of these equations is the Poisson equation in scalar potential and the other two are the convection-diffusion type DD equations in electron and hole densities.These three equations are nonlinearly coupled via the drift current and the charge distribution.The charge-density dependent recombination rate, together with the fielddependent mobility and diffusion coefficients, makes the nonlinearity even stronger.In this work, for the first time, a DG-based numerical framework is formulated and implemented to solve this coupled nonlinear system of equations.More specifically, we use the local DG (LDG) scheme [45] in cooperation with the Gummel method [46] to simulate the non-equilibrium steady-state response of semiconductor devices.To construct the (discretized) DG operator for the convection-diffusion type DD equations (linearized within the Gummel method), the LDG alternate numerical flux is used for the diffusion term [47] and the local Lax-Friedrichs flux is used for the convection term.Similarly, the discretized DG operator for the Poisson equation (linearized within the Gummel method) is constructed using the alternate numerical flux.The resulting DG-based framework is used to simulate geometrically intricate semiconductor devices with realistic models of the mobility and the recombination rate [2].Its accuracy is demonstrated by comparing the results to those obtained by the FVM and FEM solvers implemented within the commercial software package COMSOL [30].We should note here that other DG schemes, such as discontinuous Petrov Galerkin [53], hybridizable DG [48], exponential fitted DG [51], and DG with Lagrange multipliers [52] could be adopted for the DG-based framework proposed in this work.
The rest of the paper is organized as follows.Section II starts with the mathematical model where the coupled nonlinear system of Poisson and DD equations is intro-duced, then it describes the Gummel method and provides the details of the DG-based discretization.Section III demonstrates the accuracy and the applicability of the proposed framework via simulations of two realistic device examples.Finally, Section IV provides a summary and discusses possible future research directions.

A. Mathematical Model
The DD model describes the (semi-classical) transport of electrons and holes in an electric field under the drift-diffusion approximation [1], [2].It couples the Poisson equation that describes the behavior of the (static) electric potential and the two continuity equations that describe the behavior of electrons and holes.This (coupled) system of equations reads where r represents the location vector, n e (r) and n h (r) are the electron and hole densities, ϕ(r) is the electric potential, J e (r) and J h (r) are the electron and hole current densities, ε(r) is the dielectric permittivity, q is the electron charge, and R(n e , n h ) is the recombination rate.In (15) and other equations in the rest of the text, s ∈ {e, h}, and the upper and lower signs should be selected for s = e and s = h, respectively.The current densities J s (r) are given by where µ e (E) and µ h (E) are the (field-dependent) electron and hole mobilities, d e (E) = V T µ e (E) and are the electron and hole diffusion coefficients, respectively, V T = k B T /q is the thermal voltage, k B is the Boltzmann constant, T is the absolute temperature, and is the (static) electric field intensity.Inserting (3) into (2) yields The recombination rate R(n e , n h ) describes the recombination of carriers due to thermal excitation and various scattering effects.In this work, we consider the two most common processes, namely the trap assisted recombination described by the Shockley-Read-Hall (SRH) model [2] as and the three-particle band-to-band transition described by the Auger model [2] as Here The mobility models have a significant impact on the accuracy of semiconductor device simulations.Various field-and temperature-dependent models have been developed for different semiconductor materials and different device operating conditions [1], [2], [30], [49], [50].
Often, high-field mobility models, which account for the carrier velocity saturation effect, are more accurate [2], [30], [49], [50].In this work, we use the Caughey-Thomas model [2], which expresses µ e (E) and µ h (E) where E (r) is amplitude of the electric field intensity parallel to the current flow, µ 0 e and µ 0 h are the low-field electron and hole mobilities, respectively, and V sat s , β e and β h are fitting parameters obtained from experimental data.
The electric field moves the carriers through the drift term in the expressions of J e (r) and J h (r) [first term in (3)].The carrier movements change n e (r) and n h (r), which in turn affect E(r) through the Poisson equation [see (1)].Furthermore, R(n e , n h ) [in (6)] and µ e (E) and µ h (E) [in (7)] are nonlinear functions of n e (r) and n h (r), and E(r), respectively.This system can be solved using either a decoupled approach such as the Gummel method or a fully-coupled scheme such as the direct application of the Newton method [2], [22].The Gummel methods memory requirement and computational cost per iteration are less than those of the Newton method.
In addition, accuracy and stability of the solution obtained by the Gummel method are less sensitive to the initial guess [2], [22].On the other hand, the Gummel method converges slower, i.e., takes a higher number of iterations to converge to the solution [2], [22].Since the simulations of the nanostructured devices considered in this work are memory-bounded, we prefer to use the Gummel method.
The Gummel iterations operate as described next and shown in Fig. 1.To facilitate the algorithm, we first introduce the quasi-Fermi potentials [1], [2], [22] "Inverting (8) for n e (r) and n h (r), respectively, and inserting the resulting expressions into (1) yield Newton method [1], [2], [22] (see below).The Gummel method decouples the NLP equation and the DD equations (2); the nonlinearity is "maintained solely in the NLP equation and the DD equations are treated as linear systems [1], [2], [22] as shown by the description of the Gummel method below.To solve the NLP equation in ( 9), we write it as a root-finding problem The Frechet derivative of F (ϕ(r), ϕ e (r), ϕ h (r)) with respect to ϕ(r) is The root-finding problem ( 10) is solved iteratively as where subscript "t refers to the variables at iteration t.
In (12), δ t+1 ϕ (r) is obtained by solving where ϕ t (r) is the solution at iteration t (previous iteration), ϕ t e (r) and ϕ t h (r) are computed using using n t e (r) and n t h (r) in (8).At iteration t = 0, initial guesses for ϕ t (r), n t e (r) and n t h (r) are used to start the iterations.Note that, in practice, one can directly compute ϕ t+1 (r) without using the variable δ t+1 ϕ (r).This is done by adding F (ϕ t (r), ϕ t e (r), ϕ t h (r); ϕ t (r)) to both sides of (13), and using (4) and the fact that which result in the coupled system of equations in unknowns φ t+1 (r) and E t+1 (r) Here, and are known coefficients obtained from the previous iteration.
Unknowns ϕ t+1 (r) and E t+1 (r) are obtained by solving (14).Then, µ e (E t+1 ) and µ h (E t+1 ) are computed using E t+1 (r) in (7).Finally, n t+1 e (r) and n t+1 h (r) can be obtained by solving where R(n t e , n t h ) on the right hand side is computed using n t e (r) and n t h (r) (from previous iteration) in (6).Note that a "lagging technique may also be applied to R(n t e , n t h ) to take advantage of the solutions at the current iteration.This technique expresses R(n e , n h ) as a summation of functions of n t e (r) and n t h (r) and n t+1 e (r) and n t+1 h (r), and moves the functions of n t+1 e (r) and n t+1 h (r) to the left hand side of (15).More details about this technique can be found in [44].
At this stage of the iteration, ϕ t+1 (r), n t+1 e (r) and n t+1 h (r) are known; one can use these to compute ϕ t+1 e (r) and ϕ t+1 h (r) and move to the next iteration.Convergence of the iterations can be checked by either the residuals of ( 10) and (15) or by the difference between the solutions of two successive iterations.

C. DG Discretization
As explained in the previous section, at every iteration of the Gummel method, one needs to solve three linear systems of equations, namely ( 14) and ( 15) (s = e, h).This can only be done numerically for arbitrarily shaped devices.To this end, we use the LDG method [45], [47] to discretize and numerically solve these equations.We start with the description of the discretization of (14).
First, we re-write (14) in the form of the following boundary value problem In ( 16)- (19), ϕ(r) and E(r) are the unknowns to be solved for and Ω is the solution domain.Note that in LDG, E(r) is introduced as an auxiliary variable to reduce the order of the spatial derivative in (16).Here it is also a "natural unknown to be solved for within  16) and (17) with the Lagrange polynomials i (r), i = 1, . . ., N p , on element k and applying the divergence theorem to the resulting equation yield the following weak form is used in the interior of Ω.Here, averaging operators are defined as {a} = 0.5(a + +a − ) and {a} = 0.5(a + +a − ) and "jumps are defined as , where superscripts -and + refer to variables defined in element k and in its neighboring element, respectively.The vector β determines the upwinding direction of ϕ and (εE).In LDG, it is essential to choose opposite directions for ϕ and (εE), while the precise direction of each variable is not important [36], [38], [45].In this work, we choose β = n on each element surface.On boundaries of Ω, the numerical fluxes are choosen as ϕ * = f D and (εE) * = (εE) − on ∂Ω D and ϕ * = ϕ − and (εE) * = f N on ∂Ω N , respectively [47].
We expand ϕ k (r) and E ν k (r) with the same set of Lagrange polynomials i (r) where r i , i = 1, . . ., N p , denote the location of interpolating nodes, and ϕ i k and E ν,i k , ν ∈ {x, y, z}, k = 1, . . ., K, are the unknown coefficients to be solved for.
Substituting (22) and ( 23) into ( 20) and ( 21) for k = 1, . . ., K, we obtain a global matrix system Here, the global unknown vectors Φ = [ Φ1 , . . ., ΦK ] ], ν ∈ {x, y, z}.The dimension of (24) can be further reduced by substituting Ē = M −1 ( BE − Ḡ Φ) (from the second row) into the first row, which results in In ( 24) and ( 25), M g and M are mass matrices.M g is a K × K block diagonal matrix, where each N p × N p block is defined as M is also a K × K block diagonal matrix, where each block is a 3 × 3 block diagonal matrix with N p × N p identical blocks defined as ε is a diagonal matrix with entries (ε 1 , . . ., εK ), where where ∂Ω kk nν (r) i (r) j (r)dS respectively, ν ∈ {x, y, z}.The right hand side terms in (24) and ( 25) are contributed from the force term and boundary conditions and are expressed as The DD equations in (15) (within the Gummel method) are also discretized using the LDG scheme as described next.Note that, here, we only discuss the Here n(r) and q(r) are the unknowns to be solved for and Ω is the solution domain.The auxiliary variable q(r) is introduced to reduce the order of the spatial derivative.Following the same procedure used in the discretization of ( 14), we discretize the domain into nonoverlapping tetrahedrons and test equations ( 26) and (27) using Lagrange polynomials on element k.Applying the divergence theorem yield the following weak form: where n * , (dq) * , and (vn) * are numerical fluxes "connecting element k to its neighboring elements.Here, for the simplicity of notation, we have dropped the explicit dependency on r on element surfaces.For the diffusion term, the LDG alternate flux is used for the primary variable n * and the auxiliary variable (dq) * [45] Here, averages and "jumps, and the vector coefficient β are same as those defined before.For the drift term, the local Lax-Friedrichs flux is used to mimic the path of information propagation [36] (vn We note (dq) * and (vn) * are not assigned independently on ∂Ω R .
Expanding n k (r) and q ν k (r) with Lagrange polynomials i (r) where r i , i = 1, . . ., N p , denote the location of interpolating nodes, n i k and q ν,i k , ν ∈ {x, y, z} , k = 1, . . ., K are the unknown coefficients to be solved for.Substituting (32) and ( 33) into (30) and (31), we obtain a global Here, the global unknown vectors N = [ N1 , ..., NK ] T and Q 34) yields In ( 34) and ( 35), the mass matrix M , the gradient matrix Ḡ and the divergence matrix D are same as those defined before.d is a diagonal matrix with entries ( d1 , . . ., dK ), where dk = ( dx The block sparse matrix C has contribution from the third term (the volume integral) and the fourth term (the surface integral) in (30).Each block is of size N p × N p .

The volume integral term only contributes to diagonal blocks as Cvol
The surface integral term contributes to both the diagonal and off-diagonal blocks as respectively, where ν ∈ {x, y, z}, and k , ∂Ω kk , and θ k (j) are defined the same as before.
The right hand side terms in (34) are contributed from the force term and boundary conditions and are expressed as October 15, 2019 DRAFT

D. Sparse Linear Solver
The sparse linear systems ( 25) and ( 35) are constructed and solved in MATLAB.For small systems, one can use a direct solver.For large systems, when the number of unknowns is larger than 1 000 000 (when We note here that one can reuse the preconditioner throughout the Gummel iterations.Because the matrix coefficients change gradually between successive iterations, we can store the preconditioner obtained in the first iteration (t = 0) and reuse it as the preconditioner in the following few iterations.In practice, the preconditioner only needs to be updated when the convergence of the sparse iterative solver becomes slower than it is in the previous Gummel iteration.For the devices considered in this work, the number of Gummel iterations is typically less than 50 and we find the preconditioners of the initial matrices work well throughout these iterations.

III. NUMERICAL EXAMPLES
In this section, we demonstrate the accuracy and the applicability of the proposed DG-based framework by numerical experiments as detailed in the next two sections.We have simulated two practical devices and compared the results to those obtained by the COMSOL semiconductor module [30].

A. Metal-Oxide Field Effect Transistor
First, we simulate a metal-oxide semiconductor fieldeffect transistor (MOSFET).The device is illustrated in  I. To reach a relative error of 10 −2 , the SG-FVM solver uses a mesh with h = 1nm, the DG solver uses a mesh with h = 3nm, while the GLS-FEM solver requires h to be as small as 0.3nm.This sharp boundary layer of carriers is the reason why a very fine mesh is required to obtain accurate results  The SG-FVM solver requires the mesh to be admissible, which is often difficult to satisfy for 3D devices [21], [22], [27].Implementation of SG-FVM in COMSOL uses a prism mesh generated by sweeping triangles defined on surfaces (for 3D devices) [30].the boundaries where the solution changes fast, which is not possible to do using the prism mesh generated by sweeping triangles (Fig. 5).This mesh flexibility compensates for the larger number of unknowns required by the DG solver, which results from defining local expansions only connected by numerical flux.

B. Plasmonic-Enhanced Photoconductive Antenna
For the second example, we consider a plasmonicenhanced photoconductive antenna (PCA).The operation of a PCA relies on photoelectric effect: it "absorbs optical wave energy and generate terahertz (THz) shortpulse currents.Plasmonic nanostructures dramatically Fig. 6 illustrates the device structure that is optimized to enhance the plasmonic fields near the operating optical frequency [56].The semiconductor layer is LT-GaAs that is uniformly doped with a concentration of 10 16 cm −3 .
The substrate layer is semi-insulating GaAs.We should note here that it is crucial to employ the appropriate fielddependent mobility models to accurately simulate this device [57].The Caughey-Thomas model is used here.
Other material parameters same as those used in [57].within the COMSOL semiconductor module.Just like FEM, the proposed DG solver is higher-order accurate but it does not require the stabilization techniques (such as GLS and SUPG), which are used by FEM.The main drawback of the proposed method is that it requires a larger number of unknowns than FEM for the same geometry mesh.But the difference in the number of unknowns gets smaller with the increasing order of basis function expansion.Additionally, DG can account for non-conformal meshes and benefit from local h-/prefinement strategies.Indeed, we are currently working on a more "flexible version of the current DG scheme, which can account for multi-scale geometric features more efficiently by making use of these advantages.

the
Gummel method.Dirichlet and Neumann boundary conditions are enforced on surfaces ∂Ω D and ∂Ω N , and f D (r) and f N (r) are the coefficients associated with these boundary conditions, respectively.In (19), n(r) denotes the outward normal vector ∂Ω N .For the problems considered in this work, ∂Ω D represents the metal contact surfaces with f D (r) = V contact (r), where V contact (r) is the potential impressed on the contacts.The homogeneous Neumann boundary condition, i.e., f N (r) = 0, is used to truncate the simulation domain [55].To facilitate the numerical solution of the boundary value problem described by (16)-(19) (within the Gummel method), Ω is discretized into k non-overlapping tetrahedrons.The (volumetric) support of each of these elements is represented by Ω k , k = 1, . . ., K. Furthermore, let ∂Ω k denote the surface of Ω k and n(r) denote the outward unit vector normal to ∂Ω k .Testing equations ( p = (p + 1)(p + 2)(p + 3)/6 is the number of interpolating nodes, p is the order of the Lagrange polynomials and subscript ν ∈ {x, y, z} is used for identifying the components of the vectors in the Cartesian coordinate system.We note here ϕ k (r) and E k (r) denote the local solutions on element k and the global solutions on Ω are the sum of these local solutions.ϕ* and (εE) * are numerical fluxes "connecting element k to its neighboring elements.Here, the variables are defined on the interface between elements and the dependency on r is dropped for simplicity of notation/presentation.In LDG, the alternate flux, which is defined as[45]

 1 ,
y, z}.We note that ε(r) is assumed isotropic and constant in each element.Matrices Ḡ and D represent the gradient and divergence operators, respectively.For LDG, one can show that D = − ḠT[47].The gradient matrix Ḡ is a K × K block sparse matrix, where each block is of size 3N p × N p and has contribution from the volume integral term and the surface integral term in(21).The October 15dν j (r)dV, ν ∈ {x, y, z} .The surface integral term contributes to both the diagonal blocks Ḡkk and off-diagonal blocks Ḡkk , where k corresponds to the index of the elements connected to element k.Let ∂Ω kk be the interface connecting elements k and k , and let θ k (j) select the interface nodes from element k, θ k (j) =   r j ∈ Ω k , r j ∈ ∂Ω kk 0, otherwise .Then, the contributions from the surface integral term to the diagonal block and the off-diagonal blocks are Ḡsurf kk discretization of the electron DD equation (s = e) and that of the hole DD equation (s = h) only differs by the sign in front of the drift term and the values of physical parameters.To simplify the notation/presentation, we drop the subscript denoting the species (electron and hole).The electron DD equation in (15) is expressed as the following boundary value problem d(r) = d(E) and v(r) = µ(E)E(r) become known coefficients during the solution of (15) within the Gummel method.Dirichlet and Robin boundary conditions are enforced on surfaces ∂Ω D and ∂Ω R , and f D (r) and f R (r) are the coefficients associated with these boundary conditions, respectively.n(r) denotes the outward normal vector of the surface.For the problems considered in this work, represents electrode/semiconductor interfaces and, based on local charge neutrality [55], f D (r) = (C + √ C 2 + 4n i 2 )/2 and f D (r) = n 2 i /n s e for electron and hole DD equations, respectively.The homogeneous Robin boundary condition, i.e., f R (r) = 0, is used on semiconductor/insulator interfaces, indicating no carrier spills out those interfaces [55].
double precison on a computer with 128GB RAM), it is preferable to use sparse iterative solvers to reduce the memory requirement.During our numerical experiments, we have found that the generalized minimum residual (GMRES) (the MATLAB built-in function "gmres') outperforms other iterative solvers in execution time.Incomplete lower-upper (ILU) factorization is used to obtain a preconditioner for the iterative solver.The drop tolerance of the ILU is critical to keep the balance between the memory requirement and the convergence speed of the preconditioned iterative solution.A smaller drop tolerance usually results in a better preconditioner, however, it also increases the amount of fill-in, which increases the memory requirement.

Fig. 2 .Fig. 3 .
Fig. 2. The background is uniformly-doped p-type silicon and source and drain regions are uniformly-doped ntype silicon.The doping concentrations in p-and n-type regions are 10 17 cm −3 and 10 18 cm −3 , respectively.The source and drain are ideal Ohmic contacts attached to n-type regions.The gate contact is separated from the semiconductor regions by a silicon-oxide insulator layer.The dimensions of the device and the different material regions are shown in Fig. 3. Material parameters at 300K are taken from [54].Special care needs to be taken to enforce the boundary conditions [55].The DD equations are only solved in the semiconductor regions.Dirichlet boundary conditions are imposed on semiconductor-contact interfaces, where the electron and hole densities are determined from the local charge neutrality as n e = (C + √ C 2 + 4n i 2 )/2 and n h = n i 2 /n e , respectively.Homogeneous Robin boundary condition, which enforces zero net current flow, is imposed on semiconductor-insulator interfaces.Poisson equation is solved in all regions.Dirichlet boundary conditions that enforce impressed external potentials are imposed on metal contacts (the contact barrier is ignored

Figs. 4
Figs. 4 (a) and (b) and Figs. 4 (c) and (d), respectively, compare the electron density and electric field intensity distributions computed by the DG and GLS-FEM solvers on the plane z = 0 for V g = 3V and V d = 0.5V.Figs.4(a) and (b) illustrate the field-effect introduced by the gate voltage, i.e., a sharp conducting channel forms near the top interface facing the gate (y = 0.2µm).

− 3 from
this simulation.In Fig. 4 (b), the carrier density decays more slowly [compared to the result in Fig. 4(a)] away from the gate interface and suddenly drops to the Dirichlet boundary condition values at the bottom interface (y = 0).This demonstrates the unphysical smearing of the boundary (carrier) layers observed in GLS-FEM solutions.Figs. 4 (c) and (d) show the x-component of the electric field intensity distribution computed by the DG and the GLS-FEM solvers, respectively.One can clearly see that the solution computed by the GLS-FEM solver is smoother (compared to the DG solution) at the sharp corners of the gate.The unphysical effects, as demonstrated in Figs. 4 (b) and (d), result from the GLS testing, which lacks of control on local conservation law [30].

Fig. 4 .
Fig. 4. Electron density distribution computed on the plane z = 0 by (a) the DG solver (b) the GLS-FEM solver for gate voltage Vg = 3V and drain voltage V d = 0.5V.Electric field intensity distribution computed on the plane z = 0 by (c) the DG solver (d) the GLS-FEM solver for gate voltage Vg = 3V and drain voltage V d = 0.5V.

2 n FVM e 2 ,
Fig. 7 (a) shows the electron density distribution computed by the proposed DG solver.Fig. 7 (b) plots the electron density computed by the DG and the SG-FVM solvers along lines (x, y = 0.5µm, z = 0) and (x, y = 0, z = 0) versus x.The results agree well.The relative difference, which is defined as n DG e − n FVM e 2 n FVM e 2 , between the solutions obtained by the DG and the SG-FVM solvers is 0.78%.Here, . 2 denotes L2 norm and n DG e and n FVM e

Fig. 7 .
Fig. 7. (a) Electron density distribution computed by the DG solver on the plane z = 0 in the device region of the plasmonic PCA.(b) Electron density computed by the DG and the SG-FVM solvers along lines (x, y = 0.5µm, z = 0) and (x, y = 0, z = 0) versus x.
, n i is the intrinsic carrier concentration, τ e and τ h , n h ) = R SRH (n e , n h ) + R Auger (n e , n h )

TABLE I RELATIVE
ERROR IN I(V ) CURVES