CN106339562B - Simulation Method of Nonlocal Quantum Tunneling Between Energy Bands with Current Conservation Property - Google Patents

Abstract

The invention discloses a method for simulating non-local quantum tunneling between energy bands with the characteristic of current conservation, which includes: generating a linear equation system with the node value increment of the solved physical variable as a variable in the physical area of a discrete semiconductor device; Local quantum tunneling, determine the public energy interval to define the first side and the second side; weight distribution of point-to-point tunneling current density and energy integral volume generated by interpolation; coefficient matrix of linear equations is divided into ordinary differential equations caused by discretization The main matrix and the additional coefficient matrix caused by non-local quantum tunneling with current conservation characteristics are stored and stored; the coefficient matrix is solved by Gaussian elimination method. The invention guarantees the conservation of non-local quantum tunneling current and eliminates simulation errors and uncertainties caused by non-conservation; the method of combining the main matrix and the auxiliary correlation matrix can quickly solve the coefficients obtained by the linearization of nonlinear equations Linear equations to improve the efficiency of solving.

Description

Simulation Method of Nonlocal Quantum Tunneling Between Energy Bands with Current Conservation Property

technical field

The invention relates to a semiconductor device simulation method, in particular to a non-local quantum tunneling simulation method between energy bands with the characteristic of current conservation.

Background technique

Nonlocal quantum tunneling between different energy bands is a very common physical effect in compound semiconductor devices, which is common in many semiconductor devices, such as multi-junction solar cells and tunneling field effect transistors. How to accurately simulate this physical phenomenon is always concerned with semiconductor device physics and numerical mathematics. The usual method is to convert the point-to-point tunneling current density into a composite term of various continuity equations. The more commonly used is the Kane local model (E.O.Kane, "Zener tunneling in semiconductors", J.Phys.Chem.Solids, vol.12, pp.181-188, 1959), in this model, the energy band Non-local quantum tunneling is treated as a quantity that is only related to the local electric field strength, and the software that adopts this processing method is CROSSLIGHT. Hurkx et al. further added the defect-assisted tunneling enhancement effect, and expressed the effect in the form of carrier lifetime attenuation coefficient (G.A.M. Hurkx, D.B.M. Klaassen and M.P.G. Knuvers, "Anew recombination model for device simulation including tunneling", IEEE Trans. Electron Devices, vol.39, no.2, pp.331-338, 1992), the software using this model is wxAMPS. Commercial software such as Silvaco and Synopsys use the so-called nonlocal quantum tunneling model (Ieong MK, Solomon PM, Laux SE, Wong HSP, Chidambarrao D. Comparison of raised and Schottky source/drain MOSFETs using a novel tunneling contact model. International Electron Devices Meeting San Fransisco, California, USA, 1998, 733–736.), this model uses the method of carrier tunneling in quantum mechanics to deal with tunneling, and the carriers on one side have a certain probability of tunneling through the material band gap. The potential barrier to the other side, the size of the tunneling current also depends on the distribution of carriers on both sides, so that the current density caused by the non-local ship is distributed at different points in space. This model has been applied to multi-junction solar cells (M. Hermle, G. Letay, "Numerical simulation of tunnel diodes formulti-junction solar cells", Progress in Photovoltaics Research and Applications Vol.16, No.5, pp.409-418, 2008).

However, the above-mentioned non-local model still treats the point-to-point tunneling current as a compound item in the carrier continuity equation, and does not consider the conservation of the overall tunneling current on both sides, that is, the quantum tunneling current from one side The current density across to the other side should be equal. For example, in practice, it is found that the current density on both sides of the above-mentioned non-localized model will be twice as large, which will bring certain errors to the numerical analysis of the device structure.

Contents of the invention

The present invention provides a non-local quantum tunneling simulation method between energy bands with the characteristic of current conservation, which eliminates the simulation error and uncertainty caused by non-conservation, and improves the solution efficiency.

In order to achieve the above object, the present invention provides an inter-band non-local quantum tunneling simulation method with current conservation characteristics, which is characterized in that the method includes:

S1. In the physical area of discrete semiconductor devices, using partial differential equation numerical methods, such as finite volume method, finite difference method, and finite element method, etc., the generation of discrete semiconductor differential equations and boundary conditions takes the node value increment of the physical variable to be solved as a variable The linear equation system;

