Glass substrate rdl multi-field coupling prediction method with fused physical constraint neural operator

Abstract

The application discloses a glass substrate RDL multi-field coupling prediction method based on a thermal-force-fracture coupling neural operator model (TMFO), a force-optical mapping neural operator model (SOMO) and an electrodeposition multi-field coupling neural operator model (EMDO) cascade architecture, and realizes field quantity lossless transmission through a unified data interface. The TMFO embeds the heat conduction equation, the elastic mechanics balance equation and the exponential cohesive force model residual constraint, and outputs the residual stress field and the interface damage field; the SOMO adopts a U-Net and FNO hybrid architecture, applies the pupil function and the Hopkins partial coherence imaging constraint, and maps the stress field into the defocus field and the light intensity distribution; the EMDO processes the moving boundary based on the space-time Fourier neural operator, embeds the Butler-Volmer equation, the Langmuir adsorption kinetics and the Hammerstad-Jensen loss model constraint, and predicts the copper filling morphology, the roughness and the high-frequency loss. The simulation data pre-training and the MMD domain adaptive fine-tuning strategy are adopted, so that the second-level end-to-end prediction from the process parameters to the key quality indicators is realized.

Description

Technical Field

[0001] This invention relates to the technical field of semiconductor advanced packaging process simulation and artificial intelligence intersection, and in particular to a multi-field coupling prediction method for glass substrate redistribution layer (RDL) process that integrates physically constrained neural operators. Background Technology

[0002] With the rapid development of artificial intelligence, high-performance computing, and 5G communication technologies, semiconductor packaging technology is evolving towards higher integration, higher I / O density, and higher signal transmission rates. Glass substrates, due to their excellent dimensional stability, low coefficient of thermal expansion, and high-frequency signal transmission performance, have become the core material for next-generation advanced packaging substrates. The redistribution layer (RDL) fabricated on the surface of the glass substrate is a key structure for rewiring signals between the chip and package pins; its process quality directly determines the electrical performance and reliability of the package.

[0003] The manufacturing process of large-size glass substrate RDLs involves strong coupling effects across multiple physical domains, including thermo-mechanical-fracture, mechanical-optical, and electrochemical deposition. During the thermal cycling process, the mismatch in thermal expansion coefficients between the glass substrate and the copper interconnect layer generates non-uniform residual stress at the interface, inducing substrate warping and interface damage. This residual stress-induced substrate warping further causes non-uniform defocusing within the photolithographic exposure surface, affecting critical dimension (CD) accuracy. During electrochemical deposition, the transport of copper ions within the trenches, electrochemical reaction kinetics, and additive adsorption behavior collectively determine the filling morphology, while the surface roughness of the filling morphology directly affects high-frequency signal transmission loss. Therefore, the manufacturing quality of RDLs is determined by the synergistic effects of these multiple physical domains.

[0004] In existing technologies, the prediction of multi-field coupling problems in RDL processes has the following main shortcomings:

[0005] First, traditional numerical simulation methods are computationally expensive. When dealing with full-size multiphysics coupled models of large-size glass substrates, numerical simulation methods such as finite element analysis (FEA), finite difference method (FDM), and computational fluid dynamics (CFD) typically have tens of millions to hundreds of millions of degrees of freedom. The computation time for a single thermo-mechanical coupling simulation can range from several hours to tens of hours, and the joint solution cycle across physical domains can reach several days, which cannot meet the needs of rapid process iteration and online real-time prediction.

[0006] Second, the simulation tools for each physical domain are fragmented. Thermo-mechanical analysis, optical simulation, and electrochemical simulation each use different dedicated software platforms, and the data interfaces between the tools are not unified. The transfer of field quantities across physical domains relies on manual data conversion and simplification assumptions, resulting in the loss of strong coupling information between physical domains and making it impossible to establish end-to-end quantitative causal relationships from upstream process parameters to downstream quality indicators.

[0007] Third, purely data-driven machine learning methods lack sufficient physical interpretability. Deep learning methods such as convolutional neural networks (CNN) and recurrent neural networks (RNN) treat the physical system as a black box, and the prediction results are not constrained by physical conservation laws. They have poor generalization ability in process parameter ranges outside the coverage of training data, and require a large amount of actual measured and labeled data from the production line, thus limiting their practicality.

[0008] Fourth, existing neural operator methods fail to effectively handle moving boundary problems and multi-domain cascade coupling. Operator learning methods such as Fourier neural operators (FNO) are mainly aimed at steady-state or fixed-domain problems of a single physical field, lacking an effective characterization mechanism for time-varying geometric boundaries in the electrodeposition process, and have not established a unified architecture to support lossless transfer and collaborative training of cross-domain field quantities.

[0009] Fifth, the model is not adaptable to batch-to-batch material variations and production line process fluctuations. The parameters of glass substrate materials fluctuate between batches, and the electroplating solution formulation and lithography equipment status also change with the production line maintenance cycle. Existing methods lack the ability to adapt to parameter distribution shifts, requiring the re-collection of a large amount of data and retraining, which cannot meet the dynamic process control requirements of mass production environments. Summary of the Invention

[0010] The purpose of this invention is to overcome the shortcomings of the prior art and provide a glass substrate RDL multi-field coupling prediction method that integrates physical constraint neural operators.

[0011] To achieve the above objectives, the technical solution provided by this invention is as follows:

[0012] A method for multi-field coupling prediction of RDL on glass substrates incorporating physically constrained neural operators includes the following steps:

[0013] S1. Obtain the basic input data for the large-size glass substrate RDL process, including the material properties of the glass substrate and each functional layer, the characteristic parameters of the interface fracture toughness and cohesion model, the discrete sequence of the process temperature history, the photolithography process parameters, the electrodeposition process parameters, and the geometric parameters of the glass substrate and trench; and perform dimensional normalization on these basic input data, uniformly discretize them into a standard tensor format, and complete the initialization of the unified data interface.

[0014] S2. Construct a thermo-mechanical-fracture coupled neural operator (TMFO) model. The material properties parameters of the glass substrate and each functional layer and the discrete sequence data of the process temperature history are used as inputs to the thermo-mechanical-fracture coupled neural operator model. The output of the thermo-mechanical-fracture coupled neural operator model is the residual stress tensor field and the interface damage variable field.

[0015] S3. Construct a Stress-Optics Mapping Operator (SOMO) model. The residual stress tensor field output by the thermo-mechanical-fracture coupled neural operator model, the refractive index distribution in the material property parameters obtained in step S1, and the glass substrate thickness distribution in the geometric parameters of the glass substrate and the trench are used as inputs to the stress-optics mapping neural operator model after passing through a unified data interface. The output of the stress-optics mapping neural operator model is the defocus field in the lithography exposure plane and the normalized in-plane light intensity distribution.

[0016] S4. The defocus field and normalized light intensity distribution in the lithography exposure plane output in step S3 are used as additional feature channels. These are tensor-stitched and spatially aligned with the glass substrate and trench geometry parameters obtained in step S1 to obtain a composite trench morphology tensor. Then, based on the composite trench morphology tensor, the sign distance function of each trench filling interface relative to the deposition start surface at the initial moment is calculated as an additional input channel of the spatiotemporal Fourier neural operator. The sign distance function, electrodeposition process parameters, and glass substrate and trench geometry parameters are integrated together to form the complete input tensor of the spatiotemporal Fourier neural operator.

[0017] S5. Construct an Electrodeposition Multi-physics Operator (EMDO) model and input the obtained complete input tensor into the Electrodeposition Multi-physics Operator model through a unified data interface. Output the copper filling morphology evolution field, surface roughness, and high-frequency signal loss.

[0018] S6. The constructed thermo-mechanical-fracture coupled neural operator model, force-optical mapping neural operator model, and electrodeposition multi-field coupled neural operator model are combined together as the basic network architecture. The thermo-mechanical-fracture coupled neural operator model, force-optical mapping neural operator model, and electrodeposition multi-field coupled neural operator model are pre-trained using simulation datasets generated based on finite element simulation and Monte Carlo parameter sampling, respectively, to obtain a pre-trained model with general multi-field coupling representation capability.

[0019] S7. For the target production line, collect actual measured labeled data; using the pre-trained model obtained in step S6 as initialization, freeze the weights of the bottom encoders of each operator model, and adjust the weights of the top Fourier modes with a small learning rate. Selective fine-tuning of output mapping layer parameters; simultaneously, to address the data distribution shift caused by fluctuations in material properties and changes in process conditions in different batches, an MMD domain adaptive constraint is introduced into the total fine-tuning loss. The Gaussian kernel function is used to measure the distribution difference of feature embedding between the simulation source domain and the production line target domain. By minimizing the MMD loss, the intermediate feature distribution of each operator model is aligned between the two domains, enabling the pre-trained model to adapt to changes in target production line material batches and processes under the conditions of measured labeled data.

[0020] S8. Through the fine-tuned pre-trained model, the full-link multi-field coupling of residual stress tensor field, interface damage variable field, in-plane defocus field of photolithography exposure, normalized in-plane light intensity distribution, copper filling morphology evolution field, surface roughness and high-frequency signal loss can be quickly predicted.

[0021] Furthermore, the thermal-mechanical-fracture coupled neural operator model is constructed based on Fourier neural operators, including:

[0022] Let the input function be... Defined in the spatial domain The network updates the hidden layer representations iteratively layer by layer. :

[0023] ;

[0024] in, For point-by-point linear transformation, It is a Fourier integral operator whose kernel function operates in the frequency domain in a truncated Fourier mode. Parameterization:

[0025] ;

[0026] in, For learnable complex-valued weight tensors, and These are the two-dimensional fast Fourier transform and its inverse transform, respectively.

[0027] Furthermore, the force-light mapping neural operator model adopts a hybrid architecture of U-Net and Fourier neural operator network. The encoder and decoder of U-Net are responsible for capturing local stress gradients and multi-scale spatial features. Fourier neural operator networks are embedded at the skip connections of each resolution layer to supplement long-range Fourier correlation. The encoder path consists of L downsampling blocks, each of which contains two convolutional layers and one Fourier neural operator network layer. The decoder path is symmetrically set and multi-scale features are fused through bilinear upsampling and skip connections. The final output head is a 1×1 convolution, which outputs the in-focus field and the normalized in-plane light intensity distribution of the lithography exposure plane, respectively.

[0028] Furthermore, the electrodeposition multi-field coupled neural operator model adopts a spatiotemporal Fourier neural operator model, simultaneously applying Fourier integral operators in both the spatial and temporal domains. Parameterization is achieved by truncating the number of modes in the three-dimensional frequency domain to capture the spatiotemporal evolution of the filling morphology. The moving boundary processing employs level set-assisted encoding, encoding the trench filling interface position as a signed distance function. As an additional input channel, it enables the model to implicitly perceive the position of the moving boundary.

[0029] Furthermore, the physical constraint loss function of the thermo-mechanical-fracture coupled neural operator model consists of a data fitting term and a physical residual term:

[0030] ;

[0031] in, , , These are the weighting coefficients for the residual terms of the heat conduction equation, the residual terms of the elasticity equilibrium equation, and the residual terms of the cohesive force model, respectively. This refers to the corresponding data fitting loss term;

