Stochastic simulation and calibration of organometallic photoresists for extreme ultraviolet lithography

Abstract. Organometallic photoresists are being pursued as an alternative photoresist material to push the current extreme ultraviolet lithography (EUVL) to the next generation of high-NA EUVL. In order to improve the photoresist performance, an understanding of the photoresist’s response to different process conditions is required. In this endeavor, a stochastic development model is implemented, integrated into full photoresist process steps, and applied for photoresist performance investigations. The model is applied to Inpria-YA photoresist, which works mainly by the process of aggregation. Previously published modeling approaches for metal-organic photoresists assume that the development characteristics of these materials depend only on the size of the created oxo-clusters. In contrast to that, we propose a modeling approach that provides a more detailed description of the interaction among the developer, ligands, and oxo-bonds. Further, the calibration procedures conducted to extract the model parameters to match experimental data are discussed. The model approximated the experimental data with CD RMSE and LWR RMSE of 0.60 and 0.40 nm, respectively. We also investigated the impact of photoresist parameters on the process metrics, line width roughness (LWR), critical dimension (CD), dose-to-size (DtS), and exposure latitude (EL) with the calibrated model. The investigation shows that details of the interaction of photoresist and developer, especially, the so-called development critical value, have a significant impact on the LWR and DtS.


Introduction
Extreme ultraviolet lithography (EUVL) increases the resolution using a smaller wavelength (13.5 nm) for exposure compared to DUV lithography (193 nm). However, this improvement using smaller wavelength also led to several problems, such as stochastic effects. 1 These stochastic effects deteriorate the photoresist performance and limit the economically viable scaling. So far, even if the feature size is decreasing, the line width roughness (LWR) remained constant, which makes the LWR as the limiting factor.
Due to the small number of photons (large shot noise) in EUV lithography and the small absorption efficiency of chemically amplified photoresists (CARs), new photoresist materials are essential. As a result, the industry is pursuing different photoresist materials for the next generation of EUV high-volume manufacturing (HVM). Molecular organometallic resists for EUV are one of the alternative photoresist materials. They have a higher absorption (∼16 to 20 μm −1 ) compared to CARs (∼4 μm −1 ). 2 These photoresists have a core-shell structure where the core contains metal oxide molecules, and the shell contains radiation-sensitive ligands. The absorption efficiency depends mainly on the type of metal used. Inpria-YA photoresists have a tin-oxide core that enables their high photon absorption. [2][3][4] The radiation-sensitive ligands control the reaction behavior of the photoresist, preventing background condensation reaction. 5 Organometallic photoresists are designed to absorb photons that lead to cleavage of the ligands and subsequent chemical reactions to change the solubility. The cleavage of the ligands site creates an "active" site that can create a bond with another "active" site on an adjacent metal-oxide core. As this condensation reaction continues and active sites on neighbor metal-oxide cores form a bond, it creates a large networked structure (oxo-network) that is resistant to the development. 6 Previously, several modeling approaches and investigations of organometallic photoresists have been reported. [6][7][8][9][10][11] The presented models provide a good understanding of the organometallic photoresist behavior during exposure. In this paper, a modeling procedure to characterize and quantify the development process is presented together with a calibration of the model with experimental data to contribute further to understanding of the photoresists' response to different processes. Section 2 explains the modeling procedures applied for exposure, condensation, and development. The calibration of the developed model with experimental Bossung data and verification procedures of the calibration results are discussed in Sec. 3. A simulation study on the impact of photoresist parameters on the lithography metrics, LWR, critical dimension (CD), exposure latitude (EL), and dose-to-size (DtS) is presented in Sec. 4. Finally, the results and findings are summarized, and an outlook for future work is discussed.
2 Modeling Procedure

