1730 IEEE TRANSACTIONS ON ELECTRON DEVICES, VOL. 53, NO. 7, JULY 2006 Briefs Accurate Boundary Integral Calculation in Semiconductor Device Simulation Riccardo Gusmeroli and Alessandro S. Spinelli Abstract—An accurate method for the evaluation of boundary integrals in semiconductor device modeling is presented. The method is compared against more common algorithms and can achieve both high accuracy and low computational time. The range of applications is extremely wide, in- cluding contact currents and charges, carrier quantum probability fluxes, and heat fluxes. Index Terms—MOSFETs, partial differential equations (PDEs), semi- conductor device modeling. I. INTRODUCTION Evaluation of boundary integrals is a difficult task that occurs routinely in electron device simulations when, for example, current densities are calculated on the basis of an approximate solution of basic semiconductor equations [1]–[3]. In general, given a contact Γ i , fluxes to be calculated assume the following form: Φ i =  Γ i a∇φ · ndΓ i (1) where n is the unit outward normal of the domain boundary. The flux Φ i can represent different physical quantities, depending on the con- text where it is evaluated. For example, when a represents the dielectric constant and φ the electrostatic potential, Φ i equals the charge stored on the contact, which is useful for capacitance evaluation. Another notable example is when the current flowing through a contact must be evaluated in a drift-diffusion framework: In this context, a represents the electrical conductivity and φ the carrier quasi-Fermi potential. Difficulties in the numerical evaluation of (1) arise from singular- ities in spatial derivatives of the approximate solution φ h near the contact edges due to a change in the boundary condition type [4]. Accurate ad hoc methods are then required to evaluate the boundary integrals up to the required tolerance, which usually trades off accuracy with computational time. In this paper, we present a novel method for the boundary integral calculation and compare its performance in terms of accuracy and time against more common ones. The new method guarantees good accuracy and fast computational time and can thus be usefully implemented in device simulation tools. II. METHODS Calculation of Φ i is usually a postprocess task that starts once an approximate solution φ h of φ is available. In the following, we will focus on the two-dimensional (2-D) case (even if extension to the three-dimensional (3-D) case is straightforward) and assume φ h as piecewise linear on the domain mesh, choosing the (possibly stabilized) continuous Galerkin method for the discretization. Manuscript received November 11, 2005; revised April 6, 2006. The review of this brief was arranged by Editor C. McAndrew. The authors are with the Dipartimento di Elettronica e Informazione, Politecnico di Milano-IU.NET, 20133 Milano, Italy. Digital Object Identifier 10.1109/TED.2006.875806 Fig. 1. Schematic view of the simulated nMOSFET device. The simplest approach to the boundary integrals computation is to calculate the (piecewise constant) gradient of φ on the triangles adjacent to the terminal Γ i and then use a quadrature rule to approx- imate the integral (1). We will refer to this procedure as the simple method (SM). Though extremely simple and fast, the SM suffers from high inaccuracies due to the change in boundary condition type (from Dirichlet to Neumann) at the contact ends, making the result often unreliable. A better approach, which improves the accuracy of the calculation, is described in [5] for the case of contact current calculations (although extension to the general case can be inferred) and will be referred to as the Nanz method (NM). The main idea is to introduce a (globally defined) weighing function w i for each contact Γ i and then convert the surface integral (1) to a volume integral via the Green/Gauss theorem. The computational procedure requires the solution of a partial differential equation (PDE) for each w i and the evaluation of the following surface integral: Φ i =  Ω [∇w i · (a∇φ) - w i f ] dΩ (2) where Ω is the device domain. The advantage of this method in terms of accuracy comes from the fact that ∇w i can be chosen to be negligible where ∇φ suffers from accuracy problems, so that regions with more accurate solutions are more heavily weighed. The main drawback is instead the need for the solution of one PDE for each device contact, which is a computationally expensive task. In the following, these two methods are compared against a novel one, which is called the residue method (RM) and is based on the analysis in [6] and relies on properties of the chosen discretization scheme. The method is based on the expansion of φ on the basis functions omitted to satisfy the Dirichlet boundary conditions on the contacts and gives good accuracy and fast computation. We will defer the numerical details of the method to the Appendix, focusing instead on the results in the following section. III. NUMERICAL RESULTS A. Simulation Benchmark For ease of comparison, we simulated a bulk nMOSFET device, which is schematically shown in Fig. 1, using the same device simula- tor developed in [7]. For simplicity, quantum–mechanical corrections 0018-9383/$20.00 © 2006 IEEE