S2. Deal with non-local quantum tunneling, determine the public energy interval to define the first side and the second side; the first side determines the point-to-point tunneling current density of the corresponding discrete node non-local quantum tunneling in the public energy interval and the integral energy unit; the carrier quasi-Fermi level on the second side is obtained by the interpolation method corresponding to the current dispersion;

S3. Weight distribution of point-to-point tunneling current density and energy integral volume generated by interpolation;

S4. The coefficient matrix of the linear equation system is divided into the main matrix caused by the discretization of the usual differential equations and the additional coefficient matrix caused by the non-local quantum tunneling with current conservation characteristics and stored;

S5. Using the Gaussian elimination method to solve the nonlinear equation system.

The above S1 contains:

In the physical area of discrete semiconductor devices, using numerical methods for partial differential equations, such as finite volume method, finite difference method, and finite element method, discrete semiconductor differential equations and boundary conditions generate linear equations with the node value increment of the physical variable to be solved as the variable. equation set;

In the above S2, the numerical model of nonlocal quantum tunneling for node I, electrons with energy integral volume dE, and the current density generated by nonlocal quantum tunneling at grid point i is as follows:

where N is the band-edge density of states, _νth is the carrier thermal velocity, T(E) is the tunneling probability of a carrier with energy E, is the carrier quasi-Fermi potential with energy E _i on both sides, and exists with the quasi-Fermi level Relationship.

In the above S2, when dealing with non-localized quantum tunneling, first define the first side and the second side, and determine the common energy range in which non-localized quantum tunneling can occur.

In the above S2, the first side determines the point-to-point current density and integral energy unit of the corresponding discrete node nonlocal quantum tunneling in the common energy interval.

In the above S2, the carrier quasi-Fermi level of the second side is obtained by the interpolation method corresponding to the current density dispersion, and then the correction of the discrete node coefficient and its derivative part are added to the relevant point value corresponding to the coefficient matrix superior;

Among them, the current density obtained by holes is discrete as formula (2):

In formula (2), 0 and 1 mark the nodes of the first side and the second side on a grid line segment in the discrete direction, μ is the hole mobility, N is the hole band-edge state density, φ _p , E _v and V are quasi-Fermi potential, valence band edge energy and electrostatic potential of hole respectively, h is the length of line segment, χ _p is the ratio of Fermi-Dirac statistical distribution to Bose-Einstein statistical distribution;

The relationship between the quasi-Fermi potential and the tunneling energy obtained by transforming formula (2) at a certain point in a grid cell is as follows:

In the above S3, the second side performs weight distribution according to the overlapping degree of point-to-point tunneling current density and energy integral volume of the first side and the energy integral volume of the second side node, and adds the allocated energy integral volume of the first side and its derivative to to the associated node value on the second side of the coefficient matrix.

The above main matrix has a tridiagonal form and adopts a column storage form with a row bandwidth of 2*kl+ku+1.

The additional coefficient matrix storage described above contains:

define a compound variable consisting of an integer variable and a real variable to store the position and value of the additional incidence matrix element;

Define a variable-length row vector composed of the composite data to store the additional variable of each row except the main matrix element generated by the additional incidence matrix;

For the energy integral volume in the increasing direction of the encoding of the fixed side i node, if there are multiple energy integral volumes on the other side, the second side node k and the equal energy interpolation node pair There is a relationship, if then in and Additional matrix elements are introduced below the diagonal elements two columns down, with row numbers starting from to k; if node k and or Equal, the main matrix covers the additional incidence matrix, and the additional incidence matrix only stores elements that are not stored in the main matrix. For the upper double-column association form, the additional incidence matrix only stores the additional two-column matrix element values and their column indices;

For the energy integral volume in the coding reduction direction of the fixed side i node, if there are multiple energy integral volumes on the other side, the first side node k and the equal energy interpolation node pair There is a relationship, if then in and Introduces additional matrix elements above the diagonal matrix elements on both columns, row numbers from k to If node k and or Equal, the main matrix covers the additional incident matrix, and the additional incident matrix only stores the elements that are not stored in the main matrix. For the upper double-column association form, the additional incident matrix not only stores the values of the additional two columns of matrix elements and their column indices, but also stores the elements from to the matrix element position of k and initialized to 0.

The above S5 contains:

Suppose the current column index is i;

