A high-performance simulation method for simulating total dose effects on MOSFET devices

Abstract

The present invention discloses a high-performance simulation method for simulating the total dose effect of a MOSFET device, and provides the electrical information inside and outside the MOSFET device under the total dose effect. The method first divides the geometric model of the MOSFET device by a free hexahedral mesh; then, the control volume finite element method (CVFEM) is used to discretize the drift-diffusion equations, wherein the flux term of the drift-diffusion equations is interpolated by the hexahedral vector basis function, the finite difference method and the Scharfetter-Gummel (SG) method; finally, the Newton iteration method is used to solve the discrete drift-diffusion equations to obtain the potential and carrier concentration field information inside the MOSFET device, and the current-voltage curve of the MOSFET device is obtained by post-processing. The simulation method of the present invention has the advantages of high confidence, high numerical accuracy, high numerical stability and scalability. It can be used for efficient and high-precision numerical simulation of multi-physical effects such as the total dose of MOSFET devices, thereby improving the performance, robustness and reliability of integrated circuits.

Description

Technical Field

[0001] The invention relates to a high-performance simulation method for simulating the total dose effect of semiconductor devices in the chip field, in particular to an integrated circuit numerical simulation method for MOSFET device design. Background Art

[0002] Numerical simulation is a practical computer-aided design tool today. It can be used to solve the carrier drift-diffusion transport process in semiconductor devices and solve many problems that require huge calculations, complex structures, and analytical methods that are difficult to obtain accurate results. Prior art, China's invention patent "Analysis method for electric and thermal integration of PIN limiter under high-power electromagnetic pulse" (publication number: CN105718618A) discloses a numerical simulation method for PIN limiter under high-power electromagnetic pulse. This method uses the time domain spectral element method to calculate the semiconductor drift-diffusion equation group, and then obtains the internal electric field and current density of the semiconductor to solve the heat conduction equation, and then obtains the temperature distribution. However, due to the special application requirements of semiconductor devices in space environments, simulation requirements for simulating the total dose effect are proposed, but the existing technology cannot meet the requirements of high precision and high efficiency.

[0003] The numerical simulation of semiconductor devices requires high-quality meshes. In the simulation of three-dimensional structures, tetrahedral meshes have been widely used (J.Xu, Z.Ma, H.Li, Y.Song, L.Zhang and B.Lu, "A Multi-Time-Step FiniteElement Algorithm for 3-D Simulation of Coupled Drift-Diffusion Reaction Process in Total Ionizing Dose Effect," in IEEE Transactions on SemiconductorManufacturing, vol.31, no.1, pp.183-189, Feb.2018, doi:10.1109 / TSM.2017.2779058.), but the splitting error caused by the mismatch between the mesh edge and the current density direction under the tetrahedral mesh is large, and the tetrahedral mesh has high complexity. Compared with the tetrahedral mesh, the hexahedral mesh has the advantages of higher numerical accuracy and lower mesh complexity.

[0004] In view of the mixed (elliptic / parabolic / hyperbolic) properties and diffusion-convection characteristics of the carrier drift-diffusion equation in semiconductor devices, a variety of discretization schemes have emerged, such as Scharfetter-Gummel (SG), Finite Box Integration SG method (FBSG), quasi-fermi FEM and streamline upwinding Petrov Galerkin method (SUPG). However, Quasi-Fermi FEM does not meet physical conditions in some cases, and other methods are limited to using high-quality grids because the use of low-quality grids is prone to non-physical oscillations in the solution. The use of SG and CVFEM methods to discretize the drift-diffusion equation can solve this limitation.

[0005] The drift-diffusion equations for semiconductor devices also have nonlinear and strong coupling characteristics, and the decoupled Gummel iteration method (Mauri, A. and Bortolossi, A. and Novielli, G. and Sacco, R, "3D FiniteElement Modeling and Simulation of Industrial Semiconductor Devices Including Impact Ionization," in Journal of Mathematics in Industry, vol. 5, 2015, doi: 10.1186 / s13362-015-0015-z.) is usually used to obtain numerical solutions. However, its convergence rate is linear, and it is usually used to simulate low-doped semiconductor devices. It may not converge in high injection and avalanche scenarios, and its number of iterations increases with the increase of applied voltage. The fully coupled Newton iteration method has better convergence when solving this nonlinear strong coupling system, and its number of iterations is insensitive to the applied voltage. It can be used to simulate highly doped semiconductor devices. Summary of the invention

