Physics-Based Simulation of Master Template Fabrication: Integrated Modeling of Resist Coating, Electron Beam Lithography, and Reactive Ion Etching

1. Introduction

1.1. Background and Motivation

Nanoimprint lithography (NIL) has attracted extensive interest as a low-cost and high-resolution patterning technique with the ability to replicate near-10-nm features over a large area [1,2,3]. The ongoing miniaturization of semiconductor devices has commonly been framed in the context of Moore’s law. While occasionally it is proclaimed that “Moore’s law is dead,” these proclamations are only partly true. While geometric scaling of feature sizes is nearing its physical limits, the enhancement of functional density and performance per unit area continues to be possible by employing novel lithographic techniques. In this context, NIL represents another dimension of progress, enabling performance scaling that is not solely dependent on shrinking the diffraction-limited resolution of conventional optical lithography. Recent industrial assessments further indicate that NIL maturation is accelerating, particularly as its applicability expands into semiconductor, LED, and photovoltaic manufacturing, where cost-efficient, high-resolution pattern replication is increasingly required. These developments underscore that NIL is advancing not only as a laboratory-scale patterning method but also as a technology transitioning toward broader industrial deployment, supported by improvements in materials, tooling, and replication workflows [4,5].

NIL is commonly acknowledged as one of the next-generation lithography (NGL) family members, and it plays a significant role in both microelectronics and other nanomanufacturing applications. As an NGL candidate, NIL offers a compelling balance of cost efficiency and resolution, which makes it an interesting complement or alternative to other advanced lithographic technologies [4,6].

The NIL process can be generally categorized into three stages: (i) master template fabrication, (ii) replica template production, and (iii) final wafer imprinting [7]. Among these stages, the fidelity of the master template is paramount, as it dictates the precision of both replica templates and wafer-level patterns. Recent advances in master-template electron-beam lithography have demonstrated hp10-nm and sub-10-nm patterning through resist–process co-optimization and SADP-based pattern multiplication, underscoring the continued importance of high-resolution template fabrication for NIL template performance [8]. Any imperfection in the master template is unavoidably passed on to the subsequent processes, rendering its fabrication as one of the most technically challenging tasks in the development of NIL [9,10,11,12]. Specifically, the surface roughness of the master template needs to be regulated at the low-nanometer level (approximately 2–3 nm root mean square (RMS)) to prevent demolding issues or defect formation during the NIL process [13,14].