Exposure
The photoresist simulation volume is discretized into lattice cells of size δx, δy, and δz, where δx Ã δy Ã δz is the volume of the photoresist molecule (metal-oxide core with ligands). The x and y represent directions perpendicular to the exposure direction, and −z represents the direction of exposure. The intensity distribution inside the photoresist (bulk image) is simulated using the Fraunhofer lithography simulator, Dr. LiTHO. 12 Then the intensity absorbed by the photoresist molecule is computed from the bulk image according to the Beer-Lambert law, as described in the following equation: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 3 4 5 I abs;i ¼ I i Ã ð1 − e −αÃδz Þ; (1) where I i and I abs;i are the bulk image intensity and the absorbed intensity for the i th lattice cell, respectively. α is the absorption coefficient and δz represents the thickness of a single lattice cell. From the computed absorbed intensity distribution, the average number of photons absorbed ðhN p;i iÞ in the i th single lattice cell is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 2 6 3 hN p;i i ¼ I abs;i Ã D Ã δx Ã δy hc∕λ ; (2) where D is the exposure dose and ðδx Ã δyÞ is the surface area of a single lattice cell. h represents the Plank's constant (6.62607004 × 10 −34 Js), c is the speed of light (3 × 10 8 m∕s), and λ is the wavelength of the exposure (13.5 nm for EUV). The stochastic distribution of the photons in the photoresist volume follows a Poisson distribution. 13,14 Therefore, the actual number of photons absorbed by the photoresist molecules is distributed according to a Poisson distribution from the average number of absorbed photons.
The absorbed photons trigger subsequent chemical changes to modify the solubility of the photoresist. For EUVL, the chemical change is mainly initiated by the photoelectric effect where photon absorption leads to the generation of electrons. 7,15 The generated photoelectron, in turn, produces secondary electrons for further ionization of the photoresist molecules at a distance from the point of photon absorption-resulting in a blurring effect on the distribution of the absorbed photons. 16 As tracking each electron and its interaction with the photoresist molecules is computationally intensive, a simplified approach is implemented in our model to simulate these processes.
The average number of generated electrons per absorbed photon is approximated by the quantum yield. E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 6 8 7 where N e and N p represent the number of generated electrons and absorbed photons, respectively. Φ is the quantum yield, the actual number of electrons generated per absorbed photon. The actual number of electrons generated per photon is computed by generating the actual quantum yield according to the Poisson probability distribution from the average quantum yield hΦi and Eq. (3). After the computation of the generated electrons, the electron blur due to the electrons' random walk is approximated by Gaussian convolution. Such approach enables a significant reduction of computing time. The Gaussian kernel applied for the convolution is defined as in the following equation: where σ blur is the electron blur length in nm. Application of this convolution results in a (deterministic) blur of the location of the generated electrons with respect to the point of photon absorption. Typical photon and electron distributions that are computed from the bulk image (generated by Dr. LiTHO) are shown in Fig. 1.

Condensation Reaction
The photoresist changes its solubility after exposure due to physical or chemical structure changes. Several studies demonstrated that the metal-oxide cores absorb the photons to generate photoelectrons. However, the ligands act only as non-reactive spacers inhibiting the metal-oxide cores from reacting. 17,18 During exposure, the protecting ligands are cleaved by electrons. This cleavage of the ligands leads to a creation of what is called an "active site," i.e., a site where a condensation reaction [i.e., the creation of a bond (M-O-M oxo-bridge) with an active site of another adjacent metal-oxide core] can occur. 6,17 The number of active sites generated on each particular metal-oxide core follows directly from the number of electrons that land in the δx, δy, δz volume element of the metal-oxide cluster. Whether or not an active site will lead to the formation of an M-O-M bridge is simulated in a probabilistic approach, using percolation theory. 6,9,19 In the model used in Ref. 6, a metal-oxide core and a neighbor metal-oxide core are randomly selected, and if both cores have active sites, a bond is created. If there are no available active sites, another core and its neighbor are randomly selected. However, this implementation is computationally intensive, and it is not well suited for the large number of computations required for a calibration procedure and for extensive modeling studies. Therefore, our model uses a different, pseudo-random, approach that is faster without losing the stochastic nature of the process. The operation scheme is summarized in the pseudo-code in Algorithm 1. The metal-oxide cores with active sites are visited sequentially. Then one (and just one) of the six nearest neighbor metal-oxide cores is randomly selected. If this particular neighbor core also has an active site, an oxo-bond between the two cores is created, replacing the two active sites (i.e., the bond consumes the two active sites). Otherwise, the active site remains active. Then the next active site on the given metal-oxide core is selected, and another neighbor metal-oxide core is randomly selected. This process is iterated until all metal-oxide cores, and all of their active sites are visited. The pseudo-code of the implementation is presented in Algorithm 1.
A periodic boundary condition is applied for xand ydirections. If the randomly selected neighbor's index gets out of the simulation domain, a bond is created with the metal-oxide core on the opposite side.
Post exposure bake (PEB) impacts the condensation reaction and, in turn the sensitivity of the photoresist. 20 However, our investigation focuses on the modeling of the stochastic development process for a finished condensation reaction. Therefore, we do not consider explicit modeling of the PEB process.
The oxo-networks computed from the electron distribution, for a contact feature, with exposure doses of 20, 40, and 90 mJ∕cm 2 are shown in Fig. 2. For a small exposure dose, the created oxo-networks are small in size. As the dose increases, the size of the oxo-networks that are created in the exposed region increase, whereas small oxo-clusters are created in the unexposed regions. For large exposure dose, a large oxo-cluster is created in the exposed region [Figs. 2(c) and 2(d)].