[0006] In view of the defects in the prior art, the object of the present invention is to provide a high-performance simulation method for simulating the total dose effect of MOSFET devices, comprising the following steps:

[0007] Step 1: Build the geometric model of the MOSFET device and perform hexahedral meshing;

[0008] Step 2: Determine the mathematical form of the drift-diffusion equations describing the semiconductor potential, electron concentration and hole concentration field information for the segmented geometric model, wherein the drift-diffusion equations include Poisson's equation, current continuity equation and carrier drift-diffusion equation;

[0009] Step 3: Determine the simplified gate boundary conditions and ohmic contact boundary conditions of the carrier drift-diffusion equations;

[0010] Step 4: The drift-diffusion equations are discretized by the CVFEM method. During the discretization process, the flux terms of the potential and current density are interpolated using hexahedral vector basis functions, finite difference method and SG method.

[0011] Step 5: Use the fully coupled Newton iteration method to solve the discretized drift-diffusion equations;

[0012] Step 6: Continuously iterate step 5 until the solution of the drift-diffusion equations reaches the convergence condition, and obtain the potential, electron concentration and hole concentration field distribution in the semiconductor device.

[0013] Specifically, in step 2, the Poisson equation is:

[0014]

[0015] The current continuity equation is:

[0016]

[0017]

[0018] The carrier drift diffusion equation is:

[0019]

[0020]

[0021] In the formula, ε is the dielectric constant of silicon, φ is the electric potential, q is the unit charge, n and p are the electron concentration and hole concentration respectively. and are the ionized acceptor concentration and the ionized donor concentration, respectively; J n and J p are the electron current density and hole current density respectively, t is the time, R n and R p are the electron generation recombination rate and hole generation recombination rate respectively; E is the electric field, μ n and μ p are the electron mobility and hole mobility, respectively, and D n and D p are the electron diffusion coefficient and the hole diffusion coefficient, respectively.

[0022] Specifically, in step 3, the simplified gate boundary condition is:

[0023]

[0024]

[0025]

[0026] The ohmic contact boundary condition is:

[0027] φ=V app +φ eq

[0028]

[0029]

[0030] in, is the external normal vector of the hexahedron, ε ins is the dielectric constant of the gate oxide, d ins is the thickness of the gate oxide, V metal is the metal work function, V app is the applied gate voltage, Q t is the interface charge, φ eq is the built-in potential, n i is the intrinsic carrier concentration.

[0031] Specifically, step four is:

[0032] The drift-diffusion equations are expressed as:

[0033]

[0034] Where a, b, u, and h are 3×1 vectors and G is a 3×3 matrix:

[0035]

[0036]

[0037]

[0038] The simplified gate boundary conditions and ohmic contact boundary conditions are simplified as follows:

[0039]

[0040] Among them, u D is a 3×1 vector consisting of the three variables on the right side of the Ohmic contact boundary condition equation, β and γ are 3×1 vectors consisting of the three equations on the right side of the simplified gate boundary condition equation, Γ D is the Dirichlet frontier, Γ M For mixed boundaries;

[0041]

[0042]

[0043]

[0044] In step 4, the following integration is performed on the hexahedral control volume:

[0045]

[0046] Among them, C i represents any hexahedron Ω obtained by step 1 k The control volume corresponding to the vertex i of , the control volume is composed of the area surrounded by the centroids of all hexahedrons around the vertex i, the center of the face where the vertex i is located, and the midpoint of the edge with the vertex i as the endpoint; dV represents the volume integral, and dS represents the surface integral; Represents the control body C i The kth surface of j ∈Edge(Ω k ) represents the hexahedron Ω k Traverse all the edges of W j Represents a hexahedron Ω k The weight of the jth edge, (0) g j , (1) g j , (2) g j is a hexahedron Ω k The flux term of the jth edge is related to the electric field, electron current density, and hole current density, respectively, and is obtained as follows:

[0047]

[0048] In the formula, Hexahedron Ω k The two endpoints j of the jth edge 1 and j 2 The potential value, represents the length of the jth edge, The endpoints j 1 and j 2 The electron concentration, The endpoints j 1 and j 2 The hole concentration, B 1 , B 2 and V T Defined by:

[0049]

[0050]

[0051] In the formula, k b is the Boltzmann constant and T is the temperature.

[0052] Specifically, step five is:

[0053] Construct the following system equations:

[0054]

[0055]

[0056] Where N is the total number of vertices in the mesh, F i (u i ) and u i is a 3×1 vector, F(U) and U are 3N×1 vectors, where i=0,1,2,......,N-1;

[0057] The initial value of U 0 Elements in All are set to

[0058] in, and are the quasi-Fermi energies of electrons and holes, E c and E v are the conduction band bottom energy and the valence band top energy respectively;

[0059] Based on the initial value, U is iterated and the solution U is obtained through the Kth step. K Calculate the solution U of step K+1 K+1 ,as follows:

[0060] U K+1 =U K +rΔU K

[0061] Where r is the damping factor, which can be 1 or 0, ΔU K Solve by the following equation:

[0062] JΔU k =-F(U k )

[0063] Here J is a 3N×3N Jacobian matrix, In the formula, 0 in the brackets represents the first row of the corresponding vector, 1 and 2 represent the second and third rows respectively, and u j Represents the u matrix of vertex j associated with vertex i.

[0064] Specifically, in step 6, the convergence condition is: the absolute error of the potential between the solution of the previous iteration step and the solution of the current iteration step is less than 10 -10 The absolute errors of V, electron and hole concentrations are all less than 10 -5 m -3 .

[0065] Compared with the prior art, the present invention has the following significant advantages: (1) MOSFET and multi-finger MOSFET are solved using a physical model, and the electrical characteristic distribution inside the device under the total dose effect can be obtained. (2) A free hexahedral grid is used for discretization, which has the advantages of high numerical accuracy and low grid complexity. (3) The drift-diffusion equations are discretized using the CVFEM method, and the flux terms of the drift-diffusion equations are interpolated and calculated using the hexahedral vector basis function, the finite difference method, and the Scharfetter-Gummel method, which can maintain good stability under low grid quality conditions. (4) The Newton iteration method is used to solve the discretized drift-diffusion equations, which improves the convergence speed and maintains good stability under strong electric fields and high doping conditions. (5) The present invention has been tested to be able to solve a problem scale of at least 20 million degrees of freedom. BRIEF DESCRIPTION OF THE DRAWINGS

[0066] Figure 1 is the vertex V i Hexahedral control volume C i .

[0067] Figure 2 is the vertex V i The quadrilateral control volume C i .

[0068] Figure 3 It's the border And its corresponding area Ω izmn , where the first-order vector basis function W j Marked with a red arrow.

[0069] Figure 4 is the MOSFET geometric model.

[0070] Figure 5 is the electron concentration distribution of the MOSFET device.

[0071] Figure 6 is the gate voltage-drain current curve of a MOSFET device.

[0072] Figure 7 is the model configuration of a multi-finger MOSFET, where Figure 7 (a) is the geometrical dimension feature of multi-finger MOSFET. Figure 7 (b) is the physical model of the multi-finger MOSFET, and (c) is the hexahedral mesh used for calculation.

[0073] Figure 8 is the three-dimensional field distribution of the multi-finger MOSFET, where Figure 8 (a) is the potential distribution, (b) is the electron concentration n distribution, and (c) is the hole concentration h distribution.

[0074] Fig. 9 It is along Figure 8 The one-dimensional field distribution of the midline AB, where Fig. 9 (a) is the potential distribution, Fig. 9 (b) is the distribution of electron concentration n, Fig. 9 (c) is the hole concentration h distribution.

[0075] Fig.10 It is a specific flow chart of the method of the present invention. DETAILED DESCRIPTION

[0076] The present invention is further described in detail below in conjunction with the accompanying drawings. The following examples will help those skilled in the art to further understand the present invention, but are not intended to limit the present invention in any form. It should be noted that, for those of ordinary skill in the art, several changes and improvements can be made without departing from the concept of the present invention. These all belong to the protection scope of the present invention.

[0077] 1. The carrier drift-diffusion model in semiconductor devices can be described by partial differential equations such as Poisson's equation, carrier current continuity equation and carrier drift-diffusion equation.

