Method and system for simulating and optimizing performance of semiconductor multilayer heterostructure

Abstract

The application relates to the technical field of semiconductor device simulation, in particular to a semiconductor multilayer heterostructure performance simulation optimization method. The method obtains initial energy bands and carrier distribution by self-consistently solving Poisson and Schrodinger equations, reversely regulates the thickness and composition of each layer according to the target energy band gradient, and realizes energy band gradient optimization. Then, the lattice mismatch strain and dislocation density are calculated, and a strain buffer layer is inserted at the high defect interface to suppress defects. The application realizes the collaborative optimization of the energy band gradient and the lattice matching, and improves the heterostructure performance and the preparation feasibility.

Description

Technical Field

[0001] This invention relates to the field of semiconductor device simulation technology, and in particular to a method and system for performance simulation and optimization of semiconductor multilayer heterostructures. Background Technology

[0002] In the design and fabrication of semiconductor multilayer heterostructures, performance simulation optimization is a crucial step guiding material growth and device development. Existing techniques typically rely on pre-setting structural geometry and material composition parameters, using physical models such as the Poisson and Schrödinger equations to simulate and analyze the electrical properties of the heterojunction, including its band structure and carrier distribution. This forward design process aims to predict device performance under given structural parameters, providing theoretical reference for experiments.

[0003] Conventional simulation optimization methods focus on achieving a single physical objective, such as a specific bandgap or carrier concentration. Designers often determine the thickness and composition of each layer based on experience or simplified models, and then verify whether the electrical properties meet expectations through numerical calculations. If the results are unsatisfactory, parameters need to be manually adjusted and recalculated repeatedly. This process relies on the operator's expertise and trial and error, which is inefficient and makes it difficult to systematically approximate the global optimum. Especially for complex multilayer structures with huge parameter spaces, manual iteration is difficult to coordinate multiple interrelated or even conflicting performance indicators.

[0004] A more significant drawback is that existing conventional methods often separate the optimization of electrical performance from the mechanical reliability of the structure. Electrical simulations are typically based on the assumption of ideal lattice matching, or only roughly estimate the strain caused by lattice mismatch, failing to incorporate the generation and evolution of defects such as dislocations and cracks caused by strain into the core optimization loop. In actual material growth, excessive lattice mismatch stress can induce high-density dislocations. These defects can not only penetrate the active region, severely degrading carrier mobility and device lifetime, but may even lead to thin film cracking. The conventional approach is to perform mechanical reliability verification separately after the electrical design is completed, or to remedy the situation by introducing a buffer layer afterward. This sequential design mode of "performance first, reliability second" easily leads to electrically optimized structures that are mechanically infeasible, resulting in design iterations and wasted resources.

[0005] Therefore, the industry urgently needs a method that can collaboratively optimize electrical performance and mechanical reliability during the intrinsic design phase to overcome the shortcomings of traditional serial design processes, which are inefficient and prone to reliability issues. Summary of the Invention

[0006] The present invention provides a method and system for performance simulation and optimization of semiconductor multilayer heterostructures, which can solve the problems in the prior art.

[0007] A first aspect of the present invention provides a method for performance simulation optimization of semiconductor multilayer heterostructures, comprising: The initial geometric parameters, material composition distribution parameters, and target band gradient distribution of the semiconductor multilayer heterostructure are obtained. Based on the initial geometric parameters and material composition distribution parameters, the band discontinuity values ​​and carrier concentration distribution at each heterojunction interface are calculated using the Poisson equation and the Schrödinger equation in a self-consistent manner. Based on the deviation between the band discontinuity value and the target band gradient distribution, the thickness parameters and material composition parameters of each heterojunction are adjusted in reverse to ensure that the conduction band order and valence band order at the interface of each heterojunction satisfy the preset gradient continuity constraint condition, thereby obtaining the optimized geometric parameters and material composition parameters. Based on the optimized geometric parameters and material composition parameters, the lattice mismatch strain distribution and dislocation line density distribution at each heterojunction interface are calculated. When there is an interface location in the dislocation line density distribution that exceeds a preset density threshold, a strain buffer layer structure is inserted at the corresponding interface. By adjusting the thickness and composition gradient distribution of the strain buffer layer structure, the dislocation line density distribution is reduced to below the preset density threshold, and the final structural parameters after defect suppression are obtained. The final structural parameters after defect suppression are output as the fabrication process parameters of the semiconductor multilayer heterostructure.

[0008] Based on the initial geometric parameters and material composition distribution parameters, the band discontinuities and carrier concentration distributions at each heterojunction interface are calculated using the Poisson and Schrödinger equations in a self-consistent manner, including: The spatial coordinate range of each heterojunction layer is determined based on the initial geometric parameters. The dielectric constant distribution and effective mass distribution of each layer are determined based on the material composition distribution parameters. The spatial region of each heterojunction layer is discretized into a grid. A difference scheme of the Poisson equation is established on the discrete grid nodes. The Poisson equation describes the relationship between electrostatic potential and carrier charge density. An eigenvalue problem of the Schrödinger equation is established on the discrete grid nodes. The potential energy term of the Schrödinger equation is composed of the electrostatic potential and the intrinsic band structure of the material. Solving the eigenvalue problem yields the wave function distribution and intrinsic energy of each quantum state. Based on the wave function distribution and intrinsic energy, the electron concentration and hole concentration at each grid node are calculated using the Fermi-Dirac distribution function. The calculated carrier concentration is substituted into the Poisson equation to update the electrostatic potential distribution. The Schrödinger equation and the Poisson equation are solved iteratively until the maximum difference between the electrostatic potential distributions obtained from two adjacent iterations is less than the preset convergence criterion, thus obtaining the self-consistent electrostatic potential distribution and carrier concentration distribution. Based on the electrostatic potential distribution obtained after self-consistent solution, the energy values ​​of the conduction band edge and the valence band edge on both sides of each heterojunction interface are extracted. The difference between the energy values ​​of the conduction band edge is calculated to obtain the conduction band discontinuity value, and the difference between the energy values ​​of the valence band edge is calculated to obtain the valence band discontinuity value.

[0009] An eigenvalue problem of the Schrödinger equation is established on the discrete grid nodes. The potential energy term of the Schrödinger equation is composed of the electrostatic potential and the intrinsic band structure of the material. Solving the eigenvalue problem yields the wave function distribution and eigenenergy of each quantum state, including: The numerical distribution of the electrostatic potential is extracted on the discrete grid nodes, and the conduction band bottom energy and valence band top energy of the material corresponding to each grid node are read. The conduction band bottom energy and valence band top energy are used as characterization parameters of the intrinsic band structure of the material. The effective potential energy distribution of the electronic state is constructed by superimposing the numerical distribution of the electrostatic potential with the bottom energy of the conduction band, and the effective potential energy distribution of the hole state is constructed by superimposing the numerical distribution of the electrostatic potential with the top energy of the valence band. Based on the spatial coordinates of the discrete grid nodes, a discretized expression for the kinetic energy operator is constructed. The discretized expression for the kinetic energy operator uses a finite difference scheme to discretize the second spatial derivative of the wave function. Substituting the discretized expression of the kinetic energy operator and the effective potential energy distribution into the Schrödinger equation, an eigenvalue problem in matrix form is formed on the discrete grid nodes. The coefficient matrix of the eigenvalue problem is composed of kinetic energy terms and potential energy terms. The coefficient matrix of the eigenvalue problem is numerically solved to obtain the eigenvalues ​​and eigenvectors of each quantum state. The quantum states are sorted by energy level according to the eigenenergy, and the eigenenergy and wavefunction distribution of the electronic state near the bottom of the conduction band and the eigenenergy and wavefunction distribution of the hole state near the top of the valence band are extracted. The wave function distribution is normalized in space, and the spatial distribution of wave function probability density for each quantum state is calculated. The spatial distribution of wave function probability density is used to characterize the spatial localization characteristics of charge carriers in the heterojunction structure.

[0010] Based on the deviation between the band discontinuity value and the target band gradient distribution, the thickness parameters and material composition parameters of each heterojunction layer are adjusted in reverse to ensure that the conduction band order and valence band order at each heterojunction interface satisfy the preset gradient continuity constraint condition, thereby obtaining the optimized geometric parameters and material composition parameters, including: Calculate the actual band gradient value at each heterojunction interface. The actual band gradient value is the ratio of the difference in band discontinuity between two adjacent interfaces to the interface spacing. Extract the target band gradient value corresponding to each interface position from the target band gradient distribution. Calculate the deviation vector between the actual band gradient value and the target band gradient value. A sensitivity matrix is ​​established between band gradient deviation and structural parameter adjustment. The elements of the sensitivity matrix represent the influence coefficients of unit thickness change or unit composition change on the band gradient of each interface. Based on the sensitivity matrix and the deviation vector, the structural parameter adjustment that minimizes the norm of the deviation vector is solved by the least squares method. The structural parameter adjustment includes the thickness adjustment and material composition adjustment of each layer. The structural parameter adjustment is superimposed on the current thickness parameter and material composition parameter to obtain the updated thickness parameter and material composition parameter. Based on the updated thickness parameter and material composition parameter, the band discontinuity value at each heterojunction interface is recalculated. When the conduction band order and valence band order satisfy the preset gradient continuity constraint condition, the current thickness parameter and material composition parameter are output as the optimized geometric parameter and material composition parameter.