Development
After the solubility is changed for the exposed regions of the photoresist due to the condensation reaction, the photoresist is developed with a solvent. This development step creates the final photoresist profile.
The development behavior of organometallic photoresist that aggregates during exposure, to create a large oxo-network, is mainly dependent on the size of the created oxo-cluster. Small oxo-clusters can develop quickly, whereas large oxo-clusters are resistant to the developer. There are several alternatives to implement the development process, such as simple threshold, 8 percolation model, 22 or tracking of the oxo-network in contact with the photoresist. 6 The investigation of stochastic effects during the development requires a more detailed description of the  [23][24][25][26][27] demonstrated that the photoresist dissolution behavior depends on the behavior of the ligands on the metal cores. Our stochastic development model tracks the substitution of the ligand and oxo-bond sites by the solvent molecules based on balance equations. If the number of these solvent substituted sites per metal-oxide core is above a certain number and all the oxo-bonds on the core are broken, the core is washed away. This procedure is iterated until the developer time is reached. For this purpose, we employed the basic approach of the critical ionization model to track the interaction of the ligands with solvent molecules in stochastic development model, based on the implementations described in Refs. 28 and 29. The reasoning and details of the model formulations are explained below.
At the beginning of the development process, the metal-oxide cores in the photoresist have three kinds of sites, as shown in Fig. 3. 1. Ligand site. The ligand is not cleaved, i.e., it is not affected by the exposure. 2. Active site. An active site is created, but it did not form an oxo-bond during condensation. 3. Oxo-bond site. An active site created during exposure is replaced by a bond during condensation.
The illustrations in Fig. 3 place the components on a specific geometrical location in the grid cell. In the simulation, however, no such geometrical-location information is used. In each cell, we only keep track of the number of L, A, B, and S components. During the development in organic solvents, the ligands and active sites can be dissociated from the metal-oxide core and substituted by a solvent molecule. 30,31 Ligands from the developer bulk can also substitute a solvent site on a metal-oxide core again, in a reverse reaction. Oxobond sites are normally assumed to be resistant to the developer. With strong developer solutions, however, the bonds can be broken and substituted by solvent molecules. A metal-oxide core is removed from the photoresist into the developer solution, if the number of solvent substituted sites per metal-oxide core is larger than the development critical number of sites (L th ).
In order to reduce the complexity and the number of parameters in our model, ligand and active sites are assumed to be similar and will be treated identically in the development step. In this case, the active sites will be treated as ligand sites. This limits the types of sites on the metaloxide cores, at the start of the development process, to only two kinds-ligand and oxo-bond sites. For example, MOL 2 AB is treated as MOL 3 B. Using this assumption, the population balance equations for substitution reactions of the development processes at the interface between the developer and the photoresist can be defined by the following equation: 30 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 a ; 1 1 6 ; 3 4 3 (5a) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 b ; 1 1 6 ; 2 8 6 where S is the solvent and, K 1 , K 2 , and K 3 are the rate constants for the ligand substitution, the solvent substitution, and the oxo-bond substitution, respectively. M, O, L, and B represent metal, oxygen, ligand, and oxo-bond on the photoresist molecule and n and m represent the number of ligands on the given metal-oxide cores. The mechanisms of the population balance equations, Eqs. (5a) and (5b) are illustrated in Figs. 3(e) and 3(f), respectively. Note that, as we consider MO-bonds with one of the six direct neighbors of each cell, and our model only counts the number of L, A, B, and S in each cells, we need to distinguish no more than six ligands per core, even if the real resist molecules may well have more than six ligands per core. Our model describes the process of stochastic events during the development by an implementation of the Gillespie algorithm. This algorithm employs Eq. (5) and given reaction rates to generate stochastically correct, possible trajectories of the developer front. This means that the reactions are simulated explicitly as a stochastic process instead of solving the equation analytically.
In the model's implementation, the solvent molecules' reaction with ligand and oxo-bond sites of the cores are simulated by discrete events with probabilities that are governed by the reaction rates. Only the cores that are in contact with the developer are considered. To track the reactions by separate events, a time step, with a probability that a single reaction occurs in this time interval, must be first defined. This time step is computed based on the reaction rates (K 1 , K 2 , and K 3 ) and the corresponding number of sites on the cores. Then the time is updated. In each time step, substitution reactions are determined with a probability for all the cores on the development front. The reaction can be either ligand or oxo-bond or solvent substitution reaction. The probability of the selected substitution reaction depends on the reaction constants and governs the update of the sites on the cores. If the number of solvent substituted sites on a core is larger than the development critical value (L th ), and all oxo-bonds are broken, the core is removed, and its neighbor cores are added to the developer front. Then the next time step is computed, and the process is repeated. These processes are iterated until the total development time is reached. The developer front trajectory, tracked by the steps described above, creates the final photoresist profile.
Details on this implementation are given in Appendix A. The implemented development model transforms an oxo-bonds distribution, as obtained from the condensation reaction, into the final latent profile. The distribution of the oxo-bonds per core is shown in Fig. 4(a) and the developer front at different development time steps are shown in Figs. 4 For a standard implementation of the model in Python, simulation of a vertical line with 36 nm pitch and 100 nm line length, carried out on only a single-core on a PC with Intel-Core i7-4770 at 3.4 GHz (CentOS Linux 7) and 4 GB RAM, requires ∼9.5 h. The optimized implementation of the algorithm described in Appendix A reduces the computation time to only ∼100 s.