[0078] Poisson's equation can be expressed as:

[0079]

[0080] The current continuity equation can be expressed as:

[0081]

[0082] The drift-diffusion equation can be expressed as:

[0083]

[0084] Where ε is the dielectric constant of silicon, φ is the electric potential, E is the electric field, n and p are the electron concentration and hole concentration, and are the ionization acceptor concentration and the ionization donor concentration, J n and J p are the electron current density and the hole current density.

[0085] Equations (1), (2) and (3) contain the distribution laws of the electric potential, electron concentration and hole concentration field in semiconductor devices. We can describe equations (1), (2) and (3) in a unified form:

[0086]

[0087] Here a, b, u, and h(u) are 3×1 vectors, and G(u) is a 3×3 matrix

[0088]

[0089]

[0090]

[0091] 2. The carrier drift-diffusion model in semiconductor devices meets specific boundary conditions. The interface boundary conditions between the gate oxide and the semiconductor can be described by simplifying the gate boundary conditions:

[0092]

[0093]

[0094]

[0095] In order to describe the total radiation dose damage to the gate oxide layer, the equivalent interface charge Q is usually used. t Applied to the gate oxide and semiconductor interface as boundary conditions.

[0096] The ohmic contact boundary condition can be written as:

[0097] φ=V app +φ eq (7a)

[0098]

[0099]

[0100] in is the external normal vector, ε ins is the dielectric constant of the gate oxide, d ins is the thickness of the gate oxide, V metal is the metal work function, V app is the applied gate voltage, Q t is the interface charge,

[0101] The simplified gate boundary condition (6) and the ohmic contact boundary condition (7) belong to the mixed boundary condition and the Dirichlet boundary condition respectively, so they can be simplified to the following form:

[0102]

[0103] Among them, u D is a 3×1 vector consisting of the three variables on the right side of equation (7), β and γ are 3×1 vectors consisting of the three variables on the right side of equation (6), and Γ D is the Dirichlet frontier, Γ M A mixed boundary.

[0104]

[0105]

[0106]

[0107] 3. In the process of numerical solution of multi-finger MOSFET, (4) needs to be discretized, and the discretization adopts CVFEM method. i like Figure 1 Here we use a two-dimensional example to better describe the discrete process. The three-dimensional case is the same as the two-dimensional case, but the three-dimensional integral calculation is used. Two-dimensional quadrilateral control volume C i (The corresponding node V i )like Figure 2 As shown, it consists of five quadrilaterals.

[0108] In the control volume finite element method, the node V i The control body C i like Figure 2 As shown, we integrate on this control volume and then apply the divergence theorem to obtain

[0109]

[0110] Without loss of generality, we take Figure 3 For example, the boundary Can be split into 10 boundaries Therefore, the second term of (10) can be discretized as:

[0111]

[0112] The flux term G in the figure can be obtained by using vector basis functions instead of scalar basis functions, because vector basis functions are more accurate than scalar basis functions in describing the current density flux term in the case of strong drift. We define the dth behavior of G as (d) g T (d=0,1,2), Middle (d) g can be obtained by interpolating the vector basis function in its corresponding unit. Figure 7 As shown, in Middle (d) g can be expressed as:

[0113]

[0114] (d) g j is the flux value at the jth edge (j=0,1,2,3), which depends on u, but W j It varies with space, such as Figure 7 As shown. j1j2 Represents the slave node V j1 To V j2 The j-th edge of can be obtained by the finite difference method and the SG method:

[0115]

[0116] here

[0117]

[0118]

[0119] as well as

[0120]

[0121] is the Bernoulli function.

[0122] C i The integral of can be written as:

[0123]

[0124] Using unconditionally stable backward differences, (16) can be transformed into solving the following equation:

[0125] F i (u t+Δt )=0 (17)

[0126] The steady-state equation of (17) can be simplified as:

[0127] F i (u)=0 (18)

[0128] Here (1) g j and (2) g j is a nonlinear function of the electric potential φ, which makes (17) or (18) a nonlinear system. Since it is impossible to obtain analytical solutions for the electric potential, electron concentration, and hole concentration, we use the fully coupled Newton iteration method to obtain its numerical solution.

[0129] 4. After assembling the control volume equation (17) or (18) corresponding to each vertex, the following system equation can be obtained:

