Process parameter reverse generation method and system based on level set micro morphology evolution

Abstract

The present application belongs to the technical field of micro process design, and discloses a process parameter reverse generation method and system based on level set micro morphology evolution, which comprises the following steps: projecting a transient local normal vector field to a spherical coordinate system, and substituting a statistically generated micro facet normal distribution function into a rendering equation to generate a transient bidirectional scattering distribution matrix; constructing a time domain optical target evaluation function, and performing process locking based on minimum value tracking to obtain a standard judgment result and a process parameter formula; and based on the standard judgment result, performing reverse iterative optimization of the process parameter formula based on an outer loop to obtain a globally optimal physical manufacturing formula. The present application closes a loop from process parameters to stray light indicators, and can complete reverse calibration in a computer, thereby significantly saving wafer trial and error flow cost; and through time dimension optimization of optics and process, the nonlinear scattering characteristics of micro morphology in evolution are revealed and utilized, and the distortion of empirical fitting is abandoned.

Description

Technical Field

[0001] This invention relates to the field of micro-process design technology, and in particular to a method and system for reverse generation of process parameters based on the evolution of level set micro-morphology. Background Technology

[0002] In the development of modern high-end optical systems, such as AR / VR diffractive waveguides, automotive long-range lidar, and CMOS microlens arrays, wafer-level optics (WLO) has become mainstream as optical components become miniaturized and integrated. The roughness of optical surfaces at the micro- and nano-scale directly determines the system's stray light, haze, and contrast. Current research and development models forcibly combine two completely different theoretical frameworks with entirely different physical foundations between microscopic process simulation and macroscopic optical scattering analysis, specifically in the following two stages: Phase 1: Technology-Analog Processing (TCAD) based on fluid dynamics, where semiconductor process engineers use level sets to simulate morphology evolution at the microscale. The formation of microscopic surface morphology is rigorously modeled as a partial differential equation (PDE) of fluid dynamics and surface chemical reactions; the physical surface is rigorously defined as a set of isosurfaces where the function value is zero.

[0003] Phase Two: Optical scattering assessment based on phenomenological analytical models, i.e., non-sequential tracing. Optical engineers primarily rely on the two-way scattering distribution function (BSDF) when evaluating stray light in non-sequential tracing software. However, due to a lack of microscopic data, they can only use empirical phenomenological models or analytical models to statically define surface scattering. The most typical example is the ABg scattering model.

[0004] However, the aforementioned existing technologies suffer from significant drawbacks, including the physical mismatch and black-box dilemma inherent in traditional analytical models, as detailed below: (1) The collapse of the underlying physical assumptions, namely the contradiction between isotropy and anisotropy: Analytical models like ABg assume that surface roughness is completely random and isotropic. However, semiconductor processes exhibit strong directionality, producing periodic scalloping on the grating sidewalls or exposing specific crystal planes. This results in extremely strong anisotropic stray light distribution. Forcing an anisotropic semiconductor surface into a pre-defined, isotropic ABg model is a serious physical mismatch, leading to complete distortion of the simulation.

[0005] (2) The disconnect between trial-and-error chip fabrication and the lack of reverse engineering: Mathematical constants for empirical fitting ( ) and specific semiconductor equipment process parameters There are no physical equations connecting them. When optical designers provide stray light suppression specifications, foundries simply cannot translate them into specific equipment formulations. Current technology can only rely on massive DOE (Design of Experiments) for blind fabrication, and cannot discern transient optical extrema during the evolution process.

[0006] Therefore, how to provide a method and system for reverse generation of process parameters based on the evolution of micro-morphology of level sets is an urgent problem to be solved. Summary of the Invention

[0007] This invention provides a method and system for reverse generation of process parameters based on the evolution of microstructure of level sets, in order to solve the above-mentioned technical problems in the prior art.

[0008] According to a first aspect of the present invention, a method for reverse generation of process parameters based on the evolution of microstructure of level sets is provided.

[0009] In one embodiment, the method for reverse generation of process parameters based on the evolution of level set microstructure includes: The method for reverse generation of process parameters based on the evolution of micromorphology of level sets is characterized by including: The propagation velocity of the level set is decomposed, and the scalar field of the symbolic distance function is evolved and updated based on the decomposition results to obtain the updated scalar field. The first-order gradient of the space is calculated on the updated scalar field to obtain the transient local normal vector field. The transient local normal field is projected onto the spherical coordinate system, and the micro-surface element normal distribution function is statistically generated. The micro-surface element normal distribution function is then substituted into the rendering equation to generate the transient two-way scattering distribution matrix. A time-domain optical target evaluation function is constructed based on the transient bidirectional scattering distribution matrix. The time-domain optical target evaluation function is then used for process locking based on minimum value tracking to obtain the target process time. Based on the target process time, a global compliance determination and process parameter determination are performed to obtain the compliance determination result and process parameter formula. Based on the compliance determination results, the process parameter formula is optimized by reverse iteration based on the outer loop to obtain the globally optimal physical manufacturing formula.

[0010] According to a second aspect of the present invention, a reverse generation system for process parameters based on the evolution of level set microstructure is provided.

[0011] In one embodiment, the process parameter reverse generation system based on level set micromorphological evolution includes: The transient normal vector field extraction module is used to decompose the propagation velocity of the level set and evolve and update the scalar field of the sign distance function based on the decomposition results to obtain the updated scalar field. The first-order gradient of the space is calculated on the updated scalar field to obtain the transient local normal vector field. The transient bidirectional scattering distribution matrix construction module is used to project the transient local normal vector field onto the spherical coordinate system, statistically generate the micro-surface element normal distribution function, and substitute the micro-surface element normal distribution function into the rendering equation to generate the transient bidirectional scattering distribution matrix. The process reverse locking module is used to construct a time-domain optical target evaluation function based on the transient two-way scattering distribution matrix, perform process locking based on minimum value tracking on the time-domain optical target evaluation function to obtain the target process time, and perform global compliance judgment and process parameter determination based on the target process time to obtain the compliance judgment result and process parameter formula. The formula reverse iterative optimization module is used to perform reverse iterative optimization of the process parameter formula based on the compliance judgment result, so as to obtain the globally optimal physical manufacturing formula.

[0012] According to a third aspect of the present invention, a computer device is provided.

[0013] In some embodiments, the computer device includes a memory and a processor. The memory stores a computer program, and the processor executes the computer program to implement the steps of the above-described method for reverse generation of process parameters based on the evolution of level set microstructure.

[0014] According to a fourth aspect of the present invention, a computer-readable storage medium is provided.

[0015] In one embodiment, a computer program is stored on a computer-readable storage medium, and when the computer program is executed by a processor, it implements the steps of the above-described method for reverse generation of process parameters based on the evolution of microstructures of level sets.

[0016] The technical solutions provided by the embodiments of the present invention may include the following beneficial effects: 1. This invention reduces costs and increases efficiency at the physical level, eliminating blind testing in Design of Experiments (DOE), that is, it connects the process parameters at the algorithm level. p To stray light index J The closed loop can be reverse calibrated within a computer, significantly reducing wafer trial and error costs.

[0017] 2. This invention reveals and utilizes the core physical fact of the nonlinearity of scattering characteristics in the evolution of micro-morphology by optimizing the time dimension of optics and process; and obtains the golden time window that cannot be found by traditional static post-evaluation by tracking the time-domain extremum of the two-way scattering distribution function (BSDF);

[0018] 3. This invention completely eliminates the distortion of empirical fitting, that is, the optical evaluation is directly derived from the rigorous physical statistics of the normal gradient of the real process morphology, ensuring that the reverse-designed process formula has extremely high production line feasibility.

[0019] It should be understood that the above general description and the following detailed description are exemplary and explanatory only, and are not intended to limit the invention. Attached Figure Description

[0020] The accompanying drawings, which are incorporated in and form part of this specification, illustrate embodiments consistent with the invention and, together with the description, serve to explain the principles of the invention.

[0021] Figure 1 This is a flowchart illustrating a method for reverse generation of process parameters based on the evolution of microstructure of level sets, according to an exemplary embodiment. Figure 2 This is a schematic diagram of a process parameter reverse generation system based on the evolution of level set micromorphology, according to an exemplary embodiment. Figure 3 This is a schematic diagram of the structure of a computer device according to an exemplary embodiment; Figure 4 This is a flowchart illustrating the dynamic generation and reverse design of multi-scale topography and scattering properties according to an exemplary embodiment.

[0022] Figure label: 201. Transient normal vector field extraction module; 202. Transient two-way scattering distribution matrix construction module; 203. Process reverse locking module; 204. Reverse formulation iterative optimization module. Detailed Implementation

[0023] The following description and accompanying drawings fully illustrate specific embodiments described herein to enable those skilled in the art to practice them. Some portions and features of certain embodiments may be included in or replace portions and features of other embodiments. The scope of the embodiments herein includes the entire scope of the claims and all available equivalents thereof. The various embodiments described herein are presented in a progressive manner, with each embodiment focusing on its differences from other embodiments; similar or identical parts between embodiments can be referred to interchangeably.

