A semiconductor multi-physics effect simulation method based on hybrid discontinuous finite element

Abstract

The present invention discloses a semiconductor multi-physical effect simulation method based on hybridizable discontinuous finite elements. The method first constructs a high-order basis function by using a tetrahedral unit-based mesh partitioning for 2D / 3D semiconductor devices. Then, the drift-diffusion equation system of carriers in the device structure is decoupled and split into a nonlinear Poisson equation and a linearized drift-diffusion equation group; then, the hybridizable discontinuous Galerkin HDG method is used to numerically discretize the two equations respectively; the flux term in the equation is described using a stability parameter, and the Gummel iteration method is applied to solve it to obtain information such as the carrier concentration, potential and electric field distribution in the device. The method of the present invention has the advantages of supporting arbitrary high-order basis functions, high numerical accuracy, good numerical stability and support for large-scale parallelism. It can be used for efficient and high-precision numerical calculation of physical effects in complex semiconductor device structures, and for high-performance design optimization of radio frequency, microwave and millimeter wave integrated circuits.

Description

Technical Field

[0001] The present invention relates to a multi-physics effect simulation method in semiconductor device structures based on high-order hybridizable discontinuous Galerkin (HDG) in the field of semiconductor devices and integrated circuits, and is particularly used for numerical simulation of radio frequency, microwave and millimeter wave semiconductor devices such as diodes and transistors. Background Art

[0002] The physical effects in RF, microwave, and millimeter-wave semiconductor devices are complex and highly nonlinear, making reliable results difficult to obtain using analytical methods. Numerical calculation, however, is an extremely important research method. By establishing a carrier drift-diffusion model within the semiconductor device structure and selecting appropriate numerical calculation methods, accurate predictions of the physical effects within complex device structures can be achieved. Prior art, the Chinese invention patent "A Numerical Simulation Method and System for Rapid Ionization Devices" (Publication No.: CN113076669A) discloses a numerical simulation method and system for rapid ionization devices. This method discretizes the differential equations in the carrier drift-diffusion model based on the doping concentration values ​​at each grid point using the finite difference method. The method simultaneously solves the device characteristic parameters at each grid point, completing the numerical calculation of the device characteristics at different times during the device's dynamic triggering process. However, this method places high demands on the quality and quantity of the grid and is difficult to generalize to arbitrarily high-order basis functions.

[0003] In three-dimensional structural simulations, the numerical simulation of semiconductor devices requires high-quality tetrahedral and hexahedral meshes (DW Yang, W. Chen, W. S. Zhao, G. Zhu, Z. G. Zhao, J. Schutt Aine, and W. Y. Yin, “An Improved Algorithm for Drift Diffusion Transport and Its Application on LargeScale Parallel Simulation of Resistive Random Access Memory Arrays,” IEEE Access, pp. 31 273–31 285, 2019). However, for complex semiconductor device structures, mesh generation is extremely difficult, and poor mesh quality becomes a significant constraint in solving the drift-diffusion equations in semiconductor device structures. Compared to traditional methods, the use of arbitrary high-order basis functions can effectively reduce the requirements for mesh quantity and quality while maintaining numerical accuracy.

[0004] For carrier drift-diffusion models in semiconductor devices, various discretization and numerical calculation methods have emerged, including the Scharfetter-Gummel (SG), SGFinite-boxes (SG-FB) scheme, the Galerkin Least-square (GLS) method, and the Control Volume Finite Element Method with SG (CVFEM-SG). However, these methods can produce unphysical spurious solutions when solving with low-quality meshes, necessitating the use of high-quality meshes. The HDG method uses numerical fluxes to exchange information between cells, eliminating the need to impose tangential continuity conditions at the interfaces between cells. The stability of the system is controlled by the cell stability parameters. In addition, computational mathematical theory has proven its convergence on the drift-diffusion equation (NC CNguyen, J. Peraire, and B. Cockburn, “An Implicit High-order Hybridizable Discontinuous Galerkin Method for Nonlinear Convection–diffusion Equations,” Journal of Computational Physics, vol. 228, no. 23, pp. 8841–8855, 2009.), and it can be extended to arbitrary high-order basis functions. Summary of the Invention