[0130]

[0131]

[0132] Here F i and u i is a 3×1 vector (i=0,1,2,......,N-1), and F(U) and U are 3N×1 vectors.

[0133] In order to ensure the convergence of the iterative process, the initial value of the fully coupled Newton iterative method can be set as:

[0134]

[0135] here and It's quasi-Fermi energy.

[0136] At step k, the fully coupled method solves the following correlation equations:

[0137] JΔU k =-F(U k ) (twenty two)

[0138] Here J is a 3N×3N Jacobian matrix, in which the derivative terms are obtained by automatic differentiation,

[0139]

[0140] U k+1 This equation can be obtained

[0141] U k+1 =U k +rΔU k (twenty four)

[0142] 5. Through the hexahedral CVFEM method in specific implementation mode 3 and the Newton iteration method in specific implementation mode 4, the drift-diffusion equations of carriers in the semiconductor device can be solved, and then the potential, electron concentration and hole concentration field distribution in the semiconductor device can be obtained, and through post-processing, the current-voltage curve of the semiconductor device can be obtained.

[0143] 6. Establish the geometric model of the MOSFET device that needs to be solved, such as Figure 4 As shown in the figure, the geometric model is divided by a free hexahedral mesh. The doping model of the MOSFET device is set, and the acceptor concentration N a Set to 1017 cm -3 , donor concentration N b Set to 10 20 cm -3 , doping adopts Gaussian distribution, junction depth in X and Y direction is 0.2um, junction depth in Y direction is 0.1um, and junction depth in Z direction is 0.25um. Set source, drain, and bulk as ohmic contact boundary conditions, and gate as simplified gate boundary conditions. Apply 10mV voltage to the drain, 0V voltage to the source and substrate, scan the gate voltage from 0 to 2V, and set the radiation interface charge Q t 0,5×10 11 cm -3 , and 1×10 12 cm -3 The electron concentration distribution of the MOSFET device obtained during the calculation process is as follows: Figure 5 As shown in the calculation results, it can be seen that due to the total dose effect caused by the radiation interface charge, a conductive channel appears. The gate voltage-drain current curve of the MOSFET device is obtained as follows Figure 6 As shown, it can be seen from the calculation results that due to the total dose effect caused by the radiation interface charge, the threshold voltage drift and the leakage current increase occur.

[0144] According to Figure 7 Based on the geometric features shown in (a), the geometric model of the multi-finger MOSFET device to be solved is established, as shown in Figure 7(b). The geometric model is divided using a free hexahedral mesh, as shown in Figure 7 As shown in (c), the hexahedral mesh used for calculation contains 1234180 grid points and 1168464 grid cells. Therefore, 3702540 degrees of freedom need to be solved in the Newton iteration calculation. Set the doping model of the multi-finger MOSFET device, the acceptor concentration N a Set to 10 15 cm -3 , donor concentration N b Set to 5×10 20 cm -3 , doping adopts Gaussian distribution, junction depth in X and Y direction is 0.04um, junction depth in Z direction is 0.015um. Set source (Source), drain (Drain), substrate (Bulk) as ohmic contact boundary conditions, gate (Gate) as simplified gate boundary conditions. Apply 10mV voltage to drain, 0V voltage to source and substrate, 0.5V voltage to gate, set radiation interface charge Q t is 0 and 1×10 12 cm -3 . Figure 8 It is the three-dimensional distribution of electric potential, electron concentration and hole concentration. Fig. 9 It is the one-dimensional distribution of potential, electron concentration and hole concentration along line AB. It can be observed that the potential, electron concentration and hole concentration change due to the total dose effect.

[0145] The above describes the specific embodiments of the present invention. It should be understood that the present invention is not limited to the above specific embodiments, and those skilled in the art may make various changes or modifications within the scope of the claims, which will not affect the essence of the present invention.

Claims

