Element Edge Based Discretization for TCAD Device Simulation

—Technology computer-aided design (TCAD) semiconductor device simulators solve partial differential equations (PDE) using the ﬁnite volume method (FVM), or related methods. While this approach has been in use over several decades, its methods continue to be extended, and are still applicable for investigating novel devices. In this paper, we present an element edge based (EEB) FVM discretization approach suitable for capturing vector-ﬁeld effects. Drawing from a 2D approach in the literature, we have extended this method to 3D. We implemented this method in a TCAD semiconductor device simulator, which uses a generalized PDE (GPDE) approach to simulate devices with the FVM. We describe how our EEB method is compatible with the GPDE approach, allowing the modeling of vector effects using scripting. This method is applied to solve polarization effects in a 3D ferro capacitor, and a 2D ferroelectric ﬁeld-effect transistor. An example for ﬁeld-dependent mobility in a 3D MOSFET is also presented.


I. INTRODUCTION
T ECHNOLOGY computer-aided design (TCAD) semicon- ductor device simulators solve continuum partial differential equations (PDEs) of the drift-diffusion model for semiconductors on a discretized mesh [1].The finite volume method (FVM) is the most widely used approach for assembling the PDEs.In the TCAD device simulation literature, the FVM is referred as finite boxes or control volume and is closely related to finite differences [2].
The popularity of the FVM is due to the success of the Scharfetter-Gummel method (SGM) in calculating the electron and hole current densities from exponentially varying carrier densities [3], which exhibits superior numerical stability than other approaches [1], [4].While the finite element method (FEM) presents advantages in some circumstances [5], [6], the lack of an equivalent to the SGM often leads to hybrid FVM/FEM approaches when solving the device equations with an otherwise FEM solver [7]- [9].
Standard simulation models fit well into a node based approach, as volume integration and surface integration of the PDEs are accounted for in the calculation of quantities at the mesh nodes, and at the edges connecting to adjacent mesh nodes.
Advanced simulation models use vector fields and require information from node quantities off of the edge.For mobility models, this requires knowledge of electric field parallel and normal to the direction of current flow [10].For the ferroelectric effect, the vector components of electric field and polarization are considered [11], [12].Nodal approaches have been proposed for evaluation of these models, which consider all adjacent nodes to the node of interest [13], [14].These approaches consider a constant field over the entire control volume bounded by these nodes.
In contrast, an element edge based (EEB) approach has the advantage that the control volume is bounded by the sub volume along each edge of the element being considered.This makes it possible to consider phenomena occurring over smaller scales, such as for impact ionization models [15].Laux used a 2D EEB approach to model generalized mobility [10] and impact ionization [15], where vector fields are calculated for each edge of a triangular element.In [16], this discretization was compared with other approaches for impact ionization, and it was found to be less susceptible to changes in mesh size than the other methods they considered.This was attributed to the EEB approach using a smaller effective control volume for each element edge.In this paper, we adopt this approach for 2D and extend it for tetrahedral elements in 3D.
Generalized PDE (GPDE) simulators use an equation description as input from the user [4], [6], [14], [17].The goal is to allow the rapid development of new models and applications for the continuum based approach.
Many GPDE simulators do not consider vector effects, and equation assembly is often restricted to node or edge based quantities [4].Other simulators do not completely model equation sensitivities to vector effects.Some provide problem specific operators, limiting the general purpose nature of their approach.One group compared 3 TCAD simulation approaches to model the ferroelectric capacitance effect [14].They report better convergence for methods where more complete model sensitivities are considered, such as ensuring that the electric field is fully coupled with the field dependent polarization model.
In this paper, we describe an EEB approach to modeling vector-field effects in DEVSIM, an open source GPDE semiconductor device simulator, for both 2D and 3D device simulation [18].We describe the approach in Section II.In Section III, we present simulation examples for hysteresis effects in a ferroelectric capacitor in 3D and in a ferroelectric field effect transistor (FeFET) in 2D.In addition, we present simulation results for transverse electric field effects on mobility in a 3D MOSFET.