[0005] In view of the defects in the prior art, the present invention aims to provide a semiconductor multi-physics effect simulation method based on hybrid discontinuous finite element, which comprises the following steps:

[0006] Step 1: Mesh the semiconductor device and determine the physical governing equations of the semiconductor device structure, including the carrier drift-diffusion equation, Poisson's equation, and current continuity equation, which describe the electron concentration, hole concentration, potential, and electric field distribution in the device.

[0007] Step 2: Determine the Ohmic contact boundary conditions and the initial electron concentration, hole concentration, and potential analytical distribution function corresponding to the physical governing equations of the semiconductor device structure;

[0008] Step 3: The carrier drift-diffusion equation, Poisson equation, and current continuity equation are decoupled using the Gummel iteration method, splitting them into the nonlinear Poisson equation and the linearized carrier drift-diffusion equation group;

[0009] Step 4: Using the hybridized discontinuous Galerkin (HDG) method, the nonlinear Poisson equation and the carrier drift-diffusion equations are numerically discretized respectively. Mixed variables are introduced to reduce the second-order partial differential equation to two first-order equations. Stability parameters are used to process the numerical flux terms between discontinuous units to obtain the discretized nonlinear Poisson equation and the carrier drift-diffusion equations.

[0010] Step 5: Use the Gummel iteration method to solve the discretized nonlinear Poisson equation and update the potential distribution. Once convergence conditions are reached, substitute the result into the discretized carrier drift-diffusion equations to solve and update the electron and hole concentration distributions.

[0011] Step 6: Repeat step 5 until the physical field information in the solved semiconductor device structure reaches the set convergence condition, and finally obtain the electron concentration, hole concentration, potential and electric field distribution characteristics in the semiconductor device structure.

[0012] Specifically, in the first step,

[0013] The Poisson equation is:

[0014]

[0015] The carrier drift-diffusion equations are:

[0016]

[0017]

[0018] Where ε is the dielectric constant of the material in the device, φ is the potential, q is the charge, n and p are the electron concentration and hole concentration respectively, N A and N D are the ionized acceptor concentration and the ionized donor concentration, R n and R p are the electron generation recombination rate and hole generation recombination rate, D n and D p are the electron diffusion coefficient and hole diffusion coefficient, μ n and μ p are electron mobility and hole mobility respectively, and t is time.

[0019] Specifically, in the second step, the ohmic contact boundary condition is:

[0020] φ=V app +φ eq (3a)

[0021]

[0022]

[0023] Where V app is the voltage applied at the boundary, φ eq is the built-in potential, N0 is the ionization doping concentration, n i is the intrinsic carrier concentration.

[0024] Specifically, the fourth step is:

[0025] Note that Poisson's equation (1) is mathematically equivalent to the drift-free The electron drift-diffusion equation (2a) or the drift-free equation The hole drift-diffusion equation (2b) is derived from the hybridizable discontinuous finite element method. Therefore, in the following derivation of the hybridizable discontinuous finite element method, only the carrier drift-diffusion equation group (2) is derived. However, it must be noted that the same method can also be used to obtain the nonlinear Poisson equation after discretization based on the hybridizable discontinuous finite element method.

[0026] The conservative form of the carrier drift-diffusion equations is:

[0027]

[0028] in

[0029]

[0030]

[0031]

[0032] Introducing a mixed variable J to unify the above equations, we have

[0033]

[0034]

[0035] Where E is a coefficient introduced for simplicity of writing. Its physical meaning is the coefficient of the drift term of the carrier drift-diffusion equations. (1) is further transformed into two first-order partial differential equations, and we get

[0036]

[0037]

[0038]

[0039]