Determine whether there is a matrix element that is larger than the absolute value of the diagonal matrix element below the i-th column diagonal matrix element;

If so, this row is exchanged with the i-th row, and then the i-th column of Gaussian vector is calculated, which is a column vector formed by dividing the matrix elements below the diagonal matrix element by the diagonal matrix elements after the column selection pivot, as in the formula ( 4):

The elimination process is to subtract each matrix element from the product of the Gaussian vector element corresponding to the row where the matrix element is located and the column element of the row corresponding to the Gaussian vector, as shown in formula (5):

Compared with the prior art, the non-local quantum tunneling simulation method between energy bands of the present invention with current conservation characteristics has the advantage that the method proposed by the present invention can ensure the conservation of non-local quantum tunneling current and satisfy the quantum Mechanics requires that carrier tunneling is conserved in the entire space, thus eliminating the simulation errors and uncertainties caused by non-conservation;

The method proposed by the invention can quickly solve the coefficient linear equation group obtained by the linearization of the nonlinear equation group by using the method of combining the main matrix and the auxiliary correlation matrix, and improves the solution efficiency.

Description of drawings

Fig. 1 is the flow chart of the non-local quantum tunneling simulation method between energy bands with current conservation characteristics of the present invention;

Fig. 2 is the flowchart of semiconductor device simulation;

Fig. 3 is a discrete geometric relationship diagram of node variables;

Fig. 4 is a schematic diagram of nonlocal quantum tunneling;

Figure 5 is a schematic diagram of the corresponding relationship between the distribution of nodes on both sides of the tunnel in terms of energy;

Fig. 6 is a schematic diagram of the coding increase direction of the fixed side i node;

Fig. 7 is a schematic diagram of the coding reduction direction of the fixed edge i node;

Fig. 8 is the structural representation of embodiment two;

Fig. 9 is a schematic structural diagram of the third embodiment.

Detailed ways

Specific embodiments of the present invention will be further described below in conjunction with the accompanying drawings.

As shown in Figure 1, the present invention discloses Embodiment 1 of an inter-band nonlocal quantum tunneling simulation method with current conservation characteristics, and the method specifically includes the following steps:

S1. In the physical area of discrete semiconductor devices, using partial differential equation numerical methods, such as finite volume method, finite difference method, and finite element method, etc., the generation of discrete semiconductor differential equations and boundary conditions takes the node value increment of the physical variable to be solved as a variable system of linear equations.

As shown in FIG. 2 , the general process of semiconductor device simulation includes: firstly, mesh generation is performed according to the geometric shape of the device structure. Secondly, based on the partial differential equation discretization method, the Poisson equation, electron and hole continuity equations are discretized into nonlinear equations with node variables as parameters on the grid points. The discretization methods usually include finite difference, finite volume and finite element etc., the compound term in the continuity equation is usually assumed to be constant on a small discrete interval.

As shown in Fig. 3, the finite volume method is taken as an example here, and of course it is not limited to the finite volume method. Assuming that the function value on node i is fi, the space integral volume belonging to node i is defined as the interval from the midpoint of the left interval to the midpoint of the right interval, and the geometric volume of this interval is xVol=0.5*(h1+h2). At the same time, the energy integral volume belonging to node i is defined as the difference between the midpoint energy band value of the left interval and the midpoint energy band value of the right interval, and its value is EVol=0.5*(Ei+1-Ei-1). If it is assumed that f in this interval is approximated by fi, then the contribution to the whole in this interval is xVol*fi. Nonlinear equations with node variables as parameters are usually solved iteratively. The basic process is to expand the linear equations with node variable increments as coefficients through appropriate series. The solution of linear equations usually has direct elimination method And iterative solution, linear equations usually have some special structures, such as tridiagonal or block tridiagonal, can use some fast solving algorithms such as the Gaussian elimination method of the selected pivot. Continue to modify the node variable value according to the obtained variable increment until the node variable increment or the function value of the nonlinear equation system meets certain requirements. After that, the result output and subsequent processing are carried out.

S2. Deal with non-local quantum tunneling, determine the public energy interval to define the first side and the second side; the first side determines the point-to-point tunneling current density of the corresponding discrete node non-local quantum tunneling in the public energy interval and the integral energy unit; the carrier quasi-Fermi level on the second side is obtained by the interpolation method corresponding to the current dispersion.