A. Equation Assembly
1) Newton Method: The Newton method requires the evaluation of the model equations, as well as the derivatives with respect to the solution variable [1].This results in a formulation: where J is the Jacobian, ∆x is update to the solution vector, and F is referred to as the right-hand side (RHS) vector.As the method converges at each iteration, |F | → 0 and |∆x| → 0, to the limits of floating point precision.For TCAD simulation, the nonlinear nature of the PDEs requires a reasonable initial guess, and accurate derivatives, to get convergence or a good convergence rate.The PDEs considered in this paper fit in the form where C a is an EEB quantity, A a is an edge quantity, and B a is a node quantity.The label F a refers to PDE a being considered, as there are multiple simultaneously coupled PDEs in semiconductor simulation 1 .
2) Node Assembly: Fig. 1 shows a mesh node in 2D surrounded by triangular elements.The B a term is integrated over the NodeVolume.The RHS entry for this scalar quantity is then: for each node i on the simulation mesh.In this, and subsequent equations, += represents addition of the term to an existing matrix or vector entry.Similarly, −= represents subtraction.For the Jacobian: where x i is a simulation variable on the same node.The Jacobian entry is placed in the row corresponding to equation a and the column corresponding to simulation variable x.
3) Edge Assembly : Fig. 2 shows the EdgeCouple over which flux is integrated for the A a term in (2).The calculation of EdgeCouple is described in Section II-A.4.
For each edge between adjacent nodes i, j where A ai,j is the vector quantity evaluated on a mesh edge, and EdgeCouple i,j is the cross section of the mesh edge. 1 This paper does not consider our approach to time integration methods or boundary conditions.

NodeVolume EdgeCouple
Fig. 1.The 2D mesh cell with an arbitrary number of triangle elements connected to node.The volume of a node is bounded by the perpendicular bisectors of each triangle [19].
EdgeNodeVolume Fig. 2. The 2D mesh cell depicting how flux is integrated along each edge [19].The EdgeCouple is the length of the perpendicular bisectors for the 2 triangles along the edge connecting nodes n0, n1 .
The Jacobian entries are: where x k is the simulation variable on one of the two nodes on the edge.4) Element Edge Assembly: We consider the C a contribution to (2), which is EEB, and is dependent on variables on all element nodes.Fig. 3 shows a 2D representation of the volume and area over which the integration occurs.In 2D, the ElementEdgeCouple on each edge is from the triangle circumcenter to the center of the edge.Fig. 4 shows the element edge volume in 3D, where the ElementEdgeCouple is from the circumcenter of the tetrahedral element, to the circumcenters of the element triangles connected to the element edge, to the center of the element edge.In both 2D and 3D, the EdgeCouple in Section II-A.3 is calculated by summing the The mesh cell in 2D.Each triangle is subdivided into 3 sub volumes.The labels for nodes en0, en1, en2 are specific to the element edge in the sub element bounded by the solid line [19].ElementEdgeCouple for all elements connected to the edge being considered.
Similar to ( 5)-( 6), the RHS contribution for C a is: where t is the index of the element containing the edge.Similar to ( 7)-( 8), the Jacobian entries are: where x k is a simulation variable on one of the nodes on the element.In 2D, this is the 2 nodes, en0, en1 , on the element edge being considered, and a third node, en2, on a triangular element, as shown in Fig. 3.In 3D, there are 2 nodes on the each element edge, en0, en1 , as well as 2 additional nodes on the tetrahedral element en2, en3 , as shown in Fig. 4.

5) Edge Volume Assembly :
It is also possible to do a volume integration using either the edge based or EEB assembly.This is useful for models such as the density gradient model [20] or impact ionization model [15].In equations, ( 9)-( 12), the ElementNodeVolume in Fig. 3 is substituted for ElementEdgeCouple and the assembly is only done with the += operation.