[0024] The modules in the apparatus or system of this application can be implemented entirely or partially through software, hardware, or a combination thereof. These modules can be embedded in or independent of the processor in a computer device in hardware form, or stored in the memory of a computer device in software form, so that the processor can call and execute the operations corresponding to each module.

[0025] Where there is no conflict, the embodiments and features in the embodiments of the present invention can be combined with each other.

[0026] Figure 1 An embodiment of the method for reverse generation of process parameters based on the evolution of microstructure of level sets according to the present invention is shown.

[0027] In this optional embodiment, the method for reverse generation of process parameters based on the evolution of level set microstructure includes: Step S101: Decompose the propagation velocity of the level set, and evolve and update the scalar field of the symbolic distance function according to the decomposition result to obtain the updated scalar field. Calculate the first-order gradient of the space on the updated scalar field to obtain the transient local normal vector field. Step S102: Project the transient local normal field onto the spherical coordinate system, statistically generate the micro-surface element normal distribution function, and substitute the micro-surface element normal distribution function into the rendering equation to generate the transient two-way scattering distribution matrix; Step S103: Construct a time-domain optical target evaluation function based on the transient bidirectional scattering distribution matrix, perform process locking based on minimum value tracking on the time-domain optical target evaluation function to obtain the target process time, and perform global compliance judgment and process parameter determination based on the target process time to obtain compliance judgment results and process parameter formula; Step S104: Based on the compliance determination results, perform reverse iterative optimization of the process parameter formula based on the outer loop to obtain the globally optimal physical manufacturing formula.

[0028] In this optional embodiment, projecting the transient local normal field onto a spherical coordinate system, statistically generating a micro-surface element normal distribution function, and substituting the micro-surface element normal distribution function into the rendering equation to generate a transient bidirectional scattering distribution matrix includes: Normal vectors are extracted from the transient local normal vector field to obtain local physical normal vectors. These are then combined with preset macroscopic surface reference normal vectors for spherical coordinate projection to obtain micro-surface element angle parameter pairs. Based on the probability density statistics and solid angle correction of all micro-facet angle parameters, the normal distribution function of the micro-facet is obtained. Substitute the micro-surface element normal distribution function and the macro-surface reference normal vector into the rendering equation, and combine the Fresnel term and the geometric occlusion attenuation term to generate the transient bidirectional scattering distribution matrix.

[0029] In this optional embodiment, the transient local normal vector field is extracted to obtain the local physical normal vector, and then combined with the preset macroscopic surface reference normal vector for spherical coordinate projection to obtain the micro-element angle parameter pair, including: Define the macroscopic reference normal vector of the optical surface, and align the macroscopic reference normal vector to the positive direction of the vertical axis of the global coordinate system through coordinate rotation to obtain the macroscopic surface reference normal vector; Traverse every valid grid point in the narrow band of the zero isosurface in the updated scalar field, and use the central difference method to calculate and normalize the spatial first-order gradient of each valid grid point to obtain the local physical normal vector. By using inverse trigonometric functions, spherical coordinate analytical mapping is performed on each local physical normal vector and macroscopic surface reference normal vector to obtain the micro-element angle parameter pairs in spherical coordinates.

[0030] In this optional embodiment, probability density statistics and solid angle correction are performed on all micro-facet angle parameters to obtain the micro-facet normal distribution function, including: The effective angle space of spherical coordinates is divided into a two-dimensional angle discretization grid, and a two-dimensional accumulator matrix is ​​constructed. Based on all micro-facet angle parameter pairs, the two-dimensional accumulator matrix is ​​subjected to surface area weighted accumulation processing to obtain the updated two-dimensional accumulator matrix. The updated two-dimensional accumulator matrix is ​​then subjected to solid angle correction and normalization processing to obtain the probability density matrix. By using a two-dimensional Gaussian kernel to perform convolution smoothing on the probability density matrix, the normal distribution function of the micro-surface element is obtained.

[0031] In this optional embodiment, based on all micro-element angle parameter pairs, the two-dimensional accumulator matrix is ​​subjected to surface area weighted accumulation processing to obtain an updated two-dimensional accumulator matrix. Then, the updated two-dimensional accumulator matrix is ​​subjected to solid angle correction and normalization processing to obtain the probability density matrix, including: Based on all micro-facet angle parameter pairs and combined with the grid resolution in the two-dimensional accumulator matrix, calculate the two-dimensional grid index corresponding to all valid grid points within the narrow band; By using the surface area weighted accumulation method and combining the equivalent area weight of micro-facets at effective grid points, the matrix elements in the two-dimensional accumulator matrix are updated to obtain the updated two-dimensional accumulator matrix. Based on the solid angle differential formula, solid angle distortion correction is performed on the updated two-dimensional accumulator matrix to obtain the true density matrix; Based on a preset normalization constant, the true density matrix is ​​subjected to integral normalization based on energy conservation to obtain the probability density matrix.

[0032] In this optional embodiment, the expression for calculating the transient two-way scattering distribution matrix is: ; In the formula, Represents the transient two-way scattering distribution matrix; Represents the macroscopic incident ray vector; Represents the macroscopic outgoing ray vector; Represents the macroscopic surface reference normal vector; Represents the local half-angle vector of a micro-element; Represents the normal distribution function of a micro-element; Indicates the current time step in the process; Indicates Fresnel terms; This represents the geometric occlusion attenuation term.

[0033] In this optional embodiment, a time-domain optical target evaluation function is constructed based on the transient bidirectional scattering distribution matrix. The time-domain optical target evaluation function is then subjected to process locking based on minimum value tracking to obtain the target process time. Finally, a global compliance determination and process parameter determination are performed based on the target process time to obtain the compliance determination result and process parameter formula, including: In the stray light sensitive solid angle integral region, combined with the outgoing direction weight penalty function, the transient two-way scattering distribution function is weighted and integrated with respect to the solid angle corresponding to the outgoing direction to generate a time-domain optical target evaluation function. A time-domain sliding data window is constructed based on the objective function value calculated by the time-domain optical target evaluation function. The changing trend of the data in the time-domain sliding data window is judged, and the minimum value locking trigger result is determined based on the judgment result. The evolution engine is interrupted based on the trigger result, and the time point when the corresponding minimum value actually occurred is backtracked and the optimal time is locked to obtain the target process time. The smoothed trend value corresponding to the target process time is compared with the preset macroscopic optical shipment standard threshold to determine global compliance and obtain the compliance determination result. If the result of the standard assessment is that the optimization is successful, the process parameter formula will be output, and the process parameters will be set based on the process parameter formula and the target process time.

[0034] In this optional embodiment, a time-domain sliding data window is constructed based on the target function value calculated by the time-domain optical target evaluation function. The changing trend of the data within the time-domain sliding data window is judged, and the minimum value locking trigger result is determined based on the judgment result, including: By using a first-in-first-out queue of preset length, the target function values ​​of the time-domain optical target evaluation function are calculated for several time steps and stored to obtain a time-domain sliding data window; Causal low-pass filtering is used to smooth the data in the time-domain sliding data window to obtain smoothed trend values. The gradient of the optical quality index is obtained by performing discrete first-order gradient calculation on the smoothed trend value at adjacent time steps. Based on the tolerance mechanism for continuous deterioration, the gradient sign of the changing trend gradient is continuously monitored and judged, and the minimum value locking trigger result is determined according to the judgment result.

[0035] In this optional embodiment, based on the compliance determination result, the process parameter formula is optimized by reverse iterative optimization based on the outer loop to obtain the globally optimal physical manufacturing formula, including: If the compliance judgment result is unqualified, the process parameter search space is set and initial sampling is performed according to the process parameter formula to obtain several initial parameter formula samples, and an initial observation dataset is constructed based on several initial parameter formula samples. A Gaussian process regression surrogate model was trained using the initial observation dataset to obtain the predicted mean and uncertainty of the optical evaluation index corresponding to each initial parameter formulation sample. New candidate process formulations are identified based on the predicted mean and uncertainty, and new validation data points are generated based on the new candidate process formulations. The Gaussian process regression surrogate model is then iteratively trained using the new validation data points. After each iteration of training, check whether the termination condition is met. If the termination condition is met, extract the parameter set that minimizes the absolute value of the minimum optical inferiority from the evaluated initial observation dataset and output the globally optimal physical manufacturing formula.

[0036] In this optional embodiment, the process parameter search space is set and initial sampling is performed according to the process parameter formula to obtain several initial parameter formula samples, and an initial observation dataset is constructed based on the several initial parameter formula samples, including: The macroscopic process parameters in the process parameter formulation are used as outer-layer optimization variables, and the physical boundaries of the macroscopic process parameters are defined. A multidimensional parameter space is constructed based on the physical boundary. The Latin hypercube sampling method is used to perform spatial sampling in the multidimensional parameter space to generate several initial parameter formula samples. Several initial parameter formulation samples are input into the evolution engine, and the inner-layer extreme value tracking method is called to perform corresponding inner-layer optimal optics evaluation for each initial parameter formulation sample. The initial observation dataset is then constructed based on the evaluation results.