In contrast to traditional photomasks, NIL master templates need to reconcile nanoscale resolution with resilience to multiple uses. Even slight variations in resist thickness or etch profile can result in cumulative distortions after being transferred to wafers [15,16]. Consequently, computational methods that anticipate systematic fabrication errors are becoming ever more necessary to aid process optimization and optical proximity correction (OPC) [17]. In particular, OPC not only serves to correct for process-induced distortions in the final patterns but also offers a pathway to optimize residual layer thickness and to reduce non-uniform resist distribution near shot-edge regions. These corrections are particularly critical for NIL, where residual layers and local resist filling behavior significantly impact the fidelity of the imprinted structures [18]. Recent studies on confinement-induced polymer flow resistance further show that nanoscale film thickness can significantly modulate flow mobility and apparent viscosity during cavity filling, reinforcing the need for physics-aware modeling of resist transport in NIL master-template fabrication [19]. Beyond these pattern-transfer considerations, maintaining template fidelity also requires accurate monitoring of alignment errors and subtle defect formation across multi-step imprinting. Advances in optical metrology—from Mie-resonator–base[15,16d overlay sensing to FDTD-driven inspection-stack optimization—have substantially improved sensitivity to overlay shifts and line-cut defects, thereby strengthening NIL process control and complementing OPC in ensuring consistent template quality [20,21].

In the OPC perspective, the master template is not a direct replica of the design layout but an pre-distorted pattern that compensates for process-induced distortions [15,16]. By introducing shape adjustments such as edge biasing, OPC ensures that, after resist development and etching, the final wafer structures produced more closely resemble the target design. Although the technique has been widely studied in optical lithography, its systematic integration into NIL template fabrication remains relatively underexplored [17].

A number of process variability sources complicate master template fabrication. First, resist spin-coating typically yields non-uniform thickness distributions, which modify local exposure response during EBL [3]. Second, the EBL process suffers from the proximity effect, wherein electron forward scattering, backscattering, and secondary electrons cause unwanted dose accumulation and pattern blur [4,7,22]. Third, RIE introduces anisotropy and microloading effects that distort sidewall angles and create etch depth non-uniformities [3,7]. Together, these error sources degrade dimensional control and pose critical challenges to NIL master template reliability.

To overcome these difficulties, this research constructs an integrated physics-based simulation framework that connects three multiphysical models: resist spin-coating, EBL exposure and development, and anisotropic RIE transfer. Through the coupling of these steps, the simulation captures the evolution of upstream variations—such as resist thickness non-uniformity—to downstream processes. While no experimental dataset is utilized here, the results are compared to known literature trends to verify their feasibility [23,24].

The proposed approach is part of the broader objective of integrating OPC optimization into predictive simulation frameworks. The integration is expected to minimize expensive fabrication iterations, offer early feedback during mask layout design, and eventually enhance wafer-level fidelity in NIL fabrication [4] [25].

1.2. Three Models Simulation and Parameter Setting

The fabrication of high-quality and durable master templates in NIL requires a judicious balance of material selection and process parameter optimization. Substrate and resist selection are of paramount importance, as they directly dominate dimensional stability, etch durability, and the reproducibility of nanoscale features. Silicon (Si) substrates have traditionally been the industry standard due to their exceptional mechanical strength, thermal stability, and sub-nanometer surface flatness that routinely achieves an RMS roughness of less than 1 nm. These qualities afford a robust platform for high-resolution patterning, in addition to guaranteeing compatibility with following EBL and RIE steps [26,27,28,29,30]. These qualities are also essential for master templates, where surface integrity and long-term durability are requisites for large-scale NIL replication [31,32].

On this foundation, the resist material determines the achievable resolution and the fidelity of pattern transfer. In this work, we used the positive-tone resist CSAR 62 AR-P 6200.09. This chemically semi-amplified resist, based on poly(α-methylstyrene-co-methyl chloroacrylate), has seen widespread use in advanced EBL as it offers a good compromise between high sensitivity, near-10-nm resolution capability, and high etch resistance [33,34]. Importantly, CSAR 62 is a practical alternative to the more common PMMA (polymethyl methacrylate) resist while provides higher sensitivity and better line edge definition. Spin-coating at 4000 rpm, it provides uniform films with thicknesses ranging between 180–200 nm, a thickness range that optimizes fine resolution and yet provides enough robustness to function as an RIE mask. Strong resist-to-silicon adhesion also aids in dimensional control while high etching selectivity (up to 10~16:1) enables accurate transfer of nanoscale features with minimized edge roughness [35].

It is worth noting that in industrial NIL practice, the combination of Si substrate and CSAR 62 resist has become a mature and frequently applied configuration. This system ensures both compatibility with large-scale semiconductor workflows and reliability in repeated template fabrication. Accordingly, many of the parameters implemented in the following simulations—covering spin-coating thickness modeling, electron beam exposure, and dry etching behavior—were derived directly from the characteristics of CSAR 62. The alignment of simulation inputs with well-established process materials ensures that the models not only capture essential physical mechanisms but also retain strong applicability to manufacturing environments [35].

(1) 

Spin coating

Resist thickness control prior to exposure conclusively influences dose absorption, side-wall profile, and critical-dimension uniformity (CDU) in EBL. Classical studies by Emslie et al. and Meyerhofer describe the initial $h\propto {\omega }^{-{\displaystyle \frac{1}{2}}}$ thinning of a Newtonian, non-evaporating film, whereas subsequent one-dimensional models include viscosity enhancement due to solvent evaporation [36,37]. These models emphasize that spin speed ($\omega$) not only controls average film thickness but also imposes radial non-uniformities including an elevated center, a mid-radius depression, and an edge bead. Such features, if left uncorrected, will propagate through to the exposure and development steps, then perturbing local effective dose and line-width control.

For practical use, it is important to anticipate these variations in the context of OPC and mask-layout design. Although OPC has traditionally accounted for optical and stochastic blur, the addition of spin-induced topography is another way to improve across-wafer critical-dimension (CD) control. A predictive spin-coating model thus forms the foundation for combining resist topography with dose modulation in later simulation steps, and hence preserves lithographic fidelity in the 10–20 nm class of features.

(2) 

E-beam lithography (EBL)

Expanding upon the precise spin-coating simulation, the EBL process greatly benefits from better resist thickness maps. It has been widely known that resist thickness variations can affect CDU directly. For example, an IMEC/ASML study on EUV via patterning demonstrated that increasing resist thickness (e.g., from 40 to 60 nm) reduced post-development and post-etch CDU variation when the nominal dose was adjusted accordingly [38]. Thus, the previous spin-coating model not only provides nominal resist thickness but also generates spatial thickness variation normalized maps,

$$\widehat{h}(x,y)={\displaystyle \frac{h(x,y)-\overline{h}}{\overline{h}}}$$

(1)

These parameters are then implemented in the EBL simulation in terms of dose-to-clear correction factors. Such implementation ensures that the proximity effect modeling does include both electron scattering and resist geometry and does provide improved CDU prediction.

EBL has long been recognized as a versatile high-resolution lithographic technique capable of exposing features well below 10 nm [39]. However, its resolution is severely limited by the proximity effect. Forward scattering in the resist and backscattering within the substrate redistributes the deposited energy, causing feature broadening and edge distortion. If left uncorrected, these effects widen linewidths and distort feature edges [40]. Substantial efforts have gone into correcting the proximity effect, ranging from empirical dose modulation [41,42] to fast algorithmic methods such as multipole expansions and cloud acceleration for large-scale patterns [43]. Such corrections have been instrumental in allowing demanding applications such as the fabrication of suspended masks and nanofocusing optics [42,44].

In this work, EBL is simulated not in isolation but within a multi-step process (spin-coating, EBL, and RIE). The EBL module combines beam-column dynamics with a hybrid proximity kernel and resist thresholding, achieving a physics-based link between the focusing behavior of the electron column and the final developed contour. With incorporation of outputs from spin-coating and implicit beam-propagation simulation, this methodology provides a better coherent and predictive groundwork for master template fabrication compared to stand-alone exposure models.

(3) 

Dry etching

In modeling the dry etching stage of NIL, accurately capturing the evolving material interface under directional ion bombardment presents significant computational challenges. The etching process is inherently geometry-dependent and anisotropic, with profile evolution governed by both the local surface orientation and material-specific etch rates. Addressing these complexities requires numerical techniques that can robustly track moving interfaces, handle sharp material boundaries, and incorporate directional etch physics [45].

The level set method has emerged as a powerful tool for interface tracking, as it implicitly represents complex surface geometries and accommodates topological changes naturally [45,46,47,48]. To simulate the directional evolution of material interfaces during plasma etching, this work adopts a level set formulation governed by a Hamilton–Jacobi type equation, allowing interface motion to be driven by local etch rates and geometric gradients [47,48]. To resolve spatial gradients needed for normal vector estimation and interface evolution, central difference schemes are commonly used in interior regions, while forward and backward difference schemes provide stability and accuracy near domain boundaries [49,50]. Temporal advancement of the level set function can be efficiently achieved through explicit Euler integration, which offers simplicity and computational efficiency [48,51]. Additionally, angular-dependent etch rate models allow the etching velocity to be modulated by the orientation of the surface relative to the incident ion flux, enabling anisotropic behavior to be accurately reproduced [48,52,53].

By integrating these physics-informed numerical components into a unified simulation framework, this work provides a robust and scalable approach to modeling reactive ion etching. The combined methodology enables accurate prediction of etch profiles across layout-rich domains, offering significant advantages in terms of geometric flexibility, directional fidelity, and suitability for full-process NIL simulation.

1.3. Contribution and Scope of This Work

The scope of this work is to construct and verify a modular, end-to-end simulation chain that (i) presents first-principles transport and interface physics within each module, (ii) propagates stage-to-stage variations without ad-hoc post-biasing correction, and (iii) provides figures of merit in common with OPC and process control. Specifically, we focus on master-template fabrication on silicon substrate with a CSAR-62 resist stack. The three modeled stages are resist spin-coating, EBL exposure and development, and anisotropic RIE with ARDE terms.

The test pattern set is designed to be limited to isolated and dense patterns with orthogonal layout. The reason reflects standard OPC practice: verifying model fidelity on these two bounding cases is sufficient to expose density-driven bias, proximity blur, and transport-limited etch behaviors without confounding shape factors. In practical terms, if the simulation chain responses the expected trends and metric in isolated and dense patterns, this enables applying a similar set of parameters to more complex pattern layouts with greater confidence as more details are layered. When building OPC models and recipes, the proposed scope helps to distinguish between density and geometry-related corrections, reduce the number of required experiments, and shorten iteration time between dose/focus/etch settings and mask pre-distortion updates.

2. Materials and Methods

2.1. Materials and Process Parameters

In this study, the modeling of the master template was simulated using materials and process parameters widely recognized in industrial NIL processes, ensuring both practical applicability and compatibility with established manufacturing standards. The substrate material used for the master template was a standard Si substrate. Silicon substrates are most desirable due to their high mechanical strength and excellent surface flatness. Typically, the RMS surface roughness is less than 1 nm. In addition, silicon offers outstanding compatibility with both EBL and RIE processes [26,27,28,29,30]. The choice ensures dimensional stability and surface integrity during the fabrication process, while also supporting high-resolution patterning and robust etch resistance. These features comprise essential specifications for durable NIL master templates [31,32].

For the resist model, CSAR 62 AR-P 6200.09 was selected as the electron beam resist. This chemically amplified resist, based on poly(α-methylstyrene-co-methyl chloroacrylate), is designed to deliver high sensitivity and excellent resolution, which makes it suitable for advanced EBL applications. By spin coating at 4000 rpm, AR-P 6200.09 forms a uniform film with a typical thickness of 180 to 200 nm. This thickness range provides a trade-off between high-resolution pattern definition with sufficient durability for use as an RIE etch mask [33,34]. Silicon substrate and CSAR 62 AR-P 6200.09 resist compatibility is well developed. The Si surface, being smooth and inert, promotes uniform resist adhesion and development, while the resist provides high Si to resist etch selectivity, up to 10~16:1. This ensures precise pattern transfer and minimal line edge roughness during RIE. As a result, this material system facilitates reproducible and high-fidelity nanostructure definition as needed in NIL template fabrication [54].

In order to facilitate controlled OPC-oriented comparison, the test layouts included two canonical pattern classes: isolated and dense. Both isolated and dense layouts consist of axis-aligned contours arranged without overlap on a grid, with uniform inter-feature spacing throughout, as indicated in Figure 1. The isolated pattern employed pitches sufficiently large to exclude inter-feature interaction, while the dense pattern used pitches close to the practical resolution limit to emphasize proximity and etch-loading effects. Using two extreme layout cases—isolated and dense—suffices to enable a clean comparison of comparing EBL proximity behavior, RIE transfer fidelity, and downstream OPC responses. The design layouts are currently adopted by semiconductor manufacturers, thereby increasing the practical applicability and validity of the study.

2.2. Spin-Coating Thickness Model Simulation

The spin-coating step determines the initial resist thickness distribution, which has a profound impact on the fidelity of subsequent EBL processes. The thickness non-uniformity causes non-uniform exposure conditions since the local electron scattering and secondary electron yield depend on the thickness. Consequently, line edge roughness and pattern fidelity can be degraded if the resist thickness profile is not well modeled and compensated. Therefore, the inclusion of spin-coating impacts in the OPC design phase is essential to ensure robust downstream pattern transfer.

From the transport phenomenon perspective, the resist film deposited by spin-coating can be modeled as a thin Newtonian liquid layer under radially spreading conditions during rotation. The governing forces include centrifugal force, which drives outward flow; viscous shear stress, which provides resistance to radial stretching; capillary forces, which dominate near the wafer edge and generate the edge bead; and solvent evaporation, which reduces film thickness uniformly over time.

Treating the spin-coated resist as a Newtonian liquid under axisymmetry, the depth-averaged radial flux $q_r$ is obtained from the radial momentum balance in the lubrication approximation:

$$q_r={\int }_0^h{u}_{r}dz={\displaystyle \frac{h^3}{3\mu }}\left(\rho {\omega }^{2}r-{\partial }_{r}p\right)$$

(2)

where $h=h(r,t)$ is the local film thickness as a function of radial position $r$ and time $t$, the viscosity $\mu$, the density $\rho$, the spin speed $\omega$ ($rad{s}^{-1}$), and the pressure $p$. Substituting this constitutive relation into mass conservation under axial symmetry yields

$${\displaystyle \frac{\partial h}{\partial t}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r{q}_r\right)=-E$$

(3)

where $E$ represents the evaporation rate. Having symmetry with respect to angle $\theta$, only the radial direction $r$ remains; ${\displaystyle \frac{\partial }{\partial \theta }}=0$. In the wafer interior (away from the edge), curvature effects are weak so that ${\partial }_{r}p\cong 0$, so it follows the classic ODE for thickness,

$${\displaystyle \frac{\partial h}{\partial t}}=-{\displaystyle \frac{2\rho {\omega }^2}{3\mu }}{h}^3-E$$

(4)

which recovers the classical Emslie–Meyerhofer spin-speed scaling for the “core” thickness once the evaporation term is accounted for. The boundary conditions are symmetry at the wafer center ($r=0$, ${\displaystyle \frac{\partial h}{\partial r}}=0$) and no flux at the wafer edge ($r=R$), where bead buildup is accounted for by the Bornside correction. These fundamental considerations form the theoretical basis from which empirical and semi-empirical relations are obtained.

Collecting these results, the thickness profile $h\left(r,\omega \right)$, which represents the local resist thickness as a function of radial position $r$ and spin speed $\omega$, is written as the sum of a core term and a capillary edge-bead correction:

$$h\left(r,\omega \right)=h_c\left(\omega \right)\left[1+a_2{\left({\displaystyle \frac{r}{R}}\right)}^2+a_4{\left({\displaystyle \frac{r}{R}}\right)}^4\right]+h_e\left(\omega \right)exp\left[-{\displaystyle \frac{{\left(R-r\right)}^2}{2w_e{\left(\omega \right)}^2}}\right]$$

(5)

where the three functional components scale with spin speed according to:

$$h_c\left(\omega \right)=A{\left({\displaystyle \frac{\omega }{{\omega }_0}}\right)}^{-n}$$

(6)

$$h_e\left(\omega \right)=h_{e0}{\left({\displaystyle \frac{\omega }{{\omega }_0}}\right)}^{-m}$$

(7)

$$w_e\left(\omega \right)=w_{e0}{\left({\displaystyle \frac{\omega }{{\omega }_0}}\right)}^{-l}$$

(8)

In this case, ${\omega }_0$ is a reference spin speed, and $\left(n,m,l\right)$ are the scaling exponents for central thickness, edge height, and edge-bead width, respectively. This setup combines the Emslie–Meyerhofer spin-speed relation with Bornside-type edge-bead corrections, ensuring fidelity at both the global and local levels.

The first term in Eq. (5) corresponds to interior (“core”) thickness driven by the balance between centrifugal forcing and viscous dissipation, which—together with evaporation—gives rise to the classic Emslie–Meyerhofer relation for the dependence on the spin speed $h_c\left(\omega \right)\propto {\omega }^{-n}$ (with $n\approx 0.5-0.7$ depending on the type of evaporation). The polynomial term ${\displaystyle \frac{r}{R}}$ in reflects finite-radius corrections to the almost nearly uniform core. Closer to the edge of the wafer, capillarity becomes non-negligible: the curvature-induced pressure $p_{\gamma }~\gamma {\nabla }^{2}h$ introduces a capillary pressure gradient that competes with centrifugal and viscous stresses, creating a boundary-layer-like region near the edge bead. Rather than re-solving the full lubrication PDE with curvature, a Bornside semi-empirical correction is adopted with a localized envelope $h_e\left(\omega \right)exp\left[-{\displaystyle \frac{{\left(R-r\right)}^2}{2w_e{\left(\omega \right)}^2}}\right]$, where $h_e\left(\omega \right)$ is the edge-bead amplitude and $w_e\left(\omega \right)$ its characteristic width. Both parameters inherit spin-speed scaling from the underlying capillary–viscous–centrifugal balance.

For exposure coupling, the local resist thickness is normalized and included in the electron-beam development threshold as

$$T\left(x,y\right)=T_0\left[1+\alpha \widehat{h}\left(x,y\right)\right]$$

(9)

where $T_0$ is the threshold of base development, α is a sensitivity coefficient, and $\widehat{h}\left(x,y\right)={\displaystyle \frac{\left(h\left(x,y\right)-\overline{h}\right)}{\overline{t}}}$ is the spatially normalized thickness map. This couple affect ensures that thickness differences in absorption affect local dose deposition and CD variability directly.

Model validation

The validity of the proposed spin-coating thickness formulation was assessed by comparison between the simulated results and both widely adopted literature correlations and manufacturers’ specifications. Figure 2 shows the correlation between center thickness and spin speed given by the spin curve. With a spin rate of 4000 rpm, the simulated thickness matches the datasheet value of 0.20 μm of AR-P 6200.09 as reported by Allresist [55]. The proposed power-law index of $n=0.54$ in our study falls within the conventional Meyerhofer evaporation-dominated range (0.5–0.7) [37], and also agrees with the viscous flow scaling relation of the type $h\propto {\omega }^{-{\displaystyle \frac{1}{2}}}$ by Emslie–Bonner–Peck [36] that prevails at the onset of spinning. This alignment with both empirical data and classical scaling theories demonstrates that the model closely approximates the expected speed-thickness relation over a very wide range of operations.

Figure 3 presents the radial thickness profile at 4000 rpm, which has been checked against the one-dimensional solution of Bornside et al. [56]. Least-squares calibration of the polynomial coefficients provides $a_2=-0.05$ and $a_4=-0.02$, which reduces the across-wafer RMS error to 3.8% with respect to the Bornside numerical profile. The mid-radius depression and the Gaussian-type thickening at the wafer edge have been modeled successfully, and the parameterization of the edge bead peak at approximately $r\approx 46mm$ has been derived by the parameter set $\left\{h_{e0},m,{\omega }_{e0},l\right\}=\left\{0.02\mu m,0.60,2.0mm,-0.25\right\}$.The parameterization agrees with the 5–10% excess thickness at the wafer edge that has been determined by the FAB.nano application note [57].

We should note that minor deviations between the simulation outputs and experimental data are typically discernible in actual spin-coating processes due to the uncontrollable dynamics of the edge beads and the wafer–chuck interface boundary conditions. These can be corrected by Bayesian inference or Kalman filtering–based data assimilation frameworks, which allows the reconciliation of the process measurements and the physical principle-based models. In practice, the region of minimum deviation is selected for following exposure coupling; in this work, a representative window of 600 nm × 600 nm was chosen in Section 3.1 to provide reliable calibration and smooth transition to the subsequent EBL simulation stage.

2.3. E-Beam Lithography (EBL) Simulation

The resist exposure during EBL is precisely governed by the movement of charged particles within the column and their corresponding interaction of resultant character with the resist and the substrate. In this case, electron dynamics are defined through an explicit employment of the Lorentz force equation,

$${\displaystyle \frac{d}{dt}}\left(m_e\stackrel{⃑}{v}\right)=q\stackrel{⃑}{E}+q\left(\stackrel{⃑}{v}\times \stackrel{⃑}{B}\right)$$

(10)

where $m_e$ and $q$ are the electron mass and charge, respectively. In practice, the magnetic field inside the accelerating region is negligible. This approximation is justified because the accelerating gap sustains MV/m-scale electric fields while residual magnetic flux is strongly suppressed by μ-metal shielding (high-permeability alloy), making the magnetic-to-electric force ratio ${\displaystyle \frac{q\left(\stackrel{⃑}{v}\times \stackrel{⃑}{B}\right)}{q\stackrel{⃑}{E}}}={\displaystyle \frac{\left(\stackrel{⃑}{v}\times \stackrel{⃑}{B}\right)}{\stackrel{⃑}{E}}}\ll 1$. For 30-keV electrons (non-relativistic approximation$:v\approx 0.33c=9.9\times {10}^7\mathrm{m}{\mathrm{s}}^{-1}$) and typical fields $E={\displaystyle \frac{V}{L}}\approx \left(3~10\right)\times {10}^{6}V{m}^{-1}$, the equivalent magnetic field is $B_{eq}={\displaystyle \frac{E}{v}}\approx 30~100mT$, whereas stray magnetic field in the shielded gun region is far below this threshold $B_{eq}$; thus, magnetic contributions in the accelerating region can be neglected to first order. Consequently, the Lorentz equations simplify to a coupled system driven solely by the radial and axial electric fields:

$${\displaystyle \frac{d{v}_r}{dt}}={\displaystyle \frac{q}{m_e}}{E}_r$$

(11)

$${\displaystyle \frac{d{v}_z}{dt}}={\displaystyle \frac{q}{m_e}}{E}_z$$

(12)

These fields come from the structure of the electrodes: the anode accelerates electrons with an electrode voltage of $V_{anode}=30kV$, and the focusing electrode applies $V_{focus}=10kV$ along a column of length $L_{column}=10mm$. The conversion between electrode voltage $V$ and electric field components ${(E}_r,E_z)$ can be directly obtained from electrostatics mathematically represented as $E=-\nabla V$. For application in this simplified model of a column, every electrode can be represented as a localized potential source,

$$V(r,z)~{\displaystyle \frac{V_0}{\sqrt{{(z-z_0)}^2+r^2}}}$$

(13)

where $V_0$ can be either the anode voltage $V_{anode}$ or a focusing voltage $V_{focus}$. A derivative of this potential with respect to $r$ and $z$ then generates the radial and vertical components of this electric field used within this trajectory integration process. Physically, this process corresponds to Lorentz force law considering a charge on an electron to be $q=-e$. By substituting this negative charge within Eq. (11) (12), we guarantee a positive voltage on an electrode will create an attractive force upon an electron beam to reproduce the actual behavior of electron optical columns.

Trajectories are also set up at a finite source size and angular divergence ($r_{std}=10\mathrm{n}\mathrm{m}$, ${\theta }_{std}=0.6mrad$) by optimization in Section 3.2, representing the physical spread of the electron at emission. The integration procedure is based on the Dormand–Prince 5(4) method that assures accurate time-resolved tracking of focusing behavior over nanosecond scales [58].

The trajectory set also offers simultaneously a description of the beam focusing and a quantitative characterization of the effective spot size at the resist plane. The landing radius set $\left\{r_i\right\}$ is transformed into a Gaussian-equivalent width using a robust median-based estimator,

$${\sigma }_{traj}=1.4826MAD\left\{r_i\right\}$$

(14)

Here, MAD denotes the median absolute deviation. This approach creates a direct mapping between the beam blur used in exposure calculations and the related column physics that governs particle transport. A wide focusing voltage (5~25kV) adjustment in Section 3.2 is performed to minimize ${\sigma }_{traj}$​ and thus practically mimicking the autofocusing action of a real EBL system.

Having established the forward beam width, the exposure is modeled by rastering the pattern at a fixed step size $s=2nm$. Any $s\times s$ pixel tile that has an exposed area greater than a specified fraction of $\rho =0.25$ receives a charge shot,

$$Q_{shot}=D·s^2$$

(15)

where $D=120\mu C/{cm}^2$ is the nominal base dose. Superposition of Gaussian spots of width ${\sigma }_{traj}$gives us the fundamental exposure distribution $E_0(x,y)$, accounting for the influence of a finite beam size.

To simulate those scattering phenomena, the exposure map is convolved with a hybrid kernel that consists of two physically distinct components: a Gaussian $G_{\sigma }$ that models the tightly concentrated primary beam, and a squared-Lorentzian $L_2$ that accounts for the slow drop of the tails due to the backscattered electrons. The modeling expression for the hybrid kernel is represented by

$$G_{\sigma }\left(x,y\right)=A_1·exp\left(-{\displaystyle \frac{{\left(x-x_0\right)}^2+{\left(y-y_0\right)}^2}{2{\sigma }^2}}\right)$$

(16)

$$L_2\left(x,y\right)={\displaystyle \frac{A_2}{1+\frac{{\left(x-x_0\right)}^2+{\left(y-y_0\right)}^2}{{\gamma }^2}}}$$

(17)

$$K_{hybrid}\left(x,y\right)={{\omega }_{1}G}_{\sigma }\left(x,y\right)+{\omega }_2{L}_2\left(x,y\right)$$

(18)

with weights ${\omega }_1={\displaystyle \frac{A_1}{A_1+A_2}}$ and ${\omega }_2={\displaystyle \frac{A_2}{A_1+A_2}}$. In the present work $A_1=0.74,A_2=0.26$, and the Lorentzian scale is $\gamma =45nm$. The convolution,

$$D_{map}\left(x,y\right)=\left(E_0\ast K_{hybrid}\right)\left(x,y\right)={\sum }_{m=-\infty }^{\infty }{\sum }_{n=-\infty }^{\infty }{E}_0\left(m,n\right)·K_{hybrid}\left(x-m,y-n\right)$$

(19)

generates the effective dose distribution, whereby at the same time short-range forward scattering and long-range proximity effects are accounted for.

The development of the resist is modeled by a modulated global dose threshold from the previous model of spin-coating. The threshold is given by the following

$$T\left(x,y\right)=\eta D_{max}\left[1+\alpha \widehat{h}(x,y)\right]$$

(20)

where $\widehat{h}\left(x,y\right)={\displaystyle \frac{\left(h\left(x,y\right)-\overline{h}\right)}{\overline{t}}}$ is normalized thickness map after normalization, $\eta =0.56$ is the base threshold ratio, and $\alpha =0.2$ quantifies the sensitivity to thickness variation. This expresses the physical truth that slightly higher exposures are needed for thicker resist areas in order for them to fully clear. The final resist pattern is described as

$$P_{dev}\left(x,y\right)=1\left\{D_{map}\left(x,y\right)\ge T\left(x,y\right)\right\}$$

(21)

This simulation pipeline thus integrates column electrostatics, beam-size proxies from trajectories, scattering-induced proximity effects, and solubility of the resist into one physically valid model. The results—spatially resolved dose maps, proximity-corrected exposure, and developed contours—become direct inputs for following reactive ion etching simulation and NIL template pattern transfer.

2.4. Reactive Ion Etching (RIE) Simulation

To model the profile evolution during the dry etching stage, a physics-informed level set method is employed to simulate anisotropic RIE. This stage constitutes the third phase of the proposed three-stage NIL simulation framework, following EBL patterning. At this point in the process, a patterned photoresist layer has been transferred onto the substrate, and both materials are subjected to plasma etching to further define the desired topography. During RIE, vertical etching is dominated by directional ion bombardment; lateral erosion arises primarily from ion-assisted effects due to finite angular spread, with an additional minor omnidirectional chemical contribution. The simulation captures this dual-mode etching behavior by combining an angular-dependent etch rate model with a material-specific selectivity function. This allows for accurate prediction of profile evolution, accounting for the fact that photoresist masks remain largely unaffected during etching, while underlying substrate materials are etched based on local surface orientation and exposure.

The interface between photoresist and substrate is represented using a level set function $\varphi \left(x,y,z\right)$, where in Figure 4 the zero-level set $\varphi =0$ defines the material surface. In this formulation, $\varphi >0$ denotes the vacuum (air) region, while $\varphi <0$ corresponds to material regions, including both photoresist and substrate. The level set evolves in time according to the following Hamilton–Jacobi equation:

$${\displaystyle \frac{\partial \varphi }{\partial t}}+R_{etch}\left|\nabla \varphi \right|=0$$

(22)

Here, $R_{etch}$ is the total etch rate at the interface and $\left|\nabla \varphi \right|$ is the magnitude of the spatial gradient of $\varphi$, ensuring proper geometric propagation along the normal surface. To account for material-dependent etch selectivity, photoresist is assigned a finite selectivity of 0.10, yielding a vertical etch rate that is one-tenth of the substrate’s and thus only minor mask erosion. In contrast, substrate regions undergo directional and chemical etching as described by the anisotropic etch model.

The etch rate is computed as a combination of directional (ion-driven) and isotropic (chemical) components:

$$R_{etch}=R_{base}\left(\left|n_z\right|+r_{lat}\times \left(1-\left|n_z\right|\right)+f_{chem}\right)$$

(23)

where $R_{base}$ denotes the intrinsic basic material-dependent vertical etch rate. The term $\left|n_z\right|$ represents the cosine of the angle between the local surface normal and the vertical direction, thereby quantifying how aligned the surface is with the incoming ion flux. The parameter $r_{lat}$ introduces lateral etching behavior by allowing reduced etch rates on sloped or horizontal surfaces, mimicking the effect of angular broadening due to off-normal ion incidence and chemical diffusion. Finally, $f_{chem}$ accounts for the isotropic chemical etching contribution, treated here as a dimensionless scaling factor added to the angular term and modulated by $R_{base}$.

The surface normal vector $\stackrel{⃑}{n}=\left(n_x,n_y,n_z\right)$ is computed from the normalized gradient of the level set function:

$$\nabla \varphi =\left({\displaystyle \frac{\partial \varphi }{\partial x}},{\displaystyle \frac{\partial \varphi }{\partial y}},{\displaystyle \frac{\partial \varphi }{\partial z}}\right)$$

(24)

$$\stackrel{⃑}{n}={\displaystyle \frac{\nabla \varphi }{\left|\nabla \varphi \right|}}$$

(25)

To approximate the spatial derivatives of the level set function, a central difference scheme is employed at interior grid points. This method, derived from a symmetric Taylor expansion, achieves second-order accuracy, meaning that the associated truncation error decreases proportionally to $\mathcal{O}\left({∆x}^2\right)$ as the spatial grid is refined. The discrete gradient at an interior point is computed as:

$${\left.{\displaystyle \frac{\partial \varphi }{\partial x}}\right|}_i\approx {\displaystyle \frac{{\varphi }_{i+1}-{\varphi }_{i-1}}{2\Delta x}}$$

(26)

At the computational boundaries, where data from only one direction is available, first-order forward and backward difference schemes are applied. These one-sided methods yield first-order accuracy, with truncation errors that scale linearly as $\mathcal{O}\left(∆x\right)$. Specifically, the boundary derivatives are computed using:

$${\left.{\displaystyle \frac{\partial \varphi }{\partial x}}\right|}_0\approx {\displaystyle \frac{{\varphi }_1-{\varphi }_0}{\Delta x}}$$

(27)

$${\left.{\displaystyle \frac{\partial \varphi }{\partial x}}\right|}_{n-1}\approx {\displaystyle \frac{{\varphi }_{n-1}-{\varphi }_{n-2}}{\Delta x}}$$

(28)

These schemes ensure stable and reasonably accurate estimation of spatial gradients, which are essential for normal vector calculation and interface evolution in the level set framework. In the temporal domain, the level set equation is integrated using an explicit Euler scheme, which advances the interface over time based on the local etch rate and gradient magnitude:

$${\varphi }^{(t+∆t)}={\varphi }^{\left(t\right)}-R_{etch}\left|\nabla \varphi \right|∆t$$

(29)

To justify the experimentally observed speed reduction inside high-aspect-ratio (AR) trenches, the basic level-set update (Eq. (22) ~ Eq. (29)) is supplemented by Aspect Ratio Dependent Etching (ARDE) factors that modulate the vertical and lateral channels as functions of local AR. Here, we explicitly integrated ARDE into the Hamilton–Jacobi scheme and numbered the equations for reference [59]. ARDE is a transport-limited phenomenon in which the local etch rate decreases as a feature becomes deeper and narrower because of reactant depletion, by-product trapping, and partial ion-flux shadowing. We let the interfacial etch rate depend on the local aspect ratio $AR\left(x,y\right)={\displaystyle \frac{d(x,y)}{\omega \left(x,y\right)+\epsilon }}$, where $d(x,y)$ is the etched depth accumulated at the substrate plane, $\omega \left(x,y\right)$ is the local opening width estimated from the ADI via a Euclidean distance transform, and $\epsilon >0$ avoids division by zero. In this manner, the level-set evolution becomes

$${\displaystyle \frac{\partial \varphi }{\partial t}}+R_{etch}(\stackrel{⃑}{n},AR;M)\left|\nabla \varphi \right|=0$$

(30)

with material label $M\in \left\{resist,substrate\right\}$ and unit normal $\stackrel{⃑}{n}={\displaystyle \frac{\nabla \varphi }{\left|\nabla \varphi \right|}}$. For the substrate we write the anisotropic, transport-modified etch rate as

$$R_{etch}=R_{base}\left[f_{ang}\left(\left|n_z\right|\right)f_{acc}\left(AR\right)\left|n_z\right|+r_{lat}{f}_{lat}\left(AR\right)\left(n_x^2+n_y^2\right)+f_{chem}\right]$$

(31)

Eq. (31) outlines the contributions of vertical, lateral, and isotropic-chemical factors while keeping only a single material rate scale $R_{base}$ for dimensional consistency. The angular ion-flux envelope $f_{ang}\left(\left|n_z\right|\right)$, controlling the vertical projection, is modeled as the Gaussian function of the incidence angle:

$$f_{ang}\left(\left|n_z\right|\right)=exp\left[-{\left({\displaystyle \frac{arccos\left(\left|n_z\right|\right)}{\sigma }}\right)}^2\right]\mathrm{where}\sigma >0\left(radians\right)$$

(32)

A smaller $\sigma$ represents a more collimated beam and hence stronger anisotropy; conversely, a larger $\sigma$ captures a wider angular distribution. This selection is convenient to calibrate (the single parameter $\sigma$ maps to angular half-widths obtained from the instrument), numerically smooth (beneficial for stable normal estimations), and multiplicative. Therefore, it composes cleanly with other transport penalties outlined below.

Access limitation due to ARDE is encoded through a monotone penalty,

$$f_{acc}\left(AR\right)={\displaystyle \frac{1}{1+\alpha {AR}^p}},\mathrm{where}\alpha >0,p\ge 1$$

(33)

which reduces the effective vertical contribution as the feature deepens and narrows. Sidewall-dominated transport and passivation are captured by a stronger attenuation on the lateral channel,

$$f_{lat}\left(AR\right)={\left[f_{acc}\left(AR\right)\right]}^{\beta }exp\left(-\gamma AR\right),\mathrm{where}\beta >1,\gamma \ge 0$$

(34)

so that narrow, high-AR trenches experience a disproportionately larger reduction in sideways advance than in vertical descent. In Eq. (21), the lateral geometric factor is written as $n_x^2+n_y^2$, which is equivalent to ${sin}^2\theta$ by the identity ${\left|\stackrel{⃑}{n}\right|}^2=n_x^2+n_y^2+n_z^2=1$. Using $n_x^2+n_y^2$ rather than $1-\left|n_z\right|$ (or a signed $sin\theta$) ensures smoothness and non-negativity, avoids a cusp at $\left|n_z\right|=1$, and weights slopes proportionally according to their squared tilt. These features improve stability and eliminate grid-orientation bias in the three-dimensional updates. An expression for the chemical term is modeled as a dimensionless constant scaling of the base rate, $f_{chem}$, i.e., independent of AR in this implementation. When desired, geometric shadowing in the vertical channel can be incorporated by replacing $f_{ang}$ with $f_{ang}{f}_{shad}\left(AR\right)$, where $f_{shad}\left(AR\right)={\displaystyle \frac{1}{1+kAR}}$ with $k\ge 0$.

With Eq. (30) ~ Eq. (34), inward CD shrinkage arises during the evolution rather than being imposed afterward. As $AR$ grows in a trench, $f_{acc}$ lowers the vertical transport and $f_{lat}$ suppresses lateral advance even more strongly. Integrated over time in the Hamilton–Jacobi update, narrow features thus propagate more slowly both downward and sideways than wide openings. For illustration, Figure 5 shows the simulated 3D final etched topography for isolated design layout under the RIE settings.

In addition to the continuous profile evolution, the computed etch depth distribution is subsequently mapped to a binary after etch inspection (AEI) image. The conversion from a continuous depth field

to a binary mask is performed by operating with a depth-biased threshold criterion:

$$B(x,y)=\left\{\begin{array}{c}1whered(x,y)\ge a{h}_{resist}\\ 0whered(x,y)<a{h}_{resist}\end{array}\right.$$

(35)

where $h_{resist}$ denotes the nominal resist thickness and $a$ is an adjustable scaling factor (typically set to 0.8). This operation ensures that only regions with sufficient penetration across the resist are classified as etched, thus reproducing CD shrinkage observed in experimental AEI profiles.

The resulting AEI contour exhibits a measurable inward bias that is consistent with transport-limited wafer behavior. In practice, $\alpha$ (ARDE strength), $p$ (ARDE nonlinearity), $r_{lat}$ (lateral weighting), $\beta$ (lateral attenuation exponent), $\gamma$ (lateral passivation factor), and $\sigma$ (beam collimation) are calibrated against cross-sectional SEMs so that simulated depth, sidewall angle, and CD loss match process trends without the need for ad-hoc post-processing. Together, the depth-biased thresholding offers a compact yet physically consistent approximation of AEI formation, complementing the physics-based level set evolution and facilitating quantitative analysis of CD shrinkage, line edge roughness, and pattern fidelity degradation.

2.5. Pattern Fidelity Metrics to be Used

To quantify lithographic pattern fidelity produced by simulation, three mutually complementary metrics were applied: the normalized mean squared error (NMSE), the critical-dimension error (∆CD), and the line edge roughness (LER). These indicators present distinct but related views of pattern quality, covering overall image similarity, dimensional accuracy, and edge sharpness, respectively. All three metrics are widely used in OPC schemes of relevance to semiconductor lithography and are directly applicable to both EBL and RIE stages in this study.

(4) 

Normalized Mean Squared Error (NMSE)

The NMSE compares the overall pixel-wise difference between the designed binary pattern and the simulated or developed pattern. It is given by:

$$NMSE={\displaystyle \frac{\frac{1}{N}{\sum }_{x,y}{\left(P\left(x,y\right)-\widehat{P}(x,y)\right)}^2}{\frac{1}{N}{\sum }_{x,y}{P\left(x,y\right)}^2}}$$

(36)

where $N$ is the number of total pixels. This normalized version provides scale invariance and a globally defined contour fidelity measure. In semiconductors, NMSE represents how well the exposure and development process can replicate desired pattern density and hence suits best for analyzing stochastic variations in near-10-nm features.

(5) 

Critical-Dimension Error (∆CD)

The ∆CD quantifies the difference in linewidths from design to simulation on each contour feature. Rather than focusing on overall area differences, this metric directly measures the length of the shortest edge of each contour pattern printed in nanometers. The ∆CD is given by the expression:

$$\mathsf{\Delta }CD=\left|{CD}_{sim}-{CD}_{design}\right|$$

(37)

where ${CD}_{design}$ stands for the nominal target width defined in the design layout, and ${CD}_{sim}$ denotes the measured contour width after the simulation process (EBL or RIE). The design layouts used in this simulation investigation are represented in Figure 1. From a lithographic perspective, ΔCD offers a direct quantification of dimensional accuracy, which serves as a key indicator of process proximity effects and pattern transfer fidelity across exposure and etching stages. While considering OPC, minimizing ΔCD ensures effective contour replication in accordance with design intent and reduces bias arising from local pattern density or exposure-induced effects.

(6) 

Line-Edge Roughness (LER)

LER is used to find out how much edges after processing deviate from their initial straight boundaries. It is computed discrete representation by:

$$LER={\displaystyle \frac{\sum \left|∆P_{edge}\left(x\right)-∆{\widehat{P}}_{edge}\left(x\right)\right|}{\sum ∆P_{edge}\left(x\right)}}$$

(38)

Here, $∆P_{edge}\left(x\right)$ and $∆{\widehat{P}}_{edge}\left(x\right)$ are edge profiles from the design layout and simulated patterns, respectively. This technique offers a feature process-independent edge fidelity measurement that is standardized. From a semiconductor standpoint, LER plays a robust role in transistor performance, line resistance, and device variability, thus forming one of the most significant parameters both in OPC and process qualification.

In the present framework, these three metrics are applied after both the EBL and RIE simulations to benchmark pattern fidelity across different process conditions. NMSE captures global deviations in the aerial or dose map, ΔCD highlights systematic dimensional deviations of printed contours relative to the intended design geometry, and LER provides sensitivity to fine-scale roughness amplified during resist development or etching. Their combined use ensures a balanced evaluation of proximity effects, beam focusing conditions, and downstream transfer fidelity. Importantly, this triad of metrics parallels the evaluation protocols used in OPC workflows, thereby reinforcing the industrial relevance of the physics-based simulation results presented here.

3. Results

3.1. Spin-Coating Check: 600 nm × 600 nm Window

Figure 6 illustrates the simulated resist-thickness map for a 100 mm wafer spun at 4000 rpm, with a color scale spanning 0.190–0.206 μm, indicating a wafer-scale thickness variation on the order of 10–20 nm. The map exhibits the canonical features of a spin-coated film: a raised center, a shallow mid-radius depression, and a ∼7% edge bead, with a maximum center-to-edge deviation of only 0.0077 μm, thereby satisfying the uniformity requirements for advanced lithography.

To enable the following EBL, we further implemented an automated algorithm to identify the optimal 600 nm × 600 nm exposure window. The selection procedure is grounded in quantitative metrics derived from the simulated wafer-scale thickness map. First, the normalized local deviation field $\widehat{h}(x,y)$ was computed in Eq. (1), and the local gradient magnitude was obtained from finite differences,

$$\left|\nabla h\right|=\sqrt{{\left({\displaystyle \frac{\partial h}{\partial x}}\right)}^2+{\left({\displaystyle \frac{\partial h}{\partial y}}\right)}^2}$$

(39)

which captures intra-window slope and thus quantifies local planarity. Both quantities were normalized by their respective 95th-percentile values ${(p}_t,p_g)$ to avoid bias from outliers. The final scoring function was defined as

$$\mathcal{L}\left(x,y\right)={\displaystyle \frac{\left|\widehat{h}\left(x,y\right)\right|}{p_t}}+\beta {\displaystyle \frac{\left|\nabla h(x,y)\right|}{p_g}}$$

(40)

with $\beta =1.0$ assigning equal weight to thickness, uniformity and gradient suppression. Candidate points within a radial mask ($r\le R-3w_e$) were then evaluated, where $w_e$ is the edge-bead half-width. The global minimum of $\mathcal{L}\left(x,y\right)$ was chosen as the window center.

This procedure yielded an optimal window at $\left(x_0,y_0\right)\approx (25.088mm,-17.984mm)$ from the wafer center, i.e. near the 25 mm radius, consistent with the minimum-deviation zone predicted by radial analysis. The extracted local sub-map (600 nm × 600 nm, interpolated from the radial profile) shows mean thickness $h=0.196\mu m$ with relative deviation $\left|\widehat{h}\right|<8.5\times {10}^{-4}$, corresponding to <0.1 nm variation across the tile. The associated gradient magnitude is below ${1\times 10}^{-4}\mu m/mm$, confirming negligible slope.

In summary, the 600 nm × 600 nm window selection combines statistical normalization, gradient minimization, and edge-avoidance criteria to guarantee that the chosen EBL field is located in the most uniform tile of the wafer. This not only validates the spin-coating model at the sub-micron scale but also demonstrates how data assimilation principles can guide optimal region selection for high-precision exposure.

It is important to note that, in this first-stage demonstration, thickness-driven CDU is intentionally excluded from the EBL and RIE analyses. The optimal window is deliberately chosen to isolate exposure and etch physics from thickness-induced dose-to-clear variation. This design choice allows the present work to focus on establishing the physics-coupled spin → EBL → RIE framework, while the incorporation of non-uniform thickness fields—and their quantitative impact on CD and LER—will be addressed in subsequent calibration studies.

3.2. Trajectory-Derived Spot Size and Pattern Fidelity

The electron–optical blur controls the deposit of energy in the resist and as well as the extent of spatial coupling between nearby features. Although the forward scattering determines the fundamental width of each shot, long-range backscattering redistributes part of the dose over micrometer scales and is typically represented by the hybrid point spread function (PSF). In this study, the effective spot size ${\sigma }_{traj}$ is obtained from Lorentzian trajectory tracing and then used to compute two dose fields for each layout: the exposure map (the forward-broadened dose delivered by the raster) and the proximity-corrected dose calculated by convolving the exposure with an area-normalized hybrid PSF $K\left(x,y\right)$, which has the property $\iint K\left(x,y\right)dxdy=1$.

This normalization ensures field-integral conservation:

$$\sum D_{exposure}\approx \sum D_{proximity-corrected}$$

(41)

with variations by approximately 1–2% in our grids due to boundary padding. Peak values in the proximity-corrected maps are thus reduced relative to the exposure maps, not due to under-dose, but due to the fact that spatial averaging distributes the same total energy over an expanded footprint. Interestingly, this integral check avoids unit and scale ambiguity when making comparisons between layouts with varying fill factors.

To illustrate the impact of ${\sigma }_{traj}$ on these fields, three representative settings are analyzed with fixed source statistic of $r_{std}=10nm$ (see Table 1). In the lower-bound case (${\sigma }_{traj}=0.156nm$), the exposure lines are highly confined, but the proximity-corrected maps indicate weak halos around isolated features; the crosstalk between neighboring lines is kept low. With typical operating conditions (${\sigma }_{traj}=4.031nm$), the exposure footprints undergo minor expansion, and the proximity background rises to a moderate base level, especially in dense regions where local superposition is stronger; even so, the integral conservation keeps the total dose budget unchanged. This balance between forward localization and modest proximity spreading indicates that ${\sigma }_{traj}$ in the 2–5 nm range provides the most consistent match to experimental EBL performance, making the 4.031 nm case representative of optimal process fidelity. The over-focused state (${\sigma }_{traj}=39.075nm$), the exposure ridges are wide, and the proximity base level becomes virtually uniform over the feature block, signaling proximity-dominated energy delivery and, in practical terms, an enhanced risk of line bridging in dense patterns.

From an engineering perspective, the analysis indicates that ${\sigma }_{traj}$ within the low-single-nanometer range of 2~5 nm achieves a balance between localization and overlap: isolated layouts preserve contrast while dense layouts mitigate the risk of significant base level accumulation. That being said, the specific range will depend on the resist stack and substrate backscatter. The integral conservation diagnostics employed in this context offer a straightforward quantitative assessment to ensure that the correlation between simulation results and process-related dosage remains consistent.

The influence of beam spot size derived from Lorentzian trajectory tracing on pattern fidelity was quantified. Three representative conditions were selected from the parameter sweeping, all maintaining a standard deviation of the initial radial position distribution with ($r_{std}$) equal to 10 nm: a numerical lower bound (${\sigma }_{traj}=0.156nm$), a practical operating condition (${\sigma }_{traj}=4.031nm$), and an over-focused state (${\sigma }_{traj}=39.075nm$).

The corresponding exposure map and proximity-corrected dose for isolated and dense layouts are shown in Figure 7, Figure 8 and Figure 9. To avoid unit and scale ambiguity, the proximity-corrected dose is obtained by convolving the exposure map with an area-normalized hybrid PSF $\iint K\left(x,y\right)dxdy=1$. Consequently, the field integral is conserved: $\sum D_{exposure}\approx$ $\sum D_{prox}$ (within ~2% due to boundary padding), while peak values are lower because of spatial averaging.

In the numerical lower-bound condition (${\sigma }_{traj}=0.156nm$), the final contours appear continuous because the global threshold integrates sufficient local dose to connect features (Figure 7). Yet the metrics reveal poor geometric fidelity: the isolated array exhibits $NMSE=0.846$ and the dense layout shows a very high $LER=11.19$. This combination indicates overly sharp, pixel-level transitions that the resist model interprets as jagged edges—a behavior inconsistent with real EBL where electron scattering and resist chemistry impose finite blur. Hence, ${\sigma }_{traj}<1nm$ should be regarded as a numerical limit: the pattern can look connected, while the edge quality degrades and the result lacks experimental validity.

The realistic condition (${\sigma }_{traj}=4.031nm$) balances primary-beam localization with modest overlap (Figure 8). Here the ΔCD are minimal in both layouts ($\mathsf{\Delta }\mathrm{C}\mathrm{D}=0.382$ isolated; $\mathsf{\Delta }\mathrm{C}\mathrm{D}=0.382$ dense), and the dense-pattern LER drops to 3.78. Although the pixel-wise NMSE is near unity, this reflects small, distributed differences across the field rather than dimensional bias; the ΔCD and LER jointly indicate that this spot size best preserves the designed linewidth and edge continuity. Practically, this window (about 2-5 nm) aligns with reported high-resolution EBL performance and yields the most transferable predictions for NIL master fabrication.

For over-focus condition (${\sigma }_{traj}=39.075nm$), the broadened Gaussian-like spot allows for dose accumulation over neighboring shots (Figure 9). The isolated pattern greatly increases in size, leading to a strong dimensional bias ($\mathsf{\Delta }\mathrm{C}\mathrm{D}=15.93$), while the dense array tends to lump together into a thick block. Interestingly, the isolated-pattern NMSE (0.342) is smaller than that for the lower-bound condition; this reflects the fact that NMSE is a global, pixel-energy metric that might insufficiently penalize uniform over-filling with respect to fine-scale edge roughness. In contrast, the ΔCD explicitly considers linewidth expansion and thus offers a better indicator of the engineering penalty for proximity-dominated exposure in this situation.

In summary, the trajectory-coupled study shows that ${\sigma }_{traj}\approx 2~5nm$ gives best consistency with process behavior: edges remain continuous, ΔCD is low, and LER is controlled. The final ADI binary plots for spot size ${\sigma }_{traj}=4.031nm$ are presented in Figure 10. Values far below this range are numerically sharp but physically unrealistic; values well above it amplify proximity effects and erode pattern selectivity.

Table 2. Pattern fidelity metrics (NMSE, ΔCD, and LER) for isolated and dense layouts at different spot size values from the trajectory.

${\mathit{\sigma }}_{\mathit{t}\mathit{r}\mathit{a}\mathit{j}}$ (nm) entry 1	Isolated layout				Dense layout
	NMSE	ΔCD	LER	NMSE		ΔCD	LER
0.156	0.8453	5.3269	2.1952	0.6613		5.0239	11.5968
4.031	1.0242	2.3815	2.1820	1.0524		2.7319	4.8997
39.075	0.3411	18.6698	2.1443	0.6882		4.9998	11.0269

Table 2. Pattern fidelity metrics (NMSE, ΔCD, and LER) for isolated and dense layouts at different spot size values from the trajectory.

${\mathit{\sigma }}_{\mathit{t}\mathit{r}\mathit{a}\mathit{j}}$ (nm) entry 1	Isolated layout				Dense layout
	NMSE	ΔCD	LER	NMSE		ΔCD	LER
0.156	0.8453	5.3269	2.1952	0.6613		5.0239	11.5968
4.031	1.0242	2.3815	2.1820	1.0524		2.7319	4.8997
39.075	0.3411	18.6698	2.1443	0.6882		4.9998	11.0269

3.3. Consistency of PSF Parameters with Published Proximity Data

To determine whether these Gaussian–Lorentzian PSF parameters are consistent with reported proximity behavior for CSAR-62 on Si at 30–40 keV, relevant literature was surveyed (Table 3). Recent Monte Carlo and experimental studies demonstrate that chemically semi-amplified resists such as AR-P 6200 (CSAR-62) contain distinct short-, intermediate-, and long-range scattering components, which are more faithfully reproduced using hybrid Gaussian–Lorentzian or power-law–enhanced Gaussian (PLG) reconstructions than with simplified double-Gaussian models [60,61]. Reported PSF fits indicate short-range Gaussian widths of 2–5 nm—representing forward scattering and beam blur—and long-range backscattering radii on the order of 20–60 nm, even at lower accelerating voltages. Across published Monte Carlo studies, the long-range decay never collapses to only a few nanometers, confirming that tens-of-nanometers tails remain essential at 30–40 keV [54,60,61].

In this regard, for PSF parameters selected for use within this analysis: σ = 4.031 nm and γ = 45 nm, both values clearly lie within experimentally reported values for PSF parameters. This choice of γ simulates correctly at large distances for backscattered dose distributions, while any value for γ below approximately 10nm would cause dose distributions to be underestimated for dose compensation and thus any resultant proximity and CDs to remain inaccurate before OPC compensation.

To further verify physical consistency, a two-dimensional sensitivity sweep was performed by varying the short-range Gaussian width (σ) between 2–5 nm and the long-range backscattering radius (γ) between 30–60 nm. Figure 11 presents ΔCD and LER maps for isolated patterns, and Figure 12 shows the corresponding behavior for dense patterns. Across both layouts, the overall result is strongly σ-driven, reflecting the fact that σ directly controls the width of the forward-scattered energy distribution. Since forward scattering governs the shape of the deposited intensity profile at the resist surface, variations in σ modify the local gradient of the dose field, resulting in pronounced changes in LER and moderate changes in ΔCD for both isolated and dense features.

For ΔCD, the isolated and dense patterns exhibit similar magnitude ranges, but the dense configuration shows noticeably stronger γ-dependence. This arises because dense features accumulate a larger amount of long-range backscattered dose, which increases the effective background exposure and amplifies γ-sensitivity. In contrast, isolated patterns receive substantially less cumulative long-range contribution, so γ primarily shifts ΔCD only modestly. An additional feature in the isolated ΔCD map is the shallow valley located at intermediate σ values and high γ values. This reflects an intermediate regime in which forward blur (increasing with σ) and backscatter redistribution (increasing with γ) partially counterbalance each other, yielding a local minimum in ΔCD. The observed local minimum arises from a balance between σ-driven forward blur and γ-dependent backscatter accumulation within this isolated pattern geometry. Its location varies with layout and process conditions and should therefore be interpreted as a geometry-specific feature of the ΔCD landscape.

For LER, the contrast between isolated and dense patterns becomes more pronounced. LER values for dense features span a significantly wider dynamic range and increase sharply—up to approximately 8.25 in the explored domain—when σ exceeds ~4 nm. This enhanced LER sensitivity in dense layouts is physically expected: small changes in σ broaden the forward-scattered profile, which steepens or smooths local intensity transitions at neighboring feature edges, thereby magnifying edge roughness after development. Dense layouts, having multiple interacting dose gradients at closely spaced edges, exhibit stronger amplification of this σ-induced modulation, while isolated patterns lack such cumulative gradient coupling and therefore display a narrower LER range.

Taken together, the literature comparison and σ–γ sensitivity maps confirm that the adopted PSF parameters fall within a physically meaningful region and reproduce the expected scattering behavior of CSAR-62 on Si at 30–40 keV. The chosen γ preserves realistic long-range redistribution without collapsing the backscatter tail, while σ remains consistent with reported beam-blur and forward-scattering widths. These results collectively support the internal consistency of the trajectory-coupled exposure model used in this study.

3.4. Density-Dependent Shrinkage in RIE: A Loading-like Trend Predicted by the Model

RIE shows a characteristic of density dependence since local reactant transport, by-product removal, as well as the angular distribution of ion flux vary with neighborhood density. In our physics-informed level-set model (Eqs. (24)–(33) in Section 2.4), we quantify the ADI→AEI change by the area-based shrinkage $S\equiv 1-{\displaystyle \frac{A_{AEI}}{A_{ADI}}}$ with a fixed $600\times 600{nm}^2$ window (pixel size 1 nm). Under identical process parameters, the simulation yields $S=4.415\%$ for isolated features and $S=7.026\%$ for dense layouts, i.e., an additional 2.611 percentage points for the dense case.

Table 4. Quantitative summary (window $600nm\times 600nm$, $dx=1nm$).

Pattern type	${\mathit{A}}_{\mathit{A}\mathit{D}\mathit{I}}$ (px)	${\mathit{A}}_{\mathit{A}\mathit{E}\mathit{I}}$ (px)	$\mathbf{Shrinkage}(\mathit{S}\equiv 1-\frac{{\mathit{A}}_{\mathit{A}\mathit{E}\mathit{I}}}{{\mathit{A}}_{\mathit{A}\mathit{D}\mathit{I}}}$)
Isolated	20816	19897	4.415 %
Dense	57798	53737	7.026 %

Table 4. Quantitative summary (window $600nm\times 600nm$, $dx=1nm$).

Pattern type	${\mathit{A}}_{\mathit{A}\mathit{D}\mathit{I}}$ (px)	${\mathit{A}}_{\mathit{A}\mathit{E}\mathit{I}}$ (px)	$\mathbf{Shrinkage}(\mathit{S}\equiv 1-\frac{{\mathit{A}}_{\mathit{A}\mathit{E}\mathit{I}}}{{\mathit{A}}_{\mathit{A}\mathit{D}\mathit{I}}}$)
Isolated	20816	19897	4.415 %
Dense	57798	53737	7.026 %

In agreement with such quantities, the AEI plots (Figure 13) show narrower post-etch linewidths in the dense layout due to stronger bilateral retreat, more significant line-end pullback/rounding (shorter lines with rounded ends), and slightly widened gaps driven by shrinkage from both sidewalls. Such visual signatures coincide with transport-limited behavior as AR grows: the AR-dependent transport penalty suppresses vertical advance, while the lateral/passivation terms dampen sideways propagation even more strongly. Though the case here stems from simulation and therefore requires pitch-sweep SEM calibration for quantitative claims, it indicates a loading-like density trend with practical implications. In dense neighborhoods, larger shrinkage motivates more conservative OPC/recipe guard-bands, whereas isolated features can target tighter CDs without comparable lateral retreat. For this reason, process knobs that relax ARDE penalties—such as pressure/flow tuning, ion-angle collimation (small $\sigma$), or passivation balance—are expected to narrow the density gap in ΔCD.

3.5. Consistency of ARDE Parameters with Published RIE Behavior

The ARDE parameters adopted in this work reflect well-established plasma-etching behavior observed in silicon RIE. As summarized in Table 5, the vertical etch rate is modulated through the transport-access function $f_{acc}\left(AR\right)$, whereas the lateral channel is subjected to the compounded attenuation $f_{lat}\left(AR\right)={\left[f_{acc}\left(AR\right)\right]}^{\beta }exp\left(-\gamma AR\right)$. This structure mirrors physical differences between downward ion-driven removal and sideways chemically dominated erosion. Classical findings by Coburn and Winters [62] showed that strongly directional ions sustain vertical etching even in deep trenches, while lateral advance rapidly becomes transport-limited due to reliance on chemical pathways and off-angle ions. Additional studies by Gottscho et al. [63], Oehrlein [64], and Tserepi and Gogolides [65] established that reactant depletion, by-product stagnation, and polymer-layer growth selectively suppress lateral motion as the aspect ratio increases. More recent measurements further report that the ratio commonly increases to between 2 and 5 times the vertical suppression once AR exceeds 2–3 for Cl₂, HBr/O₂, and fluorocarbon-based RIE processes [66,67].

These experimentally observed trends justify the multiplicative lateral-penalty structure in Eq. (34). The shared transport-limited reduction is captured by $f_{acc}\left(AR\right)$, while the additional power term with β > 1 increases the curvature of the lateral attenuation, producing the approximate relation

$${\displaystyle \frac{f_{lat}\left(AR\right)}{f_{acc}\left(AR\right)}}\approx {AR}^{-(\beta -1)}$$

(42)

which aligns with the reported lateral suppression being between two and five times stronger than bottom etching. Meanwhile, the exponential factor $exp\left(-\gamma AR\right)$ accounts for polymer-driven sidewall passivation and by-product accumulation, phenomena widely observed in HBr/O₂ and fluorocarbon chemistries. The chosen values β = 1.5 and γ = 2.6 therefore arise naturally from plasma-etch physics rather than empirical curve fitting. Taken together, these considerations demonstrate that the ARDE behavior encoded in the present model is consistent with measured silicon RIE characteristics and supports realistic predictions of CD shrinkage and sidewall-angle evolution.

To clarify how the ARDE formulation maps to pattern-scale consequences, Figure 14 plots the normalized CD shrinkage as a function of aspect ratio for four representative ARDE parameter sets:

Case L1 (α = 4.5, p = 2.5, β = 1.0, γ = 0) – lateral suppression absent; CD remains nearly constant with AR, reflecting purely vertical transport-limited behavior.

Case L2 (α = 4.5, p = 2.5, β = 1.5, γ = 2.6) – the baseline parameter set used in this study; CD shrinkage increases monotonically with AR, matching experimentally observed behavior for Si RIE under Cl₂ and HBr/O₂ chemistries.

Case L3 (α = 4.5, p = 2.5, β = 2.0, γ = 4.0) – strong lateral attenuation; CD collapses rapidly, consistent with aggressive fluorocarbon passivation conditions.

Case L4 (α = 8, p = 3, β = 1.2, γ = 1.0) – enhanced vertical ARDE combined with weaker lateral suppression; CD shrinkage becomes smoother and more distributed across AR.

Across these cases, the curvature and rate of CD loss track directly with the magnitude of the lateral ARDE penalty, demonstrating that the multiplicative term $f_{lat}\left(AR\right)={\left[f_{acc}\left(AR\right)\right]}^{\beta }exp\left(-\gamma AR\right)$ governs how rapidly the CD narrows as AR increases. This figure therefore provides direct visual evidence that the selected ARDE formulation reproduces physically meaningful trends and responds predictably to parameter variations, thereby addressing the Reviewer’s concern regarding the physical justification of the model structure.

4. Discussion

Currently available commercial lithography tools like PROLITH, Sentaurus Lithography, and typical OPC/RET simulators can model high-fidelity physics for optical or e-beam exposure but remain at the single-process level without the ability to model a full NIL master-fabrication flow. These platforms lack explicit couplings between spin-coating thickness variation, proximity kernels, development thresholds, and ARDE-driven etch anisotropy. Moreover, NIL master-template fidelity will not be a focused simulation requirement of available software today because the NIL technique itself remains a developing field in the context of the larger OPC/MPC community. A formal model of the calibration needs of the NIL technique has yet to be defined. Our study directly addresses this gap by assembling a physics-linked spin coating → EBL → RIE propagation chain tailored for NIL master fabrication. By tracking how resist-thickness variations affect exposure thresholds, how trajectory-dependent spot size redistributes energy, and how ARDE limits etch advance, the model quantifies stage-to-stage accumulation of variation—a capability unavailable in standard tools.

Within the OPC realm, NIL is considered emerging, and therefore MPC models, calibration recipes, and even standard test protocols are still in early development. This study assembled a first-stage NIL master-template model system that links spin coating, exposure, and etch, and used two canonical layouts—isolated and dense—as boundary conditions to check whether the coupled physics simulate density-driven trends without confounding from geometry. This choice follows long-standing OPC practice: if physics correctly distinguish between the isolated and dense conditions, then the same calibrated parameters can be extended to more complex shapes with confidence.

Spin-coating. The spin-coating model combines Emslie–Meyerhofer scaling with a Bornside-type edge correction so that wafer-scale thickness and edge-bead behavior translate directly into thickness-dependent development thresholds in the exposure step. This connection is significant in NIL because residual thickness and local slopes can subtly bias dose-to-clear. To verify the model at the field scale used for exposure, a 600 nm × 600 nm window was selected after minimizing normalized thickness deviation and local gradient. The obtained sub-map reveals an average thickness of $h=0.196\mu m$ with $\left|\widehat{h}\right|<8.5\times {10}^{-4}$, corresponding to thickness variation < 0.1 nm in the tile; associated gradient magnitude is below ${1\times 10}^{-4}\mu m/mm$, indicating negligible deviation. These checks justify treating the exposure window as effectively uniform in thickness, while keeping the global thickness modulation available to the EBL threshold model. As this study aims to establish a physics-linked workflow rather than a full CDU prediction, thickness-induced CD variation across the wafer is not modeled in this stage. The EBL and RIE analyses therefore operate on an intentionally uniform tile to isolate proximity and ARDE mechanisms. Future extensions will lift this restriction by propagating spatially varying thickness fields into dose-to-clear modulation and etch response.

EBL exposure and development. Instead of just scattering via Monte-Carlo, the EBL module combines column electrostatics and trajectory-derived spot size to a hybrid Gaussian–Lorentzian proximity kernel, with thickness-modulated development thresholds. Parameter sweeping for focusing and energy identifies a low single-nanometer spot-size window ($\approx 2-5nm$) that simultaneously optimizes dimensional accuracy and edge continuity. Approaching extreme under-blur may sharpen pixelwise appearance but degrades edge statistics and violates physically realistic scattering; approaching over-blur raises proximity-driven bias and potentially leads to line bridging in dense arrays. Side-by-side isolated and dense layout comparisons of exposure and proximity-corrected maps visualize this transition and provide a compact diagnostic for recipe tuning.

RIE transfer. On the mask etching simulation, the level-set approach is ARDE-informed. Therefore, the result is naturally a reduction in vertical motion and horizontal motion, in line with the expectation of transportation limitations and shadowing of the ion current. When run with the same process conditions, the result is a shrink with characteristics of the loading effect, stronger in dense design layouts. For illustration, the model yields $S=4.415\%$ for isolated patterns and $S=7.026\%$ for dense patterns, an extra 2.611 percentage points in the dense condition. This captures the CD movement inward, as well as sidewall recession, as if reactant supply and by-product extrusion were impeded.

These findings, taken in combination, tell a coherent story of stage-to-stage propagation. The spin-coating thickness profiles set the conditions for exposure, then the trajectory-correlated spot size governs how proximity redistributes energy, with the ARDE-informed interface motion determining how much of the previous variation preserves into etched topography. From a manufacturing perspective, two implications follow. First, the approach acts as a virtual design of experiment (DOE) for screening dose, focus, and etch combinations, and for quantifying isolated-versus-dense biases with OPC-aligned metrics (ΔCD, LER, NMSE), which shortens mask design iteration cycles. Second, the functional form of the approach facilitates data assimilation: limited thickness tiles, AEI binaries, or cross-sectional data can be used to recalibrate the proximity and ARDE factors so that the predictions track tool and lot drift over time. From an application perspective in NIL OPC, this approach also supports MPC for constrained optimization, in which the sensitivities and constraints are defined by simulators, and the optimized parameters satisfy CD/LER thresholds and sidewall criteria under equipment limitations. These directions are realistic next steps since each module can be characterized in terms of meaningful, process-quantifiable parameters that correspond directly to process-control variables across other coating, exposure, and etch domains.

5. Conclusions and Outlook

This paper establishes the first stage of a three-stage NIL framework that extends over the master template, the replica template, and the wafer. The present contribution is a physics-based, sequential model system for the master template fabrication that serves as a knowledge platform. It specifies governing assumptions, interfaces between modules, and OPC-aligned metric definitions that can be shared by others. Parameter values and correction formulas remain process dependent and tool specific, so the study is intended as a baseline guideline for subsequent calibration rather than a universal recipe.

In this regard, the reported metrics should be interpreted as demonstrating the internal consistency of the physics-linked spin coating → EBL → RIE propagation chain, rather than as calibrated quantitative predictions. The framework is designed as a conceptual and methodological foundation, establishing how cross-stage interactions can be represented and analyzed. Its full quantitative power will be realized when experimentally measured CDs, profiles, and cross-sectional data are incorporated to calibrate proximity kernels, development thresholds, and ARDE parameters. Once such measurements are available, the same structure can be extended toward NIL-specific OPC/MPC workflows and integrated process-control strategies.

With these limitations acknowledged, the next stage of development will focus on extending the framework to replica-template and wafer-level modeling. Future studies will include numerical and physical simulations for imprint mechanics process, resist flow and solidification, demolding behavior, and following transfer. Cross-stage calibration will maintain the same set of interface contracts and metrics, while validation will rely on targeted metrology such as AEI or SEM contours, selected cross-sections, and thickness tiles.

Moving forward, MPC is necessary for NIL OPC. The present framework provides the predictive basis for mask-side correction under equipment limitations and will be used to guide constrained optimization that respects CD and LER targets as well as sidewall criteria. Uncertainty estimates from parameter posteriors can then be used in robust MPC algorithms to ensure coverage under tool drift and lot-to-lot shifts. In this context, two practical directions are identified as promising in the coming future. First, a grey-box model will combine physics principles with future process data to develop surrogates that balance model interpretability with data-driven updates. Second, a learning solution will explore reinforcement learning to optimize dose, focus, and mask-etch settings under constraint handling and safety filters. Both solutions will leverage GPU acceleration for kernel-parameter calibration and broader parameter sweeps, together with fast ARDE-informed level-set evolution in order to ensure that process optimization remains computationally manageable.

Appendix A.1. Spin-Coating Model Setting Parameters and Variables

Table A1. Setting parameters.

Parameter Name	Symbol	Value	Unit	Definition
Wafer radius	$R$	50.0	mm	Radius of the wafer used in thickness simulation.
Reference spin speed	${\omega }_0$	4000	rpm	Spin speed at which reference coefficients are defined.
Central thickness coefficient	$A$	0.20	µm	$\mathrm{Reference}\mathrm{resist}\mathrm{thickness}\mathrm{at}\mathrm{wafer}\mathrm{center}\mathrm{for}\omega ={\omega }_0$.
Central thickness exponent	$n$	0.54	-	Power-law scaling exponent for center thickness vs. spin speed.
Edge-bead height coefficient	$h_{e0}$	0.02	µm	$\mathrm{Edge}-\mathrm{bead}\mathrm{height}\mathrm{at}\omega ={\omega }_0$.
Edge-bead height exponent	$m$	0.60	-	Power-law scaling exponent for edge-bead height vs. spin speed.
Edge-bead width coefficient	$w_{e0}$	2.0	mm	$\mathrm{Edge}-\mathrm{bead}\mathrm{Gaussian}\mathrm{width}\mathrm{at}\omega ={\omega }_0$.
Edge-bead width exponent	$l$	0.25	-	Power-law scaling exponent for edge-bead width vs. spin speed.
Radial correction coefficient	$a_2$	-0.05	-	Second-order radial correction in thickness profile.
Radial correction coefficient	$a_4$	-0.02	-	Fourth-order radial correction in thickness profile.
Wafer grid resolution	-	2000	-	Number of discretization points for wafer thickness map.
Local tile size	-	600	nm	Side length of square window for local resist thickness extraction.
Normalization coefficient	$\alpha$	0.2	-	Sensitivity factor coupling thickness variation into dose threshold.

Table A1. Setting parameters.

Parameter Name	Symbol	Value	Unit	Definition
Wafer radius	$R$	50.0	mm	Radius of the wafer used in thickness simulation.
Reference spin speed	${\omega }_0$	4000	rpm	Spin speed at which reference coefficients are defined.
Central thickness coefficient	$A$	0.20	µm	$\mathrm{Reference}\mathrm{resist}\mathrm{thickness}\mathrm{at}\mathrm{wafer}\mathrm{center}\mathrm{for}\omega ={\omega }_0$.
Central thickness exponent	$n$	0.54	-	Power-law scaling exponent for center thickness vs. spin speed.
Edge-bead height coefficient	$h_{e0}$	0.02	µm	$\mathrm{Edge}-\mathrm{bead}\mathrm{height}\mathrm{at}\omega ={\omega }_0$.
Edge-bead height exponent	$m$	0.60	-	Power-law scaling exponent for edge-bead height vs. spin speed.
Edge-bead width coefficient	$w_{e0}$	2.0	mm	$\mathrm{Edge}-\mathrm{bead}\mathrm{Gaussian}\mathrm{width}\mathrm{at}\omega ={\omega }_0$.
Edge-bead width exponent	$l$	0.25	-	Power-law scaling exponent for edge-bead width vs. spin speed.
Radial correction coefficient	$a_2$	-0.05	-	Second-order radial correction in thickness profile.
Radial correction coefficient	$a_4$	-0.02	-	Fourth-order radial correction in thickness profile.
Wafer grid resolution	-	2000	-	Number of discretization points for wafer thickness map.
Local tile size	-	600	nm	Side length of square window for local resist thickness extraction.
Normalization coefficient	$\alpha$	0.2	-	Sensitivity factor coupling thickness variation into dose threshold.

Table A2. Setting variables.

Variable Name	Symbol	Unit	Definition
Radial coordinate	$r$	mm	Distance from wafer center.
Normalized radius	$\overline{r}={\displaystyle \frac{r}{R}}$	-	Dimensionless radial position.
Thickness profile	$h(r,\omega )$	µm	Total resist thickness as a function of radius and spin speed.
Central thickness	$h_c\left(\omega \right)$	µm	Thickness contribution at wafer center from Emslie–Meyerhofer relation.
Edge-bead height	$h_e\left(\omega \right)$	µm	Edge-bead contribution as function of spin speed.
Edge-bead width	$w_e\left(\omega \right)$	µm	Effective Gaussian width of edge-bead.
Thickness map (global)	$h(x,y)$	µm	Interpolated wafer-scale thickness distribution.
Mean thickness	$\overline{h}$	µm	Average wafer thickness used for normalization.
Normalized thickness variation	$\widehat{h}(x,y)$	-	Local deviation of thickness relative to wafer mean.

Table A2. Setting variables.

Variable Name	Symbol	Unit	Definition
Radial coordinate	$r$	mm	Distance from wafer center.
Normalized radius	$\overline{r}={\displaystyle \frac{r}{R}}$	-	Dimensionless radial position.
Thickness profile	$h(r,\omega )$	µm	Total resist thickness as a function of radius and spin speed.
Central thickness	$h_c\left(\omega \right)$	µm	Thickness contribution at wafer center from Emslie–Meyerhofer relation.
Edge-bead height	$h_e\left(\omega \right)$	µm	Edge-bead contribution as function of spin speed.
Edge-bead width	$w_e\left(\omega \right)$	µm	Effective Gaussian width of edge-bead.
Thickness map (global)	$h(x,y)$	µm	Interpolated wafer-scale thickness distribution.
Mean thickness	$\overline{h}$	µm	Average wafer thickness used for normalization.
Normalized thickness variation	$\widehat{h}(x,y)$	-	Local deviation of thickness relative to wafer mean.

Appendix A.2. EBL Model Setting Parameters and Variables

Table A3. Setting parameters.

Parameter Name	Symbol	Value	Unit	Definition
Beam energy	$E_{beam}$	25	keV	Electron kinetic energy.
Step size	$s$	2	nm	Spatial interval between beam shots.
Beam current	$I$	0.2	nA	Column current.
Dose	$D$	120.0	µC/cm²	Nominal base dose.
Focus voltage	$V_{focus}$	-5	kV	$\mathrm{Lens}\mathrm{setting}\mathrm{optimized}\mathrm{for}{\sigma }_{traj}$.
Anode voltage	$V_{anode}$	+30	kV	Acceleration potential.
Gate voltage	$V_{gate}$	+3.5	kV	Extraction potential controlling beam emission and entry into the accelerating region.
Column length	$L_{column}$	10.0	mm	Distance from the e-beam source to the wafer.
Number of electrons	$N_e$	500	-	Samples for trajectory statistics.
Initial radial spread	$r_{std}$	10	nm	Effective source size (1σ radial spread at entrance plane).
Initial angular spread	${\theta }_{std}$	0.6	mrad	Effective beam divergence (1σ angular spread at entrance plane).
Gaussian weight	$A_1$	0.74	-	Forward scattering weight
Lorentzian weight	$A_2$	0.26	-	Backscattering weight.
Lorentzian scale	$\gamma$	45	nm	Backscattering range parameter
Occupancy threshold	$\rho$	0.25	-	Tile fill criterion for shot placement
Threshold ratio	$\eta$	0.56	-	Dose ratio for development.
Thickness sensitivity	$\alpha$	0.2	-	$\mathrm{Correction}\mathrm{factor}\mathrm{from}\widehat{h}(x,y)$

Table A3. Setting parameters.

Parameter Name	Symbol	Value	Unit	Definition
Beam energy	$E_{beam}$	25	keV	Electron kinetic energy.
Step size	$s$	2	nm	Spatial interval between beam shots.
Beam current	$I$	0.2	nA	Column current.
Dose	$D$	120.0	µC/cm²	Nominal base dose.
Focus voltage	$V_{focus}$	-5	kV	$\mathrm{Lens}\mathrm{setting}\mathrm{optimized}\mathrm{for}{\sigma }_{traj}$.
Anode voltage	$V_{anode}$	+30	kV	Acceleration potential.
Gate voltage	$V_{gate}$	+3.5	kV	Extraction potential controlling beam emission and entry into the accelerating region.
Column length	$L_{column}$	10.0	mm	Distance from the e-beam source to the wafer.
Number of electrons	$N_e$	500	-	Samples for trajectory statistics.
Initial radial spread	$r_{std}$	10	nm	Effective source size (1σ radial spread at entrance plane).
Initial angular spread	${\theta }_{std}$	0.6	mrad	Effective beam divergence (1σ angular spread at entrance plane).
Gaussian weight	$A_1$	0.74	-	Forward scattering weight
Lorentzian weight	$A_2$	0.26	-	Backscattering weight.
Lorentzian scale	$\gamma$	45	nm	Backscattering range parameter
Occupancy threshold	$\rho$	0.25	-	Tile fill criterion for shot placement
Threshold ratio	$\eta$	0.56	-	Dose ratio for development.
Thickness sensitivity	$\alpha$	0.2	-	$\mathrm{Correction}\mathrm{factor}\mathrm{from}\widehat{h}(x,y)$

Table A4. Setting variables.

Variable Name	Symbol	Unit	Definition
Spot size from trajectory	${\sigma }_{traj}$	nm	Beam spot size computed based on electron trajectories at the wafer plane.
Electric field (axial)	$E_z(x,y)$	V/mm	Axial electric field in column.
Electric field (radial)	$E_r(x,y)$	V/mm	Radial electric field in column.
Primary exposure map	$E_0(x,y)$	µC/cm²	Shot-based Gaussian sum.
Electron trajectory	$r\left(t\right),z\left(t\right)$	mm	Radial and axial position vs. time.
Proximity dose map	$D_{map}(x,y)$	µC/cm²	Exposure after convolution
Hybrid kernel	$K_{hybrid}(x,y)$	-	$\mathrm{Weighted}\mathrm{sum}\mathrm{of}{G}_{\sigma }$$\mathrm{and}{L}_2$.
Thickness map (normalized)	$\widehat{h}\left(x,y\right)$	-	From the spin-coating model.
Threshold map	T(x,y)	µC/cm²	Development criterion.
Final resist pattern	$P_{dev}\left(x,y\right)$	-	Binary developed contour.

Table A4. Setting variables.

Variable Name	Symbol	Unit	Definition
Spot size from trajectory	${\sigma }_{traj}$	nm	Beam spot size computed based on electron trajectories at the wafer plane.
Electric field (axial)	$E_z(x,y)$	V/mm	Axial electric field in column.
Electric field (radial)	$E_r(x,y)$	V/mm	Radial electric field in column.
Primary exposure map	$E_0(x,y)$	µC/cm²	Shot-based Gaussian sum.
Electron trajectory	$r\left(t\right),z\left(t\right)$	mm	Radial and axial position vs. time.
Proximity dose map	$D_{map}(x,y)$	µC/cm²	Exposure after convolution
Hybrid kernel	$K_{hybrid}(x,y)$	-	$\mathrm{Weighted}\mathrm{sum}\mathrm{of}{G}_{\sigma }$$\mathrm{and}{L}_2$.
Thickness map (normalized)	$\widehat{h}\left(x,y\right)$	-	From the spin-coating model.
Threshold map	T(x,y)	µC/cm²	Development criterion.
Final resist pattern	$P_{dev}\left(x,y\right)$	-	Binary developed contour.

Appendix A.3. RIE Model Setting Parameters and Variables

Table A5. Setting parameters.

Parameter Name	Symbol	Value	Unit	Definition
Grid size (x, y, z)	$n_x,n_y,n_z$	$600\times 600\times 360$	-	Number of grid points in each spatial direction
Spatial resolution	$∆x$	1.0	nm	Grid spacing in all directions
Time step	$∆t$	0.2	s	Integration step size for level-set evolution
Total simulation time	$t_{total}$	40.0	s	Duration of RIE process
Vertical etch rate (Si)	$V_{vert,Si}$$(R_{base}$)	2.8	nm/s	Vertical etch rate of silicon substrate
Vertical etch rate (Resist)	$V_{vert,R}$	0.28	nm/s	$\mathrm{Vertical}\mathrm{etch}\mathrm{rate}\mathrm{of}\mathrm{CSAR}62\mathrm{resist}(V_{vert,Si}$$\times S_{mask}$)
Mask selectivity	$S_{mask}$	0.10	-	Resist-to-silicon etch rate ratio
Lateral-rate scaling (ion-assisted)	$r_{lat}$	0.03	-	$V_{lat}=r_{lat}{R}_{base}\sqrt{n_x^2+n_y^2}$
Ion angular width	${\theta }_{ion}$	2.5	degree	$\mathrm{Gaussian}\mathrm{angular}\mathrm{envelope};\sigma ={\theta }_{ion}{\displaystyle \frac{\pi }{180}}$
Chemical factor (Resist)	$f_{chem,R}$	0.06	-	Chemical etching contribution for resist
Chemical factor (Si)	$f_{chem,Si}$	0.12	-	Chemical etching contribution for silicon substrate
ARDE strength	$\alpha$	4.5	-	$f_{acc}\left(AR\right)={\displaystyle \frac{1}{1+\alpha {AR}^p}}$
ARDE nonlinearity	$p$	2.5	-
Shadowing factor	$k$	1.6	-	$f_{shad}\left(AR\right)={\displaystyle \frac{1}{1+kAR}}$
Lateral passivation factor	$\gamma$	2.6	-	$f_{lat}\left(AR\right)={\left[f_{acc}\left(AR\right)\right]}^{\beta }exp\left(-\gamma AR\right)$
Lateral attenuation exponent	$\beta$	1.5	-
AR regularizer	$\epsilon$	0.001	nm	$AR\left(x,y\right)={\displaystyle \frac{d(x,y)}{\omega \left(x,y\right)+\epsilon }}$

Table A5. Setting parameters.

Parameter Name	Symbol	Value	Unit	Definition
Grid size (x, y, z)	$n_x,n_y,n_z$	$600\times 600\times 360$	-	Number of grid points in each spatial direction
Spatial resolution	$∆x$	1.0	nm	Grid spacing in all directions
Time step	$∆t$	0.2	s	Integration step size for level-set evolution
Total simulation time	$t_{total}$	40.0	s	Duration of RIE process
Vertical etch rate (Si)	$V_{vert,Si}$$(R_{base}$)	2.8	nm/s	Vertical etch rate of silicon substrate
Vertical etch rate (Resist)	$V_{vert,R}$	0.28	nm/s	$\mathrm{Vertical}\mathrm{etch}\mathrm{rate}\mathrm{of}\mathrm{CSAR}62\mathrm{resist}(V_{vert,Si}$$\times S_{mask}$)
Mask selectivity	$S_{mask}$	0.10	-	Resist-to-silicon etch rate ratio
Lateral-rate scaling (ion-assisted)	$r_{lat}$	0.03	-	$V_{lat}=r_{lat}{R}_{base}\sqrt{n_x^2+n_y^2}$
Ion angular width	${\theta }_{ion}$	2.5	degree	$\mathrm{Gaussian}\mathrm{angular}\mathrm{envelope};\sigma ={\theta }_{ion}{\displaystyle \frac{\pi }{180}}$
Chemical factor (Resist)	$f_{chem,R}$	0.06	-	Chemical etching contribution for resist
Chemical factor (Si)	$f_{chem,Si}$	0.12	-	Chemical etching contribution for silicon substrate
ARDE strength	$\alpha$	4.5	-	$f_{acc}\left(AR\right)={\displaystyle \frac{1}{1+\alpha {AR}^p}}$
ARDE nonlinearity	$p$	2.5	-
Shadowing factor	$k$	1.6	-	$f_{shad}\left(AR\right)={\displaystyle \frac{1}{1+kAR}}$
Lateral passivation factor	$\gamma$	2.6	-	$f_{lat}\left(AR\right)={\left[f_{acc}\left(AR\right)\right]}^{\beta }exp\left(-\gamma AR\right)$
Lateral attenuation exponent	$\beta$	1.5	-
AR regularizer	$\epsilon$	0.001	nm	$AR\left(x,y\right)={\displaystyle \frac{d(x,y)}{\omega \left(x,y\right)+\epsilon }}$

Table A6. Setting variables.

Variable Name	Symbol	Unit	Definition
Level-set function	$\varphi (x,y,z)$	nm	$\mathrm{Signed}\mathrm{distance}\mathrm{function}\mathrm{representing}\mathrm{material}\mathrm{interface}(\varphi =0$$\mathrm{at}\mathrm{boundary};\varphi >0$$\mathrm{in}\mathrm{air};\varphi <0$ in material)
Normal vector (x,y,z)	$n_x,n_y,n_z$	-	Normal direction of surface obtained from level-set gradient
Vertical velocity term	$V_{vert}$	nm/s	Ion-driven anisotropic etch velocity along z-axis
Lateral velocity term	$V_{lat}$	nm/s	Etch velocity component parallel to x–y plane
Chemical velocity term	$V_{chem}$	nm/s	Isotropic etch velocity contribution
Total etch rate	$V$$(R_{etch}$)	nm/s	$\mathrm{Combined}\mathrm{etch}\mathrm{rate}\mathrm{at}\mathrm{each}\mathrm{surface}\mathrm{point}\text{: }V=V_{vert}+V_{lat}+V_{chem}$
Etch depth map	$d(x,y)$	nm	Final etched depth distribution relative to initial substrate surface
Opening width	$\omega \left(x,y\right)$	nm	From Euclidean distance transform of the open region

Table A6. Setting variables.

Variable Name	Symbol	Unit	Definition
Level-set function	$\varphi (x,y,z)$	nm	$\mathrm{Signed}\mathrm{distance}\mathrm{function}\mathrm{representing}\mathrm{material}\mathrm{interface}(\varphi =0$$\mathrm{at}\mathrm{boundary};\varphi >0$$\mathrm{in}\mathrm{air};\varphi <0$ in material)
Normal vector (x,y,z)	$n_x,n_y,n_z$	-	Normal direction of surface obtained from level-set gradient
Vertical velocity term	$V_{vert}$	nm/s	Ion-driven anisotropic etch velocity along z-axis
Lateral velocity term	$V_{lat}$	nm/s	Etch velocity component parallel to x–y plane
Chemical velocity term	$V_{chem}$	nm/s	Isotropic etch velocity contribution
Total etch rate	$V$$(R_{etch}$)	nm/s	$\mathrm{Combined}\mathrm{etch}\mathrm{rate}\mathrm{at}\mathrm{each}\mathrm{surface}\mathrm{point}\text{: }V=V_{vert}+V_{lat}+V_{chem}$
Etch depth map	$d(x,y)$	nm	Final etched depth distribution relative to initial substrate surface
Opening width	$\omega \left(x,y\right)$	nm	From Euclidean distance transform of the open region