[0011] Establish a sensitivity matrix between band gradient bias and structural parameter adjustments, including: A set of structural parameters that need to be controlled is determined, the set of structural parameters including the thickness parameters of each heterojunction layer and the material composition parameters of each layer. A preset perturbation increment is applied to each structural parameter in the set of structural parameters, the perturbation increment being a preset percentage of the current value of the structural parameter. After applying the aforementioned perturbation increment, the band discontinuity values ​​at each heterojunction interface are recalculated using the self-consistent solution method of the Poisson equation and the Schrödinger equation, thus obtaining the distribution of the perturbed band discontinuity values. The actual band gradient value after the disturbance is calculated based on the distribution of the band discontinuities after the disturbance. The difference between the actual band gradient value after the disturbance and the actual band gradient value before the disturbance is obtained to obtain the change in band gradient. Calculate the ratio of the band gradient change to the perturbation increment to obtain the sensitivity coefficient of the structural parameter to the band gradient of each interface. Arrange the sensitivity coefficients according to the interface number and the structural parameter number to construct a sensitivity matrix. The row index of the sensitivity matrix corresponds to the interface position, the column index corresponds to the structural parameter, and the matrix element is the sensitivity coefficient of the corresponding structural parameter to the band gradient of the corresponding interface.

[0012] Based on the optimized geometric parameters and material composition parameters, the lattice mismatch strain distribution and dislocation line density distribution at each heterojunction interface are calculated, including: Based on the optimized material composition parameters, the lattice constant and elastic modulus tensor of each heterojunction layer are obtained by querying the material database. The difference in lattice constants between the materials on both sides of each heterojunction interface is calculated. The ratio of the difference in lattice constants to the lattice constant of the substrate material is defined as the lattice mismatch degree. Based on the lattice mismatch and the elastic modulus tensor of each layer, the in-plane strain component and vertical strain component of each heterojunction layer are calculated using the elastic mechanical continuum model to obtain the strain tensor distribution at each heterojunction interface, and the strain tensor distribution is converted into the lattice mismatch strain distribution. Based on the lattice mismatch strain distribution and the optimized geometric parameters, the Matthews-Blakeslee critical thickness theory is used to determine whether each heterojunction interface exceeds the critical thickness. For heterojunction interfaces that exceed the critical thickness, the equilibrium spacing of mismatched dislocations is calculated based on the lattice mismatch strain distribution. The equilibrium spacing is jointly determined by the strain relaxation condition and the dislocation elastic energy balance condition. The reciprocal of the equilibrium spacing is used as the dislocation line density at the interface. The dislocation line density distribution is obtained by summing the dislocation line densities of all interfaces.

[0013] A strain buffer layer structure is inserted at the corresponding interface. By adjusting the thickness and composition gradient distribution of the strain buffer layer structure, the dislocation line density distribution is reduced to below the preset density threshold, resulting in the final structural parameters after defect suppression, including: Identify the interface location in the dislocation line density distribution that exceeds a preset density threshold, insert a strain buffer layer structure at the interface location, set the initial thickness of the strain buffer layer structure to a preset ratio of the thickness of adjacent layers, and divide the strain buffer layer structure into several sub-layers along the growth direction. The material composition of each sublayer is set to vary in a piecewise manner between the material components on both sides of the interface. The initial component distribution is determined by linear interpolation, and the component gradient value of each sublayer is used as the parameter to be optimized. The corresponding lattice constant is calculated based on the material composition of each sublayer. The lattice mismatch strain of each sublayer is calculated based on the difference between the lattice constant and the lattice constant of the adjacent layer. The lattice mismatch strain of each sublayer is substituted into the elastic strain energy formula to calculate the spatial distribution of strain energy density inside the strain buffer layer. The strain energy density characterizes the driving force of dislocation nucleation. Based on the spatial distribution of strain energy density and strain relaxation theory, the dislocation line density update value at the interface position after inserting the strain buffer layer is calculated. When the dislocation line density update value decreases to below the preset density threshold, the current total thickness of the strain buffer layer and the material composition distribution of each sublayer are recorded. The total thickness of the strain buffer layer and the material composition distribution of each sublayer are integrated into the optimized geometric parameters and material composition parameters. The number of heterojunction layers and the position coordinates of each layer are updated to obtain the final structural parameters after defect suppression.

[0014] A second aspect of the present invention provides a performance simulation and optimization system for semiconductor multilayer heterostructures, comprising: The parameter acquisition unit is used to acquire the initial geometric parameters, material composition distribution parameters, and target band gradient distribution of the semiconductor multilayer heterostructure. Based on the initial geometric parameters and material composition distribution parameters, the band discontinuity values ​​and carrier concentration distribution at each heterojunction interface are calculated by self-consistently solving the Poisson equation and the Schrödinger equation. The parameter optimization unit is used to adjust the thickness parameters and material composition parameters of each heterojunction layer in reverse according to the deviation between the band discontinuity value and the target band gradient distribution, so that the conduction band order and valence band order at the interface of each heterojunction meet the preset gradient continuity constraint conditions, and obtain the optimized geometric parameters and material composition parameters. The strain calculation unit is used to calculate the lattice mismatch strain distribution and dislocation line density distribution at each heterojunction interface based on the optimized geometric parameters and material composition parameters. The defect suppression unit is used to insert a strain buffer layer structure at the corresponding interface when there is an interface position in the dislocation line density distribution that exceeds a preset density threshold. By adjusting the thickness and composition gradient distribution of the strain buffer layer structure, the dislocation line density distribution is reduced to below the preset density threshold, and the final structural parameters after defect suppression are obtained. The final structural parameters after defect suppression are then output as the fabrication process parameters for the semiconductor multilayer heterostructure.

[0015] A third aspect of the present invention provides an electronic device, comprising: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0016] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0017] This method enables the synergistic optimization design of band structures and defects in semiconductor multilayer heterostructures. By self-consistently solving the Poisson and Schrödinger equations, the method accurately obtains the band discontinuities and carrier distribution under the initial structure, providing an accurate physical basis for subsequent optimization. Based on the target band gradient distribution, the thickness and material composition of each layer are controlled inversely to ensure that the conduction and valence band orders satisfy the preset gradient continuity constraints, thereby effectively controlling the transport behavior of carriers and optimizing the electrical performance of the device.

[0018] After obtaining the optimized geometric and compositional parameters, the lattice mismatch strain and dislocation linear density distribution were further calculated, enabling a comprehensive evaluation from electrical performance to structural reliability. When interfaces with excessive dislocation density were identified, the interfacial stress was effectively relaxed by intelligently inserting and optimizing the strain buffer layer and adjusting its thickness and compositional gradient, thus suppressing the dislocation density below the threshold and significantly improving the structural quality and interfacial stability of the heterojunction.

[0019] This method closely integrates band engineering and defect engineering, simultaneously solving two core challenges in high-performance heterostructure design: band modulation and defect suppression, through a closed-loop feedback optimization process. The final output fabrication process parameters possess both optimized band gradient and low defect density, providing a direct and reliable design basis for fabricating high-performance, high-reliability semiconductor optoelectronic and electronic devices. Attached Figure Description

[0020] Figure 1 A flowchart illustrating the performance simulation and optimization method for semiconductor multilayer heterostructures; Figure 2 A flowchart illustrating the process for determining the optimized geometric parameters and material composition parameters. Detailed Implementation

[0021] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0022] The technical solution of the present invention will be described in detail below with reference to specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described again in some embodiments.

[0023] Figure 1 This is a flowchart illustrating the performance simulation and optimization method for semiconductor multilayer heterostructures according to an embodiment of the present invention, as shown below. Figure 1 As shown, the performance simulation and optimization method for semiconductor multilayer heterostructures includes: The initial geometric parameters, material composition distribution parameters, and target band gradient distribution of the semiconductor multilayer heterostructure are obtained. Based on the initial geometric parameters and material composition distribution parameters, the band discontinuity values ​​and carrier concentration distribution at each heterojunction interface are calculated using the Poisson equation and the Schrödinger equation in a self-consistent manner. Based on the deviation between the band discontinuity value and the target band gradient distribution, the thickness parameters and material composition parameters of each heterojunction are adjusted in reverse to ensure that the conduction band order and valence band order at the interface of each heterojunction satisfy the preset gradient continuity constraint condition, thereby obtaining the optimized geometric parameters and material composition parameters. Based on the optimized geometric parameters and material composition parameters, the lattice mismatch strain distribution and dislocation line density distribution at each heterojunction interface are calculated. When there is an interface location in the dislocation line density distribution that exceeds a preset density threshold, a strain buffer layer structure is inserted at the corresponding interface. By adjusting the thickness and composition gradient distribution of the strain buffer layer structure, the dislocation line density distribution is reduced to below the preset density threshold, and the final structural parameters after defect suppression are obtained. The final structural parameters after defect suppression are output as the fabrication process parameters of the semiconductor multilayer heterostructure.

[0024] In one optional implementation, based on the initial geometric parameters and material composition distribution parameters, the band discontinuities and carrier concentration distributions at each heterojunction interface are calculated using a self-consistent solution of the Poisson and Schrödinger equations, including: The spatial coordinate range of each heterojunction layer is determined based on the initial geometric parameters. The dielectric constant distribution and effective mass distribution of each layer are determined based on the material composition distribution parameters. The spatial region of each heterojunction layer is discretized into a grid. A difference scheme of the Poisson equation is established on the discrete grid nodes. The Poisson equation describes the relationship between electrostatic potential and carrier charge density. An eigenvalue problem of the Schrödinger equation is established on the discrete grid nodes. The potential energy term of the Schrödinger equation is composed of the electrostatic potential and the intrinsic band structure of the material. Solving the eigenvalue problem yields the wave function distribution and intrinsic energy of each quantum state. Based on the wave function distribution and intrinsic energy, the electron concentration and hole concentration at each grid node are calculated using the Fermi-Dirac distribution function. The calculated carrier concentration is substituted into the Poisson equation to update the electrostatic potential distribution. The Schrödinger equation and the Poisson equation are solved iteratively until the maximum difference between the electrostatic potential distributions obtained from two adjacent iterations is less than the preset convergence criterion, thus obtaining the self-consistent electrostatic potential distribution and carrier concentration distribution. Based on the electrostatic potential distribution obtained after self-consistent solution, the energy values ​​of the conduction band edge and the valence band edge on both sides of each heterojunction interface are extracted. The difference between the energy values ​​of the conduction band edge is calculated to obtain the conduction band discontinuity value, and the difference between the energy values ​​of the valence band edge is calculated to obtain the valence band discontinuity value.