S2.1. To deal with non-local quantum tunneling, first define the first side and the second side, and determine the common energy interval where non-local quantum tunneling can occur.

Numerical model of nonlocal quantum tunneling For node I, electrons with energy integral volume dE, the current density generated by nonlocal quantum tunneling at grid point i is as follows:

where N is the band-edge density of states, _νth is the carrier thermal velocity, T(E) is the tunneling probability of a carrier with energy E, is the carrier quasi-Fermi potential with energy E _i on both sides, and exists with the quasi-Fermi level Relationship. Since the current density is an integral function of carrier energy, according to our description above, it can also be converted into the form of a space integral function.

As shown in Figure 4, it is a schematic diagram of nonlocal quantum tunneling, and the lines at the bottom indicate the distribution of grid points on both sides of the material. The basic physical process of nonlocal quantum tunneling also marks the distribution of grid points on both sides of the emission quantum tunneling. It can be seen from Figure 3 that if the discrete nodes on the first side (left) are used as fixed points, due to The energies of the discrete nodes on both sides may not coincide, and this is impossible in actual operation, because the carrier energy is also one of the variables that need to be solved, but it is certain that for each node on a certain side where tunneling occurs i (corresponding to the tunneling energy Eqt), the energy of node i falls on the other side to and In the grid cell of the end point, there are equal energy interpolation node pairs on the other side Therefore, interpolation is required to obtain the quasi-Fermi potential of carriers at the equal energy point on the second side (right).

S2.2. The first side determines the point-to-point current density and integral energy unit of the corresponding discrete node nonlocal quantum tunneling in the public energy interval.

S2.3. The carrier quasi-Fermi energy level on the second side is obtained by the interpolation method corresponding to the current density dispersion, and then the correction of the discrete node coefficient and its derivative part are added to the relevant point value corresponding to the coefficient matrix superior.

The invention proposes an interpolation method for obtaining quasi-Fermi energy levels of carriers at equal energy points according to the current density discrete method, which is obtained by transforming the discrete format of the carrier current density equation. Usually the current density is assumed to be fixed or linearly changed between two nodes, and the band edge value is also considered to be linearly changed. This is the Sharfetter-Gummel discretization mechanism. The current density obtained by holes is discretized as follows: 2):

In formula (2), 0 and 1 mark the nodes of the first side and the second side on a grid line segment in the discrete direction, μ is the hole mobility, N is the hole band-edge state density, φ _p , E _v and V are the quasi-Fermi potential, valence band edge energy and electrostatic potential of holes, respectively, h is the length of the line segment, and χ _p is the ratio of the Fermi-Dirac statistical distribution to the Bose-Einstein statistical distribution.

If the current density adopts the Sharfetter-Gummel discretization mechanism, the change relationship of the quasi-Fermi potential of the carrier is actually established. In a certain grid cell obtained by transforming the formula (2) (grid points 0 and 1 on both sides) The relationship between the quasi-Fermi potential at which quantum tunneling occurs at a certain point and the tunneling energy is shown in formula (3):

It can be seen from formula (3) that the carrier quasi-Fermi potential is neither constant nor linearly changing, but presents a distribution related to the difference of the grid endpoint values.

S3. The point-to-point tunneling current density and energy integral volume generated by the interpolation are weighted. The second side performs weight distribution according to the overlapping degree of point-to-point tunneling current density and energy integral volume of the first side and the energy integral volume of the second side node, and adds the distributed first side energy integral volume and its derivative to the coefficient matrix On the value of the relevant node on both sides.

The invention proposes a weight distribution method of point-to-point tunneling current and energy integral volume with the characteristic of current conservation. In order to ensure that the tunneling current densities on both sides are equal, Σ _i J[E _i ]EVol _i =Σ _k J[E _k ]EVol _k is required. To satisfy the above formula, only point-to-point tunneling current density and energy integral volume on one side can be fixed. According to The overlapping degree of the energy integral volume of the nodes on both sides is carried out on the other side for weight distribution.

As shown in Figure 5, if the tunneling current density on the left is fixed, there is a situation where the energy integral volumes of multiple nodes on the right fall within the energy integral volume of a node on the left, that is, J[E _i ]EVol _i and Multiple nodes EVol _k on the right overlap, so the weight of the tunneling current density of the left i node assigned to each k node on the right is J[E _i ]EVol _i ∩EVol _k . According to the previous description, the point-to-point tunneling current density J[E _i ] of the node i on the left contains the interpolation point on the right The carrier quasi-Fermi potential of , such that the right node k and the node There is a connection.