Illustrations of the Complete Process Flow
A schematics that illustrates the working principle of the model is presented in Fig. 5. The schematics shows the photoresist molecules and the effect of the different processes on the sites of the metal-oxide cores; the activation of the ligands during exposure, random creation of oxo-bonds during condensation, and the evolution of the developer front during development.
Representative simulation results for the modeling process flow of the organometallic photoresist are shown in Fig. 6. The process starts with the computation of the bulk image [ Fig. 6(a)] based on the exposure conditions, such as the source and mask specifications. The average number of photons absorbed in the photoresist is computed from the bulk image by specifying the exposure dose and and the absorption coefficient (α) based on Eq. (2). The average number of photons absorbed in the photoresist is redistributed according to the Poisson probability distribution. This redistribution process generates a stochastic distribution of the photons [ Fig. 6(b)].
The electrons generated from the absorbed photons are computed from the stochastic distribution of the photons and the quantum yield (Φ), based on Eq. (3). To describe the electron blur effect, the distribution of the electrons is convoluted with the Gaussian distribution with a standard deviation equal to the blur length (σ blur ), as described in Eq. (4). The final electron distribution is shown in Fig. 6(c). Afterward, the electrons' distribution is related to the activation of ligands on the metal-oxide cores, and oxo-bonds are created between activated sites of neighboring cores by applying the percolation model (see Algorithm 1). This procedure results in the distribution of oxo-bonds, displayed in Fig. 6(d). Finally, the development step is simulated based on the distribution of oxo-bonds, applying the model parameters, the development critical value (L th ), the rate constants (K 1 , K 2 , and K 3 ), and the development time (see Algorithm 2 in Appendix A). The final result of the development simulation, the developer arrival time (DArT), is shown in Fig. 6(e). Table 1 summarizes the model parameters and sources of variability in each process step, together with the implementation procedures that are applied to capture these variabilities. The variability that is expressed by the LWR is a combined result of stochastics in the photon absorption, the electron generation, the bond creation, and the interaction of the developer molecules with the ligands and bonds of the photoresist molecules.

Calibration of Model with Experimental Data
The parameter values for the developed photoresist model, summarized in Table 1, are extracted by a model calibration with experimental data. The calibration procedure is conducted using lines-and-space (L/S) Bossung data (CD and LWR values) with exposure dose values ranging The bulk image; the stochastic distributions of (b) absorbed photons; (c) generated electrons; (d) oxo-bonds per core; and (e) the DArT. The experimental data contain 758 measurement points for vertical L/S patterns with 5 mask CD and pitch variations, as shown in table in Fig. 7(a) for a center-slit position. For our model calibration, only 73 data points are selected to keep the computation time to a minimum but also include enough data points that can define the model behavior, simultaneously. The calibration data points are shown in red X marks in Fig. 10. All measurement data points are used for verification of the calibrated parameter values.
The experimental data are from CD-SEM measurements conducted for 16 scans (or frames). During multiple scans of the CD-SEM measurement, the electron beam shrinks the photoresist lines. Our calibration is based on CD-SEM data from the single-frame and 16th-frame. These data exhibit a CD shrinkage in the range of 2.5 to 6 nm between single frame and 16th-frame. Single-frame SEM images have a non-negligible but small shrinkage compared to the 16thframe image. In our calibration and verification, we used the CD data from the single-frame measurement. However, the single-frame measurements are noisy and cannot provide a reliable LWR measurement. Therefore, we used the LWR data from the 16th-frame SEM image for the calibration. The effect of the SEM on the roughness due to the photoresist lines shrinkage is (a) (b) Fig. 7 (a) Feature types and (b) source used for the measurement. The feature types are specified as P "pitch" V "mask CD," with the specified mask CD and pitch. "V" stands for vertical line.