[0025] After obtaining the initial geometric parameters and material composition distribution parameters of the semiconductor multilayer heterostructure, it is necessary to accurately calculate the band discontinuities and carrier concentration distributions at each heterojunction interface through self-consistent solutions of the Poisson and Schrödinger equations. This calculation process is the core step in achieving accurate bandgap engineering design, and can accurately reflect the influence of quantum confinement effects and carrier redistribution on the band structure.

[0026] When determining the spatial coordinate range of each heterojunction layer based on initial geometric parameters, a unified spatial coordinate system needs to be established. This system is based on the growth direction. Axis direction, set the first The starting coordinates of the layered heterojunction are Thickness is Then the spatial coordinate range of this layer is For those containing The multi-layered heterojunction structure has a total spatial range of [missing information]. This coordinate definition method clearly describes the spatial relationships between each layer, facilitating subsequent physical field calculations. When determining the dielectric constant distribution and effective mass distribution of each layer based on material composition distribution parameters, the influence of material composition on these physical parameters needs to be considered. For ternary compound semiconductor materials, such as... Its relative permittivity It can be expressed as a linear interpolation relationship. ,in It is an aluminum component. and These are the dielectric constants of the binary compounds. Effective electron mass. Similar interpolation relationships are followed, but for some material systems, a quadratic nonlinear interpolation formula may be needed to more accurately describe component dependence. Determining the effective mass of a hole is more complex, requiring differentiation between heavy and light holes and separate interpolation calculations.

[0027] When discretizing the spatial regions of each heterojunction layer, a one-dimensional non-uniform meshing scheme is adopted. Near the heterojunction interface and in regions with drastic changes in doping concentration, a smaller mesh spacing is used to improve computational accuracy; a typical mesh spacing can be set to 0.1 nm to 0.5 nm. In regions with relatively uniform material parameters, the mesh spacing can be appropriately increased to 1 nm to 2 nm to reduce the total number of computational nodes and improve computational efficiency. Assuming that the discretized result is... There are 1 grid node, and the node coordinates are as follows: ,in .

[0028] When establishing a difference scheme for the Poisson equation on discrete grid nodes, the central difference method is used for spatial derivative discretization. The Poisson equation describes the relationship between electrostatic potential and carrier charge density, and its standard form is: ,in It is the electrostatic potential. For dielectric constant distribution, This represents the charge density distribution. Charge density includes the concentration of ionized donors. ionization acceptor concentration electron concentration With hole concentration The contribution is expressed as ,in This represents the elementary charge. For nodes... The Poisson equation at point A can be written in the difference scheme as follows: ,in , , For nodes and The average dielectric constant between.

[0029] When establishing the eigenvalue problem of the Schrödinger equation on discrete grid nodes, the finite difference method or the finite element method is used for discretization. The Schrödinger equation under the one-dimensional effective mass approximation is as follows: ,in To reduce Planck's constant, For effective mass distribution, For wave function, For intrinsic energy, This is the effective potential energy. Effective potential energy is derived from electrostatic potential energy. Together with the intrinsic band edge positions of the material, it forms the effective potential energy for conduction band electrons. ,in The energy at the bottom of the conduction band. For a hole in the valence band, the effective potential energy is... ,in The valence band top energy is used. By using finite difference discretization, the Schrödinger equation is transformed into a matrix eigenvalue problem. ,in For Hamiltonian matrix, Let be the wave function vector. Numerical linear algebra methods, such as the Lanczos algorithm or the Jacobi-Davidson algorithm, are used to solve this eigenvalue problem to obtain the wave function distribution and eigenenergy of each quantum state.

[0030] When calculating the electron and hole concentrations at each grid node based on the wave function distribution and intrinsic energy, the occupancy probability of quantum states needs to be considered. Under thermal equilibrium conditions, the energy is... The occupancy probability of a quantum state is given by the Fermi-Dirac distribution function, expressed as: ,in For Fermi level, Boltzmann constant, Temperature. Node The electron concentration at that location is calculated as follows: The summation iterates through all conduction band quantum states. The hole concentration is calculated as follows: The summation iterates through all valence band quantum states. In actual calculations, the Fermi level... This is determined by the condition of electrical neutrality, which requires that the total charge of the entire structure be zero.

[0031] When the calculated carrier concentration is substituted into the Poisson equation to update the electrostatic potential distribution, a self-consistent iterative loop is formed. The specific steps are as follows: given an initial electrostatic potential distribution... Solving the Schrödinger equation yields the wave function and eigenenergy, and the carrier concentration is then calculated. and Substituting this into the Poisson equation and solving for the updated electrostatic potential distribution yields the updated electrostatic potential distribution. Iteratively solving the Schrödinger equation and the Poisson equation yields a sequence. To accelerate convergence, a hybrid update scheme can be adopted, namely... ,in The electrostatic potential obtained directly from this iteration, This is a mixed parameter, typically ranging from 0.1 to 0.5.

[0032] When determining whether a self-consistent iteration has converged, the maximum difference criterion is used. This involves calculating the maximum difference in electrostatic potential distributions obtained from two consecutive iterations. ,when If the value is less than the preset convergence criterion, the self-consistent solution is considered convergent. A typical value for the convergence criterion is... Electron Volt Electron volts.

[0033] An eigenvalue problem of the Schrödinger equation is established on the discrete grid nodes. The potential energy term of the Schrödinger equation is composed of the electrostatic potential and the intrinsic band structure of the material. Solving the eigenvalue problem yields the wave function distribution and eigenenergy of each quantum state, including: The numerical distribution of the electrostatic potential is extracted on the discrete grid nodes, and the conduction band bottom energy and valence band top energy of the material corresponding to each grid node are read. The conduction band bottom energy and valence band top energy are used as characterization parameters of the intrinsic band structure of the material. The effective potential energy distribution of the electronic state is constructed by superimposing the numerical distribution of the electrostatic potential with the bottom energy of the conduction band, and the effective potential energy distribution of the hole state is constructed by superimposing the numerical distribution of the electrostatic potential with the top energy of the valence band. Based on the spatial coordinates of the discrete grid nodes, a discretized expression for the kinetic energy operator is constructed. The discretized expression for the kinetic energy operator uses a finite difference scheme to discretize the second spatial derivative of the wave function. Substituting the discretized expression of the kinetic energy operator and the effective potential energy distribution into the Schrödinger equation, an eigenvalue problem in matrix form is formed on the discrete grid nodes. The coefficient matrix of the eigenvalue problem is composed of kinetic energy terms and potential energy terms. The coefficient matrix of the eigenvalue problem is numerically solved to obtain the eigenvalues ​​and eigenvectors of each quantum state. The quantum states are sorted by energy level according to the eigenenergy, and the eigenenergy and wavefunction distribution of the electronic state near the bottom of the conduction band and the eigenenergy and wavefunction distribution of the hole state near the top of the valence band are extracted. The wave function distribution is normalized in space, and the spatial distribution of wave function probability density for each quantum state is calculated. The spatial distribution of wave function probability density is used to characterize the spatial localization characteristics of charge carriers in the heterojunction structure.

[0034] After completing the self-consistent iterative solution of the Poisson equation, the numerical distribution of electrostatic potential is extracted at the discrete grid nodes. This numerical distribution of electrostatic potential reflects the space charge effect caused by carrier redistribution in the heterostructure. Specifically, during the extraction process, the converged electrostatic potential values ​​at each node are read sequentially according to their spatial coordinates, forming a one-dimensional or multi-dimensional potential energy array. Simultaneously with the extraction of electrostatic potential values, the conduction band bottom energy and valence band top energy of the corresponding material at each grid node are read. These band parameters reflect the intrinsic electronic level structure characteristics of the material; the conduction band bottom energy represents the minimum allowable energy of the electronic state, and the valence band top energy represents the maximum allowable energy of the hole state. For multilayer heterostructures, there are differences in the conduction band bottom energy and valence band top energy between different material layers, forming band order at the heterojunction interface. During the reading process, based on the material layer type of each grid node, the corresponding conduction band bottom energy and valence band top energy values ​​are retrieved from the material parameter database, and a correspondence is established between them and the spatial position of that node.

[0035] To construct the effective potential energy distribution of electronic states in the Schrödinger equation, the numerical distribution of the electrostatic potential is superimposed with the conduction band bottom energy. This superposition process is performed independently at each grid node. Specifically, the electrostatic potential value of the node is added to the conduction band bottom energy of the corresponding material to obtain the effective potential energy of the electronic state at that location. This effective potential energy comprehensively considers the combined effects of the material's intrinsic band structure and space charge effect, truly reflecting the potential field of electron motion in the heterostructure. A similar method is used to construct the effective potential energy distribution of hole states, superimposing the numerical distribution of the electrostatic potential with the valence band top energy. Since holes carry a positive charge, their behavior in the electrostatic potential is opposite to that of electrons. Therefore, when constructing the effective potential energy of hole states, the sign difference of the electrostatic potential effect needs to be considered. In practice, the effective potential energy distribution of hole states at that location is obtained by subtracting the absolute value of the valence band top energy of the corresponding material from the electrostatic potential value of each node.