[0040] Where Ω is the solution area, Γ D is the ohmic contact boundary, Γ Nis the insulating interface, n is the unit external normal vector, g n and g p are the electron and hole concentrations at the ohmic contact boundary, respectively;

[0041] First, the entire solution area Ω is discretized into k sub-units, and a local sub-unit Ω is considered. h , using the estimated value of the numerical flux trace on the cell boundary To solve for the original variable q h and mixed variable J h ,Right now

[0042]

[0043]

[0044]

[0045]

[0046] and are the estimated values ​​of the numerical flux traces of electron and hole concentrations on the cell boundary respectively; the system of equations (8) is called a local problem because the equations are only applied to the local subcell Ω h Internal solution;

[0047] Next, conservative conditions are imposed at the global boundaries to ensure the continuity of variables and fluxes, i.e.

[0048]

[0049]

[0050]

[0051] The system of equations (9) is called a global problem because the equations are solved at the global boundaries Γ and Γ N Solve above, where Γ is the internal boundary, is the sum of the transitions of units i and j; is the numerical flux about the mixing variable, in the local subunit Ω h Inside, build and relationship, there is

[0052]

[0053] Where τ is the stability parameter defined in each subunit;

[0054] Next, we use the Galerkin method to obtain the weak form of the system of equations (8) and (9), for the variable J h ,qh , Introduce the discrete finite element space R h , Q h , To solve, the three discrete finite element spaces are the spaces to which the three variables belong respectively, in the local subunit Ω h Internal, local problem search (J h ,q h )∈R h ×Q h ,satisfy

[0055]

[0056]

[0057] Where R, W are variables J h ,q h The test function of , D is the coefficient of the diffusion term of the carrier drift-diffusion equations; for global problems, find satisfy

[0058]

[0059] Where U is a variable The test function;

[0060] In the above formulas, all subscripts h represent the corresponding variables in the local subunit Ω h The corresponding variables on ;

[0061] Based on the weak forms (11) and (12), the carrier drift-diffusion equations based on the HDG method are given in the local subunit Ω h The discrete matrix equation form is:

[0062]

[0063] The discretized form of the global problem is:

[0064]

[0065] The above A, B, CF, G, H, I are based on (11) and (12) about single variable q or J or The coefficient matrix, F q , F J , It is an incentive item.

[0066] Specifically, the fifth and sixth steps can be solved using the conventional Gummel iteration method.

[0067] Compared with the prior art, the present invention has the following significant advantages:

[0068] (1) High spatial discretization flexibility: The HDG method supports unstructured and non-conformal grids in spatial discretization. Due to the independence between units, each unit can use grids of different shapes and types. The meshing is more flexible than the traditional finite element method (FEM), and it is easier to fit complex structures while maintaining high efficiency and high accuracy.

[0069] (2) High flexibility in time discretization: The HDG method supports explicit and implicit time formats in time discretization, and can adopt local time steps according to the density of the grid in different regions.

[0070] (3) Support for arbitrary high-order basis functions; due to the high degree of independence between units, each unit can be processed using a different basis function, further reducing the requirements for mesh generation.

[0071] (4) Support for large-scale parallelization; the HDG method can be used for large-scale parallel computing problems of physical fields in multi-scale and complex semiconductor device structures. This is mainly reflected in the fact that multi-scale problems can be decomposed into multiple sub-problems, and different sub-problems can be calculated in parallel separately. BRIEF DESCRIPTION OF THE DRAWINGS

[0072] Figure 1 This is a flow chart for semiconductor device simulation based on high-order hybridized discontinuous Galerkin (HDG).

[0073] Figure 2 (a) is the geometric model of the transistor. Figure 2 (b) is the mesh used in the simulation.

[0074] Figure 3 is the physical field distribution of the transistor device, where Figure 3 (a) is the electric potential, Figure 3 (b) is the electron concentration n distribution, Figure 3 (c) is the hole concentration p distribution.

[0075] Figure 4 (a) is the 3D geometric model of the image sensor. Figure 4 (b) is a two-dimensional schematic diagram of the external voltage excitation interface.