[0032] For the residual terms of the heat conduction equation:

[0033] ;

[0034] This is the predicted value of the intermediate temperature field in the network. For the heat source term, constraints ensure that the predicted temperature field satisfies the transient heat conduction equation;

[0035] For the residual terms of the equilibrium equations in elasticity:

[0036] ;

[0037] It is a fourth-order elastic tensor. It is a volume force vector. For strain tensor, These are the predicted values ​​for the displacement field. The thermal strain tensor is constrained to ensure that the stress field satisfies mechanical equilibrium and constitutive relations; For reference temperature;

[0038] For the cohesive model residuals:

[0039] ;

[0040] This is the predicted value of the interface normal traction force. for Type fracture toughness, For feature separation length, For the Heaviside function, These are the separation values ​​for damage initiation and complete failure, respectively, and the constraints ensure that the damage variables satisfy the irreversible evolution condition. and These are the predicted values ​​of the interface normal and tangential separation, respectively. This is the predicted value of the derivative of the damage variable with respect to time; The effective separation quantity prediction value is obtained by combining the normal and tangential separation quantities; This represents the highest effective separation rate in history.

[0041] Furthermore, the physical constraint loss function of the force-light mapping neural operator model is as follows:

[0042] ;

[0043] in, and These are the weighting coefficients for the pupil function constraint term and the partially coherent illumination model constraint term, respectively. This refers to the corresponding data fitting loss term;

[0044] Constraints for the pupil function:

[0045] ;

[0046] This is a predicted value for the in-plane light intensity distribution; For the incident field, ensure that the predicted light intensity is consistent with the diffraction optical model derived from the predicted defocus amount; For the pupil function:

[0047] ;

[0048] The amplitude function of the pupil; Wavefront phase error introduced by defocusing:

[0049] ;

[0050] For the exposure wavelength, Numerical aperture, To normalize the pupil coordinates, This refers to the defocus amount;

[0051] Constraints for the partially coherent illumination model:

[0052] ;

[0053] .

[0054] Furthermore, the physical constraint loss function of the electrodeposition multi-field coupled neural operator model is as follows:

[0055] ;

[0056] in, , , , , These are the weighting coefficients for the constraint terms of the steady-state electrolyte electric field equation, the mass transfer equation, the electrochemical kinetics, the additive adsorption equation, and the end-to-end mapping constraint term of the morphology-roughness-high-frequency loss, respectively. This refers to the corresponding data fitting loss term;

[0057] Constraint terms for the steady-state electrolyte electric field equation:

[0058] ;

[0059] Electrolyte conductivity; The predicted value of the interface normal current density is constrained to ensure that the potential field satisfies the steady-state electric field equation and the continuity of the interface current within the electrolyte domain. This is the predicted value for the electrolyte potential; For the electrolyte domain; This represents the boundary of the electrodeposition interface.

[0060] For the mass transfer equation constraint terms:

[0061] ;

[0062] The diffusion coefficient of copper ions is . For the electrolyte flow velocity field, constraints ensure that the concentration field satisfies the mass conservation equation; This is a predicted value for copper ion concentration; This represents the total electrodeposition time.

[0063] Electrochemical kinetic constraint terms:

[0064] ;

[0065] For exchange current density; The anode-cathode transfer coefficient; This is an overpotential. This is the predicted value of the solid-state potential. This is the equilibrium potential for copper electrodeposition. The constant is Faraday's constant, and this constraint ensures the consistency of electrochemical kinetics at the interface; It is the ideal gas constant;

[0066] Constraint terms for the additive adsorption equation:

[0067] ;

[0068] The adsorption and desorption rate constants are constrained to ensure that the coupling relationship between the additive adsorption behavior and the local current density is satisfied. This is the predicted value for the surface coverage of the additive; This refers to the concentration of additives in the electrolyte. Reference current density;

[0069] For the end-to-end mapping constraint terms of topography-roughness-high frequency loss:

[0070] ;

[0071] This is a predicted value for high-frequency signal loss. The width of the interconnect; Surface resistance, The permeability of free space, The conductivity of copper; For skin depth; To establish the characteristic impedance of the transmission line, a complete physical causal chain was established from the predicted filling morphology through the intermediate roughness to the high-frequency loss, achieving an end-to-end interpretable mapping; the surface roughness was derived from the power spectral density (PSD) of the filling morphology.

[0072] ;

[0073] and These are the lower and upper frequency limits for the power spectral density integral, respectively.

[0074] Furthermore, in step S7, the total loss is fine-tuned. for:

[0075] ;

[0076] in, For mission losses; These are the weighting coefficients for the adaptive constraints in the MMD domain; For MMD domain adaptive constraints:

[0077] ;

[0078] This represents the number of samples in the source domain. The number of samples in the target domain; For the source domain Intermediate feature embedding of each sample; For the target domain Intermediate feature embedding of each sample; For the regenerating nucleus Hilbert space, The kernel mapping function uses a Gaussian kernel. , This is the bandwidth parameter of the Gaussian kernel function; Let be the norm in the regenerated kernel Hilbert space.

[0079] Furthermore, the unified data interface specification is as follows: the inputs and outputs of the thermal-mechanical-fracture coupled neural operator model, the force-optical mapping neural operator model, and the electrodeposition multi-field coupled neural operator model are all expressed in the form of discretized physical field tensors, and a unified approach is adopted. or The tensor format, in which For batch dimension, The number of physical field channels. For spatial resolution, This represents the number of time steps.

[0080] Compared with existing technologies, the principles and advantages of this technical solution are as follows:

[0081] 1. By adopting a surrogate model based on Fourier neural operators to replace the traditional full-scale full-process finite element simulation, the time for a single end-to-end prediction and inference is reduced from several hours to several days to seconds, meeting the engineering requirements for process parameter scanning, process window optimization, and online real-time prediction.

[0082] 2. By using a unified data interface, the lossless transfer of field quantities between the thermal-mechanical-fracture coupled neural operator model, the force-optical mapping neural operator model, and the electrodeposition multi-field coupled neural operator model is realized, forming an end-to-end quantitative causal relationship from upstream process parameters to downstream quality indicators, thus avoiding the cumulative errors introduced by the traditional fragmented simulation process.

[0083] 3. By embedding physical constraints such as heat conduction equation, elasticity equilibrium equation, cohesion model, pupil function, partially coherent imaging theory, electric field equation, mass transfer equation, Butler-Volmer equation, Langmuir adsorption kinetics and Hammerstad-Jensen loss model into the loss function, the prediction results strictly satisfy the corresponding physical laws and still have good generalization ability in the process parameter range outside the training data coverage.

[0084] 4. By implicitly embedding the time-varying geometric boundary into the input channel of the spatiotemporal Fourier neural operator through level set-assisted encoding, the network can perceive the dynamic evolution of the trench filling interface during electrodeposition, avoiding the dependence of traditional methods on unstructured meshes and effectively handling the moving boundary problem.

[0085] 5. By using transfer learning and maximum mean difference domain adaptive methods, it can quickly adapt to the fluctuations in material properties and changes in process conditions of different batches with only a small amount of production line test data, thus meeting the dynamic process control requirements in mass production environments. Attached Figure Description

[0086] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0087] Figure 1 This is a flowchart illustrating the principle of the glass substrate RDL multi-field coupling prediction method that integrates physical constraint neural operators according to the present invention. Detailed Implementation

[0088] The present invention will be further described below with reference to specific embodiments:

[0089] like Figure 1 As shown in this embodiment, the glass substrate RDL multi-field coupling prediction method integrating physically constrained neural operators includes the following steps:

[0090] S1. Obtain the basic input data for the large-size glass substrate RDL process, including the material properties (thermal conductivity) of the glass substrate and each functional layer. Volumetric heat capacity Coefficient of thermal expansion Lamé coefficient shear modulus With refractive index distribution ), interfacial fracture toughness and cohesive force model characteristic parameters ( ), complete discrete sequence of process temperature history Photolithography process parameters (exposure wavelength) Numerical aperture Partial coherence factor Pupil amplitude function Electrodeposition process parameters (current density timing) Copper ion concentration in electrolyte Additive concentration Electrodeposition temperature and the geometric parameters of the glass substrate and trench (width) ,depth Aspect Ratio ); and then perform dimensional normalization on these basic input data, uniformly discretizing them into ); or Standard tensor format, complete the Unified Data Interface (UDI) initialization.

[0091] S2. Construct a thermo-mechanical-fracture coupled neural operator model to obtain the material properties of the glass substrate and each functional layer. Discrete sequence data of process temperature history As the input to this thermo-mechanical-fracture coupled neural operator model, the output of this model is the residual stress tensor field. and interface damage variable field ;

[0092] In this step, the thermal-mechanical-fracture coupled neural operator model is constructed based on Fourier neural operators, including:

[0093] Let the input function be... Defined in the spatial domain The network updates the hidden layer representations iteratively layer by layer. :

[0094] ;

[0095] in, For point-by-point linear transformation, It is a Fourier integral operator whose kernel function operates in the frequency domain in a truncated Fourier mode. Parameterization:

[0096] ;

[0097] in, For learnable complex-valued weight tensors, and These are the two-dimensional fast Fourier transform and its inverse transform, respectively.

[0098] The physical constraint loss function of the thermo-mechanical-fracture coupled neural operator model consists of a data fitting term and a physical residual term:

[0099] ;

[0100] in, , , These are the weighting coefficients for the residual terms of the heat conduction equation, the residual terms of the elasticity equilibrium equation, and the residual terms of the cohesive force model, respectively. This refers to the corresponding data fitting loss term;

[0101] For the residual terms of the heat conduction equation:

[0102] ;

[0103] This is the predicted value of the intermediate temperature field in the network. For the heat source term, constraints ensure that the predicted temperature field satisfies the transient heat conduction equation;

[0104] For the residual terms of the elasticity equilibrium equation:

[0105] ;

[0106] It is a fourth-order elastic tensor. It is a volume force vector. For strain tensor, These are the predicted values ​​for the displacement field. The thermal strain tensor is constrained to ensure that the stress field satisfies mechanical equilibrium and constitutive relations; For reference temperature;

[0107] For the cohesive model residuals:

[0108] ;

[0109] This is the predicted value of the interface normal traction force. for Type fracture toughness, For feature separation length, For the Heaviside function, These are the separation values ​​for damage initiation and complete failure, respectively, and the constraints ensure that the damage variables satisfy the irreversible evolution condition. and These are the predicted values ​​of the interface normal and tangential separation, respectively. This is the predicted value of the derivative of the damage variable with respect to time; The effective separation quantity prediction value is obtained by combining the normal and tangential separation quantities; This represents the highest effective separation rate in history.