[0036] Based on the spatial coordinates of the discrete grid nodes, a discretized expression for the kinetic energy operator is constructed. The kinetic energy operator in the Schrödinger equation contains a second-order spatial derivative term of the wave function in continuous space. To perform numerical solutions on the discrete grid, a finite difference scheme is used to discretize and approximate this second-order derivative. Specifically, a central difference scheme is employed, for the position... At a grid node, the second derivative of its wavefunction can be approximated as: ,in Represents a node wave function value at that point This represents the spatial spacing between grid nodes. The complete discretized expression of the kinetic energy operator also needs to consider the spatial variation characteristics of the effective mass. In heterogeneous structures, the effective mass of charge carriers differs between different material layers, thus the kinetic energy operator exhibits discontinuities at the interface. To accurately handle this discontinuity, a harmonic or arithmetic mean of the effective mass is used for interpolation at the interface grid nodes to ensure the numerical stability of the discretized expression of the kinetic energy operator at the interface.

[0037] Substituting the discretized expression of the kinetic energy operator and the effective potential energy distribution into the Schrödinger equation, a matrix-form eigenvalue problem is formed on the discrete grid nodes. This matrix-form eigenvalue problem can be expressed as: ,in For Hamiltonian matrix, For wave function vectors, Eigenenergy. Hamiltonian matrix. During the construction process, the diagonal elements contain the main diagonal part of the kinetic energy term and the effective potential energy of the corresponding node, while the off-diagonal elements originate from the coupling terms between adjacent nodes in the discretized expression of the kinetic energy operator. For those containing A one-dimensional heterogeneous structure with grid nodes, the Hamiltonian matrix is... The tridiagonal matrix has a sparsity characteristic that facilitates the use of efficient numerical algorithms for solving it. During the matrix assembly process, appropriate boundary conditions need to be applied to the boundary nodes, usually Dirichlet boundary conditions, which are zero wave functions at the structural boundaries. This corresponds to the physical picture of charge carriers being confined inside the heterostructure.

[0038] The coefficient matrix of the eigenvalue problem is numerically solved to obtain the eigenvalues ​​and eigenvectors corresponding to each quantum state. Since the Hamiltonian matrix possesses Hermitian properties, its eigenvalues ​​are all real numbers and its eigenvectors are mutually orthogonal, allowing for the use of efficient algorithms specifically designed for Hermitian matrices. In practice, for small-scale systems, a fully diagonalized method can be used, diagonalizing the Hamiltonian matrix through Jacobian rotation or Haushold transformation to directly obtain all eigenvalues ​​and eigenvectors. For large-scale systems, since only a finite number of quantum states near the band edges need to be considered, a partially diagonalized method, such as the Lanczos or Davidson algorithm, can be used to calculate only the lowest energy level electronic states near the bottom of the conduction band and the highest energy level hole states near the top of the valence band, significantly reducing computational complexity. The obtained eigenvalues ​​correspond to the energy levels of each quantum state, and the eigenvectors correspond to the wavefunction values ​​of each quantum state at discrete grid nodes.

[0039] The quantum states are ordered by energy level based on their intrinsic energies. For electronic states, they are arranged in ascending order of intrinsic energy, with the lowest energy state labeled as the ground state, followed by the first excited state, the second excited state, and so on. For hole states, since the valence band apex corresponds to the highest allowed state, they are arranged in descending order of intrinsic energy, with the highest energy state labeled as the heavy hole ground state, followed by the heavy hole excited states. In semiconductor heterostructures, due to the quantum confinement effect, the allowed energy levels of charge carriers exhibit a discrete distribution, and the spacing between adjacent energy levels is closely related to the geometry of the quantum well or quantum barrier. After sorting, the intrinsic energies and wavefunction distributions of electronic states near the conduction band apex are extracted, typically selecting the 5 to 10 lowest energy electronic states for subsequent analysis. Simultaneously, the intrinsic energies and wavefunction distributions of hole states near the valence band apex are extracted, including heavy hole and light hole states, which play a dominant role in the transition processes of optoelectronic devices.

[0040] The wavefunction is distributed spatially and then normalized. Wavefunction normalization is a fundamental requirement in quantum mechanics, ensuring that the probability of finding a charge carrier in the entire space is 1. Specifically, in the normalization process, the integral of the modulus square of the wavefunction vector of each quantum state over all grid nodes is calculated. This integral is approximated by summation over the discrete grid, i.e. The normalization coefficient is the reciprocal square root of the integral value. Multiplying this normalization coefficient by each element of the original wavefunction vector yields the normalized wavefunction distribution. After normalization, the spatial distribution of the wavefunction probability density for each quantum state is calculated. This probability density is defined as the square root of the wavefunction modulus. The normalized wavefunction modulus squared can be directly calculated at each grid node. The spatial distribution of the wavefunction probability density intuitively reflects the spatial localization characteristics of charge carriers in the heterostructure. In the quantum well region, the wavefunction probability density exhibits a clear peak distribution, indicating that the charge carriers are strongly confined in the potential well; in the barrier region, the wavefunction probability density decays rapidly, exhibiting exponentially decaying tunneling tail wave characteristics. By analyzing the spatial distribution characteristics of the wavefunctions of different quantum states, the difference between bound states and continuous states can be identified, providing a wavefunction basis for subsequent calculation of charge carrier concentration and analysis of transport characteristics.

[0041] Figure 2 A flowchart illustrating the process for determining the optimized geometric parameters and material composition parameters.

[0042] Based on the deviation between the band discontinuity value and the target band gradient distribution, the thickness parameters and material composition parameters of each heterojunction layer are adjusted in reverse to ensure that the conduction band order and valence band order at each heterojunction interface satisfy the preset gradient continuity constraint condition, thereby obtaining the optimized geometric parameters and material composition parameters, including: Calculate the actual band gradient value at each heterojunction interface. The actual band gradient value is the ratio of the difference in band discontinuity between two adjacent interfaces to the interface spacing. Extract the target band gradient value corresponding to each interface position from the target band gradient distribution. Calculate the deviation vector between the actual band gradient value and the target band gradient value. A sensitivity matrix is ​​established between band gradient deviation and structural parameter adjustment. The elements of the sensitivity matrix represent the influence coefficients of unit thickness change or unit composition change on the band gradient of each interface. Based on the sensitivity matrix and the deviation vector, the structural parameter adjustment that minimizes the norm of the deviation vector is solved by the least squares method. The structural parameter adjustment includes the thickness adjustment and material composition adjustment of each layer. The structural parameter adjustment is superimposed on the current thickness parameter and material composition parameter to obtain the updated thickness parameter and material composition parameter. Based on the updated thickness parameter and material composition parameter, the band discontinuity value at each heterojunction interface is recalculated. When the conduction band order and valence band order satisfy the preset gradient continuity constraint condition, the current thickness parameter and material composition parameter are output as the optimized geometric parameter and material composition parameter.

[0043] After obtaining the band discontinuities at the interfaces of various heterojunctions in a semiconductor multilayer heterostructure, the band gradient needs to be accurately evaluated and controlled. The band gradient reflects the distribution characteristics of the transport barrier of charge carriers in the heterojunction and directly affects the electrical performance and optical response characteristics of the device.