[0076] Figure 5 is the electron concentration n distribution of the image sensor device. DETAILED DESCRIPTION

[0077] The present invention will be described in further detail below with reference to the accompanying drawings. The following examples will help those skilled in the art further understand the present invention but are not intended to limit the present invention in any way. It should be noted that those skilled in the art may make various variations and modifications without departing from the scope of the present invention. Such variations and modifications are all within the scope of protection of the present invention.

[0078] 1. The carrier drift-diffusion equation system in the semiconductor device structure can be described by second-order partial differential equations such as Poisson's equation and carrier drift-diffusion equation.

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

[0080]

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

[0082]

[0083]

[0084] Where ε is the dielectric constant of the material in the device, φ is the potential, q is the charge, n and p are the electron concentration and hole concentration respectively, N A and N D are the ionized acceptor concentration and the ionized donor concentration, R n and R p are the electron generation recombination rate and hole generation recombination rate, D n and D p are the electron diffusion coefficient and hole diffusion coefficient, μ n and μ p are electron mobility and hole mobility respectively, and t is time.

[0085] The carrier drift-diffusion equation system in the semiconductor device structure satisfies specific boundary conditions. The ohmic contact boundary condition can be written as:

[0086] φ=V app +φ eq (3a)

[0087]

[0088]

[0089] Where V app is the voltage applied at the boundary, φ eq is the built-in potential, N0 is the ionization doping concentration, n i is the intrinsic carrier concentration.

[0090] 2. In the proposed method, the Gummel iterative algorithm is used to decouple the carrier drift-diffusion equation and split it into two parts: the nonlinear Poisson equation and the linearized drift-diffusion equation group. First, the nonlinear Poisson equation is solved and the potential distribution is updated. After reaching the convergence condition, the result is substituted into the linearized drift-diffusion equation group to solve and update the electron concentration and hole concentration distribution, and then substituted into the nonlinear Poisson equation until convergence. The method flow chart is as follows: Figure 1 shown.

[0091] The method of numerically discretizing the nonlinear Poisson equation and the carrier drift-diffusion equation group using the HDG method is shown in steps 3-6.

[0092] 3. The HDG method is used to numerically discretize the nonlinear Poisson equation and the carrier drift-diffusion equations. The conservative form is:

[0093]

[0094] in

[0095]

[0096]

[0097]

[0098] Introducing a mixed variable J to unify the above equations, we have

[0099]

[0100]

[0101] Further transforming (4) into a system of first-order partial differential equations, we can obtain

[0102]

[0103]

[0104]

[0105] Where Ω is the solution area, Γ D is the ohmic contact boundary, Γ N is the insulating interface, n is the unit external normal vector, g n and g p are the electron and hole concentrations at the ohmic contact boundary, respectively.

[0106] 4. The HDG method solves the equation system (7) in two steps. First, in the local subunit Ω h Inside, use the estimated value of the numerical flux "trace" on the cell boundary To solve for the original variable q h and mixed variable J h ,Right now

[0107]

[0108]

[0109]

[0110]

[0111] The system of equations (8) is called a local problem; in the second step, conservative conditions are imposed at the global boundary to ensure the continuity of variables and fluxes, that is,

[0112]

[0113]

[0114]

[0115] The system of equations (9) is called the global problem, where Γ is the internal boundary, is the sum of the transitions of units i and j, for example It represents the sum of the n·r value in cell i and the n·r value in cell j. is the numerical flux about the mixing variable, in the local subunit Ω h Inside, there

[0116]

[0117] Where τ is a stability parameter defined in each subunit, which is a non-negative constant on each boundary of the unit; the subscript h indicates that the corresponding parameter is in the local subunit Ω h The corresponding parameters on the boundary.

[0118] 5. HDG method for variable J h ,q h , Discretize the finite element space using R h , Q h , To solve the system of equations (8) and (9), in the local subunit, the local problem is to find (J h ,q h )∈R h ×Q h ,satisfy