[0110] S3. Construct a force-light mapping neural operator model, and obtain the residual stress tensor field output from the thermo-force-fracture coupled neural operator model. and the refractive index distribution in the material properties obtained in step S1 Glass substrate thickness distribution in the geometric parameters of glass substrate and trench After passing through a unified data interface, the input to this force-light mapping neural operator model is used. This model employs a hybrid architecture of U-Net and Fourier neural operator networks. The U-Net encoder and decoder are responsible for capturing local stress gradients and multi-scale spatial features. Fourier neural operator networks are embedded at the skip connections of each resolution layer to supplement long-range Fourier correlations. The encoder path consists of L downsampling blocks, each containing two convolutional layers and one Fourier neural operator network layer. The decoder path is symmetrically set, fusing multi-scale features through bilinear upsampling and skip connections. The final output head is a 1×1 convolution, outputting the defocus field within the lithographic exposure plane. and normalized in-plane light intensity distribution .

[0111] The physical constraint loss function of the force-light mapping neural operator model is as follows:

[0112] ;

[0113] in, and These are the weighting coefficients for the pupil function constraint term and the partially coherent illumination model constraint term, respectively. This refers to the corresponding data fitting loss term;

[0114] Constraints for the pupil function:

[0115] ;

[0116] This is a predicted value for the in-plane light intensity distribution; For the incident field, ensure that the predicted light intensity is consistent with the diffraction optical model derived from the predicted defocus amount; For the pupil function:

[0117] ;

[0118] The pupil amplitude function; Wavefront phase error introduced by defocusing:

[0119] ;

[0120] For the exposure wavelength, Numerical aperture, To normalize the pupil coordinates, This refers to the defocus amount;

[0121] Constraints for the partially coherent illumination model:

[0122] ;

[0123] .