[0044] For including The multilayer structure of the layered heterojunction has a total number of interfaces. The first. For the first Each interface ( The actual band gradient value is calculated by differencing the band discontinuities at adjacent interfaces. (Conduction band gradient) It can be represented as ,in Indicates the first The conduction band discontinuity value of each interface. Indicates the first The interface and the first Spatial spacing between interfaces. Valence band gradient. Using the same calculation method, for ,in Indicates the first The valence band discontinuity of each interface. Interface spacing. By the Layer and first The thickness of the layers is determined by both layers. Under the single-interface approximation, it can be taken as the average of the thicknesses of the two layers or the distance between their geometric centers.

[0045] The design values ​​corresponding to each interface position are extracted from the target band gradient distribution. The target band gradient distribution is usually pre-defined according to the device application requirements; for example, a uniform ionization rate gradient needs to be established in avalanche photodetectors, while a specific subband transition energy distribution needs to be achieved in quantum cascade lasers. For the first... A single interface extracts the target conduction band gradient from the target distribution. Gradient of target price band .

[0046] A bias vector is constructed to quantify the difference between the actual band gradient and the target value. The conduction band gradient bias vector is defined. ,in Similarly, define the valence band gradient deviation vector. ,in The elements of these two deviation vectors reflect the degree of deviation of the band gradient at each interface position. Positive values ​​indicate that the actual gradient is greater than the target gradient, and negative values ​​indicate that the actual gradient is less than the target gradient.

[0047] Establish the sensitivity relationship between band gradient deviation and structural parameter adjustments. Structural parameters include the thickness parameters and material composition parameters of each layer. For an N-layer structure, the thickness parameter vector is... The material composition parameter vector is In ternary or quaternary compound semiconductors, the component parameters may be multidimensional, requiring the establishment of sensitivity relationships for each component separately.

[0048] The sensitivity matrix characterizes the effect of a unit change in structural parameters on the band gradient. The sensitivity matrix of the conduction band gradient to thickness is also relevant. elements Indicates the first When the layer thickness changes by one unit length, the first The change in the gradient of the interface conduction band. This element can be calculated using the numerical perturbation method, and the i-th... The layer thickness increases by a tiny amount Recalculating the Poisson and Schrödinger equations yields new band discontinuities, which in turn lead to the calculation of new band gradients. Similarly, the sensitivity matrix of the conduction band gradient to the component. elements Characterizing the first Changes in layer composition affect the first The influence of the interface conduction band gradient is calculated as follows: ,in For component perturbations. Sensitivity matrix of valence band gradient. and The same calculation method is used.

[0049] The thickness sensitivity matrix and the component sensitivity matrix are combined into a comprehensive sensitivity matrix. For conduction band gradients, the comprehensive sensitivity matrix... Its dimensions are For the valence band gradient, the comprehensive sensitivity matrix... Same dimensions. Define a vector of structural parameter adjustments. This includes thickness and composition adjustments for all layers.

[0050] A system of linear equations is established based on the sensitivity matrix and the bias vector. The adjustment relationship of the conduction band gradient can be expressed as follows: This equation represents the elimination of conduction band gradient bias by adjusting structural parameters. The adjustment relationship for the valence band gradient is as follows: Since both the conduction band and valence band need to satisfy the gradient requirement simultaneously, the two sets of equations are combined into an augmented system. .

[0051] The least squares method is used to solve the overdetermined system of equations. The cost function is defined as: The cost function represents the sum of squares of the adjusted residual biases. Finding the extreme value of the cost function yields the normal equation for the optimal adjustment. ,in , The normal equation is solved using matrix factorization to obtain the structural parameter adjustment vector. .

[0052] In practical calculations, it may be necessary to impose physical constraints on the adjustment amount. The thickness adjustment amount must meet the following requirements. ,in The minimum achievable thickness is typically a few monolayer atoms. The composition adjustment amounts must meet certain requirements. Ensure that the components are within physically permissible limits. When the unconstrained least squares solution violates these constraints, a constrained optimization algorithm must be used, such as adding a penalty term or employing the projected gradient method. The adjustments obtained from the solution are then added to the current structural parameters.

[0053] Establish a sensitivity matrix between band gradient bias and structural parameter adjustments, including: A set of structural parameters that need to be controlled is determined, the set of structural parameters including the thickness parameters of each heterojunction layer and the material composition parameters of each layer. A preset perturbation increment is applied to each structural parameter in the set of structural parameters, the perturbation increment being a preset percentage of the current value of the structural parameter. After applying the aforementioned perturbation increment, the band discontinuity values ​​at each heterojunction interface are recalculated using the self-consistent solution method of the Poisson equation and the Schrödinger equation, thus obtaining the distribution of the perturbed band discontinuity values. The actual band gradient value after the disturbance is calculated based on the distribution of the band discontinuities after the disturbance. The difference between the actual band gradient value after the disturbance and the actual band gradient value before the disturbance is obtained to obtain the change in band gradient. Calculate the ratio of the band gradient change to the perturbation increment to obtain the sensitivity coefficient of the structural parameter to the band gradient of each interface. Arrange the sensitivity coefficients according to the interface number and the structural parameter number to construct a sensitivity matrix. The row index of the sensitivity matrix corresponds to the interface position, the column index corresponds to the structural parameter, and the matrix element is the sensitivity coefficient of the corresponding structural parameter to the band gradient of the corresponding interface.

[0054] In practical band gradient manipulation, it is necessary to clarify the degree to which changes in various structural parameters affect the band gradient distribution. A set of structural parameters to be manipulated needs to be determined, encompassing all independently adjustable physical quantities in the multilayer heterostructure. Specifically, the set of structural parameters includes the thickness parameters and material composition parameters of each heterojunction layer. For structures containing... The multilayer structure of a layered heterojunction, with the set of thickness parameters denoted as . ,in Indicates the first The thickness of the layer, measured in nanometers. The set of material composition parameters is denoted as... ,in Indicates the first The material composition ratio of the layer, for ternary compound semiconductor materials such as , The values ​​range from 0 to 1. These structural parameters collectively determine the band structure's band distribution characteristics and are key variables for realizing band structure engineering design.

[0055] For each structural parameter in the set of structural parameters, a preset perturbation increment is applied. The selection of the perturbation increment follows the numerical accuracy requirements of the finite difference method; it needs to be small enough to ensure the effectiveness of the linear approximation, but not so small that numerical calculation errors dominate. Typically, the perturbation increment is set as a preset percentage of the current value of the structural parameter, with a typical range of 0.5% to 2%. Taking the thickness parameter as an example, suppose the current thickness of a certain layer is... For nanometers, if the perturbation percentage is 1%, then the perturbation increment is... Nanoscale. Regarding material composition parameters, if the current composition is... Using the same percentage of disturbance, the disturbance increment is... The perturbation method employs a single-parameter, one-by-one perturbation approach, where only one structural parameter is perturbed at a time, while other parameters remain unchanged, thus isolating the influence of individual parameters. This method ensures that the physical meaning of the sensitivity coefficient is clear, facilitating the formulation of subsequent parameter control strategies.

[0056] After applying the perturbation increment, the electronic band distribution of the entire multilayer heterostructure needs to be recalculated. A self-consistent solution method using the Poisson and Schrödinger equations is employed, which can simultaneously consider the electrostatic potential distribution and quantum confinement effects. The Poisson equation describes the relationship between the space charge distribution and the electrostatic potential. ,in It is the electrostatic potential. Where is the dielectric constant, and Electron and hole concentrations, respectively. Let and represent the concentrations of ionized donors and acceptors, respectively. The Schrödinger equation describes the quantum state distribution of charge carriers. ,in For effective quality, and The first The wave functions and energy levels of each quantum state are determined. These two coupled equations are solved iteratively until the carrier concentration distribution and electrostatic potential distribution are self-consistent, yielding a convergent band structure. Under the perturbed structural parameters, the complete self-consistent solution process is re-executed to calculate the band discontinuities at each heterojunction interface, obtaining the perturbed band discontinuity distribution.

[0057] The actual band gradient value after perturbation is calculated based on the distribution of band discontinuities after perturbation. For two adjacent interfaces... and The conduction band gradient is defined as follows: ,in Indicates the first Conduction band discontinuities at the interface. Valence band and energy band gradients. A similar definition is used. Before applying the perturbation increment, the actual band gradient values ​​at each interface have already been obtained, and are denoted as the actual band gradient values ​​before the perturbation. After applying the perturbation, the recalculated band gradient value is denoted as... The change in band gradient is obtained by subtracting the actual band gradient value after perturbation from the actual band gradient value before perturbation. This change directly reflects the degree to which the perturbation of specific structural parameters affects the energy band gradient of each interface, and is the basic data for constructing the sensitivity matrix.

[0058] The ratio of the band gradient change to the perturbation increment is calculated to obtain the sensitivity coefficient of the structural parameters to the band gradients at each interface. This is based on the thickness parameter. For example, its application to the first The sensitivity coefficient of the conduction band gradient of an interface is defined as follows: The unit is the rate of change of energy gradient per nanometer. Similarly, material composition parameters... For the The sensitivity coefficient of the conduction band gradient of the interface is This coefficient characterizes the efficiency with which compositional changes modulate the band gradient. Complete sensitivity information is obtained by applying perturbations to each structural parameter and calculating the corresponding sensitivity coefficients. The magnitude of these sensitivity coefficients reflects the effectiveness of parameter manipulation; larger absolute values ​​indicate that the parameter has a stronger ability to modulate the band gradient, and these parameters should be prioritized for adjustment in subsequent optimization processes.

[0059] Arrange the sensitivity coefficients according to the interface number and structural parameter number to construct a sensitivity matrix. The sensitivity matrix is ​​a two-dimensional array, where the row index corresponds to the interface position and the column index corresponds to the structural parameter. For parameters containing... The multi-layered heterogeneous structure has a total of There are 1 heterojunction interface, and the total number of structural parameters is . Therefore, the dimension of the sensitivity matrix is ​​. The first matrix Line number Column elements represent the first The structural parameter affects the first The sensitivity coefficients of each interface can be gradient-based. The specific matrix element arrangement is as follows: [Front...] The sensitivity of the thickness parameters for each layer is listed below. The columns correspond to the sensitivity of material composition parameters for each layer. Matrix elements. The physical meaning is that when the first When the structural parameter changes by a unit, the first... The band gradient of each interface occurs accordingly. The magnitude of the change. The construction of the sensitivity matrix transforms the multi-parameter optimization problem into a linear algebra problem, enabling the rapid determination of the required parameter adjustment scheme using matrix operations.

[0060] In practical applications, the sensitivity matrix also needs to consider numerical stability and physical rationality. Some structural parameters have minimal impact on the band gradient at long-distance interfaces, with corresponding sensitivity coefficients close to zero or within the numerical error range. These elements can be filtered out by setting thresholds. Setting sensitivity coefficients with absolute values ​​less than a preset threshold to zero simplifies subsequent matrix operations and improves computational efficiency. The threshold selection needs to be determined based on the physical properties of the specific material system, typically set to 0.1% to 1% of the absolute value of the maximum sensitivity coefficient. The condition number of the sensitivity matrix is ​​an important indicator of its numerical properties. An excessively large condition number indicates that the matrix is ​​close to singular, leading to overly sensitive parameter adjustment schemes to input errors. When the condition number exceeds a set upper limit, regularization methods or dimensionality reduction techniques are needed to improve the numerical properties of the matrix and ensure the robustness of the optimization algorithm. The parameter control model established through the sensitivity matrix enables precise adjustment of the band gradient distribution, providing a quantitative design tool for performance optimization of semiconductor multilayer heterostructures.

[0061] Based on the optimized geometric parameters and material composition parameters, the lattice mismatch strain distribution and dislocation line density distribution at each heterojunction interface are calculated, including: Based on the optimized material composition parameters, the lattice constant and elastic modulus tensor of each heterojunction layer are obtained by querying the material database. The difference in lattice constants between the materials on both sides of each heterojunction interface is calculated. The ratio of the difference in lattice constants to the lattice constant of the substrate material is defined as the lattice mismatch degree. Based on the lattice mismatch and the elastic modulus tensor of each layer, the in-plane strain component and vertical strain component of each heterojunction layer are calculated using the elastic mechanical continuum model to obtain the strain tensor distribution at each heterojunction interface, and the strain tensor distribution is converted into the lattice mismatch strain distribution. Based on the lattice mismatch strain distribution and the optimized geometric parameters, the Matthews-Blakeslee critical thickness theory is used to determine whether each heterojunction interface exceeds the critical thickness. For heterojunction interfaces that exceed the critical thickness, the equilibrium spacing of mismatched dislocations is calculated based on the lattice mismatch strain distribution. The equilibrium spacing is jointly determined by the strain relaxation condition and the dislocation elastic energy balance condition. The reciprocal of the equilibrium spacing is used as the dislocation line density at the interface. The dislocation line density distribution is obtained by summing the dislocation line densities of all interfaces.

[0062] Based on the optimized geometric and material composition parameters, the corresponding values ​​are input into the material database retrieval system. This database stores complete lattice constant data and sixth-order elastic modulus tensor data for common semiconductor materials at 300K. For ternary or quaternary compound semiconductors, their lattice constants are calculated using linear interpolation based on Vegard's law. (Al...) Ga Taking As material as an example, its lattice constant It is possible We obtained, among which It is an aluminum component. and These are the lattice constants of the binary compounds. For the elastic modulus tensor, the anisotropic elastic constants are also obtained using linear interpolation. , and .

[0063] After extracting the lattice constant values ​​of the materials on both sides of each heterojunction interface, the difference in lattice constant between adjacent layers is calculated. Taking the first... Layer and First Taking the heterojunction interface formed by the layers as an example, the difference in lattice constants for ,in and The first The lattice constants of the i+1th and i-th layers. In practical calculations, the substrate material is used as a reference, and the lattice mismatch is defined. The lattice constants of each layer and the substrate are given. The ratio of the difference to the substrate lattice constant, i.e. This dimensionless parameter directly reflects the degree of lattice mismatch between the epitaxial layer and the substrate; a positive value represents tensile strain, and a negative value represents compressive strain.

[0064] After obtaining the lattice mismatch of each layer, an elastic continuum model is established for strain analysis. It is assumed that in-plane lattice matching is maintained at the heterojunction interface, and the epitaxial layer generates biaxial strain in a direction parallel to the interface. For cubic semiconductor materials, the in-plane strain components... Equal to lattice mismatch According to the Poisson effect, the strain component in the vertical direction It is possible Calculate, where and For the first The elastic constants of the layer material. Constructing a complete strain tensor requires combining the in-plane biaxial strain and the perpendicular strain component into a symmetric second-order tensor. For heterostructures grown along the

[001] crystal orientation, the strain tensor... The non-zero components are diagonal elements and .

[0065] The strain tensors of each heterojunction layer are arranged along the growth direction to form the strain tensor distribution of the entire multilayer structure. This distribution exhibits a step change characteristic with the interface location as the dividing point. In the extremely thin regions on both sides of the interface, the actual strain state deviates from the prediction of the continuum model due to atomic-scale interface reconstruction, but in regions several atomic layers away from the interface, the continuum model has sufficient accuracy. The strain tensor distribution data is stored as an array structure indexed by interface location, with each array element containing the values ​​of the three principal strain components of the corresponding layer. To facilitate subsequent defect analysis, the strain tensor distribution is converted into a lattice mismatch strain distribution based on lattice mismatch degree, retaining the sign and numerical information of the in-plane strain components at each interface.

[0066] Based on lattice mismatch strain distribution data and optimized geometric parameters, it is determined whether the interfaces of each heterojunction exceed the critical thickness. The Matthews-Blakeslee theory gives the critical thickness of the mismatched epitaxial layer. With lattice mismatch The relationship is ,in The Burgers vector of the edge dislocation is the magnitude. For the first The Poisson's ratio of the layer material. This transcendental equation requires a numerical iterative method to solve. Initial settings. The estimated value is substituted into the right side of the equation to calculate the new critical thickness value. This iteration is repeated until the relative error between the two calculations is less than 0.1%. For zincblende structure semiconductors grown in the

[001] crystal orientation, the Burgers vector size... Pick Poisson's ratio pass calculate.

[0067] The calculated critical thickness value is compared one by one with the optimized actual thickness of each layer. If the... Actual thickness of the layer Less than critical thickness If the thickness is zero, the layer is in a pseudocrystalline growth state, and the lattice mismatch strain is entirely borne by elastic deformation, without generating mismatch dislocations. The dislocation line density at this interface is recorded as zero. If the actual thickness exceeds the critical thickness, the epitaxial layer will release some strain energy by forming mismatch dislocations, requiring further calculation of the equilibrium spacing and dislocation line density of the mismatch dislocations.

[0068] For heterojunction interfaces exceeding the critical thickness, mismatched dislocations are distributed in an array at the interface. The equilibrium spacing of the mismatched dislocations... It is jointly determined by the strain relaxation condition and the dislocation elastic energy equilibrium condition. The strain relaxation condition requires that the plastic strain introduced by the dislocation array compensates for part of the lattice mismatch, and the residual strain after relaxation... satisfy The dislocation elastic energy equilibrium condition is based on the self-energy expression of a dislocation line of unit length, considering the elastic interaction between dislocations and the equilibrium spacing. It is possible The approximate calculation shows that the formula comprehensively considers the effects of layer thickness and strain state on dislocation nucleation.

[0069] After calculating the equilibrium spacing of mismatched dislocations at each supercritical interface, its reciprocal is used as the dislocation linear density at that interface. ,Right now The unit is per centimeter. This parameter characterizes the number of dislocation lines passing through a unit length interface, directly reflecting the structural defect density at the interface. For interfaces not exceeding the critical thickness, the dislocation line density remains zero. The dislocation line density values ​​of all heterojunction interfaces are arranged according to the interface number or spatial position to form a dislocation line density distribution dataset. This distribution reveals the spatial distribution law of defects in multilayer heterostructures, providing a quantitative basis for determining the subsequent insertion position of strain buffer layers. In actual calculations, the orientation factor of dislocation lines also needs to be considered. For heterojunctions grown on

[001] , the main formation along the direction of the dislocation line is the orientation of the dislocation line. For 60° mixed dislocations in the direction of displacement, the dislocation line density calculation needs to consider the contributions of dislocations in two equivalent directions. The final line density is the value calculated in a single direction. After calculating the dislocation linear density distribution, the results are output in the form of a data table or curve, annotating key parameters such as the location of each interface, lattice mismatch, actual thickness, critical thickness, and dislocation linear density. This quantitative analysis provides a direct criterion for judging the structural quality; the interface locations where the dislocation linear density exceeds the preset threshold are the target areas where a strain buffer layer needs to be inserted.

[0070] A strain buffer layer structure is inserted at the corresponding interface. By adjusting the thickness and composition gradient distribution of the strain buffer layer structure, the dislocation line density distribution is reduced to below the preset density threshold, resulting in the final structural parameters after defect suppression, including: Identify the interface location in the dislocation line density distribution that exceeds a preset density threshold, insert a strain buffer layer structure at the interface location, set the initial thickness of the strain buffer layer structure to a preset ratio of the thickness of adjacent layers, and divide the strain buffer layer structure into several sub-layers along the growth direction. The material composition of each sublayer is set to vary in a piecewise manner between the material components on both sides of the interface. The initial component distribution is determined by linear interpolation, and the component gradient value of each sublayer is used as the parameter to be optimized. The corresponding lattice constant is calculated based on the material composition of each sublayer. The lattice mismatch strain of each sublayer is calculated based on the difference between the lattice constant and the lattice constant of the adjacent layer. The lattice mismatch strain of each sublayer is substituted into the elastic strain energy formula to calculate the spatial distribution of strain energy density inside the strain buffer layer. The strain energy density characterizes the driving force of dislocation nucleation. Based on the spatial distribution of strain energy density and strain relaxation theory, the dislocation line density update value at the interface position after inserting the strain buffer layer is calculated. When the dislocation line density update value decreases to below the preset density threshold, the current total thickness of the strain buffer layer and the material composition distribution of each sublayer are recorded. The total thickness of the strain buffer layer and the material composition distribution of each sublayer are integrated into the optimized geometric parameters and material composition parameters. The number of heterojunction layers and the position coordinates of each layer are updated to obtain the final structural parameters after defect suppression.

[0071] After identifying interfaces where the dislocation linear density exceeds a preset density threshold, strain buffer layers need to be introduced at these high-strain interface locations to suppress defects. By scanning the dislocation linear density distribution data across the entire heterostructure, the dislocation linear density value at each interface location is compared with the preset density threshold. This preset density threshold is typically set based on device performance requirements. When the dislocation linear density at a certain interface location exceeds this threshold, the spatial coordinates of that interface and the material parameters of the two adjacent layers are recorded, including material composition, lattice constant, and elastic constants.

[0072] For the identified high-strain interfaces, a strain buffer layer structure is inserted at the location. The initial thickness of the strain buffer layer needs to be set by comprehensively considering strain relaxation efficiency and material growth feasibility. Typically, the initial thickness is set as a preset ratio of the thickness of the thinner of the two adjacent layers, and this ratio is determined based on the lattice mismatch degree of the material system. For lattice mismatch less than 1%, the preset ratio can be set to 0.3 to 0.5 times; for lattice mismatch between 1% and 3%, the preset ratio increases to 0.5 to 0.8 times; when the lattice mismatch exceeds 3%, the preset ratio needs to reach 0.8 to 1.2 times. Assuming the two adjacent layers have thicknesses of... AlGaAs layer and thickness The GaAs layer, and If the lattice mismatch is approximately 0.15%, then the initial thickness of the strain buffer layer is... Set as .

[0073] After determining the total thickness of the strain buffer layer, it needs to be spatially discretized along the growth direction. The strain buffer layer is then divided into several sub-layers, the number of which is... The choice of sublayers directly affects the precision of composition gradient control. Too few sublayers will result in overly steep composition changes, failing to effectively suppress strain accumulation; too many sublayers will increase computational complexity and be difficult to control precisely in actual fabrication. The number of sublayers is usually determined based on the total thickness and typical growth precision, with the thickness of each sublayer preferably controlled between 2 nm and 5 nm. If the total thickness of the strain buffer layer is 30 nm, it can be divided into 8 to 12 equal sublayers, each with a thickness of approximately 2.5 nm to 3.75 nm.

[0074] The design of the material composition of each sublayer is the core of realizing the function of the strain buffer layer. The role of the strain buffer layer is to achieve a smooth transition from one lattice constant to another in space; therefore, the material composition of each sublayer should change according to a specific pattern between the material compositions on both sides of the interface. For ternary compound semiconductors, such as... Component parameters The value range is determined by the component values ​​of the two adjacent layers. and The initial component distribution is determined using a linear interpolation method, i.e., the first... Components of each sublayer Calculated as ,in The value range is 1 to This linear interpolation provides a reasonable initial component distribution, but it is not the optimal solution. The component gradient values ​​of each sublayer need to be adjusted as parameters to be optimized.

[0075] Based on the material composition parameters of each sublayer, the corresponding lattice constants can be calculated. For ternary compounds, the lattice constants typically follow Vegard's law, meaning the lattice constant of the mixture is a linear combination of the lattice constants of its components. System, lattice constant It can be represented as ,in The lattice constant of AlAs is approximately 0.56605 nm. The lattice constant of GaAs is approximately 0.56533 nm. After obtaining the lattice constant of each sublayer, the lattice mismatch strain of each sublayer is calculated by comparing it with the lattice constants of adjacent layers. Lattice mismatch strain of each sublayer relative to its adjacent layers Defined as ,in is the lattice constant of the adjacent layer.

[0076] Lattice mismatch strain generates elastic strain energy within a material, which is the primary driving force for dislocation nucleation and propagation. Substituting the lattice mismatch strain of each sublayer into the elastic strain energy formula allows calculation of the spatial distribution of strain energy density within the strain buffer layer. Elastic strain energy density It is proportional to the square of the strain, and can be expressed as ,in For effective Young's modulus, For lattice mismatch strain, The thickness is denoted as ...

[0077] Based on the spatial distribution of strain energy density and combined with the strain relaxation theory model, the updated dislocation line density at the interface after inserting a strain buffer layer can be calculated. Strain relaxation theory describes the nucleation, multiplication, and movement of dislocations driven by a given strain energy. Critical thickness theory gives the maximum thickness at which a material can maintain quasi-isostrain growth without relaxation under a specific lattice mismatch. When the actual thickness exceeds the critical thickness, some strain energy is released by introducing dislocations, thus increasing the dislocation line density. This is related to excess strain energy and the dislocation nucleation barrier of the material. By numerically solving the strain relaxation dynamics equation, the updated value of the dislocation linear density after introducing a strain buffer layer can be obtained.

[0078] In calculating the updated dislocation linear density, the structural parameters of the strain buffer layer need to be iteratively adjusted. If the initially designed strain buffer layer fails to reduce the dislocation linear density below the preset density threshold, the composition gradient distribution of the sublayers needs to be adjusted or the total thickness of the buffer layer needs to be increased. The adjustment of the composition gradient can be achieved using an optimization algorithm, with the dislocation linear density as the objective function and the composition values ​​of each sublayer as optimization variables, performing numerical optimization under the condition of satisfying physical constraints. Optimization algorithms can include gradient descent, genetic algorithms, or simulated annealing. Through multiple iterative calculations, when the updated dislocation linear density value first decreases below the preset density threshold, the optimization iteration stops, and the current total thickness of the strain buffer layer and the material composition distribution of each sublayer are recorded.

[0079] After obtaining the optimized strain buffer layer structural parameters, they need to be integrated into the original heterostructure parameters. The total thickness of the strain buffer layer, as well as the material composition distribution and thickness information of each sublayer, are added to the optimized geometric and material composition parameter data. Due to the introduction of the new strain buffer layer, the number of layers in the original heterostructure needs to be increased accordingly, and the spatial coordinates of each layer also need to be recalculated. Assume the original structure is at the interface... The position was inserted containing The strain buffer layer consists of several sub-layers, with a total thickness of [missing information]. Then the position coordinates of all layers after this interface need to be translated as a whole. The updated structural parameters include detailed information about the strain buffer layer, forming the final structural parameters after defect suppression. These parameters can be directly used to guide the epitaxial growth process to achieve the fabrication of low-defect-density semiconductor multilayer heterostructures.

[0080] A second aspect of the present invention provides a performance simulation and optimization system for semiconductor multilayer heterostructures, comprising: The parameter acquisition unit is used to acquire the initial geometric parameters, material composition distribution parameters, and target band gradient distribution of the semiconductor multilayer heterostructure. Based on the initial geometric parameters and material composition distribution parameters, the band discontinuity values ​​and carrier concentration distribution at each heterojunction interface are calculated by self-consistently solving the Poisson equation and the Schrödinger equation. The parameter optimization unit is used to adjust the thickness parameters and material composition parameters of each heterojunction layer in reverse according to the deviation between the band discontinuity value and the target band gradient distribution, so that the conduction band order and valence band order at the interface of each heterojunction meet the preset gradient continuity constraint conditions, and obtain the optimized geometric parameters and material composition parameters. The strain calculation unit is used to calculate the lattice mismatch strain distribution and dislocation line density distribution at each heterojunction interface based on the optimized geometric parameters and material composition parameters. The defect suppression unit is used to insert a strain buffer layer structure at the corresponding interface when there is an interface position in the dislocation line density distribution that exceeds a preset density threshold. By adjusting the thickness and composition gradient distribution of the strain buffer layer structure, the dislocation line density distribution is reduced to below the preset density threshold, and the final structural parameters after defect suppression are obtained. The final structural parameters after defect suppression are then output as the fabrication process parameters for the semiconductor multilayer heterostructure.

[0081] A third aspect of the present invention provides an electronic device, comprising: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0082] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0083] This invention can be a method, apparatus, system, and / or computer program product. The computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for performing various aspects of the invention.