1. A high-performance simulation method for simulating the total dose effect of MOSFET devices, characterized in that: The steps include: Step 1: Build the geometric model of the MOSFET device and perform hexahedral meshing; Step 2: Determine the mathematical form of the drift-diffusion equations describing the semiconductor potential, electron concentration and hole concentration field information for the segmented geometric model. The drift-diffusion equations include the Poisson equation, the current continuity equation and the carrier drift-diffusion equation. The Poisson equation is: The current continuity equation is: The carrier drift diffusion equation is: In the formula, ε is the dielectric constant of silicon, φ is the electric potential, q is the unit charge, n and p are the electron concentration and hole concentration respectively. and are the ionized acceptor concentration and the ionized donor concentration, respectively; J n and J p are the electron current density and hole current density respectively, t is the time, R n and R p are the electron generation recombination rate and hole generation recombination rate respectively; E is the electric field, μ n and μ p are the electron mobility and hole mobility, respectively, and D n and D p are the electron diffusion coefficient and the hole diffusion coefficient, respectively; Step 3: Determine the simplified gate boundary conditions and ohmic contact boundary conditions of the carrier drift-diffusion equations. The simplified gate boundary conditions are: The ohmic contact boundary condition is: φ=V app +φ eq in, is the external normal vector of the hexahedron, ε ins is the dielectric constant of the gate oxide, d ins is the thickness of the gate oxide, V metal is the metal work function, V app is the applied gate voltage, Q t is the interface charge, φ eq is the built-in potential, n i is the intrinsic carrier concentration; Step 4: The drift-diffusion equations are discretized by the CVFEM method. During the discretization process, the flux terms of potential and current density are interpolated using hexahedral vector basis functions, finite difference method and SG method, as follows: The drift-diffusion equations are expressed as: Where a, b, u, and h are 3×1 vectors and G is a 3×3 matrix: The simplified gate boundary conditions and ohmic contact boundary conditions are simplified as follows: Among them, u D is a 3×1 vector consisting of the three variables on the right side of the Ohmic contact boundary condition equation, β and γ are 3×1 vectors consisting of the three equations on the right side of the simplified gate boundary condition equation, Γ D is the Dirichlet frontier, Γ M For mixed boundaries; In step 4, the following integration is performed on the hexahedral control volume: Among them, C i represents any hexahedron Ω obtained by step 1 k The control volume corresponding to the vertex i of , the control volume is composed of the area enclosed by the centroids of all hexahedrons around the vertex i, the center of the face where the vertex i is located, and the midpoint of the edge with the vertex i as the endpoint; dV represents the volume integral, and dS represents the surface integral; Represents the control body C i The kth surface of j ∈Edge(Ω k ) represents the hexahedron Ω k Traverse all the edges of W j Represents a hexahedron Ω k The weight of the jth edge, (0) g j , (1) g j , (2) g j is a hexahedron Ω k The flux term of the jth edge is related to the electric field, electron current density, and hole current density, respectively, and is obtained as follows: In the formula, Hexahedron Ω k The potential values ​​of the two endpoints j1 and j2 of the jth edge are: represents the length of the jth edge, are the electron concentrations at endpoints j1 and j2, are the hole concentrations at endpoints j1 and j2, B1, B2, and V T Defined by: In the formula, k b is the Boltzmann constant, T is the temperature; Step 5: Use the fully coupled Newton iteration method to solve the discretized drift-diffusion equations, specifically: Construct the following system equations: Where N is the total number of vertices in the mesh, F i (u i ) and u i is a 3×1 vector, F(U) and U are 3N×1 vectors, where i=0,1,2,......,N-1; The initial value of U 0 Elements in All are set to in, and are the quasi-Fermi energies of electrons and holes, E c and E v are the conduction band bottom energy and the valence band top energy respectively; Based on the initial value, U is iterated and the solution U is obtained through the Kth step. K Calculate the solution U of step K+1 K+1 ,as follows: U K+1 =U K +rΔU K Where r is the damping factor, which can be 1 or 0, ΔU K Solve by the following equation: JΔU k =-F(U k ) Here J is a 3N×3N Jacobian matrix, In the formula, 0 in the brackets represents the first row of the corresponding vector, 1 and 2 represent the second and third rows respectively, and u j The u matrix representing vertex j associated with vertex i; Step 6: Continuously iterate and solve step 5 until the solution of the drift-diffusion equations reaches the convergence condition, and obtain the potential, electron concentration and hole concentration field distribution in the semiconductor device. The convergence condition is: the absolute error of the potential between the solution of the previous iteration step and the solution of the current iteration step is less than 10 -10 The absolute errors of V, electron and hole concentrations are all less than 10 -5 m -3 .

Images