S4. The coefficient matrix of the linear equation system is divided into the main matrix caused by the discretization of the general differential equation and the additional coefficient matrix caused by the non-local quantum tunneling with the characteristic of current conservation, and stored.

The storage of the coefficient matrix is divided into two parts, one part is caused by the discretization of the usual differential equations, called the main matrix, usually in the form of a tridiagonal, and the other part is caused by the nonlocal quantum tunneling with the characteristic of current conservation, called is an additional coefficient matrix, which usually has a row form or a column form, but does not have both forms at the same time. In order to improve the speed, the solution of the coefficient matrix adopts the modified Gaussian elimination method with column selection.

Main matrix storage: Since the nonlocal quantum tunneling between energy bands generally only occurs in a small device area, most of the coefficient matrices of linear equations are just the main matrix caused by the discretization of the carrier continuity differential equation, which has three pairs In the form of corners, the column storage form with a row bandwidth of 2*kl+ku+1 is adopted, which is convenient for the Gaussian elimination method of column selection pivot.

Additional coefficient matrix storage:

Defines a compound variable consisting of an integer variable and a real variable to store the location and value of the additional incidence matrix element.

Define a row vector of variable length consisting of the composite data to store additional variables for each row in addition to the main matrix elements produced by the additional incidence matrix.

For the correlation matrix that exists in the upper triangle and the lower triangle, a relatively different storage method is adopted according to the basic principle of the Gaussian elimination method of selecting the pivot.

As shown in Figure 6, for the energy integral volume in the increasing direction of the encoding of the node i on the fixed side, if there are multiple energy integral volumes on the other side, the second side node k and the equal energy interpolation node pair There is a relationship, if then in and Additional matrix elements are introduced below the diagonal elements two columns down, with row numbers starting from to k; if node k and or Equal, the main matrix covers the additional incidence matrix, and the additional incidence matrix only stores elements that are not stored in the main matrix. For the upper double-column association form, the additional incidence matrix only stores the additional two-column matrix element values and their column indices;

As shown in Figure 7, for the energy integral volume in the lower direction of the encoding of the node i on the fixed side, if there are multiple energy integral volumes on the other side, the first side node k and the equal energy interpolation node pair There is a relationship, if then in and Introduces additional matrix elements above the diagonal matrix elements on both columns, row numbers from k to If node k and or Equal, the main matrix covers the additional incident matrix, and the additional incident matrix only stores the elements that are not stored in the main matrix. For the upper double-column association form, the additional incident matrix not only stores the values of the additional two columns of matrix elements and their column indices, but also stores the elements from to the matrix element position of k and initialized to 0.

Considering the non-overlapping nature of the energy integral volumes of the left nodes, the elements of the additional incidence matrix can only have the row form or the column form, but not both forms at the same time.

S5. Using the Gaussian elimination method to solve the system of linear equations.

The invention discloses a Gaussian element elimination method for solving a coefficient matrix. The method combines a main matrix and an additional matrix, adopts a method of column selection of principal elements, quickly solves a coefficient equation, and saves the calculation times of directly using an iterative algorithm. The basic operation flow of the Gaussian elimination method based on the selected pivot:

Suppose the current column index is i.

Determine whether there is a matrix element that is larger than the absolute value of the diagonal matrix element below the i-th column diagonal matrix element;

If so, this row is exchanged with the i-th row, and then the i-th column of Gaussian vector is calculated, which is a column vector formed by dividing the matrix elements below the diagonal matrix element by the diagonal matrix elements after the column selection pivot, as in the formula ( 4):

The elimination process is to subtract each matrix element from the product of the Gaussian vector element corresponding to the row where the matrix element is located and the column element of the row corresponding to the Gaussian vector, as shown in formula (5):

According to above-mentioned idea, the specific way of this embodiment is as follows:

a), for the upper double-column association form, no matter the column selects the pivot or the process of eliminating elements, it will not produce additional additional matrix elements, and the method of the present invention directly exchanges and eliminates elements for this situation;