[0037] In this optional embodiment, the process of locking in a new candidate process formulation based on the predicted mean and uncertainty, generating new validation data points based on the new candidate process formulation, and iteratively training a Gaussian process regression surrogate model using the new validation data points includes: The predicted mean and uncertainty are input into the acquisition function, and new candidate process formulations are locked in by maximizing the acquisition function. The new candidate process formulation is sent to the physical evolution engine, and the new candidate process formulation is subjected to expensive inner-layer verification to obtain new verification data points. The newly validated data points were added to the initial observation dataset, and the Gaussian process regression surrogate model was retrained and corrected.

[0038] Figure 2An embodiment of the process parameter reverse generation system based on the evolution of level set micromorphology of the present invention is shown.

[0039] In this optional embodiment, the process parameter reverse generation system based on level set micromorphological evolution includes: The transient normal vector field extraction module 201 is used to decompose the propagation velocity of the level set, and to evolve and update the scalar field of the sign distance function according to the decomposition result to obtain the updated scalar field. The first-order gradient of the space is calculated on the updated scalar field to obtain the transient local normal vector field. The transient bidirectional scattering distribution matrix construction module 202 is used to project the transient local normal vector field onto the spherical coordinate system, statistically generate the micro-surface element normal distribution function, and substitute the micro-surface element normal distribution function into the rendering equation to generate the transient bidirectional scattering distribution matrix. The process reverse locking module 203 is used to construct a time-domain optical target evaluation function based on the transient two-way scattering distribution matrix, perform process locking based on minimum value tracking on the time-domain optical target evaluation function, obtain the target process time, and perform global compliance judgment and process parameter determination based on the target process time to obtain compliance judgment results and process parameter formula. The formula reverse iterative optimization module 204 is used to perform reverse iterative optimization of the process parameter formula based on the compliance judgment result to obtain the globally optimal physical manufacturing formula.

[0040] To facilitate understanding of the above technical solutions of the present invention, the following further explains the above technical solutions of the present invention from the perspective of architecture and principle, as follows: It should be noted that the partial differential equations (PDEs) in Phase 1, based on fluid dynamics-driven process evolution (TCAD), are as follows: ; in, Signed Distance Function (SDF): A scalar field function defined in a three-dimensional continuous space, whose absolute value represents the shortest geometric distance from a point in space to a real physical surface; its sign is used to distinguish material properties, for example, Inside the solid material, The external etching gas is used. The physical surface is strictly defined as the set of isosurfaces where the function value is zero. : Three-dimensional spatial coordinate vector, i.e. . Process evolution time, in seconds or minutes. Total Propagation Velocity: Characterized by the velocity of the material surface at a given process formulation at a given point in space. The growth or etching rate along the normal direction. : Process Parameter Set. A vector containing the underlying physical variables that control the operation of the equipment, such as etching gas flow rate ratio, cavity gas pressure, radio frequency (RF) bias power, etc. Spatial gradient operator. This represents the gradient magnitude of the SDF scalar field. According to the mathematical definition of SDF, its magnitude is always equal to 1 in an ideal case.

[0041] In Phase Two, the ABg scattering model, based on a phenomenological analytical model, is used to derive the classic ABg inverse power-law formula: ; in, Two-way scattering distribution function. Defined as the ratio of the emitted radiance in one direction to the incident radiance in the other direction. and : These are the cosine vectors of the projection directions of the scattered light and the specular reflected light onto the macroscopic surface, respectively. It is used to characterize the degree to which scattered light deviates from the ideal specular reflection direction. : Purely empirical fitting constants. To control the amplitude coefficient of the overall scattering intensity; To control the transition parameters in the roll-off region of the small-angle scattering curve; The slope exponent is used to control the rate of energy decay at large angle scattering.

[0042] This invention aims to overcome the isotropic assumptions of traditional optical analytical models and establish a bridge for inverse design, directly connecting macroscopic optical stray light indicators to microscopic process formulation parameters. Its fundamental logic is: at each time step of the level set evolution, the true transient physical normal is extracted directly from the underlying hydrodynamic evolution, reconstructing dynamic BSDF data containing real process imprints. By constructing an objective function characterizing optical degradation and continuously tracking the extreme points of this function in the time stream, the optimal process stop time and optimal formulation parameters can be precisely and inversely determined.

[0043] The formula derivation and reverse optimization mechanism of this invention specifically includes: Step 1: Anisotropic evolution and transient normal extraction based on level sets: In advancing the time step At that time, the propagation speed of the level set Detailed decomposition: ; in, : Rate term for isotropic chemical reactions. Characterized by events such as radical-dominated pure chemical corrosion, this rate is equal in magnitude in all directions. : Anisotropic ion bombardment rate term. Characterizes physical sputtering caused by a directional ion stream accelerated by an electric field in reactive ion etching (RIE), the angle between this rate and the surface normal and the ion stream direction. Strong correlation is the core physical source of anisotropic features such as scallop patterns. : Random perturbation noise term. Used to simulate the generation of random surface roughness caused by plasma concentration fluctuations or local microcrystalline defects in materials.

[0044] At each time step By directly calculating the first-order spatial gradient of the current implicit SDF scalar field, the transient local normal vector field of the zero isosurface (physical surface) can be obtained extremely quickly. : ; Step 2: Real-time construction of physical-level BSDF based on micro-surface element theory: normal vector within the narrow band Projecting onto a spherical coordinate system, a probability density matrix is ​​statistically generated, and after normalization, it is used as the normal distribution function of the micro-surface element. Subsequently, the transient BSDF matrix is ​​generated analytically by substituting the Cook-Torrance BRDF rendering equation: ; in, and : These are the macroscopic incident ray vector and the macroscopic outgoing (observed) ray vector, respectively. : Macroscopic surface reference normal vector, such as the vertical normal of an ideal plane. Half-vector of a micro-element. Mathematically defined as... Physically, this represents the ability to redirect light from Perfect mirror reflection The normal direction of a tiny surface. : Normal Distribution Function (NDF) of micro-area elements. The bridge variable in this invention is derived from the SDF field statistics of the previous stage of the process, and characterizes... At any moment, the normal direction points precisely. The probability area proportion of a tiny surface element in the overall surface. Fresnel term. Analytically calculated from the complex refractive index of the underlying material, it determines the proportion of energy reflected by light on the micro-surface element. : Geometry Masking / Shadowing Term. Corrects the energy attenuation caused by microscopic protrusions blocking incident light or preventing reflected light from being reflected.

[0045] The precise projection of the narrowband normal vector into the spherical coordinate system includes: reducing the dimension of the implicit field normal vector in the three-dimensional Cartesian coordinate system to a two-dimensional angle space, providing basic data for subsequent statistical analysis of the distribution of micro-surface element normals. The specific steps are as follows: 1.1 Define the macroscopic reference coordinate system, that is, set the macroscopic reference normal vector of the optical surface as... To simplify calculations, coordinate rotation is used to... Align with the positive Z-axis direction of the global coordinate system, that is: At this point, the macroscopic extension surface of the optical surface is parallel to the XOY plane of the global coordinate system.

[0046] 1.2 Extract and normalize the local normal vectors, i.e., traverse the narrow band of the zero isosurface in the level set evolution field, i.e., the signed distance function satisfies Each valid grid point in the grid point set k The spatial first-order gradient at this point is obtained using the central difference method. Then, normalization is performed to obtain the local physical normal vector. : ; in, The three orthogonal components of the local normal vector in Cartesian coordinates satisfy the normalization constraint. ; Narrowband half-width threshold, typically taken as 1-2 times the spatial grid size.

[0047] 1.3 Perform analytical mapping of spherical coordinates, that is, for each normalized local normal vector The zenith angle and azimuth angle parameters are mapped to spherical coordinates using inverse trigonometric functions. .

[0048] The mapping formula is as follows: ; .in, The angle between the local normal vector and the macroscopic reference normal (Z-axis) is used to characterize the tilt of the micro-element, with a theoretical range of [value missing]. If due to the chamfered structure ,Right now In reflection BSDF calculations, it can be removed or its normal can be flipped. The angle between the projection of the local normal vector onto the XOY plane and the X-axis is used to characterize the orientation of the micro-element. The bivariate arctangent function, compared to the conventional... According to , The sign of the output is automatically determined by the quadrant to ensure the correct azimuth angle. Accurately fall ,or A complete closed-loop interval is obtained to avoid quadrant singularity errors.

[0049] The specific statistics and generation of the probability density matrix (NDF) include: based on the angular parameters of all micro-facets. By using two-dimensional histogram statistics and solid angle correction, a discrete matrix of the micro-surface element normal distribution function (NDF) that satisfies physical conservation is generated. The specific steps are as follows: 2.1 Binning the two-dimensional angular discrete mesh involves dividing the effective angular space of the spherical coordinates into a mesh, constructing a mesh of size . Two-dimensional accumulator matrix The grid resolution is defined as follows: Zenith angle resolution: Corresponding grid index Azimuth resolution: Corresponding grid index .