B. Element Edge Based Fields
1) Calculating the Vector Field: To calculate a vector field on an edge of a tetrahedral element, we begin by calculating the components along each mesh edge in the device region.For the case of electric field, this is where ψ i and ψ j are the potentials at nodes i, j , and EdgeLength i,j is the distance between them, as shown in Fig. 2. The derivatives with respect to the node potentials are: In Fig. 5, we show how these scalar fields are used to calculate E t0 , the vector field based on all edges connect to node 0. The electric field on each edge is assumed to be where E t0 is the field over the entire element, and s i,j is the unit vector along the edge connecting nodes i, j .
For tetrahedra with nodes 0, 1, 2, 3 , we use the x, y, z components of the unit vectors to calculate the electric field for all edges connected to node 0 which can be written more compactly as so that E t0 is found by solving this small linear system.The derivatives with respect to the node potentials are then In 3D, we treat the field for the edge between nodes 0, 1 using an average: where E t1 is calculated by considering node 1 as the node of interest in ( 16)-( 21).The derivative with respect to the node potential is then: where k refers to one of the 4 nodes on the element.We assume that the average employed in ( 22)-( 23) is appropriate for creating vectors from edge based models.It can be shown for electric field that E t0 = E t1 , so the average is valid.For cases like current density, which does not have this property, this assumption appears valid for the simulations in Section III-C.
In 2D, we use an edge average from [10], [15], based on a weighting the length of perpendicular bisectors of the triangular elements.Our simulator also provides the ability to customize the weighting.
Once the vector fields are attained, they are used to calculate the components for the element edge assembly in ( 9)- (12), by first converting to a scalar value, and then adding to J and F .Scalar models defined as an edge or element edge model, can also be used in calculations with the vector field components.This will be described in Section II-C, where the approach is described for the polarization vector in ferroelectric materials, as well as for the transverse electric field used for mobility models.
2) Computation: In practice, the S t0 term in ( 17) is calculated once for each node of a tetrahedral element and is assigned as S t0 or S t1 , based on the element edge of interest.Since only geometric information is stored, the matrix factorization is done once at the beginning of the simulation, and is solved for all EEB field calculations across all iterations of the simulation.It is possible to take advantage of high performance math routines, so that the matrix factorization and back substitution is done in an efficient manner [21].
3) Limitations on Element Type : In our approach, we have restricted the method to triangular elements in 2D and tetrahedral elements in 3D.The original finite box methods restricted themselves to box and cube elements [1], and were later expanded [2].
In 3D, it is necessary that the centers of the tetrahedra circumspheres are inside the element, so that the volume and surface integrals are calculated correctly [2].In 2D, the circumcenter of triangular elements should also be inside the circumcircle, or equivalently not obtuse [1].While there are methods to accommodate elements violating this condition, their presence often leads to less accurate results and convergence difficulties.
4) Derivative Accuracy: Good convergence behavior is often dependent on having accurate derivatives for the simulation models [1].We demonstrate this with the simulation example in Section III-A.Accurate derivatives are also important smallsignal analysis, and the impedance field method [22]- [24].

C. Generalized PDE Approach
In DEVSIM, models are implemented as symbolic expressions, and they are assembled and solved [17], [19].Symbolic differentiation is employed to get the required derivatives.The simulator accepts equations in the form of: Node models on nodes (Section II-A.2) Edge models on edges (Section II-A.3)Element edge models on element edges (Section II-A.4)The GPDE equations are input by calling commands from a Python script [25].The simulator is a module loaded by the interpreter, mostly written in C++ [17].
1) Ferroelectric Model: A ferroelectric model may be implemented directly, using scripting.For a ferroelectric insulator, the equation is: where P is the polarization.We first start with the edge model for electric field and its derivatives EF, EF:psi@n0, EF:psi@n1 where EF is evaluated using an expression for (13) and the rest using (14) for nodes 0, 1 , respectively.The user defines the model and its derivatives using the symbolic math engine in the simulator.
During evaluation, indexes en0, en1 reference nodes on the element edge and en2, en3 reference the other element nodes on the tetrahedron.The simulator tracks the element edge indexes on a per element edge basis, so that the same PDE expression will work for each element edge being considered.
The simulator then evaluates the symbolic expression for (24) using (PX+Permittivity * EF_x) * unitx + (PY+Permittivity * EF_y) * unity + (PZ+Permittivity * EF_z) * unitz where unitx, unity, and unitz are the unit vector components along the element edge and Permittivity is a parameter.The terms PX, PY, and PZ are components which computed in using expressions involving the current and previous solutions for the electric field.We utilized a widely accepted empirical function of hyperbolic tangent to simulate the hysteresis loop [11], [12], [26].Scripting was used to generate these terms for all of the x, y, and z components and the derivatives of the polarization vector and electric fields.
2) Field Dependent Mobility: To calculate the surface mobility for the simulations in Section III-C, we first calculate the electric field parallel to current flow as: where j ti,j is the current density vector calculated using where j 0,1 , j 0,2 , and j 0,3 are edge current densities using a suitable low field mobility, and then averaged onto the element edge using the same average as (22).The vector j t0,1 is used to find the direction of current flow, and is not the final current density.The transverse electric field is then calculated using: The expressions for |E ⊥ ti,j | are then used to calculate the surface components of the mobility model [27].The bulk and surface mobility components are then averaged together, and the velocity saturation model is applied to get the EEB mobility.The EEB current densities are then integrated using ( 9)-( 12), having been calculated using the SGM applied on each element edge.