[0084] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for performance simulation and optimization of semiconductor multilayer heterostructures, characterized in that, include: The initial geometric parameters, material composition distribution parameters, and target band gradient distribution of the semiconductor multilayer heterostructure are obtained. Based on the initial geometric parameters and material composition distribution parameters, the band discontinuity values ​​and carrier concentration distribution at each heterojunction interface are calculated using the Poisson equation and the Schrödinger equation in a self-consistent manner. Based on the deviation between the band discontinuity value and the target band gradient distribution, the thickness parameters and material composition parameters of each heterojunction are adjusted in reverse to ensure that the conduction band order and valence band order at the interface of each heterojunction satisfy the preset gradient continuity constraint condition, thereby obtaining the optimized geometric parameters and material composition parameters. Based on the optimized geometric parameters and material composition parameters, the lattice mismatch strain distribution and dislocation line density distribution at each heterojunction interface are calculated. When there is an interface location in the dislocation line density distribution that exceeds a preset density threshold, a strain buffer layer structure is inserted at the corresponding interface. By adjusting the thickness and composition gradient distribution of the strain buffer layer structure, the dislocation line density distribution is reduced to below the preset density threshold, and the final structural parameters after defect suppression are obtained. The final structural parameters after defect suppression are output as the fabrication process parameters of the semiconductor multilayer heterostructure.