b), for the lower double-column association form, no matter the column selects the pivot or the process of eliminating elements, additional additional matrix elements will be generated, and the method of the present invention first generates the row index for this situation. to k+kl column index to k and each matrix element is a sub-matrix of the composite data declared in S4, directly exchange and eliminate elements in this sub-matrix;

The invention obtains the conservation of the left and right quantum tunneling currents by fixing the point-to-point tunneling current on one side and the energy integration unit, and satisfies the requirements of quantum mechanical tunneling. This method can be applied to the numerical analysis of devices containing tunnel junctions, such as multi-junction high-efficiency solar cells, electronic power devices, and the like.

As shown in Figure 8, it is the second embodiment of the present invention, taking the GaInP/GaAs double-junction solar cell as an example, which includes a GaAs buffer layer 81, an AlGaAs back field 82, a GaAs active layer 83, an AlInP window layer 84, a GaInP n++ Doped layer 85 , AlGaAs p++ doped layer 86 , AlGaInP back field 87 , GaInP active layer 88 , AlGaInP window layer 89 , and GaAs cap layer 810 . The structure is grown on an n-type GaAs substrate using low-pressure metal-organic chemical vapor deposition equipment, and interband non-localized quantum tunneling occurs between the GaInP n++ doped layer 85 and the AlGaAs p++ doped layer 86. Practice has proved that the tunneling characteristics It directly restricts the performance of the entire device, and accurate simulation of tunneling characteristics can significantly guide the design and development of actual devices.

As shown in FIG. 9, it is the third embodiment of the present invention, taking the tunneling field effect transistor (TFET) as an example, which includes a buffer layer 91, an n-type doped layer 92, an n++ doped layer 93, and a p++ doped layer 94 . The structure is prepared by ion implantation or epitaxial growth, and interband non-localized quantum tunneling occurs between the p++ doped layer 94 and the n++ doped layer 93 . Accurate simulation of tunneling characteristics can play a guiding role in the design and development of actual device structures

Although the content of the present invention has been described in detail through the above preferred embodiments, it should be understood that the above description should not be considered as limiting the present invention. Various modifications and alterations to the present invention will become apparent to those skilled in the art upon reading the above disclosure. Therefore, the protection scope of the present invention should be defined by the appended claims.

Claims (10)