[0124] S4, the in-focus field of the lithographic exposure plane output in step S3 and normalized in-plane light intensity distribution As an additional feature channel, it is tensor-stitched and spatially aligned with the geometric parameters of the glass substrate and trench obtained in step S1 to obtain the composite trench morphology tensor; then, based on the composite trench morphology tensor, the signed distance function of each trench filling interface relative to the deposition initiation surface at the initial moment is calculated. As an additional input channel of the spacetime Fourier neural operator (ST-FNO), and the symbolic distance function is used. With electrodeposition process parameters ( , , , The geometric parameters of the glass substrate and trench are integrated together to form the complete input tensor of the spatiotemporal Fourier neural operator, with a unified format. ;

[0125] S5. Construct a multi-field coupled neural operator model for electrodeposition, and input the obtained complete input tensor into the model through a unified data interface to output the copper filling morphology evolution field. Surface roughness and and high-frequency signal loss ;Depend on Extract final fill rate Distribution of cavities.

[0126] In this step, the electrodeposition multi-field coupled neural operator model adopts a spatiotemporal Fourier neural operator model, simultaneously applying Fourier integral operators in both the spatial and temporal domains, and through the three-dimensional frequency domain... The truncated pattern number is parameterized to capture the spatiotemporal evolution of the filling morphology; the moving boundary processing uses level set-assisted coding to encode the trench filling interface position as a signed distance function. As an additional input channel, it enables the model to implicitly perceive the position of the moving boundary.

[0127] The physical constraint loss function of the electrodeposition multi-field coupled neural operator model is as follows:

[0128] ;

[0129] in, , , , , These are the weighting coefficients for the constraint terms of the steady-state electrolyte electric field equation, the mass transfer equation, the electrochemical kinetics, the additive adsorption equation, and the end-to-end mapping constraint term of the morphology-roughness-high-frequency loss, respectively. This refers to the corresponding data fitting loss term;

[0130] Constraint terms for the steady-state electrolyte electric field equation:

[0131] ;

[0132] Electrolyte conductivity; The predicted value of the interface normal current density is constrained to ensure that the potential field satisfies the steady-state electric field equation and the continuity of the interface current within the electrolyte domain. This is the predicted value for the electrolyte potential; For the electrolyte domain; This represents the boundary of the electrodeposition interface.

[0133] For the mass transfer equation constraint terms:

[0134] ;

[0135] The diffusion coefficient of copper ions is . For the electrolyte flow velocity field, constraints ensure that the concentration field satisfies the mass conservation equation; This is a predicted value for copper ion concentration; This represents the total electrodeposition time.

[0136] Electrochemical kinetic constraint terms:

[0137] ;

[0138] For exchange current density; The anode-cathode transfer coefficient; This is an overpotential. This is the predicted value of the solid-state potential. This is the equilibrium potential for copper electrodeposition. The constant is Faraday's constant, and this constraint ensures the consistency of electrochemical kinetics at the interface; It is the ideal gas constant;

[0139] Constraint terms for the additive adsorption equation:

[0140] ;

[0141] The adsorption and desorption rate constants are constrained to ensure that the coupling relationship between the additive adsorption behavior and the local current density is satisfied. This is the predicted value for the surface coverage of the additive; The concentration of additives in the electrolyte; Reference current density;

[0142] For the end-to-end mapping constraint terms of topography-roughness-high frequency loss:

[0143] ;

[0144] This is a predicted value for high-frequency signal loss. The width of the interconnect; Surface resistance, The permeability of free space, The conductivity of copper; For skin depth; To establish the characteristic impedance of the transmission line, a complete physical causal chain was established from the predicted filling morphology through the intermediate roughness to the high-frequency loss, achieving an end-to-end interpretable mapping; the surface roughness was derived from the power spectral density (PSD) of the filling morphology.

[0145] ;

[0146] and These are the lower and upper frequency limits for the power spectral density integral, respectively.

[0147] S6. The constructed thermo-mechanical-fracture coupled neural operator model, force-optical mapping neural operator model, and electrodeposition multi-field coupled neural operator model are combined together as the basic network architecture. The thermo-mechanical-fracture coupled neural operator model, force-optical mapping neural operator model, and electrodeposition multi-field coupled neural operator model are pre-trained using a large-scale simulation dataset generated based on finite element simulation (FEM) and Monte Carlo parameter sampling, respectively, to obtain a pre-trained model with general multi-field coupling representation capability.

[0148] S7. For the target production line, collect a small amount of measured labeled data (including warpage measurement, linewidth CD measurement, resistivity test, and high-frequency S-parameter test data); using the pre-trained model obtained in step S6 as initialization, freeze the weights of the bottom encoder of each operator model (to prevent catastrophic forgetting), and apply a small learning rate to the weights of the top Fourier mode. Selective fine-tuning of output mapping layer parameters is performed. Meanwhile, for the data distribution shift caused by changes in material properties (dielectric constant, thermal expansion coefficient batch fluctuations) and process conditions (electroplating solution formula changes, lithography machine differences) of different batches, an MMD (Maximum Mean Discrepancy) domain adaptive constraint is introduced into the total fine-tuning loss. The Gaussian kernel function is used to measure the distribution difference of feature embedding between the simulation source domain (simulation data) and the production line target domain (production line measured data). By minimizing the MMD loss, the intermediate feature distribution of each operator model is aligned between the two domains, so that the pre-trained model can adapt to the changes in target production line material batches and processes under the conditions of measured labeled data.

[0149] Specifically, fine-tuning the total loss for:

[0150] ;

[0151] in, For mission losses; These are the weighting coefficients for the adaptive constraints in the MMD domain; For MMD domain adaptive constraints:

[0152] ;

[0153] This represents the number of samples in the source domain. The number of samples in the target domain; For the source domain Intermediate feature embedding of each sample; For the target domain Intermediate feature embedding of each sample; For the regenerating nucleus Hilbert space, The kernel mapping function uses a Gaussian kernel. , This is the bandwidth parameter of the Gaussian kernel function; Let be the norm in the regenerated kernel Hilbert space.

[0154] S8. Through the fine-tuned pre-trained model, the full-link multi-field coupling of residual stress tensor field, interface damage variable field, in-plane defocus field of photolithography exposure, normalized in-plane light intensity distribution, copper filling morphology evolution field, surface roughness and high-frequency signal loss can be quickly predicted.

[0155] In this invention,

[0156] Redistribution Layer (RDL): A fine-pitch metal interconnect layer fabricated on the surface of a packaging substrate using photolithography and electroplating processes, used to redistribute signals between chip pads and package pins. This invention specifically refers to submicron-scale copper interconnect structures fabricated on large-size glass substrates, the quality of which is determined by the synergistic effects of multiple physical domains, including thermal field, force field, optical imaging, and electrochemical deposition.

[0157] Fourier Neural Operator (FNO): A neural network architecture based on frequency domain integral operators, which models the mapping between function spaces by parameterizing convolution kernels in a truncated Fourier mode space with learnable complex-valued weight tensors. Compared to traditional convolutional neural networks, FNO is resolution-independent and can efficiently capture long-range spatial correlations of physical fields in the global frequency domain, making it suitable for fast solution surrogate modeling of partial differential equations. This invention uses FNO as the core backbone network of a thermo-mechanical-fracture coupled neural operator.

[0158] Cohesive Z-M (CZM) is a fracture mechanics model that describes the initiation and propagation of interfacial damage in materials. It characterizes the mechanical response of the interface during normal and tangential separation processes through the traction force-separation relationship. This invention employs an exponential cohesive Z-M model, using type I fracture toughness, characteristic separation length, and irreversible damage evolution conditions as physical constraints, embedded in the loss function of a thermo-mechanical-fracture coupled neural operator, to predict the interfacial damage variable field of large-size glass substrates during thermal cycling.

[0159] Spatiotemporal Fourier Neural Operator (ST-FNO): An extended form of neural operator that simultaneously applies Fourier integral operators in both the spatial and temporal domains. By truncating the spatial and temporal mode numbers in the three-dimensional frequency domain, it achieves efficient modeling of the evolution of time-varying physical fields. This invention uses ST-FNO as the core network of a multi-field coupled neural operator for electrodeposition to handle the dynamic evolution of moving boundaries within the trench during copper filling, while simultaneously capturing the spatiotemporal coupling characteristics of the filling morphology.

[0160] The Butler-Volmer equation is a fundamental kinetic equation describing the relationship between current density and overpotential in electrochemical interfacial reactions. It characterizes the charge transfer process using parameters such as the anodic transfer coefficient, cathode transfer coefficient, and exchange current density. This invention embeds the Butler-Volmer equation residuals as constraint terms into the loss function of the electrodeposition multi-field coupled neural operator, ensuring that the predicted interfacial current density and overpotential distribution satisfy the fundamental physical laws of electrochemical kinetics, thereby guaranteeing the physical self-consistency of the copper filling morphology prediction results.

[0161] The Hammerstad-Jensen loss model is a semi-empirical model for predicting the impact of conductor surface roughness on high-frequency transmission loss. It corrects the surface resistance of an ideal smooth conductor by introducing a roughness correction factor, using the ratio of the root mean square surface roughness Rq to the skin depth as a key parameter. This invention embeds this model as a physical constraint into an electrodeposition multi-field coupled neural operator, establishing an end-to-end interpretable mapping link from predicting the filling morphology, extracting roughness through power spectral density, and then calculating the high-frequency signal loss alpha-loss(f,x).

[0162] The above-described embodiments are merely preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Therefore, any changes made in accordance with the shape and principle of the present invention should be covered within the protection scope of the present invention.

Claims

1. A multi-field coupling prediction method for RDL on glass substrates that integrates physically constrained neural operators, characterized in that, Includes the following steps: S1. Obtain the basic input data for the large-size glass substrate RDL process, including the material properties of the glass substrate and each functional layer, the characteristic parameters of the interface fracture toughness and cohesion model, the discrete sequence of the process temperature history, the photolithography process parameters, the electrodeposition process parameters, and the geometric parameters of the glass substrate and trench; and perform dimensional normalization on these basic input data, uniformly discretize them into a standard tensor format, and complete the initialization of the unified data interface. S2. Construct a thermo-mechanical-fracture coupled neural operator model. The material property parameters and discrete sequence data of the process temperature history of the glass substrate and each functional layer are used as inputs to the thermo-mechanical-fracture coupled neural operator model. The output of the thermo-mechanical-fracture coupled neural operator model is the residual stress tensor field and the interface damage variable field. S3. Construct a force-light mapping neural operator model. The residual stress tensor field output by the thermo-mechanical-fracture coupled neural operator model, the refractive index distribution in the material property parameters obtained in step S1, and the glass substrate thickness distribution in the geometric parameters of the glass substrate and the trench are used as inputs to the force-light mapping neural operator model after passing through a unified data interface. The output of the force-light mapping neural operator model is the defocus field in the lithography exposure plane and the normalized in-plane light intensity distribution. S4. The defocus field and normalized light intensity distribution in the lithography exposure plane output in step S3 are used as additional feature channels. These are tensor-stitched and spatially aligned with the glass substrate and trench geometry parameters obtained in step S1 to obtain a composite trench morphology tensor. Then, based on the composite trench morphology tensor, the sign distance function of each trench filling interface relative to the deposition start surface at the initial moment is calculated as an additional input channel of the spatiotemporal Fourier neural operator. The sign distance function, electrodeposition process parameters, and glass substrate and trench geometry parameters are integrated together to form the complete input tensor of the spatiotemporal Fourier neural operator. S5. Construct an electrodeposition multi-field coupled neural operator model, and input the obtained complete input tensor into the electrodeposition multi-field coupled neural operator model through a unified data interface to output the copper filling morphology evolution field, surface roughness and high-frequency signal loss. S6. The constructed thermo-mechanical-fracture coupled neural operator model, force-optical mapping neural operator model, and electrodeposition multi-field coupled neural operator model are combined together as the basic network architecture. The thermo-mechanical-fracture coupled neural operator model, force-optical mapping neural operator model, and electrodeposition multi-field coupled neural operator model are pre-trained using simulation datasets generated based on finite element simulation and Monte Carlo parameter sampling, respectively, to obtain a pre-trained model with general multi-field coupling representation capability. S7. Collect measured and labeled data for the target production line; Using the pre-trained model obtained in step S6 as initialization, freeze the weights of the bottom encoders of each operator model, and apply a small learning rate to the weights of the top Fourier modes. Selective fine-tuning of output mapping layer parameters; simultaneously, to address the data distribution shift caused by fluctuations in material properties and changes in process conditions in different batches, an MMD domain adaptive constraint is introduced into the total fine-tuning loss. The Gaussian kernel function is used to measure the distribution difference of feature embedding between the simulation source domain and the production line target domain. By minimizing the MMD loss, the intermediate feature distribution of each operator model is aligned between the two domains, enabling the pre-trained model to adapt to changes in target production line material batches and processes under the conditions of measured labeled data. S8. Through the fine-tuned pre-trained model, the full-link multi-field coupling of residual stress tensor field, interface damage variable field, in-plane defocus field of photolithography exposure, normalized in-plane light intensity distribution, copper filling morphology evolution field, surface roughness and high-frequency signal loss can be quickly predicted.

2. The glass substrate RDL multi-field coupling prediction method integrating physically constrained neural operators according to claim 1, characterized in that, The thermo-mechanical-fracture coupled neural operator model is constructed based on Fourier neural operators and includes: Let the input function be... Defined in the spatial domain The network updates the hidden layer representations iteratively layer by layer. : ； in, For point-by-point linear transformation, It is a Fourier integral operator whose kernel function operates in the frequency domain in a truncated Fourier mode. Parameterization: ； in, For learnable complex-valued weight tensors, and These are the two-dimensional fast Fourier transform and its inverse transform, respectively.

3. The glass substrate RDL multi-field coupling prediction method integrating physically constrained neural operators according to claim 1, characterized in that, The force-light mapping neural operator model adopts a hybrid architecture of U-Net and Fourier neural operator network. The encoder and decoder of U-Net are responsible for capturing local stress gradients and multi-scale spatial features. Fourier neural operator network is embedded at the skip connections of each resolution layer to supplement long-range Fourier correlation. The encoder path consists of L downsampling blocks, each of which contains two convolutional layers and one Fourier neural operator network layer. The decoder path is symmetrically set, and multi-scale features are fused through bilinear upsampling and skip connections; the final output head is a 1×1 convolution, which outputs the in-focus field and the normalized in-plane light intensity distribution of the lithographic exposure plane, respectively.

4. The glass substrate RDL multi-field coupling prediction method integrating physically constrained neural operators according to claim 1, characterized in that, The electrodeposition multi-field coupled neural operator model employs a spatiotemporal Fourier neural operator model, simultaneously applying Fourier integral operators in both the spatial and temporal domains. Parameterization is achieved by truncating the number of modes in the three-dimensional frequency domain to capture the spatiotemporal evolution of the filling morphology. Moving boundary processing utilizes level set-assisted encoding, coding the trench filling interface position as a signed distance function. As an additional input channel, it enables the model to implicitly perceive the position of the moving boundary.

5. The glass substrate RDL multi-field coupling prediction method integrating physically constrained neural operators according to claim 2, characterized in that, The physical constraint loss function of the thermo-mechanical-fracture coupled neural operator model consists of a data fitting term and a physical residual term: ； in, , , These are the weighting coefficients for the residual terms of the heat conduction equation, the residual terms of the elasticity equilibrium equation, and the residual terms of the cohesive force model, respectively. This refers to the corresponding data fitting loss term; For the residual terms of the heat conduction equation: ； This is the predicted value of the intermediate temperature field in the network. For the heat source term, constraints ensure that the predicted temperature field satisfies the transient heat conduction equation; For the residual terms of the equilibrium equations in elasticity: ； It is a fourth-order elastic tensor. It is a volume force vector. For strain tensor, These are the predicted values ​​of the displacement field. The thermal strain tensor is constrained to ensure that the stress field satisfies mechanical equilibrium and constitutive relations; For reference temperature; For the cohesive model residuals: ； This is the predicted value of the interface normal traction force. for Type fracture toughness, For feature separation length, For the Heaviside function, These are the separation values ​​for damage initiation and complete failure, respectively, and the constraints ensure that the damage variables satisfy the irreversible evolution condition. and These are the predicted values ​​of the interface normal and tangential separation, respectively. This is the predicted value of the derivative of the damage variable with respect to time; The effective separation quantity prediction value is obtained by combining the normal and tangential separation quantities; This represents the highest effective separation rate in history.

6. The glass substrate RDL multi-field coupling prediction method integrating physically constrained neural operators according to claim 3, characterized in that, The physical constraint loss function of the force-light mapping neural operator model is as follows: ； in, and These are the weighting coefficients for the pupil function constraint term and the partially coherent illumination model constraint term, respectively. This refers to the corresponding data fitting loss term; Constraints for the pupil function: ； This is a predicted value for the in-plane light intensity distribution; For the incident field, ensure that the predicted light intensity is consistent with the diffraction optical model derived from the predicted defocus amount; For the pupil function: ； The amplitude function of the pupil; Wavefront phase error introduced by defocusing: ； For the exposure wavelength, Numerical aperture, To normalize the pupil coordinates, This refers to the defocusing amount; Constraints for the partially coherent illumination model: ； For the first Incident light field of a discrete point source; For the first The weights of discrete point sources.

7. The glass substrate RDL multi-field coupling prediction method integrating physically constrained neural operators according to claim 1, characterized in that, The physical constraint loss function of the electrodeposition multi-field coupled neural operator model is as follows: ； in, , , , , These are the weighting coefficients for the constraint terms of the steady-state electrolyte electric field equation, the mass transfer equation, the electrochemical kinetics, the additive adsorption equation, and the end-to-end mapping constraint term of the morphology-roughness-high-frequency loss, respectively. This refers to the corresponding data fitting loss term; Constraint terms for the steady-state electrolyte electric field equation: ； Electrolyte conductivity; The predicted value of the interface normal current density is constrained to ensure that the potential field satisfies the steady-state electric field equation and the continuity of the interface current within the electrolyte domain. This is the predicted value for the electrolyte potential; For the electrolyte domain; This represents the boundary of the electrodeposition interface. For the mass transfer equation constraint terms: ； The diffusion coefficient of copper ions is . For the electrolyte flow velocity field, constraints ensure that the concentration field satisfies the mass conservation equation; This is a predicted value for copper ion concentration; This represents the total electrodeposition time. Electrochemical kinetic constraint terms: ； For exchange current density; The anode-cathode transfer coefficient; This is an overpotential. This is the predicted value of the solid-state potential. This is the equilibrium potential for copper electrodeposition. The constant is Faraday's constant, and this constraint ensures the consistency of electrochemical kinetics at the interface; It is the ideal gas constant; Constraint terms for the additive adsorption equation: ； The adsorption and desorption rate constants are constrained to ensure that the coupling relationship between the additive adsorption behavior and the local current density is satisfied. This is the predicted value for the surface coverage of the additive; This refers to the concentration of additives in the electrolyte. Reference current density; For the end-to-end mapping constraint terms of topography-roughness-high frequency loss: ； This is a predicted value for high-frequency signal loss. The width of the interconnect; Surface resistance, The permeability of free space, The conductivity of copper; For skin depth; To establish the characteristic impedance of the transmission line, a complete physical causal chain was established from the predicted filling morphology through the intermediate roughness to the high-frequency loss, achieving an end-to-end interpretable mapping; the surface roughness was derived from the power spectral density (PSD) of the filling morphology. ； and These are the lower and upper frequency limits for the power spectral density integral, respectively.

8. The glass substrate RDL multi-field coupling prediction method integrating physically constrained neural operators according to claim 1, characterized in that, In step S7, fine-tune the total loss. for: ； in, For mission losses; These are the weighting coefficients for the adaptive constraints in the MMD domain; For MMD domain adaptive constraints: ； This represents the number of samples in the source domain. The number of samples in the target domain; For the source domain Intermediate feature embedding of each sample; For the target domain Intermediate feature embedding of each sample; For the regenerating nucleus Hilbert space, The kernel mapping function uses a Gaussian kernel. , The bandwidth parameter of the Gaussian kernel function; Let be the norm in the regenerated kernel Hilbert space.

9. The glass substrate RDL multi-field coupling prediction method integrating physically constrained neural operators according to claim 1, characterized in that, The unified data interface specification is as follows: the inputs and outputs of the thermal-mechanical-fracture coupled neural operator model, the force-optical mapping neural operator model, and the electrodeposition multi-field coupled neural operator model are all expressed in the form of discretized physical field tensors, and a unified approach is adopted. or The tensor format, in which For batch dimension, The number of physical field channels. For spatial resolution, This represents the number of time steps.

Images