2. The method according to claim 1, characterized in that, Based on the initial geometric parameters and material composition distribution parameters, the band discontinuities and carrier concentration distributions at each heterojunction interface are calculated using the Poisson and Schrödinger equations in a self-consistent manner, including: The spatial coordinate range of each heterojunction layer is determined based on the initial geometric parameters. The dielectric constant distribution and effective mass distribution of each layer are determined based on the material composition distribution parameters. The spatial region of each heterojunction layer is discretized into a grid. A difference scheme of the Poisson equation is established on the discrete grid nodes. The Poisson equation describes the relationship between electrostatic potential and carrier charge density. An eigenvalue problem of the Schrödinger equation is established on the discrete grid nodes. The potential energy term of the Schrödinger equation is composed of the electrostatic potential and the intrinsic band structure of the material. Solving the eigenvalue problem yields the wave function distribution and intrinsic energy of each quantum state. Based on the wave function distribution and intrinsic energy, the electron concentration and hole concentration at each grid node are calculated using the Fermi-Dirac distribution function. The calculated carrier concentration is substituted into the Poisson equation to update the electrostatic potential distribution. The Schrödinger equation and the Poisson equation are solved iteratively until the maximum difference between the electrostatic potential distributions obtained from two adjacent iterations is less than the preset convergence criterion, thus obtaining the self-consistent electrostatic potential distribution and carrier concentration distribution. Based on the electrostatic potential distribution obtained after self-consistent solution, the energy values ​​of the conduction band edge and the valence band edge on both sides of each heterojunction interface are extracted. The difference between the energy values ​​of the conduction band edge is calculated to obtain the conduction band discontinuity value, and the difference between the energy values ​​of the valence band edge is calculated to obtain the valence band discontinuity value.