[0119]

[0120]

[0121] Where R, W are variables J h ,q h For global problems, find satisfy

[0122]

[0123] Where U is a variable The test function.

[0124] 6. Based on (11) and (12), the discretization form of the HDG method in the local subunit becomes:

[0125]

[0126] The discretized form of the global problem is:

[0127]

[0128] The above A, B, C, E, F, G, H are coefficient matrices based on (11) and (12), F q , F J , is the excitation term. The specific implementation process of the HDG method is: in each local subunit, assemble and solve (12), and use the estimated Represents original variables and mixed variables Then substitute into (13); obviously, (13) can be further expressed as only the unknown quantity The matrix equation system for :

[0129]

[0130] It is the global matrix equation of the HDG method, K is the global coefficient matrix, and F is the excitation term. It can be found that The unknowns are only on the unit boundary. Therefore, the solution scale of the HDG method, especially when high-order basis functions are used, is much smaller than the unknowns based on the unit cell in the traditional method.

[0131] 7. Through the discrete format of the Gummel iterative method in Specific Implementation 2 and the HDG method in Specific Implementation 6, the nonlinear Poisson equation and the carrier drift-diffusion equation group of the semiconductor device structure can be iteratively solved to obtain the electron concentration, hole concentration, potential and electric field distribution characteristics in the semiconductor device.

[0132] 8. Two numerical examples are given below:

[0133] (1) Semiconductor transistor model: establish the transistor geometry model to be solved and set the doping concentration distribution, such as Figure 2 (a) is shown. The geometric model is divided into uniform free tetrahedral meshes and the calculation is performed using the 4th order basis function, as shown in Figure 2 (b) As shown. The doping model of the transistor is set according to the device manual, and the carrier doping adopts Gaussian distribution. The emitter (E), base (B), and collector (C) of the transistor are set to ohmic contact boundary conditions. The physical field distribution in the transistor obtained during the calculation process is shown as follows Figure 3 As shown, Figure 3 (a) is the electron concentration n, Figure 3 (b) is the hole concentration p distribution, Figure 3 (c) is the electric potential distribution. The calculation results show that the use of high-order basis functions can greatly reduce the requirements for the number and quality of meshes while ensuring the correctness of the results. In particular, in areas where the field changes drastically, only a single mesh unit can capture the details of the relevant physical quantities.

[0134] (2) Image sensor model, according to Figure 4 (a) shows the three-dimensional geometric features, and establishes the geometric model of the image sensor to be solved. The geometric model is divided using a free tetrahedron mesh. The carrier doping model in the image sensor is set to a Gaussian distribution. The output electrode (Vout), selection electrode (Vsl), fixed electrode (Vdd), and set electrode (Vset) of the image sensor are set to ohmic contact boundary conditions as interfaces for external voltage input, as shown in Figure 1. Figure 4 (b) shown. Figure 5 The calculated three-dimensional distribution characteristics of electron concentration are consistent with the accuracy of the third-party commercial software COMSOL.