Simulation Parameters and Process Conditions
The process conditions and simulation parameters in Table 2  Pythmea, 34 a python multi-objective evolutionary algorithm from Dr. LiTHO, is used for the calibration of the model to the experimental data. The model parameters are calibrated for CD and the LWR simultaneously using a multi-objective optimizer. The fitness of the calibration was determined, for both LWR and CD, from a root-mean-square error (RMSE) of simulation results from the experimental data. After the model is calibrated, the verification runs were conducted with 300 nm line length for 758 data points for a calibration result selected from the Paretooptimal solutions.
The photoresist model parameters, absorption coefficient, quantum yield, electron blur length, and development model parameters, were varied during the calibration (Table 4). Because the simulation zero-focus position deviates from the exposure tool zero-focus position with an unknown offset, a focus offset parameter is included in the calibration. Additionally, metrology offset parameters, to compensate for systematic deviations of CD and LWR measurements, are also included in the calibration.

Model Options, Calibration Results, and Verification of the Models
Calibration was conducted under three different assumptions.
1. Model 1. In this model, it was assumed that the developer does not affect the oxo-bonds created during the condensation reaction. For this procedure, only Eq. (5a) is considered and the metal-oxide cores are developed if the number of solvent substitute sites is above the critical value, irrespective of the number of oxo-bonds on the core. This assumption is valid for a solvent with a small dielectric constant or low polarity. 35 The CD RMSE and LWR RMSE values for the calibration and the verification, for one of the selected solution from the Pareto front, are summarized in Table 3. 2. Model 2. This model includes the substitution reaction for the oxo-bonds. Yeh et al. 23 demonstrated that developers with high dielectric constants, such as ethanol, break oxo-bonds, and damage the oxo-networks created during the condensation reaction. AFM images of the patterned lines, patterned using ethanol as a developer, exhibit a rough surface with apparent damage to the oxo-cluster. In contrast, patterned lines developed in 1-propanol have a smooth surface. Even though the patterned lines developed in 1-propanol have a smooth profile, it does not necessarily mean the breaking of the oxo-bonds does not occur. Therefore, in this model, both balance equations [Eqs. (5a) and (5b)] are considered. Moreover, it is enforced that every oxo-bond on a photoresist core should be broken before the core is considered as developed. This condition has to be satisfied even if the number of solvent substituted sites on the core is above the development critical value. The addition of this assumption increased the sensitivity of the development model with respect to number of oxo-bonds per core. As summarized in Table 3, model 2 fits the experimental data better than model 1. Figure 8 compares the RMSE data that were obtained in the verification of the models with and without breaking of the oxo-bonds. 3. Model 2 with mask bias. An additional calibration parameter, mask bias, is applied. The verification RMSE for CD exhibited an unreasonably high value for P60V27 compared to the other feature types for both models. The occurrence of a large RMSE of CD values at a single pitch is hard to explain by chemical effects in the photoresist. As a result, an additional calibration parameter, a mask bias only for P60V27 feature, is applied and a new calibration run was conducted. The results show a further improvement of the RMSE  value for the CD and LWR (Fig. 8). This calibration run results in good accuracy for the CD and LWR (results shown in Table 3).
The verification results of the three calibration runs are summarized in Table 3. The model approximates the experimental data with CD RMSE of 0.60 nm and LWR RMSE of 0.40 nm for the verification run. These results are comparable to the observed results of well-established CAR photoresist models.
The relative errors of all data points, shown in Fig. 9, were calculated from the verification errors normalized by the corresponding experimental CD or LWR values. The model can approximate the experimental data with 30% maximum error for LWR and 6% maximum error for the CD. Uncertainties in CD computation for profiles with higher roughness can lead to minor deviations on the final CD and LWR values. The parameter values for the "best" calibration solution selected from the Pareto front are summarized in Table 4.
The Bossung data for verification of the calibration parameters for the smallest pitch (P36V18) and the largest pitch (P70V27) are shown in Fig. 10. For P70V27, the simulated LWR and CD errors vary from 0 to 1.3 nm and 0 to 1.0 nm, respectively. Similarly, for P36V18, simulated LWR errors vary from 0 to 2.0 nm and simulated CD errors vary from 0 to 1.14 nm. The experimental data for feature type P36V18 contain outlier data that could not be approximated by the model [marked by blue circles in Fig. 10(b)]. Excluding the outlier data,

Relation Between Photoresist Parameters and the Lithography Metrics
This section studies the impact of the parameters of the implemented model on important lithography metrics CD, LWR, dose-to-size, and exposure latitude. To perform a qualitative analysis of observed dependencies and to derive general tendencies, we extend previously published work 36 on the application of correlation analysis (CA) for photoresists parameters. The generation of appropriate datasets for CA requires a preprocessing to provide an efficient sampling of the data space and to avoid defects in the simulated photoresist profiles. To focus our investigations on the impact of the photoresist parameters and to include knowledge from previously published work on the scaling of lithography metrics, 37,38 some of the computed lithography metrics are normalized by factors that describe the impact of imaging parameters.