3. The method according to claim 2, characterized in that, An eigenvalue problem of the Schrödinger equation is established on the discrete grid nodes. The potential energy term of the Schrödinger equation is composed of the electrostatic potential and the intrinsic band structure of the material. Solving the eigenvalue problem yields the wave function distribution and eigenenergy of each quantum state, including: The numerical distribution of the electrostatic potential is extracted on the discrete grid nodes, and the conduction band bottom energy and valence band top energy of the material corresponding to each grid node are read. The conduction band bottom energy and valence band top energy are used as characterization parameters of the intrinsic band structure of the material. The effective potential energy distribution of the electronic state is constructed by superimposing the numerical distribution of the electrostatic potential with the bottom energy of the conduction band, and the effective potential energy distribution of the hole state is constructed by superimposing the numerical distribution of the electrostatic potential with the top energy of the valence band. Based on the spatial coordinates of the discrete grid nodes, a discretized expression for the kinetic energy operator is constructed. The discretized expression for the kinetic energy operator uses a finite difference scheme to discretize the second spatial derivative of the wave function. Substituting the discretized expression of the kinetic energy operator and the effective potential energy distribution into the Schrödinger equation, an eigenvalue problem in matrix form is formed on the discrete grid nodes. The coefficient matrix of the eigenvalue problem is composed of kinetic energy terms and potential energy terms. The coefficient matrix of the eigenvalue problem is numerically solved to obtain the eigenvalues ​​and eigenvectors of each quantum state. The quantum states are sorted by energy level according to the eigenenergy, and the eigenenergy and wavefunction distribution of the electronic state near the bottom of the conduction band and the eigenenergy and wavefunction distribution of the hole state near the top of the valence band are extracted. The wave function distribution is normalized in space, and the spatial distribution of wave function probability density for each quantum state is calculated. The spatial distribution of wave function probability density is used to characterize the spatial localization characteristics of charge carriers in the heterojunction structure.

4. The method according to claim 1, characterized in that, Based on the deviation between the band discontinuity value and the target band gradient distribution, the thickness parameters and material composition parameters of each heterojunction layer are adjusted in reverse to ensure that the conduction band order and valence band order at each heterojunction interface satisfy the preset gradient continuity constraint condition, thereby obtaining the optimized geometric parameters and material composition parameters, including: Calculate the actual band gradient value at each heterojunction interface. The actual band gradient value is the ratio of the difference in band discontinuity between two adjacent interfaces to the interface spacing. Extract the target band gradient value corresponding to each interface position from the target band gradient distribution. Calculate the deviation vector between the actual band gradient value and the target band gradient value. A sensitivity matrix is ​​established between band gradient deviation and structural parameter adjustment. The elements of the sensitivity matrix represent the influence coefficients of unit thickness change or unit composition change on the band gradient of each interface. Based on the sensitivity matrix and the deviation vector, the structural parameter adjustment that minimizes the norm of the deviation vector is solved by the least squares method. The structural parameter adjustment includes the thickness adjustment and material composition adjustment of each layer. The structural parameter adjustment is superimposed on the current thickness parameter and material composition parameter to obtain the updated thickness parameter and material composition parameter. Based on the updated thickness parameter and material composition parameter, the band discontinuity value at each heterojunction interface is recalculated. When the conduction band order and valence band order satisfy the preset gradient continuity constraint condition, the current thickness parameter and material composition parameter are output as the optimized geometric parameter and material composition parameter.

5. The method according to claim 4, characterized in that, Establish a sensitivity matrix between band gradient bias and structural parameter adjustments, including: A set of structural parameters that need to be controlled is determined, the set of structural parameters including the thickness parameters of each heterojunction layer and the material composition parameters of each layer. A preset perturbation increment is applied to each structural parameter in the set of structural parameters, the perturbation increment being a preset percentage of the current value of the structural parameter. After applying the aforementioned perturbation increment, the band discontinuity values ​​at each heterojunction interface are recalculated using the self-consistent solution method of the Poisson equation and the Schrödinger equation, thus obtaining the distribution of the perturbed band discontinuity values. The actual band gradient value after the disturbance is calculated based on the distribution of the band discontinuities after the disturbance. The difference between the actual band gradient value after the disturbance and the actual band gradient value before the disturbance is obtained to obtain the change in band gradient. Calculate the ratio of the band gradient change to the perturbation increment to obtain the sensitivity coefficient of the structural parameter to the band gradient of each interface. Arrange the sensitivity coefficients according to the interface number and the structural parameter number to construct a sensitivity matrix. The row index of the sensitivity matrix corresponds to the interface position, the column index corresponds to the structural parameter, and the matrix element is the sensitivity coefficient of the corresponding structural parameter to the band gradient of the corresponding interface.

6. The method according to claim 1, characterized in that, Based on the optimized geometric parameters and material composition parameters, the lattice mismatch strain distribution and dislocation line density distribution at each heterojunction interface are calculated, including: Based on the optimized material composition parameters, the lattice constant and elastic modulus tensor of each heterojunction layer are obtained by querying the material database. The difference in lattice constants between the materials on both sides of each heterojunction interface is calculated. The ratio of the difference in lattice constants to the lattice constant of the substrate material is defined as the lattice mismatch degree. Based on the lattice mismatch and the elastic modulus tensor of each layer, the in-plane strain component and vertical strain component of each heterojunction layer are calculated using the elastic mechanical continuum model to obtain the strain tensor distribution at each heterojunction interface, and the strain tensor distribution is converted into the lattice mismatch strain distribution. Based on the lattice mismatch strain distribution and the optimized geometric parameters, the Matthews-Blakeslee critical thickness theory is used to determine whether each heterojunction interface exceeds the critical thickness. For heterojunction interfaces that exceed the critical thickness, the equilibrium spacing of mismatched dislocations is calculated based on the lattice mismatch strain distribution. The equilibrium spacing is jointly determined by the strain relaxation condition and the dislocation elastic energy balance condition. The reciprocal of the equilibrium spacing is used as the dislocation line density at the interface. The dislocation line density distribution is obtained by summing the dislocation line densities of all interfaces.

7. The method according to claim 1, characterized in that, A strain buffer layer structure is inserted at the corresponding interface. By adjusting the thickness and composition gradient distribution of the strain buffer layer structure, the dislocation line density distribution is reduced to below the preset density threshold, resulting in the final structural parameters after defect suppression, including: Identify the interface location in the dislocation line density distribution that exceeds a preset density threshold, insert a strain buffer layer structure at the interface location, set the initial thickness of the strain buffer layer structure to a preset ratio of the thickness of adjacent layers, and divide the strain buffer layer structure into several sub-layers along the growth direction. The material composition of each sublayer is set to vary in a piecewise manner between the material components on both sides of the interface. The initial component distribution is determined by linear interpolation, and the component gradient value of each sublayer is used as the parameter to be optimized. The corresponding lattice constant is calculated based on the material composition of each sublayer. The lattice mismatch strain of each sublayer is calculated based on the difference between the lattice constant and the lattice constant of the adjacent layer. The lattice mismatch strain of each sublayer is substituted into the elastic strain energy formula to calculate the spatial distribution of strain energy density inside the strain buffer layer. The strain energy density characterizes the driving force of dislocation nucleation. Based on the spatial distribution of strain energy density and strain relaxation theory, the dislocation line density update value at the interface position after inserting the strain buffer layer is calculated. When the dislocation line density update value decreases to below the preset density threshold, the current total thickness of the strain buffer layer and the material composition distribution of each sublayer are recorded. The total thickness of the strain buffer layer and the material composition distribution of each sublayer are integrated into the optimized geometric parameters and material composition parameters. The number of heterojunction layers and the position coordinates of each layer are updated to obtain the final structural parameters after defect suppression.

8. A semiconductor multilayer heterostructure performance simulation and optimization system, used to implement the method as described in any one of claims 1-N, characterized in that, include: The parameter acquisition unit is used to acquire the initial geometric parameters, material composition distribution parameters, and target band gradient distribution of the semiconductor multilayer heterostructure. Based on the initial geometric parameters and material composition distribution parameters, the band discontinuity values ​​and carrier concentration distribution at each heterojunction interface are calculated by self-consistently solving the Poisson equation and the Schrödinger equation. The parameter optimization unit is used to adjust the thickness parameters and material composition parameters of each heterojunction layer in reverse according to the deviation between the band discontinuity value and the target band gradient distribution, so that the conduction band order and valence band order at the interface of each heterojunction meet the preset gradient continuity constraint conditions, and obtain the optimized geometric parameters and material composition parameters. The strain calculation unit is used to calculate the lattice mismatch strain distribution and dislocation line density distribution at each heterojunction interface based on the optimized geometric parameters and material composition parameters. The defect suppression unit is used to insert a strain buffer layer structure at the corresponding interface when there is an interface position in the dislocation line density distribution that exceeds a preset density threshold. By adjusting the thickness and composition gradient distribution of the strain buffer layer structure, the dislocation line density distribution is reduced to below the preset density threshold, and the final structural parameters after defect suppression are obtained. The final structural parameters after defect suppression are then output as the fabrication process parameters for the semiconductor multilayer heterostructure.

9. An electronic device, characterized in that, include: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the method according to any one of claims 1 to N.

10. A computer-readable storage medium having computer program instructions stored thereon, characterized in that, When the computer program instructions are executed by the processor, they implement the method described in any one of claims 1 to N.

Images