1. A method for simulating nonlocal quantum tunneling between energy bands with current conservation characteristics, characterized in that the method comprises: S1. In the physical area of discrete semiconductor devices, using partial differential equation numerical methods, such as finite volume method, finite difference method, and finite element method, etc., the generation of discrete semiconductor differential equations and boundary conditions takes the node value increment of the physical variable to be solved as a variable The linear equation system; S2. Deal with non-local quantum tunneling, determine the public energy interval to define the first side and the second side; the first side determines the point-to-point tunneling current density of the corresponding discrete node non-local quantum tunneling in the public energy interval and the integral energy unit; the carrier quasi-Fermi level on the second side is obtained by the interpolation method corresponding to the current dispersion; S3. Weight distribution of point-to-point tunneling current density and energy integral volume generated by interpolation; S4. The coefficient matrix of the linear equation system is divided into the main matrix caused by the discretization of the usual differential equations and the additional coefficient matrix caused by the non-local quantum tunneling with current conservation characteristics and stored; S5. Using the Gaussian elimination method to solve the system of linear equations. 2. the nonlocal quantum tunneling simulation method between energy bands with electric current conservation characteristic as claimed in claim 1, is characterized in that, described S1 comprises: Discrete the physical area of semiconductor devices, and generate grids according to the geometry of semiconductor devices; Discretize semiconductor differential equations and boundary conditions according to the usual finite volume method, finite difference method and finite element method; Simultaneously generate a system of linear equations with the increment of the nodal value of the physical variable being solved for as a variable. 3. the non-local quantum tunneling simulation method between energy bands with electric current conservation characteristic as claimed in claim 1, is characterized in that, in described S2, the numerical model of non-local quantum tunneling is for node I, energy The current density generated by electrons with integral volume dE and nonlocal quantum tunneling at grid point i is as follows: where N is the band-edge density of states, _νth is the carrier thermal velocity, T(E) is the tunneling probability of a carrier with energy E, is the carrier quasi-Fermi potential with energy E _i on both sides, and exists with the quasi-Fermi level Relationship. 4. the method for simulating nonlocal quantum tunneling between energy bands with current conservation characteristics as claimed in claim 1, characterized in that, in said S2, when dealing with nonlocal quantum tunneling, first define the first side and The second side is to determine the public energy interval in which the non-local quantum tunneling can occur. 5. the method for simulating non-local quantum tunneling between energy bands with current conservation characteristics as claimed in claim 4, characterized in that, in said S2, the first side determines the corresponding discrete nodes in the public energy interval Point-to-point current density and integrated energy units for nonlocal quantum tunneling. 6. the non-local quantum tunneling simulation method between energy bands with current conservation characteristic as claimed in claim 4, is characterized in that, in described S2, the carrier quasi-Fermi level of the second side adopts current density The interpolation method corresponding to the discrete is obtained, and then the correction to the discrete node coefficient and its derivative part are added to the relevant point value corresponding to the coefficient matrix; Among them, the current density obtained by holes is discrete as formula (2): In formula (2), 0 and 1 mark the nodes of the first side and the second side on a grid line segment in the discrete direction, μ is the hole mobility, N is the hole band-edge state density, φ _p , E _v and V are quasi-Fermi potential, valence band edge energy and electrostatic potential of hole respectively, h is the length of line segment, χ _p is the ratio of Fermi-Dirac statistical distribution to Bose-Einstein statistical distribution; The relationship between the quasi-Fermi potential and the tunneling energy obtained by transforming formula (2) at a certain point in a grid cell is as follows: 7. The method for simulating non-local quantum tunneling between energy bands with current conservation characteristics as claimed in claim 4, wherein in said S3, the second side is based on the point-to-point tunneling current density and energy integral of the first side The overlapping degree of the volume and the energy integral volume of the second edge node is weighted, and the assigned energy integral volume of the first edge and its derivative are added to the relevant node value of the second edge of the coefficient matrix. 8. the method for simulating non-local quantum tunneling between energy bands with current conservation characteristics as claimed in claim 1, wherein the main matrix has a tridiagonal form, and the row bandwidth is 2*kl+ku +1 for the column store form. 9. the nonlocal quantum tunneling simulation method between energy bands with electric current conservation characteristic as claimed in claim 4, is characterized in that, described additional coefficient matrix storage comprises: define a compound variable consisting of an integer variable and a real variable to store the position and value of the additional incidence matrix element; Define a variable-length row vector composed of the composite data to store the additional variable of each row except the main matrix element generated by the additional incidence matrix; For the energy integral volume in the increasing direction of the encoding of the fixed side i node, if there are multiple energy integral volumes on the other side, the second side node k and the equal energy interpolation node pair There is a relationship, if then in and Additional matrix elements are introduced below the diagonal elements two columns down, with row numbers starting from to k; if node k and or Equal, the main matrix covers the additional incidence matrix, and the additional incidence matrix only stores elements that are not stored in the main matrix. For the upper double-column association form, the additional incidence matrix only stores the additional two-column matrix element values and their column indices; For the energy integral volume in the coding reduction direction of the fixed side i node, if there are multiple energy integral volumes on the other side, the first side node k and the equal energy interpolation node pair There is a relationship, if then in and Introduces additional matrix elements above the diagonal matrix elements on both columns, row numbers from k to If node k and or Equal, the main matrix covers the additional incident matrix, and the additional incident matrix only stores the elements that are not stored in the main matrix. For the upper double-column association form, the additional incident matrix not only stores the values of the additional two columns of matrix elements and their column indices, but also stores the elements from to the matrix element position of k and initialized to 0. 10. the non-localized quantum tunneling simulation method between energy bands with electric current conservation characteristic as claimed in claim 9, is characterized in that, described S5 comprises: Suppose the current column index is i; Determine whether there is a matrix element that is larger than the absolute value of the diagonal matrix element below the i-th column diagonal matrix element; If so, exchange the row where the matrix element whose absolute value is larger than the diagonal matrix element is located with the i-th row, and then calculate the i-th column Gaussian vector, which is the pair after dividing the matrix element below the diagonal matrix element by the column selection pivot element A column vector composed of corner matrix elements, such as formula (4): The elimination process is to subtract each matrix element from the product of the Gaussian vector element corresponding to the row where the matrix element is located and the column element of the row corresponding to the Gaussian vector, as shown in formula (5):