Methodology
First, the parameter values of the datasets have to be randomly generated. For this purpose, parameter values are selected from the given bounds with Latin hypercube sampling (LHS) instead of random sampling to cover the entire parameter space with a minimum number of datasets. 39 Each dataset contains values of process settings (mask CD, focus, and pitch) and photoresist parameters (α, σ blur , Φ, L th , K 1 , K 2 , K 3 , and th) summarized in Table 5. 600 samples are generated for the 11 parameters to ensure a good convergence for the CA. The pitch, unlike the other parameters, is not directly generated. Instead, the duty ratio is sampled; and the pitch is computed from the mask CD and the duty ratio. The specifications of the lower and upper bounds used for LHS, summarized in Table 5, are determined by the effect of the parameter values on the final results of the simulation, to avoid defects on simulated profiles. In addition, the bounds are also limited to keep the monotonic relation of the parameter with the metrics because a non-monotonic data can suppress the correct results of the correlation. The relation of the LWR with the blur length is monotonic for a small range of blur length values. Careful selection of the parameter range and the application of the Spearman's rank correlation (see below) helps to address the characteristics behavior of blur in our CA.
After the generation of parameters, simulations are performed for the computation of the metrics. First, the CD values, corresponding to the randomly generated process settings (mask CD, pitch, and focus), are computed with the calibrated parameters of model 2 from Table 4 with the nominal dose of 52 mJ∕cm 2 . In the following simulation with variable photoresist parameters, the CD values, obtained from the calibrated model 2 parameters, are used as target-CD for the computation of DtS and EL. LWR was extracted at the computed DtS. Finally, the CD was computed with the nominal dose of 52 mJ∕cm 2 .
The final step in the preprocessing of the data before CA addresses the simultaneous variation of process settings (mask CD, pitch, and focus) and photoresist parameters and their impact on the obtained results. In theory, it would be possible to vary only the photoresist parameters and apply the CA. However, the study of the photoresist response for fixed exposure conditions will not provide representative data for a wide range of applications. In order to include a wide range of process variations, the process settings are varied as well, but their dominating impact on certain lithography metrics has to be considered. This is done by appropriate normalization techniques and application of scaling rules for lithography metrics that have been described by several authors. 37,38,40 For example, it is known that the LWR decreases with increasing normalized image log-slope (NILS) and dose. Therefore, we normalize the variation of exposure conditions by application of appropriate scaling rules.
These impacts of the process settings are normalized in two steps.
1. The impact of the investigated process settings on the image quality and lithography metrics is typically quantified by image log-slope (ILS) and NILS. Here we employ analytical dependencies from Refs. 37, 38, and 40, to account for the impact of image quality and exposure dose on LWR and EL. Specifically, we normalize LWR according to Ligand substitution rate (s −1 ) In the case of CD, we use the change in CD (ΔCD) instead of the absolute value, deviation of the simulated CD from the target-CD. For DtS, we use the data without any normalization. 2. The impact of some photoresist parameters on the metrics is dependent on the process settings, and these indirect contributions of the process settings should be normalized. For example, the effect of blur length on the metrics depends on the pitch. In order to remove the pitch dependency and treat only the blur length impact, the modulation transfer function (MTF), defined in Refs 38 and 40, is fed to the CA instead of the blur length values. MTF, defined as MTF ¼ e −2ð πσ blur P Þ 2 , is the Fourier transform of the Gaussian kernel (Eq. 4) that we applied for the approximation of the electron blur effect. 41 However, the inverse relation 1∕MTF is applied because the increase in blur length values corresponds to a decrease in MTF values. Otherwise, the sign of the correlation coefficients for the blur length will be inverted.  Figure 11 demonstrates the impact of these postprocessing procedures. It can be seen that the application of the normalization increases the correlation between the photoresist parameters and the lithography metrics.
Finally, the combined impact of the parameters on the lithography metrics has to be considered in the CA of the lithography metrics. To decorrelate the combined impacts of the parameters on the lithography metrics, we compute the semipartial correlation coefficients. 36 To consider the monotonic, but non-linear relation between some of the parameters and the lithography metrics, Spearman's rank correlation coefficient is computed instead of Pearson's linear correlation coefficients. To address both non-linear dependencies and decorrelate the impact of the photoresist parameters, we calculate semipartial rank correlation coefficients (SRCCs). 39,42 SRCC values are in the range from −1 to 1. Negative SRCC values mean the increase in the photoresist parameter value leads to a reduction of the lithography metrics, whereas positive SRCC values mean the increase in the photoresist parameter value leads to an increase in the metrics. On the other hand, an absolute value of SRCC above 0.2 indicates that the photoresist parameter's variation significantly impacts the metric for the defined parameter space. Otherwise, the impact is insignificant.