A. Ferrocapacitor Simulation
TetGen [28] was used to generate a tetrahedral mesh of a cube and contacts on the top and bottom.Not all of the elements met the requirements mentioned in Section II-B.3.A script was written to mark bad elements for refinement, and it was possible to achieve a mesh meeting the simulator's requirements.
The mesh and physics are loaded via scripts into the simulator.Fig. 6 shows the compensated charge versus the applied voltage for a ferroelectric capacitor.From the initial zero bias condition, the voltage is ramped over the range between ±15 V .The hysteresis in the charge due to the polarization is apparent in the figure .As a test of the accuracy of our approach, we repeated the simulations with derivative terms with respect to nodes en2, en3 set to 0. This emulates the condition in GPDE simulators which evaluate vector effects over all adjacent nodes, but may restrict their derivative information on edges [4], [14].The primary convergence criteria in our simulator is relative error, which is calculated at each node i as: where m is the iteration number, and ψ m i is the new solution value after the update.The relative error is taken at the node with the largest value of r.Fig. 7 shows relative error versus iteration number for one of the bias points.With all derivatives accounted for, it takes 2 iterations to get a relative error near 10 −15 , which is the limit of accuracy for double precision floating point arithmetic.This number of iterations is expected, since the equation derivatives are constant with respect to bias.When the off edge derivatives are neglected, the iteration count is higher to to attain a relative error within 10 −10 .

B. FeFET Simulation
We simulated a 2D p-type FeFET with a top contact/bottom gate architecture.Ohmic contacts were set at the source and drain electrodes.A triangular mesh was generated using Gmsh [29].Extended precision math mode was enabled in the equation assembly to prevent numerical noise in the off state current results.Fig. 8 shows the drain current versus gate voltage.The simulation details and an analysis of the switching processes of FeFET with respect to its structure will be presented in [30].It is important to note that when the derivative off the edge on node en2 was set to 0, the simulation failed before the sweep of biases could finish.Therefore, the full set of derivatives we implemented using the EEB approach was necessary to run the simulations.

C. MOSFET Simulation
We simulated a 3D MOSFET, based on a doping profile and 2D structure for a 90 nm n-channel MOSFET archived in [31].For the simulations, we created a 3D mesh using Cubit [32].The structure is shown in Fig. 9.
Using the scripting interface, we implemented the Poisson, electron-continuity, and hole-continuity equations, with SRH recombination [1].We implemented the mobility model from [27], with surface mobility components dependent on the electric field normal to current flow.
After refinement, the bulk region had 270,693 tetrahedra and 50,446 nodes.Solving this region with the oxide regions, nitride regions, and a polysilicon gate electrode, the simulation matrix had 166,225 rows for the fully coupled drift-diffusion simulation.
We then compared the mobility models using these three cases: • Concentration dependent bulk mobility model • Bulk mobility with velocity saturation • Bulk and surface mobility model with velocity saturation where the normal electric field for surface mobility was implemented using (27).
Figure 10 show a contour plot of the velocity saturated mobility with both the bulk and surface models enabled.Figure 11 shows I DS versus V GS for the three cases.The I DS versus V GS dependence is shown in Fig. 12.Consistent with the theory, the E ⊥ dependent mobility results in a lowering of the drain current.2) is surrounded by oxide (5) and two nitride regions (3,4).The bulk region (1) has a 120 nm drawn gate length.The source and drain contacts are both 50 nm underneath the nitride regions.A body contact was placed on the bottom of the 60 nm silicon region.The oxide thickness is 4.9 nm and the device is 25nm thick.
Fig.3.The mesh cell in 2D.Each triangle is subdivided into 3 sub volumes.The labels for nodes en0, en1, en2 are specific to the element edge in the sub element bounded by the solid line[19].

Fig. 4 .
Fig. 4. The element cell in 3D.The tetrahedron is subdivided into 6 subvolumes along each edge.The element edge of interest has the nodes labeled en0, en1 .The other nodes are labeled en2, en3 .

Fig. 5 .
Fig.5.Depiction on how vector field is calculated from the scalar fields calculated along the 3 edges connected to node 0.

Fig. 7 .
Fig. 7.The relative error versus iteration number.The solid line is for derivatives for all nodes on the element.The dashed line is when only nodes on the edge are considered.

Fig. 9 .
Fig. 9.The 90 nm 3D MOSFET.The polysilicon gate (2) is surrounded by oxide(5) and two nitride regions(3,4).The bulk region (1) has a 120 nm drawn gate length.The source and drain contacts are both 50 nm underneath the nitride regions.A body contact was placed on the bottom of the 60 nm silicon region.The oxide thickness is 4.9 nm and the device is 25nm thick.

Fig. 10 .Fig. 12 .
Fig. 10.Contour plot of the bulk and surface electron mobility with velocity saturation for a bias of V GS = 1 V, V DS = 1 V and V SB = 0.