[0135] 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 semiconductor multi-physics effect simulation method based on hybrid discontinuous finite element, characterized in that: The steps include: Step 1: Mesh the semiconductor device and determine the physical governing equations of the semiconductor device structure, including the carrier drift-diffusion equation, Poisson's equation, and current continuity equation, which describe the electron concentration, hole concentration, potential, and electric field distribution in the device. Step 2: Determine the Ohmic contact boundary conditions and the initial electron concentration, hole concentration, and potential analytical distribution function corresponding to the physical governing equations of the semiconductor device structure; Step 3: The carrier drift-diffusion equation, Poisson equation, and current continuity equation are decoupled using the Gummel iteration method, splitting them into the nonlinear Poisson equation and the linearized carrier drift-diffusion equation group; Step 4: Using the hybridized discontinuous Galerkin HDG method, the nonlinear Poisson equation and the carrier drift-diffusion equations are numerically discretized respectively. Mixed variables are introduced to reduce the second-order partial differential equation to two first-order equations. Stability parameters are used to process the numerical flux terms between discontinuous units to obtain the discretized nonlinear Poisson equation and the carrier drift-diffusion equations. Specifically: The conservative form of the carrier drift-diffusion equations is: in φ is the electric potential, n and p are the electron concentration and hole concentration respectively, R n and R p are the electron generation recombination rate and hole generation recombination rate, D n and D p are the electron diffusion coefficient and hole diffusion coefficient, μ n and μ p are electron mobility and hole mobility respectively, and a mixed variable J is introduced to unify the above equations, Where E is the coefficient of the drift term of the carrier drift-diffusion equations; (4) is further transformed into two first-order partial differential equations, and we get Where Ω is the solution area, Γ D is the ohmic contact boundary, Γ N is the insulating interface, n is the unit external normal vector, g n and g p are the electron and hole concentrations at the ohmic contact boundary, respectively; First, the entire solution area Ω is discretized into k sub-units, and a local sub-unit Ω is considered. h , using the estimated value of the numerical flux trace on the cell boundary To solve for the original variable q h and mixed variable J h ,Right now and are the estimated values ​​of the numerical flux traces of electron and hole concentrations at the cell boundaries, respectively; The system of equations (8a) to (8d) is called a local problem. The equations are only valid in the local subunit Ω. h Internal solution; Next, conservative conditions are imposed at the global boundaries to ensure the continuity of variables and fluxes, i.e. The system of equations (9a) to (9c) is called the global problem. N Solve above, where Γ is the internal boundary; is the sum of the transitions of units i and j; is the numerical flux about the mixing variable, in the local subunit Ω h Inside, build and relationship, there is Where τ is the stability parameter defined in each subunit; Next, the Galerkin method is used to obtain the weak form of the equation systems (8a) to (8d) and (9a) to (9c), with respect to the variable J h ,q h , Introduce the discrete finite element space R h , Q h , To solve, in the local subunit Ω h Internal, local problem search (J h ,q h )∈R h ×Q h ,satisfy Where R, W are variables J h ,q h The test function of , D is the coefficient of the diffusion term of the carrier drift-diffusion equations; for global problems, find satisfy Where U is a variable The test function; In the above formulas, all subscripts h represent the corresponding variables in the local subunit Ω h The corresponding variables on ; Based on the weak forms (11a) to (11b) and (12), the carrier drift-diffusion equations based on the HDG method are given in the local subunit Ω h The discrete matrix equation form is: The discretized form of the global problem is: The above A, B, C, F, G, H, I are based on (11a) to (11b) and (12) about the single variable q or J or The coefficient matrix of is the incentive item; The nonlinear Poisson equation discretized based on the hybrid discontinuous finite element method is obtained using the same method; Step 5: Use the Gummel iteration method to solve the discretized nonlinear Poisson equation and update the potential distribution. Once convergence conditions are reached, substitute the result into the discretized carrier drift-diffusion equations to solve and update the electron and hole concentration distributions. Step 6: Repeat step 5 until the physical field information in the solved semiconductor device structure reaches the set convergence condition, and finally obtain the electron concentration, hole concentration, potential and electric field distribution characteristics in the semiconductor device structure.

2. The semiconductor multi-physics effect simulation method based on hybrid discontinuous finite element according to claim 1, characterized in that: In the first step, The Poisson equation is: The carrier drift-diffusion equations are: Where ε is the dielectric constant of the material in the device, q is the charge, N A and N D are the ionization acceptor concentration and the ionization donor concentration, respectively, and t is the time.

3. The semiconductor multi-physics effect simulation method based on hybrid discontinuous finite element according to claim 2, characterized in that: In the second step, the ohmic contact boundary conditions are: φ=V app +φ eq (3a) Where V app is the voltage applied at the boundary, φ eq is the built-in potential, N0 is the ionization doping concentration, n i is the intrinsic carrier concentration.

Images