Results of the Analysis
The results of our CA, shown in Fig. 12, provide several expected findings that are consistent with the literature. As shown in Fig. 12(a), the quantum yield (Φ) has a significant effect on the reduction of the LWR. The increase in the blur length (σ blur ) and the photoresist thickness (th) leads to an increase in LWR. The CD of lines increases with increasing absorption coefficient (α) and quantum yield [ Fig. 12(b)]. In turn, the increase in the CD for the increase in these parameters means a decrease in the DtS. EL is mainly affected by the blur length; the increase in the blur length leads to a reduction of the EL. The findings on the development parameters are less obvious. The development critical value (L th ) exhibits a significant correlation with LWR and DtS and increases the CD of the lines. The correlation of the other development parameters with the process metrics is insignificant. The development critical value has a linear effect on the lithography metrics-its value decides the number of oxo-bonds per core that create the edges of the profiles. Its variation leads to a variation in the shift of the edge position. The other development process parameters impact the path of the development, 14 but their impact on the CD and LWR is less seen in the final results of the lithography metrics variation. Nevertheless, the inclusion of these parameters and the corresponding path is important to obtain a fitness of the model calibration. The full understanding of the impact of these parameters on the process metrics needs further investigation.
Finally, we have to emphasize that the CA results are dependent on the parameter bounds chosen for the investigations. Consideration of different ranges of values can lead to deviations from the presented results. 36 Nevertheless, the observed tendencies especially for the significant impact of the development critical value on LWR are observed for all reasonable choices of parameters ranges.

Conclusion and Summary
A development model for organometallic photoresists that tracks the developer's effect on the ligand and oxo-bond sites of the metal-oxide cores for negative-tone development was proposed and implemented in a complete imaging and resist process simulation flow. The model approximates the experimental Bossung data, of lines and spaces with different features and pitches, with CD RMSE of 0.60 nm and LWR RMSE of 0.4 nm. Notably, the application of multi-objective optimization and the detailed description of the dissolution process by the model resulted in good approximations of both CD and LWR simultaneously. The model also demonstrated that the development behavior of the investigated organometallic photoresist can be approximated based on the number of created oxo-bonds on the cores. The calibration results show that, for Inpria-YA photoresists, three ligands on the photoresist core should be substituted by solvent molecules before the photoresist core is assumed to be developed.
The correlation analysis confirms that both exposure and development photoresist parameters have a significant impact on the lithography metric. The most important impact of the development parameters is seen for the development critical value that impacts the LWR, CD, and DtS. The increase in the development critical value leads to a reduction in LWR and DtS, and an increase in the CD of the lines. The increase in quantum yield has a significant influence on the reduction of the LWR and the DtS. However, the quantum yield has no significant effect on the EL. The increase in the quantum yield also leads to an increase in the CD.
Application of the model to other feature types including contact holes and line ends could provide further insights on its predictivity. In addition, calibration of the models including the LER in addition to the LWR and CD data can improve the model and give better understanding of the development process. Further investigation and extensions of the photoresist parameter CA could help to separate and optimize the impact of exposure and development effects on the trade-off among resolution, LWR, and sensitivity.