[0050] 2.2 Surface area weighted accumulation (constructing a frequency histogram), i.e., traversing all valid points within the narrow band. Calculate its corresponding two-dimensional grid index: ; In the formula, This is the floor operator. The minimum value of the azimuth angle interval, for example, the interval is hour, The interval is hour, .

[0051] To achieve high physical fidelity, surface area weighted accumulation is used instead of traditional simple counting accumulation, and matrix elements are... The update rules are as follows: ; In the formula, Grid points The equivalent area weight of the micro-element at the location; in a mesh-based scalar field, it is usually taken as ,in The volume of a single voxel can accurately characterize the implicit physical surface area at that point.

[0052] 2.3 Solid Angle Correction is the core correction step in spherical coordinate statistics. On a sphere, near the pole ( The physical solid angle of the grid is much smaller than that near the equator ( The mesh needs to be based on the solid angle differential formula. For accumulator matrix Corrections are performed to eliminate statistical artifacts caused by pole convergence, yielding the true density matrix. : ; In the formula, : No. The center angle value of each zenith angle grid is calculated using the following formula: .

[0053] 2.4 Energy conservation and normalization generate the final probability density matrix, which satisfies the normal distribution function in the Cook-Torrance optical rendering equations. The law of conservation of energy states that the sum of the areas of all micro-elements projected onto the macroscopic surface equals the total macroscopic area. This requires applying the matrix... After rigorous integral normalization, the final probability density matrix elements are obtained. The calculation formula is: ; Wherein, the normalization constant Solve using the following formula: ; In the formula, The final output of the two-dimensional probability density matrix elements is the normal distribution function of the micro-surface elements. Discrete equivalent representation; : Area reduction factor of the projection of a micro-surface element onto a macroscopic physical surface.

[0054] 2.5 Gaussian smoothing filtering, i.e., the probability density matrix generated is affected by the discreteness of the micro-grid. There may be high-frequency numerical noise. This can be smoothed by convolving the matrix with a two-dimensional Gaussian kernel to ensure the continuity of the BSDF data and avoid non-physical spot noise during ray tracing. ; In the formula, Standard deviation is Two-dimensional Gaussian smoothing kernel, .

[0055] Step 3: Construction of the objective function and inverse locking of process time, i.e., extreme value tracking: Define specific optical quality indicators that users care about As the objective function. For example, to suppress a certain solid angle region. For stray light within the area, construct a weighted integral function: ; In the formula, : Temporal optical objective function. Represents the objective function evaluated at the stage of process advancement. At any given moment, the overall quality of the optical surface is compromised. Stray light sensitive solid angle integration region. For example, all directions of space outside 2° of the main beam. : Outgoing direction weighting penalty function. Used to apply a higher mathematical penalty weight to stray light affecting specific sensitive areas within the field of view, such as directly in front of AR glasses. : The fixed incident vector of the system's main working ray. :Exit direction The corresponding solid angle differential element.

[0056] Due to the trade-off between roughness generation and sidewall smoothing effect, The curve exhibits non-monotonic characteristics in the time domain. The system automatically searches for the optimal process stop time. ; Minimum pursuit of objective function and optimal process time (i.e., target process time) The specific locking mechanism includes: while the morphology is evolving, the system synchronously executes the following online extreme value determination and locking method in the background or a separate monitoring thread to determine the optimal process stop time. The specific steps are as follows: 3.1 Constructing a time-domain sliding data window, i.e., a sliding data window that moves with time steps. As the system progresses, it will not retain all historical data, but will instead maintain a fixed-length... First-In-First-Out (FIFO) queue Used to store the most recent The objective function evaluation results at each time step: ; This window can greatly save memory consumption and make extreme value determination depend only on the current local evolution trend.

[0057] 3.2 Applying Causal Low-pass Filtering to Suppress High-Frequency False Peaks: This addresses the issue of high-frequency spurious peaks caused by discrete transitions at the mesh interface or the introduction of Monte Carlo noise. The curve will inevitably be accompanied by high-frequency numerical oscillations (glitch). To prevent the algorithm from getting trapped in local false minima prematurely, the system uses a window... The data within the data is smoothed causally, for example, by using moving average filtering or Savitzky-Golay filtering, to obtain smoothed trend values. : ; in, To smooth the window length, .

[0058] 3.3 Calculation of Discrete First-Order Gradient (Trend Gradient Calculation): At each time step, the system calculates the first-order time difference of the smoothed objective function, i.e., the gradient, to quantify the rate and direction of change of the current optical quality indicators. ; Its judgment logic: when When the sidewall smoothing effect is dominant, stray light is improving; when When this occurs, it indicates that roughness or overall morphology degradation is dominant, and stray light is deteriorating.

[0059] 3.4 Confirmation of the true valley value based on the tolerance for continuous deterioration (Patience Window), i.e., in If evolution stops immediately the moment the value just exceeds 0, it might just be a slightly large fluctuation, possibly caused by errors or other factors. The system introduces a tolerance parameter for continuous deterioration. The system continuously monitors the gradient sign. Only when the gradient changes from negative to positive, and subsequently... Within a given time step, the gradient remains positive, meaning that the optical performance exhibits a deterministic and irreversible deterioration trend. Only then does the system formally trigger the minimum locking signal. ; At this point, the system determines that the true physical minimum, i.e., the bottom of the valley, has appeared, and its location is... Before taking the step.

[0060] 3.5 System Termination and Reverse output: Once the minimum value locking signal is triggered, the system immediately executes the following circuit breaker mechanism: Interrupt the TCAD evolution engine: stop consuming expensive PDE solution computing power. Backtrack and lock the optimal time: output the true minimum occurrence point as the target process time. ; Global compliance check: Check locked minimum values Is it below the preset macroscopic optical shipment standard threshold? .like Optimization successful, directly output current process parameter formula. and optimal processing time To the contract manufacturer. If This indicates that even the best-case scenario under the current process formulation is unsatisfactory. The system transmits this signal to an outer optimizer, such as a Bayesian optimization module, to modify the process formulation. For example, adjusting the etching gas ratio and clearing the window queue. This will restart a new round of evolutionary tracking.

[0061] Step 4: Outer Circulation and Global Process Formulation Reverse optimization includes: setting up a set of process parameters As an outer optimization variable. Solving the global optimization problem: ; The system's final output is the underlying physical manufacturing formula that maximizes macroscopic optical performance. .

[0062] Outer layer global process formulation ( The specific implementation of reverse optimization includes: the system constructs an outer layer process parameter set The independent variable is the minimum value locked in the inner layer. This is a black-box evaluation mapping for the response values. To find the globally optimal formulation within a minimal number of TCAD evolutions, the system employs a Bayesian optimization control algorithm based on a Gaussian process (GP). The specific execution steps are as follows: 4.1 Define the process parameter search space and initialization sampling, i.e., define the boundary: set the macroscopic process parameter set. The physical boundaries. For example, suppose... ,in That is, the etching gas ratio boundary is , RF power boundary is Spatial sampling: Using the Latin hypercube sampling (LHS) method, initial samples are uniformly generated in the multidimensional parameter space. For example, 5, parameter formula samples Initial evaluation: The initial formulations were input into the TCAD engine, and the inner-layer extreme value tracking method (i.e., the aforementioned sliding window method) was used to obtain the corresponding optimal optical evaluation results for the inner layer. Construct the initial observation dataset .

[0063] 4.2 Constructing a global surrogate model, i.e., utilizing the current dataset The system trains a Gaussian process regression model (GPR) in memory. The Gaussian process does not directly output a definite predicted value, but rather outputs the objective function across the entire parameter space. The probability distribution on the surface. For any unknown process formulation. The GPR model can predict the mean value of the optical evaluation index corresponding to the formulation. And uncertainty, i.e., standard deviation .

[0064] It uses a rapidly computed alternative model to replace time-consuming real TCAD simulations to predict the degree of optical degradation, i.e., error, that any combination of process parameters may produce.

[0065] A Gaussian process (GP) is a set of random variables, any finite number of which follow a joint normal distribution. In process formulation... p With optical indicators The construction and training steps of the mapping are as follows: 1. Define the prior distribution, which is the assumption of the objective function. It follows a Gaussian process with zero (or a constant) mean: ; Define the prior mean function .

[0066] 2. Selection of covariance function (kernel function), i.e., kernel function This study characterizes the similarities between different process formulations. To accommodate the nonlinear characteristics of morphology variations in semiconductor processes, this invention prioritizes the Matérn5 / 2 kernel function. ; In the formula, + The distance is Euclidean. The signal variance; This is the characteristic length scale.

[0067] 3. Construct the observation likelihood function, that is, take into account the numerical noise present in the TCAD simulation and BSDF extraction process. Actual observed values Compared with the true value The relationship is .

[0068] 4. Training process: Hyperparameter optimization, that is, the essence of training GPR is to find the dataset that best interprets the existing dataset. hyperparameters The specific steps are as follows: Construct the marginal log-likelihood (MLL) function: ; in, It is based on the current hyperparameter calculations. The covariance matrix of order. Perform optimization: use the L-BFGS-B algorithm to... Perform iterative optimization to maximize the MLL.

[0069] The specific prediction steps of the GPR model are as follows: After the model training is completed, for any set of unknown process formulations to be evaluated... GPR can use Bayesian inference to provide prediction results.

[0070] 1. Construct a joint distribution, which is based on the definition of a Gaussian process, using existing observations... The function value at the point to be predicted Follows a joint normal distribution: ; 2. Calculate the conditional posterior distribution, which is to derive the conditional posterior distribution using the conditional independence property of the multivariate normal distribution. In the given The posterior distribution is still a Gaussian distribution: .

[0071] 3. Calculate the predicted mean As an estimated value for optical indicators: ; This mean is a weighted linear combination of known experimental results, with the weights determined by the correlation between the new and old formulations.

[0072] 4. Calculate the uncertainty That is, the predictive standard deviation: ; It measures the reliability of the prediction. If Located near existing data points, Smaller; if It is in an unexplored parameter vacuum region. Larger. In the outer-layer optimization process, uncertainty can not only tell the system which formula looks good (i.e., with a small mean), but also which area is not yet clear (i.e., with a large variance), thereby driving the acquisition function to achieve a dynamic balance between exploration and utilization.

[0073] 4.3 Maximizing the Acquisition Function to Propose the Next Candidate Proposal, which is the core decision-making mechanism of the optimization engine. The system determines which process formulation to try next by maximizing an Acquisition Function (AF). This invention uses Expected Improvement (EI) or Upper Confidence Bound (UCB) as the acquisition function: ; The system's extremely fast search in mathematical space makes The biggest new formula .

[0074] Decision-making logic game (Exploitation vs. Exploration): Exploitation: The collection function tends to select the predicted mean. A formulation region with very low (i.e., good predicted optical performance). Exploration: The acquisition function also tends to select for uncertainties. A very high (i.e., the system has not yet explored this area and may be hiding better results) recipe area.

[0075] 4.4 Perform inner-layer expensive validation (Nested Evaluation), i.e., lock in new candidate formulations. Then, the system sends it to the underlying TCAD physics evolution engine: initiates level set topology evolution; dynamically extracts normals and generates temporal BSDF; performs online sliding window determination to lock the optimal process stop time under the new formulation. and the corresponding true optical inferiority minimum value .

[0076] 4.5 Posterior Update and Convergence Check, i.e., data update: updating the newly validated data points. Add to observation dataset In the process, the surrogate model is retrained and refined. With each iteration, the surrogate model's fit to the physical laws becomes increasingly accurate. Termination condition determination: After each iteration, the exit condition is checked: if found... The shipment volume is already below the preset extremely high standard threshold, for example, stray light energy. Alternatively, the number of outer optimization loops may have reached the preset maximum computing power budget limit.

[0077] 4.6 Output the final optimal output, which means that when the termination condition is met, the system extracts the optimal recipe from the evaluated dataset. The parameter set with the smallest absolute value is used to formally output the globally optimal physical manufacturing formula to the user: optimal machine bottom-level parameters: Strictly matched optimal process time: Expected optimal optical performance: .

[0078] The specific modifications to the Gaussian Process Regression (GPR) model, i.e., the posterior update steps, include: when the underlying TCAD engine completes a new candidate formulation... The evolutionary assessment was performed, and the true minimum value was returned. Then, the system triggers the posterior correction mechanism of the GPR model. This is to avoid recalculating from scratch every time new data is added. For matrix inversion with low complexity, this invention employs an incremental expansion algorithm based on Cholesky decomposition (Incremental CholeskyUpdate), the specific steps of which are as follows: 1. Dimensional augmentation of the observation dataset, which involves appending newly validated data points to the historical observation matrix. Assume the original dataset contains... There are samples, and the parameter matrix is... The observation vector is Amplified The order dataset is represented as: ; 2. Block Matrix Expansion: This method eliminates the need to recalculate the covariance between old data points. Instead, it uses a kernel function (Matérn 5 / 2) to calculate only the correlation vector between new and old data points. and the variance of the new data points themselves. The new covariance matrix It is constructed losslessly as a block matrix: ; In the formula, For the old The order covariance matrix; for The cross covariance vector of order 1; It is a scalar.

[0079] 3. Dimensionality Reduction and Accelerated Update Based on Cholesky Decomposition: In GPR prediction and likelihood assessment, the core bottleneck of computation is finding the inverse of the covariance matrix. This invention maintains the Cholesky lower triangular decomposition factors of the old matrix in memory at all times. That is, satisfying When the matrix is ​​expanded to At that time, the system directly obtains the new factorization factor through the following analytical derivation. : ; Using matrix multiplication equations, the system solves for the unknown vector extremely quickly through forward substitution. and scalar : ; ; This step reduces the computational time complexity of model correction from the standard Recalculate and optimize to dimensionality reduction In the process of optimizing massive amounts of process parameters, it greatly reduces the update latency of the proxy model.

[0080] 4. Hyperparameter Re-optimization: Since newly added data points may change the system's perception of global physical trends, the system utilizes the aforementioned rapidly obtained hyperparameters for dynamic re-optimization. The marginal log-likelihood function (MLL) is recalculated. The L-BFGS-B optimizer is invoked, starting with the old hyperparameters (Warm Start), and performs a small number of fine-tuning iterations to update the hyperparameters of the kernel function. If the newly added optical measurement parameters experience a sharp, non-linear decrease, the characteristic length-scale parameter of the kernel function will... The value will be automatically reduced, allowing the surrogate model to tolerate more dramatic local gradient changes in subsequent predictions.

[0081] 5. Refreshing the posterior probability field and entering the next round of the game, i.e., after the hyperparameter calibration is completed, the predicted mean surface of the entire Gaussian process. and uncertainty surface Instant refresh. Originally in The huge uncertainty collapses into the range of observation noise, and the system then drives the acquisition function, such as EI, to find the next target to explore.

[0082] like Figure 4As shown, this invention can be divided into three decoupled modules in its implementation. For example, it uses C++ as the underlying high-performance non-sequential ray tracing calculation engine, and seamlessly bridges it with a dynamic dynamics script written in Python through cross-language memory binding technology such as pybind11. The functions of each module are as follows: Microscopic process evolution engine (module A): receives parameters. Perform high-frequency time-domain iteration of the level set PDE and maintain Field data. Meshless scattering generator (Module B): Resides in memory. Extracts narrowband gradients, performs statistics, and generates data in real time using rendering equations. The matrix is ​​overwritten directly into the material property block of the underlying C++ optical engine using zero-copy data via shared tensor memory. The inverse evaluation and tracking controller (module C) receives BSDF data from module B and calculates the objective function. When an extreme value is detected, the evolution of interrupt module A is halted, and the parameters are updated. And then start the next iteration. Figure 4 Detailed Structural Description: The forward evolution of the microscopic level set morphology demonstrates the standard four steps of the level set algorithm: initializing the 3D SDF mesh, calculating the evolutionary velocity field, solving the PDE to update the SDF field. Arrows indicate process completion. The narrowband BSDF dynamic mapping path demonstrates how this invention achieves co-evolutionary coupling from physical evolution to real-time scattering parameter generation. Its branches include extracting narrowband mesh point data and interpolating to calculate the normal vector of the physical surface. By constructing a local BSDF matrix using micro-facet theory, the system automatically outputs real-time scattering property packets. The optical index extremum tracking and inverse calibration demonstrate the inverse design path of this invention. Stray light indices within a specific angular range are input. In the flow of time t Weighted integral function for continuous computation of BSDF J ( t ).if J ( t If the maximum value (or minimum value) is reached, the system will automatically reverse and lock the optimal process time t. opt And update the process parameter p. The reverse calibration closed loop is marked by checkmarks on the reverse inference design path, indicating the efficient path of the present invention. It abandons explicit mesh reconstruction and multiple reflections, and completely opens up the complete co-evolution and reverse design closed loop from physical evolution to real-time scattering parameter generation and reverse parameter calibration.

[0083] (1) Alternatives to the micromorphological evolution engine: In addition to using the level set method, this invention can also use the phase-field method or Monte Carlo atomic-level reaction simulation. The spatial gradient of the continuous order parameter in the phase-field method, or the discrete surface voxel model generated by Monte Carlo simulation, can also be used as the original input, extracted, and constructed as the normal distribution function (NDF), thereby generating... The theoretical background of the phase-field method: Unlike level sets, which trace definite boundaries, i.e., zero isosurfaces, the phase-field method introduces a continuous order parameter. To describe the state of a material. For example, to specify... Represents the solid phase, i.e., wafer material. The phase represents the etching gas, while the interface is a diffuse layer of minimal thickness with a smooth, continuous transition between 0 and 1. The phase-field method is more natural and stable than level sets when dealing with extremely complex topological fractures and mergers. Specific steps: 1. Construct the free energy functional and governing equations, i.e., define the total free energy of the system. This incorporates local chemical energy and interfacial gradient energy. Based on the principle of minimizing system free energy, we construct the Allen-Cahn equation, which governs morphological evolution, for non-conservative fields such as etching; or the Cahn-Hilliard equation, which governs conservative fields such as surface atomic migration smoothing. ; In the formula, The kinetic mobility coefficient is determined by specific process parameters. Decide; 1. To simulate the thermal noise term of process fluctuations. 2. Solve the partial differential equations to advance the evolution, i.e., discretize the above equations on a three-dimensional spatial grid and advance them over time to obtain the results at each time step. Continuous order parametric field 3. Direct extraction of implicit normal vectors, i.e., without generating explicit boundaries, the normal vectors of the physical surface are obtained. Can be directly derived from the order parameter The spatial gradient is precisely defined. The normal vector within the transition zone, i.e., the physical interface, is: ; Mapping to the BSDF generation flow, that is, the extracted normal vectors mentioned above The data is directly input into the two-dimensional angle probability statistics module in the main process of this invention to generate the corresponding NDF and dynamic BSDF.

[0084] The theoretical background of Monte Carlo atomic-level simulation: The Monte Carlo (MC) method simulates morphology evolution by tracking the random collision and reaction probabilities of massive discrete particles, such as ions, free radicals, polymer molecules, and discrete voxels on the wafer surface. Specific execution steps: 1. Construct a discrete voxel mesh and reaction probability model, that is, divide the three-dimensional simulation space into a unit cell-sized array of discrete voxels. Based on the underlying process formulation... 1. Define the reaction probabilities when different gas particles come into contact with surface voxels, such as adsorption probability, etching escape probability, and specular reflection probability. 2. Perform random particle ray tracing and evolution, i.e., at each time step... A massive number of random particles are launched from above the simulation domain. The particles travel in a straight line, and when they collide with a solid voxel, a random number generator (RNG) is invoked to compare the result with the aforementioned probability model, determining whether to erode the voxel (etching) or add a new voxel (deposition). 3. Surface smoothing and normal vector reconstruction: Since the microscopic surface generated by the MC method is a stepped surface composed of discrete cubes, i.e., voxels, the gradient cannot be directly calculated. The system uses Principal Component Analysis (PCA) or Local Plane Fitting at the interface. Centered on a surface voxel, the system takes the neighboring voxel point cloud within a specified search radius and fits it using the least squares method to obtain a local microplane. The normal vector of this microplane is the local normal vector. 4. Mapping to the BSDF generation stream involves taking the surface area weighted set of the extracted discrete normal vectors and feeding it into the two-dimensional angle probability statistics module to generate the corresponding NDF and dynamic BSDF.

[0085] The physical definition of the normal distribution function (NDF) is as follows: In microfacet theory, any rough optical surface is microscopically assumed to be composed of countless perfectly smooth microfacets, like mirrors. The physical meaning is: describing the precise direction of the local normal vector on a unit macroscopic surface area. Direction, i.e., the probability density of micro-element in the direction of a half-angle vector; or the proportion of the total area.

[0086] The specific structure and content of NDF in this invention: In traditional computer graphics or optical software, such as Zemax's ABg model, These are all analytical mathematical formulas based on arbitrary assumptions, such as Gaussian distribution, GGX distribution, or Beckmann distribution. These formulas assume that the roughness is absolutely symmetric and regular. However, in this invention, It is an empirical discrete matrix driven by real physical evolution data, and its specific contents include: 1. Data structure, namely a two-dimensional floating-point matrix. Row index Zenith angle corresponding to micro-element Column index Azimuth of the corresponding micro-element 2. Energy conservation constraint (Normalization Constraint): All elements in the matrix must satisfy strict mathematical normalization conditions, i.e., all micro-surface elements must be normalized to the macroscopic physical surface, i.e., the reference normal. The sum of the areas of the projections must equal 1: ; 3. The decisive role of NDF in optical calculations: when a virtual ray from a non-sequential ray tracing engine hits the surface, the system will determine the direction of the incident light based on the path of the ray. and observation direction Calculate the half-angle vector Subsequently, the system directly looks up the probability density value of the corresponding direction in the NDF matrix. . The larger the value, the more surface area there is that can just barely block light from... Reflected to .this When combined with Fresnel terms (reflectivity) and geometric occlusion terms, the measured stray light BSDF energy distribution is directly output.

[0087] The replacement of the outer-layer reverse optimization algorithm, that is, the outer-layer global process formulation. In the inverse optimization closed loop, in addition to using gradient-based adjoint state method or Bayesian optimization, heuristic global search algorithms such as genetic algorithms (GA) and simulated annealing can also be used.

[0088] Detailed implementation of the outer-layer global process formulation reverse optimization algorithm: In this invention, the core task of outer-layer optimization is to establish a set of microscopic semiconductor process parameters. (Input) and macroscopic optical stray light objective function The extreme value mapping relationship between (outputs).

[0089] I. Adjoint State Method Based on Gradient: For scenarios where the influence of process parameters on morphology has a continuous analytical property, this invention employs the adjoint state method to efficiently extract the accurate gradient of the objective function with respect to process parameters, avoiding the high computational cost of massive repetitive simulations required by the traditional finite difference method (FDM). Algorithm Input and Output Data Mapping: The input data is the initial set of microscopic process parameters. For example, the etching rate coefficient distribution, the initial level set symbol distance field The output data is the exact gradient vector of the objective function with respect to the process parameters. Updated optimal process parameters Algorithm architecture and execution steps: 1. Forward pass, i.e., input... Solve the level set evolution partial differential equation (PDE). , obtain from arrive The complete morphological history trajectory is obtained, and the final optical objective function is calculated. 2. Construction and Backward Pass of the Adjoint Equation: This involves constructing the adjoint partial differential equation (Adjoint PDE) corresponding to the forward PDE. Adjoint variables are introduced. , the objective function The variational function as the adjoint equation in The final value condition at time [time]. From Integrating in reverse time to Solve for the adjoint field 3. Gradient extraction and update, i.e., the sensitivity gradient of the objective function to process parameters is updated from the positive field. and accompanying field The inner product integral can be directly calculated: Subsequently, L-BFGS (Constrained Quasi-Newton Method) was used to update the process parameters along the negative gradient direction. .

[0090] II. Bayesian Optimization (BO): When process parameters contain discontinuous categorical variables, such as gas type or PDEs with strong nonlinearity making gradient analysis difficult, Gaussian process-driven Bayesian optimization is employed. Algorithm Structure: 1. Surrogate Model: This uses Gaussian Process Regression (GPR) as a black-box mapping architecture. GPR does not output a single predicted value but rather the posterior probability distribution in the parameter space. 2. Kernel Function and Connectivity: A Matérn 5 / 2 kernel function with Automatic Relevance Decision (ARD) is selected. ,in The weighted Euclidean distance between parameter features. This kernel function exhibits local second-order non-differentiability compared to the conventional RBF kernel allowance function, perfectly matching the non-smooth abrupt changes in microstructure caused by gas switching in the etching process. 3. Acquisition Function: This function uses Expected Improvement (EI) to activate the exploration mechanism, balancing the mean and variance. Training Process: 1. Training Data Source: 20 initial formulations are generated within the process parameter boundaries using Latin hypercube sampling (LHS). The data is fed into the TCAD engine to obtain the corresponding real optical specifications. This constitutes the initial training set. 2. Data preprocessing: Preprocessing process parameters, such as RF power. W, gas flow rate sccm performs min-max scaling to... Space, eliminating the influence of dimensions; for 3. Loss Function: The loss function for training the GPR model is defined as the Negative Marginal Log-Likelihood (NMLL). 4. Optimization Algorithm Selection: The L-BFGS-B algorithm is used to minimize the NMLL in order to solve for the hyperparameters of the Gaussian process, such as the length scale. and signal variance Parameter settings: 1. Hyperparameter adjustment range: The length-scale boundary of the Matérn kernel is limited to... The initial setting for observation noise variance is: It also allows for fine-tuning during training to absorb the high-frequency truncation error introduced by the TCAD Monte Carlo mesh. 2. Optimization strategy: Set the maximum expensive TCAD validation budget to 100 runs. For each new... And obtain the truth Then, update the dataset. And retrain GPR.

[0091] III. Implementation of Genetic Algorithms (GA): For scenarios with numerous local minima in the parameter space, a genetic algorithm is used for global jump-heuristic search. Algorithm Structure: 1. Chromosome Encoding: Real-parameter encoding is employed. Machine physical parameters, such as etching duration, gas ratio, and electrode voltage, are directly mapped to chromosome gene sequences. This avoids the Hamming Cliff problem when optimizing binary encoding in a continuous space. 2. Fitness Function: Optical objective function The smaller the better, therefore the fitness function is defined as follows: ,in Preventing division by zero. 3. Genetic operator connections: 3.1 Selection layer: Tournament Selection; 3.2 Crossover layer: Simulated Binary Crossover (SBX), preserving the continuity of superior genes from the parents; 3.3 Mutation layer: Polynomial Mutation, responsible for escaping local minima. Training (evolution) process details: 1. Data sources and calculations: Initialized scale is... A random population. The genes of each individual are decoded into a process formula and input into the TCAD-optics co-simulation engine to calculate the corresponding stray light integral in parallel. This yields the population fitness distribution. 2. Evolutionary Iteration Mechanism: Individuals with low fitness (i.e., high stray light) are eliminated, and high-fitness individuals are used to generate a new generation of population through SBX and mutation. An elite retention strategy is enforced to ensure that the best process formulation ever observed unconditionally enters the next generation. Parameter Settings: 1. Population Size: Selection criteria: In terms of coverage Within the dimensional process parameter space, it is limited by the number of cores in the TCAD parallel computing cluster. 2. Crossover rate ( ) and the rate of variation ( ): Set it to 0.8 to ensure convergence speed; Set as , For the gene dimension, ensure that each mutation perturbs an average of one process parameter. The distribution exponent for SBX is set to 20, and the distribution exponent for multinomial mutation is set to 20.

[0092] IV. Simulated Annealing (SA) Implementation: For high-dimensional discrete process spaces, such as optical thin films with 1-20 layers and SiO2 or Ta2O5 materials, the target space exhibits a highly stepped and discontinuous structure, making traditional optimization prone to failure. Therefore, the simulated annealing algorithm is employed. Algorithm Structure: 1. State Transition Architecture (Markov Chain): The algorithm state represents the current process formulation. System energy Strictly mapped to the optical penalty function under this formulation 2. Neighborhood Search Network: Generates new formulations by applying random perturbations near the current process parameters. For example, randomly changing the thickness or material of a thin film. 3. Activation Criterion (Metropolis Criterion): The receiver probability function is defined as... If the new formulation leads to a reduction in stray light, that is... probability 100% acceptable; if stray light deteriorates, then... Then, based on probability Accepting this inferior solution grants the algorithm the ability to overcome local barriers. Training annealing process: 1. Preprocessing: Starting from the current empirical recipe (Base Recipe) as the initial solution, calculate the initial system energy. 2. Iterative calculation: at the current temperature Next, execute Sub-Markov chain Monte Carlo (MCMC) perturbation. Each perturbation invokes an optical simulation engine to calculate the new energy and updates the current optimal formula using the Metropolis criterion. 3. Cooling Operation: After completing the isothermal iteration, reduce the system temperature according to the cooling schedule. This limits the probability of the system accepting inferior formulas. Parameter settings: 1. Initial temperature ( ): The design is based on ensuring that the probability of accepting substandard formulas in the initial stage reaches over 80%. Average energy fluctuations are typically assessed through trial random walks. ,set up 2. Cooling Schedule: Employs an exponential cooling strategy. The attenuation coefficient This ensures the annealing process is slow enough to approach the physically lowest global energy state, i.e., optimal optical performance. 3. Number of inner cycles. : Set as parameter dimension It is 10 times larger, ensuring that the neighborhood space is fully explored under each temperature gradient.

[0093] (2) Substitution of optical target evaluation index, i.e., the objective function of reverse design Besides the stray light integral at a specific solid angle, it can also be replaced by any quantitative optical metric derived from BSDF. For example, Total Integrated Scatter (TIS), the modulation transfer function (MTF) drop-off of the imaging system, or the signal-to-noise ratio (SNR) at a specific diffraction order.

[0094] This invention relates to dynamic BSDF generation and evaluation based on implicit scalar field evolution. During the level set evolution of the process microstructure, gradient differentiation is continuously and directly performed on the zero isosurface neighborhood of the scalar field to extract transient physical normals. Two-dimensional probability histograms are then used to statistically generate dynamic BSDF data matrices corresponding to each time step in real time. A process inverse design mechanism based on time-domain BSDF objective function extremum tracking is employed, defining the time-domain objective function according to macroscopic optical indicators such as stray light or haze. The system monitors the minimum inflection point of the objective function in the morphology-time evolution sequence and automatically reverse-engineers the optimal process operation time, such as the optimal etching stop point. Based on a two-layer optimization loop for closed-loop optimization of the underlying manufacturing parameters, an automated optimization architecture with inner and outer layers is constructed. The inner layer tracks the temporal extrema under a specific formulation, while the outer layer uses underlying manufacturing parameters, such as gas ratios and cycle time, as variables to solve for the global optimal bound of the objective function through an optimization algorithm, thereby directly outputting a semiconductor manufacturing formulation that satisfies macroscopic optics requirements.

[0095] This invention is applicable to 3D effect simulation of advanced photomasks, etching process tolerance analysis of slanted gratings in AR / VR waveguide glasses, evaluation of metasurface fabrication morphology degradation, and joint optimization of deep trench isolation process and optical crosstalk in CMOS image sensors (CIS).

[0096] In one embodiment, a computer device is provided, which may be a server, and its internal structure diagram may be as follows: Figure 3As shown, the computer device includes a processor, memory, and a network interface connected via a system bus. The processor provides computing and control capabilities. The memory includes a non-volatile storage medium and internal memory. The non-volatile storage medium stores an operating system, computer programs, and a database. The internal memory provides an environment for the operation of the operating system and computer programs in the non-volatile storage medium. The database stores static and dynamic information data. The network interface communicates with external terminals via a network connection. When the computer program is executed by the processor, it implements the steps in the above-described embodiment of the reverse generation method for process parameters based on the evolution of level set microstructures.

[0097] Those skilled in the art will understand that Figure 3 The structure shown is merely a block diagram of a portion of the structure related to the present invention and does not constitute a limitation on the computer device to which the present invention is applied. A specific computer device may include more or fewer components than those shown in the figure, or combine certain components, or have different component arrangements.

[0098] In addition, the present invention also provides a computer device, including a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement the steps in the above-described embodiments of the method for reverse generation of process parameters based on the evolution of microstructures of level sets.

[0099] In addition, the present invention also provides a computer-readable storage medium storing a computer program thereon, which, when executed by a processor, implements the steps in the above-described embodiments of the method for reverse generation of process parameters based on the evolution of microstructures of level sets.

[0100] Those skilled in the art will understand that implementing all or part of the processes in the reverse generation method for process parameters based on the evolution of level set microstructures in the above embodiments can be accomplished by a computer program instructing related hardware. The computer program can be stored in a non-volatile computer-readable storage medium. When executed, the computer program can include the processes of the embodiments of the reverse generation method for process parameters based on the evolution of level set microstructures described above. Any references to memory, storage, databases, or other media used in the embodiments provided by this invention can include at least non-volatile and volatile memory. Non-volatile memory can include read-only memory (ROM), magnetic tape, floppy disk, flash memory, or optical storage, etc. Volatile memory can include random access memory (RAM) or external cache memory. By way of illustration and not limitation, RAM can be in various forms, such as static random access memory (SRAM) or dynamic random access memory (DRAM), etc.

[0101] This invention is not limited to the structures described above and shown in the accompanying drawings, and various modifications and changes can be made without departing from its scope. The scope of this invention is limited only by the appended claims.

Claims

1. A method for reverse generation of process parameters based on the evolution of microstructure of level sets, characterized in that, include: The propagation velocity of the level set is decomposed, and the scalar field of the symbolic distance function is evolved and updated based on the decomposition results to obtain the updated scalar field. The first-order gradient of the space is calculated on the updated scalar field to obtain the transient local normal vector field. The transient local normal field is projected onto the spherical coordinate system, and the micro-surface element normal distribution function is statistically generated. The micro-surface element normal distribution function is then substituted into the rendering equation to generate the transient two-way scattering distribution matrix. A time-domain optical target evaluation function is constructed based on the transient bidirectional scattering distribution matrix. The time-domain optical target evaluation function is then used for process locking based on minimum value tracking to obtain the target process time. Based on the target process time, a global compliance determination and process parameter determination are performed to obtain the compliance determination result and process parameter formula. Based on the compliance determination results, the process parameter formula is optimized by reverse iteration based on the outer loop to obtain the globally optimal physical manufacturing formula.

2. The method for reverse generation of process parameters based on the evolution of microstructure of level sets according to claim 1, characterized in that, The step of projecting the transient local normal vector field onto a spherical coordinate system, statistically generating a micro-surface element normal distribution function, and substituting the micro-surface element normal distribution function into the rendering equation to generate a transient bidirectional scattering distribution matrix includes: Normal vectors are extracted from the transient local normal vector field to obtain local physical normal vectors. These are then combined with preset macroscopic surface reference normal vectors for spherical coordinate projection to obtain micro-surface element angle parameter pairs. Based on the probability density statistics and solid angle correction of all micro-facet angle parameters, the normal distribution function of the micro-facet is obtained. Substitute the micro-surface element normal distribution function and the macro-surface reference normal vector into the rendering equation, and combine the Fresnel term and the geometric occlusion attenuation term to generate the transient bidirectional scattering distribution matrix.

3. The method for reverse generation of process parameters based on the evolution of microstructure of level sets according to claim 2, characterized in that, The process of extracting the normal vector from the transient local normal field to obtain the local physical normal vector, and then combining it with the preset macroscopic surface reference normal vector to perform spherical coordinate projection to obtain the micro-element angle parameter pair includes: Define the macroscopic reference normal vector of the optical surface, and align the macroscopic reference normal vector to the positive direction of the vertical axis of the global coordinate system through coordinate rotation to obtain the macroscopic surface reference normal vector; Traverse every valid grid point in the narrow band of the zero isosurface in the updated scalar field, and use the central difference method to calculate and normalize the spatial first-order gradient of each valid grid point to obtain the local physical normal vector. By using inverse trigonometric functions, spherical coordinate analytical mapping is performed on each local physical normal vector and macroscopic surface reference normal vector to obtain the micro-element angle parameter pairs in spherical coordinates.

4. The method for reverse generation of process parameters based on the evolution of microstructure of level sets according to claim 2, characterized in that, The step of performing probability density statistics and solid angle correction on all micro-facet angle parameters to obtain the micro-facet normal distribution function includes: The effective angle space of spherical coordinates is divided into a two-dimensional angle discretization grid, and a two-dimensional accumulator matrix is ​​constructed. Based on all micro-facet angle parameter pairs, the two-dimensional accumulator matrix is ​​subjected to surface area weighted accumulation processing to obtain the updated two-dimensional accumulator matrix. The updated two-dimensional accumulator matrix is ​​then subjected to solid angle correction and normalization processing to obtain the probability density matrix. By using a two-dimensional Gaussian kernel to perform convolution smoothing on the probability density matrix, the normal distribution function of the micro-surface element is obtained.

5. The method for reverse generation of process parameters based on the evolution of microstructure of level sets according to claim 4, characterized in that, The process involves performing surface area-weighted accumulation on the two-dimensional accumulator matrix based on all micro-element angle parameter pairs to obtain an updated two-dimensional accumulator matrix. Then, the updated two-dimensional accumulator matrix undergoes solid angle correction and normalization to obtain the probability density matrix, which includes: Based on all micro-facet angle parameter pairs and combined with the grid resolution in the two-dimensional accumulator matrix, calculate the two-dimensional grid index corresponding to all valid grid points within the narrow band; By using the surface area weighted accumulation method and combining the equivalent area weight of micro-facets at effective grid points, the matrix elements in the two-dimensional accumulator matrix are updated to obtain the updated two-dimensional accumulator matrix. Based on the solid angle differential formula, solid angle distortion correction is performed on the updated two-dimensional accumulator matrix to obtain the true density matrix; Based on a preset normalization constant, the true density matrix is ​​subjected to integral normalization based on energy conservation to obtain the probability density matrix.

6. The method for reverse generation of process parameters based on the evolution of microstructure of level sets according to claim 2, characterized in that, The expression for calculating the transient two-way scattering distribution matrix is ​​as follows: ； In the formula, Represents the transient two-way scattering distribution matrix; Represents the macroscopic incident ray vector; Represents the macroscopic outgoing ray vector; Represents the macroscopic surface reference normal vector; Represents the local half-angle vector of a micro-element; Represents the normal distribution function of a micro-element; Indicates the current time step in the process; Indicates Fresnel terms; This represents the geometric occlusion attenuation term.

7. The method for reverse generation of process parameters based on the evolution of microstructure of level sets according to claim 1, characterized in that, The process involves constructing a time-domain optical target evaluation function based on the transient bidirectional scattering distribution matrix, performing process locking based on minimum value tracking on the time-domain optical target evaluation function to obtain the target process time, and then performing global compliance determination and process parameter determination based on the target process time to obtain the compliance determination result and process parameter formula, including: In the stray light sensitive solid angle integral region, combined with the outgoing direction weight penalty function, the transient two-way scattering distribution function is weighted and integrated with respect to the solid angle corresponding to the outgoing direction to generate a time-domain optical target evaluation function. A time-domain sliding data window is constructed based on the objective function value calculated by the time-domain optical target evaluation function. The changing trend of the data in the time-domain sliding data window is judged, and the minimum value locking trigger result is determined based on the judgment result. The evolution engine is interrupted based on the trigger result, and the time point when the corresponding minimum value actually occurred is backtracked and the optimal time is locked to obtain the target process time. The smoothed trend value corresponding to the target process time is compared with the preset macroscopic optical shipment standard threshold to determine global compliance and obtain the compliance determination result. If the result of the standard assessment is that the optimization is successful, the process parameter formula will be output, and the process parameters will be set based on the process parameter formula and the target process time.

8. The method for reverse generation of process parameters based on the evolution of microstructure of level sets according to claim 7, characterized in that, The process of constructing a time-domain sliding data window based on the target function value calculated from the time-domain optical target evaluation function, judging the changing trend of the data in the time-domain sliding data window, and determining the minimum value locking trigger result based on the judgment result includes: By using a first-in-first-out queue of preset length, the target function values ​​of the time-domain optical target evaluation function are calculated for several time steps and stored to obtain a time-domain sliding data window; Causal low-pass filtering is used to smooth the data in the time-domain sliding data window to obtain smoothed trend values. The gradient of the optical quality index is obtained by performing discrete first-order gradient calculation on the smoothed trend value at adjacent time steps. Based on the tolerance mechanism for continuous deterioration, the gradient sign of the changing trend gradient is continuously monitored and judged, and the minimum value locking trigger result is determined according to the judgment result.

9. The method for reverse generation of process parameters based on the evolution of microstructure of level sets according to claim 1, characterized in that, The process parameter formulation is optimized through reverse iterative optimization based on the outer loop, based on the compliance determination results, to obtain the globally optimal physical manufacturing formulation, including: If the compliance judgment result is unqualified, the process parameter search space is set and initial sampling is performed according to the process parameter formula to obtain several initial parameter formula samples, and an initial observation dataset is constructed based on several initial parameter formula samples. A Gaussian process regression surrogate model was trained using the initial observation dataset to obtain the predicted mean and uncertainty of the optical evaluation index corresponding to each initial parameter formulation sample. New candidate process formulations are identified based on the predicted mean and uncertainty, and new validation data points are generated based on the new candidate process formulations. The Gaussian process regression surrogate model is then iteratively trained using the new validation data points. After each iteration of training, check whether the termination condition is met. If the termination condition is met, extract the parameter set that minimizes the absolute value of the minimum optical inferiority from the evaluated initial observation dataset and output the globally optimal physical manufacturing formula.

10. The method for reverse generation of process parameters based on the evolution of microstructure of level sets according to claim 9, characterized in that, The process parameter search space setting and initial sampling process based on the process parameter formula to obtain several initial parameter formula samples, and the construction of an initial observation dataset based on these initial parameter formula samples, includes: The macroscopic process parameters in the process parameter formulation are used as outer-layer optimization variables, and the physical boundaries of the macroscopic process parameters are defined. A multidimensional parameter space is constructed based on the physical boundary. The Latin hypercube sampling method is used to perform spatial sampling in the multidimensional parameter space to generate several initial parameter formula samples. Several initial parameter formulation samples are input into the evolution engine, and the inner-layer extreme value tracking method is called to perform corresponding inner-layer optimal optics evaluation for each initial parameter formulation sample. The initial observation dataset is then constructed based on the evaluation results.

11. The method for reverse generation of process parameters based on the evolution of microstructure of level sets according to claim 9, characterized in that, The process of identifying new candidate process formulations based on predicted mean and uncertainty, generating new validation data points based on these new candidate formulations, and iteratively training a Gaussian process regression surrogate model using these new validation data points includes: The predicted mean and uncertainty are input into the acquisition function, and new candidate process formulations are locked in by maximizing the acquisition function. The new candidate process formulation is sent to the physical evolution engine, and the new candidate process formulation is subjected to expensive inner-layer verification to obtain new verification data points. The newly validated data points were added to the initial observation dataset, and the Gaussian process regression surrogate model was retrained and corrected.

12. A process parameter reverse generation system based on the evolution of level set microstructure, characterized in that, The process parameter reverse generation system based on the evolution of level set microstructure includes: The transient normal vector field extraction module is used to decompose the propagation velocity of the level set and evolve and update the scalar field of the sign distance function based on the decomposition results to obtain the updated scalar field. The first-order gradient of the space is calculated on the updated scalar field to obtain the transient local normal vector field. The transient bidirectional scattering distribution matrix construction module is used to project the transient local normal vector field onto the spherical coordinate system, statistically generate the micro-surface element normal distribution function, and substitute the micro-surface element normal distribution function into the rendering equation to generate the transient bidirectional scattering distribution matrix. The process reverse locking module is used to construct a time-domain optical target evaluation function based on the transient two-way scattering distribution matrix, perform process locking based on minimum value tracking on the time-domain optical target evaluation function to obtain the target process time, and perform global compliance judgment and process parameter determination based on the target process time to obtain the compliance judgment result and process parameter formula. The formula reverse iterative optimization module is used to perform reverse iterative optimization of the process parameter formula based on the compliance judgment result, so as to obtain the globally optimal physical manufacturing formula.

13. A computer device comprising a memory and a processor, wherein the memory stores a computer program, characterized in that, When the processor executes the computer program, it implements the steps of the method for reverse generation of process parameters based on the evolution of microstructure of level sets as described in any one of claims 1 to 11.

14. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, it implements the steps of the method for reverse generation of process parameters based on the evolution of micromorphology of level sets as described in any one of claims 1 to 11.

Images