Global-in-time Spectral Duhamel-Picard Method for Post-Exposure Bake Simulation and Physics-Informed Neural Operators for EUV Lithography

The gist: The proposed method presents a global-in-time spectral Duhamel–Picard method and a physics-informed neural operator (PINO) that enable accelerated, time-marching-free prediction of concentration fields at arbitrary times within the bake interval.

Governing Equations and Solution Method

The PEB process is governed by four coupled concentration fields: acid A, quencher Q, protected group P, and deprotected group D The governing equations are formulated as a coupled partial differential equation (PDE) system, and the system is solved by explicit time marching, which scales as ∆tmax ∝ ∆x 2/D Conventional numerical methods focus on improving the accuracy and efficiency of time marching.

Duhamel Formulation for PEB Reaction-Diffusion System

The governing equations are recast in a Duhamel form that incorporates linear diffusion through the heat semigroup The nonlinear acid and quencher trajectories are then represented by global polynomial expansions, whose temporal coefficient fields are determined using a combination of Anderson-accelerated Picard iteration and Gaussian quadrature. For the coupled A and Q trajectories, Anderson-accelerated Picard iteration is used to solve the coupled nonlinear Volterra integral equations.

Global-in-Time Picard Iteration with Polynomials and Gaussian Quadrature

The Duhamel time integral at an arbitrary time t ∈ [0, tPEB] is evaluated using an N-point Gaussian quadrature rule The complete trajectory is then represented as u(x, t) ≈ K ∑ −1 j=0 c u j(x) Φ j (η(t)), t ∈ [0, tPEB], where the temporal coefficient fields are the spatially resolved temporal coefficient fields. The iteration is terminated when the prescribed residual tolerance is reached, at which point the complete trajectories are determined by the temporal coefficient fields c A(x) and c Q(x).

Mixed-transform Factorized Fourier Neural Operator

The global-in-time coefficient representation is carried over to the PINO, which predicts temporal coefficient fields rather than concentration fields at individual time locations. The mixed-transform f-FNO architecture replaces a fully coupled three-dimensional spectral kernel with an axial operator followed by factorized lateral operators. For training, the mixed-transform f-FNO predicts the temporal coefficient fields c u theta of the nonlinear correction over the full bake interval [0, tPEB].

Loss Functions for Mixed-transform f-FNO PINO

The physics-informed objective is a weighted sum: L = λPDE,AL A PDE + λPDE,QL Q PDE + λint,AL A int + λint,QL Q int + λICLIC + Lend + Lcon. The differential residual loss is defined by Lu PDE = 1/Nt u 2 0 ∑ Nt n=1 Ru [u theta](⋅, t n) squared, u ∈ A, Q. The integral-consistency losses are defined by Lu int = 1/Npair u 2 0 ∑ (t i,t j)∈P ‖Iu [u theta](⋅;t i, t j) squared, u ∈ A, Q.

Results and Discussion

The global-in-time spectral Duhamel–Picard method reaches a self-consistency level of approximately 10−7 within 10 iterations. The wall-clock performance shows that the global-in-time calculation takes about 80 s to reach a convergence tolerance of 10−5. The PINO predictions show good agreement with the reference Picard solutions both for the acid and protected groups. The PINO preserves the principal three-dimensional post-development topology across randomly sampled mask layouts and PEB times while providing a modest runtime benefit.

Conclusion

The work presents a global-in-time spectral Duhamel–Picard method and a PINO that accurately and efficiently evaluates the acid, quencher, protected-group, and deprotected-group concentration trajectories over the entire PEB interval. The three-dimensional profiles predicted by the PINO closely match those obtained using the global-in-time Duhamel–Picard method at z-planes closer to the resist surface. The PINO preserves the principal three-dimensional post-development topology across randomly sampled mask layouts and PEB times while reducing the runtime relative to the global-in-time Duhamel–Picard method. The final thresholding at T = tdev may further translate small arrival-time differences into finite contour displacements.

Acknowledgments

The authors express gratitude to Seungjin Lee for valuable discussions during the initial stages of this work. This work has been enabled in part by the NanoIC pilot line. The acquisition and operation are jointly funded by the Chips Joint Undertaking, through the European Union’s Digital Europe (101183266) and Horizon Europe programs (101183277), as well as by the participating states Belgium (Flanders), France, Germany, Finland, Ireland and Romania.

References

[1] Chris A. Mack. Fundamental Principles of Optical Lithography: The Science of Microfabrication

[2] Hiroshi Fukuda and Shinji Okazaki. Analysis of chemical amplification resist systems using a kinetic model and numerical simulation

[3] Taiki Kimura, Tetsuaki Matsunawa, and Shoji Mimotogi. Compact resist model using single convolution kernel

[4] Kyoungyoon Park, Dongyeul Lee, Yangwoo Heo, Pillsoo Kang, Useong Kim, Kyoil Koo, Seongjeon Choi, Seung-Hune Yang, and Seongtae Jeong. Development of an opc compact model integrating underlayer effects in metal oxide resists

[5] Yuri Granik and Gurdaman S. Khaira. Resist shrinkage during development: rigorous simulations and first compact model for opc

[6] Pei-Chun Lin and Charlie Chung-Ping Chen. High speed and flexible peb 3d diffusion simulation based on sylvester equation

[7] Bengt Fornberg. Generation of finite difference formulas on arbitrarily spaced grids

[8] Qipan Wang, Xiaohan Gao, Yibo Lin, Runsheng Wang, and Ru Huang. Deepeb: A neural partial differential equation solver for post exposure baking simulation in lithography

[9] Gangmin Cho, Taeyoung Kim, Seohyun Kim, and Youngsoo Shin. A fast and accurate peb simulation through recurrent neural network

[10] Alok Dutt, Leslie Greengard, and Vladimir Rokhlin. Spectral deferred correction methods for ordinary differential equations

[11] Peter Linz. Analytical and Numerical Methods for Volterra Equations

[12] Lawrence C. Evans. Partial Differential Equations, volume 19 of Graduate Studies in Mathematics

[13] Amnon Pazy. Semigroups of Linear Operators and Applications to Partial Differential Equations

[14] David Gottlieb and Steven A. Orszag.

1. No sequential time marching is required for PEB simulation, as Full trajectories are represented by finite polynomial expansions determined globally over the bake interval via Anderson-accelerated Picard iteration rather than sequential time marching. This allows for accelerated, time-marching-free prediction of concentration fields at arbitrary times within the bake interval.

2. PINO training can be performed without precomputed solutions, as the method is trained without precomputed solution trajectories, instead learning a mapping from initial fields to temporal coefficient fields. The PINO can then predict the temporal coefficient fields and reconstruct the complete trajectory of each field using a truncated Chebyshev or Legendre polynomial expansion.

3. The system can predict concentration profiles at any arbitrary time within the bake interval, as the reconstructed solution is defined by u(x, t) ≈ K ∑ −1 j=0 c u j(x) Φ j (η(t)), t ∈ [0, tPEB], where coefficients are determined by the converged temporal coefficient fields c A(x) and c Q(x).

4. The PINO can provide a modest runtime benefit over the global-in-time method, achieving an inference time of about 80% of the time required for the reference Picard solution, which is useful for high-throughput optimization.

5. The system can accurately predict post-development characteristics by evaluating CD errors at various planes, as predicted profiles are compared to experimental data, although errors increase deeper into the resist due to accumulation of differences in the acid concentration history.