Appendix A: Implementation of the Development Model
A discrete stochastic development model is implemented by simulating each reaction explicitly using the Gillespie algorithm. 28,43 The evolution of the development process through time can be tracked by computing the reaction of the solvent molecules with the sites of the metal-oxide cores that are in contact with the developer, based on Eq. (5).
The probability that a reaction occurs in time interval ½t; t þ δt is given by R i Ã δt, where i is a reaction type and R i is a reaction rate. At the start of the simulation, the developer front contains only the cores at the top of the photoresist as the cores in contact with the developer constitute the developer front. 28 The reaction rates of the population balance equations in Eq. (5) are defined based on the number of ligand, oxo-bond, and solvent sites on the metal-oxide cores.
These reactions rates determine not only the "stochastic state" of the developer front but also the next reaction. For our case, the state of the developer front is computed by the following equation: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 6 ; 6 9 9 where N L , N B , and N S are the total number of ligand sites, oxo-bond sites, and solvent sites on the developer front, respectively. These total number of sites are the sums of the respective number of sites on metal-oxide cores that are on the developer front. The cell for the next reaction and the reaction type are randomly selected based on reactivity rates, R L , R B , and R S . First, a core i is randomly selected from the cores on the developer front. Then the procedure based on Von Neumann rejection from Ref. 28, instead of the original Gillespie algorithm procedure, is applied for acceptance of the randomly selected core. The Von Neumann rejection is adapted for our application as shown in the following equation: where r 0 is generated from a uniformly distributed random numbers. For the acceptance of the randomly selected core, the sum of the reaction rates of the core, excluding the oxo-bond substitution reaction rate, is compared with the maximum of the ligand or the solvent substitution rate. This is because a metal-oxide core without an oxo-bond is expected to be developed-probability of reaction is 1-while a core that created oxo-bonds on all of the possible sites is resistant to the developer-probability of reaction is 0. As a result, the probability that the core reacts is computed from the solvent and ligand substitution rates compared to the maximum value of the ligand substitution rate.
If the selected core is accepted, the reaction type is randomly selected and the number of sites on the core (i) is updated, as described in Eq. (10): where i represents the selected core on the developer front. r 1 is generated from the uniformly distributed random numbers and N j;i is the number of sites of type j (ligand or oxo-bond or solvent sites) on core i selected from the developer front. The time step for the reaction, in the Gillespie algorithm, is calculated from the state of the system using Eq. (11). Then the time will advance from t to t þ δt: 28 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 1 1 6 ; 1 5 4 δt ¼ − lnðr 2 Þ∕R tot ; where r 2 is generated from a uniformly distributed random numbers. If the number of solvent substituted sites per core is larger than the critical number of sites (L th ), and all the oxo-bonds on the core are broken, the core is developed. If the core meets the development criteria, the developer front is updated by removing the core and inserting the neighbor cores to the developer front.
These steps-random selection of the core and the reaction type, update of the sites on the core, and update of the developer front-are executed for each time step. These processes are iterated until the development time is reached.
Similar to the previous process steps of exposure and condensation, a periodic boundary condition is applied for x and y directions.

Computational Optimization Approach
The exact stochastic simulation approach for the Gillespie algorithm (the naïve implementation), discussed above, is computationally intensive. 44 Randomly selecting a single core and the reaction type for the selected core requires several iterations until all the cores on the developer front are visited. Also due to the large R tot , because the developer front contains several cores, the update time (δt) is small. Until the development time is reached, the simulation requires several iterations. Furthermore, after each iteration, the reaction rate is recalculated, and the developer front is updated. The developer front update is also a computationally intensive process. As a result, especially for simulation of profiles with long lines, as it is required for proper computation of the LWR, the development model becomes inefficient. Therefore, to reduce the computation time, optimization of the processes is required.
In order to reduce the number of updates, the tau-leaping approximation procedure 44,45 is included the implementation. In the tau-leap method, the number of sites updated at a time is randomly selected based on Poisson probability, instead of executing a single reaction for each iteration. The generated random number of reactions is used to update the sites on the core based on the selected reaction type. Due to the Poisson distribution's unbounded nature, the maximum numbers of reactions are limited to the number of available sites on the core corresponding to the reaction type, to avoid negative values for the sites on the core. 46 Note that, in our implementation, a single reaction is executed if the Poisson generated number of reactions is zero.
However, the tau-leaping method-based optimization is still insufficient for the model to be fast enough to be applied for calibration. For a single metal-oxide core, six reaction sites are available with three possible reactions types. Each core requires several reaction steps-at a minimum, reactions equal to the critical value-before it is developed. In addition, due to the large size of the developer front, several iterations are still required to change the state of the developer front. As a result, an aggressive optimization approach is required.
The development process is a weakly coupled reaction network. 47,48 This means that the development behavior of a core on the developer front is affected only by the immediate neighbor cores, irrespective of the developer front's size. In addition, a reaction on a metal-oxide core affects the neighbor cores or changes the state of the developer front significantly if the reaction leads to the development of the core. Due to these reasons, the metal-oxide cores on the developer front can be treated to be spatially independent. In other words, the developer front can be subdivided into smaller systems where the reactions types are probabilistically selected for each subsystem. 49 As a result, the developer front is partitioned into subsystems of single cores. In this assumption, each core can be treated as a separate system, and the reaction type is sampled independently for each core. Therefore, all the cores can simultaneously undergo substitution by a solvent or by a ligand.
In the current implementation, the cores on the developer front are visited sequentially, and depending on the reaction rates of the corresponding core, the reaction type is randomly selected. The next reaction's time step is computed separately for each core from its total reaction rate based on Eq. (11). The computed update time is multiplied by the Poisson generated number of reactions, to account for the number of reactions generated based on Poisson probability. Then based on the selected reaction type, the sites on the current core are updated. After all the cores are visited, the reaction rates are recalculated, and the developer front is updated. These processes are iterated until the development time is reached. The pseudo-code for the optimized procedure is presented in Algorithm 2.