Methods and systems for parameter-sensitive and orthogonal gauge design for lithography calibration

[0001] The present application claims priority to U.S. Patent Appln. No. 61/113,004, filed on Nov. 10, 2008, and U.S. Patent Appln. No. 61/140,812, filed Dec. 24, 2008, the full contents of each being incorporated herein by reference.

[0002] The present invention relates generally to designing gauge patterns for calibration associated with a lithography process, and more specifically to computationally efficient methods for designing calibration pattern sets wherein individual patterns have significantly different responses to different parameter variations and are also extremely sensitive to parameter variations, and thus robust against random measurement errors in calibration.

[0003] Lithographic apparatuses are used, for example, in the manufacture of integrated circuits (ICs). In such a case, the mask contains a circuit pattern corresponding to an individual layer of the IC, and this pattern is imaged onto a target portion (e.g. comprising one or more dies) on a substrate (silicon wafer) that has been coated with a layer of radiation-sensitive material (resist). In general, a single wafer will contain a whole network of adjacent target portions that are successively irradiated via the projection system, one at a time. In one type of lithographic projection apparatus, each target portion is irradiated by exposing the entire mask pattern onto the target portion in one go; such an apparatus is commonly referred to as a wafer stepper. In an alternative apparatus, commonly referred to as a step-and-scan apparatus, each target portion is irradiated by progressively scanning the mask pattern under the projection beam in a given reference direction (the “scanning” direction) while synchronously scanning the substrate table parallel or anti-parallel to this direction. Since, in general, the projection system will have a magnification factor M (generally <1), the speed V at which the substrate table is scanned will be a factor M times that at which the mask table is scanned. More information with regard to lithographic devices as described herein can be gleaned, for example, from U.S. Pat. No. 6,046,792.

[0004] In a manufacturing process using a lithographic projection apparatus, a mask pattern is imaged onto a substrate that is at least partially covered by a layer of radiation-sensitive material (e.g. resist). Prior to this imaging step, the substrate typically undergoes various procedures, such as priming, resist coating and a soft bake. After exposure, the substrate may be subjected to other procedures, such as a post-exposure bake (PEB), development, a hard bake and measurement/inspection of the imaged features. This array of procedures is used as a basis to pattern an individual layer of a device, e.g., an IC. Such a patterned layer then undergoes various processes such as etching, ion-implantation (e.g. doping), metallization, oxidation, chemo-mechanical polishing, etc., all intended to finish off an individual layer. If several layers are required, then the whole procedure, or a variant thereof, will have to be repeated for each new layer. Eventually, an array of devices will be present on the substrate (i.e. wafer). These devices are then separated from one another by a technique such as dicing or sawing, whence the individual devices can be mounted on a carrier, connected to pins, etc.

[0005] For the sake of simplicity, the projection system is sometimes hereinafter referred to as the “lens”; however, this term should be broadly interpreted as encompassing various types of projection systems, including refractive optics, reflective optics, and catadioptric systems, for example. The radiation system also typically includes components operating according to any of these design types for directing, shaping or controlling the projection beam of radiation, and such components may also be referred to below, collectively or singularly, as a “lens”. Further, the lithographic apparatus may be of a type having two or more substrate tables (and/or two or more mask tables). In such “multiple stage” devices the additional tables can be used in parallel, and/or preparatory steps are carried out on one or more tables while one or more other tables are being used for exposures. Twin stage lithographic apparatus are described, for example, in U.S. Pat. No. 5,969,441.

[0006] The photolithographic masks referred to above comprise geometric patterns corresponding to the circuit components to be integrated onto a silicon wafer. The patterns used to create such masks are generated utilizing computer-aided design (CAD) programs, this process often being referred to as electronic design automation (EDA). Most CAD programs follow a set of predetermined design rules in order to create functional masks. These rules are set by processing and design limitations. For example, design rules define the space tolerance between circuit devices (such as gates, capacitors, etc.) or interconnect lines, so as to ensure that the circuit devices or lines do not interact with one another in an undesirable way. The design rule limitations are typically referred to as “critical dimensions” (CD). A critical dimension of a circuit can be defined as the smallest width of a line or hole or the smallest space between two lines or two holes. Thus, the CD determines the overall size and density of the designed circuit. Of course, one of the goals in integrated circuit fabrication is to faithfully reproduce the original circuit design on the wafer via the mask.

[0007] As noted, microlithography is a central step in the manufacturing of semiconductor integrated circuits, where patterns formed on semiconductor wafer substrates define the functional elements of semiconductor devices, such as microprocessors, memory chips etc. Similar lithographic techniques are also used in the formation of flat panel displays, micro-electro mechanical systems (MEMS) and other devices.

[0008] As semiconductor manufacturing processes continue to advance, the dimensions of circuit elements have continually been reduced while the amount of functional elements, such as transistors, per device has been steadily increasing over decades, following a trend commonly referred to as “Moore's law”. At the current state of technology, critical layers of leading-edge devices are manufactured using optical lithographic projection systems known as scanners that project a mask image onto a substrate using illumination from a deep-ultraviolet laser light source, creating individual circuit features having dimensions well below 100 nm, i.e. less than half the wavelength of the projection light.

[0009] This process, in which features with dimensions smaller than the classical resolution limit of an optical projection system are printed, is commonly known as low-k₁ lithography, according to the resolution formula CD=k₁×λ/NA, where λ is the wavelength of radiation employed (currently in most cases 248 nm or 193 nm), NA is the numerical aperture of the projection optics, CD is the “critical dimension”—generally the smallest feature size printed—and k₁ is an empirical resolution factor. In general, the smaller k₁, the more difficult it becomes to reproduce a pattern on the wafer that resembles the shape and dimensions planned by a circuit designer in order to achieve particular electrical functionality and performance. To overcome these difficulties, sophisticated fine-tuning steps are applied to the projection system as well as to the mask design. These include, for example, but not limited to, optimization of NA and optical coherence settings, customized illumination schemes, use of phase shifting masks, optical proximity correction in the mask layout, or other methods generally defined as “resolution enhancement techniques” (RET).

[0010] As one important example, optical proximity correction (OPC, sometimes also referred to as optical and process correction) addresses the fact that the final size and placement of a printed feature on the wafer will not simply be a function of the size and placement of the corresponding feature on the mask. It is noted that the terms “mask” and “reticle” are utilized interchangeably herein. For the small feature sizes and high feature densities present on typical circuit designs, the position of a particular edge of a given feature will be influenced to a certain extent by the presence or absence of other adjacent features. These proximity effects arise from minute amounts of light coupled from one feature to another. Similarly, proximity effects may arise from diffusion and other chemical effects during post-exposure bake (PEB), resist development, and etching that generally follow lithographic exposure.

[0011] In order to ensure that the features are generated on a semiconductor substrate in accordance with the requirements of the given target circuit design, proximity effects need to be predicted utilizing sophisticated numerical models, and corrections or pre-distortions need to be applied to the design of the mask before successful manufacturing of high-end devices becomes possible. The article “Full-Chip Lithography Simulation and Design Analysis-How OPC Is Changing IC Design”, C. Spence, Proc. SPIE, Vol. 5751, pp 1-14 (2005) provides an overview of current “model-based” optical proximity correction processes. In a typical high-end design, almost every feature edge requires some modification in order to achieve printed patterns that come sufficiently close to the target design. These modifications may include shifting or biasing of edge positions or line widths as well as application of “assist” features that are not intended to print themselves, but will affect the properties of an associated primary feature.

[0012] The application of model-based OPC to a target design requires good process models and considerable computational resources, given the many millions of features typically present in a chip design. However, applying OPC is generally not an exact science, but an iterative process that does not always resolve all possible weaknesses on a layout. Therefore, post-OPC designs, i.e. mask layouts after application of all pattern modifications by OPC and any other RET's, need to be verified by design inspection, i.e. intensive full-chip simulation using calibrated numerical process models, in order to minimize the possibility of design flaws being built into the manufacturing of a mask set. This is driven by the enormous cost of making high-end mask sets, which run in the multi-million dollar range, as well as by the impact on turn-around time by reworking or repairing actual masks once they have been manufactured.

[0013] Both OPC and full-chip RET verification may be based on numerical modeling (i.e. computational lithography) systems and methods as described, for example in, U.S. Pat. No. 7,003,758 and an article titled “Optimized Hardware and Software For Fast, Full Chip Simulation”, by Y. Cao et al., Proc. SPIE, Vol. 5754, 405 (2005).

[0014] As mentioned above, both OPC and RET require robust models that describe the lithography process precisely. Calibration procedures for such lithography models are thus required to achieve models that are valid, robust and accurate across the process window. Currently, calibration is done by actually printing a certain number of 1-dimensional and/or 2-dimensional gauge patterns on a wafer and performing measurements on the printed patterns. More specifically, those 1-dimensional gauge patterns are line-space patterns with varying pitch and CD, and the 2-dimensional gauge patterns typically include line-ends, contacts, and randomly selected SRAM (Static Random Access Memory) patterns. These patterns are then imaged onto a wafer and resulting wafer CDs or contact hole (also known as a via or through-chip via) energy are measured. The original gauge patterns and their wafer measurements are then used jointly to determine the model parameters which minimize the difference between model predictions and wafer measurements.

[0015] A model calibration process as described above and used in the prior art is illustrated in FIG. 3. In the prior art model calibration (FIG. 3), the process begins with a design layout 302, which can include gauges and other test patterns, and can also include OPC and RET features. Next, the design is used to generate a mask layout in step 304, which can be in a standard format such as GDSII or OASIS. Then two separate paths are taken, for simulation and measurement.

[0016] In a simulation path, the mask layout and a model 306 are used to create a simulated resist image in step 308. The model 306 provides a model of the lithographic process for use in computational lithography, and the calibration process aims to make the model 306 as accurate as possible, so that computational lithography results are likewise accurate. The simulated resist image is then used to determine predicted critical dimensions (CDs), contours, etc. in step 310.

[0017] In a measurement path, the mask layout 304 is used to form a physical mask (i.e. reticle), which is then imaged onto a wafer in step 312. The lithographic process (e.g. NA, focus, dose, illumination source, etc.) used to image the wafer is the same as that intended to be captured in model 306. Measurements (e.g. using metrology tools, etc.) are then performed on the actual imaged wafer in step 314, which yields measured CDs, contours, etc.

[0018] A comparison is made in step 316 between the measurements from step 314 and the predictions from step 310. If the comparison determines that the predictions match the measurements within a predetermined error threshold, the model is considered to be successfully calibrated in step 318. Otherwise, changes are made to the model 306, and steps 308, 310 and 316 are repeated until the predictions generated using the model 306 match the measurements within a predetermined threshold.

[0019] The inventors have noted that the design of gauge patterns such as those included in design layout 302 can greatly affect the accuracy of the model 306 and/or the time needed to successfully complete the calibration process. Unfortunately, the conventional art does not include a systematic study on how to determine the type or design of gauge patterns to be used for calibration. For example, there is no theoretical guidance on the choice of pitch and CD for the line-space patterns or the number of gauges. In current practice, the selection of gauge patterns is rather arbitrary—they are often chosen from experience or randomly chosen from the real circuit patterns. Such gauge patterns are often incomplete or super-complete or both for calibration. For example, none of the chosen gauge patterns will effectively discriminate between certain of the model parameters, thus it may be difficult to determine the parameter values due to measurement inaccuracies. On the other hand, many patterns can yield very similar responses to different parameter variations, thus some of them are redundant and wafer measurements on these redundant patterns waste resources.

[0020] The present invention relates to computationally efficient techniques for designing gauge patterns for calibrating a lithographic process model for use in a simulation process. According to some aspects, the present invention relates to methods of designing gauge patterns that are extremely sensitive to parameter variations, and thus robust against random measurement errors in calibration of a lithographic process utilized to image a target design having a plurality of features. In some embodiments, the process includes identifying a most sensitive line width/pitch combination with optimal assist feature placement which leads to most sensitive CD changes against parameter variations. Another process includes designing gauges which have two main features. The difference between the two main features' CDs is extremely sensitive against parameter variation, and thus robust against random measurement error and any measurement error in bias. In yet another process, patterns are designed that lead to most sensitive intensity.

[0021] According to further aspects, the invention includes methods for designing gauges which minimize the above-described degeneracy, and thus maximize pattern coverage for model calibration. More specifically, the present invention relates to methods of designing gauge patterns that achieve complete coverage of parameter variations with minimum number of gauges and corresponding measurements in the calibration of a lithographic process model utilized to simulate imaging of a target design having a plurality of features. According to some aspects, a method according to the invention includes transforming a model parametric space into a new space (based on CD sensitivity or Delta TCCs), then iteratively identifying the direction that is most orthogonal to existing gauges' CD sensitivities in this new space, and determining a most sensitive line width/pitch combination with optimal assist feature placement which leads to most sensitive CD changes along that direction in a model parametric space.

[0035] The present invention will now be described in detail with reference to the drawings, which are provided as illustrative examples of the invention so as to enable those skilled in the art to practice the invention. Notably, the figures and examples below are not meant to limit the scope of the present invention to a single embodiment, but other embodiments are possible by way of interchange of some or all of the described or illustrated elements. Moreover, where certain elements of the present invention can be partially or fully implemented using known components, only those portions of such known components that are necessary for an understanding of the present invention will be described, and detailed descriptions of other portions of such known components will be omitted so as not to obscure the invention. Embodiments described as being implemented in software should not be limited thereto, but can include embodiments implemented in hardware, or combinations of software and hardware, and vice-versa, as will be apparent to those skilled in the art, unless otherwise specified herein. In the present specification, an embodiment showing a singular component should not be considered limiting; rather, the invention is intended to encompass other embodiments including a plurality of the same component, and vice-versa, unless explicitly stated otherwise herein. Moreover, applicants do not intend for any term in the specification or claims to be ascribed an uncommon or special meaning unless explicitly set forth as such. Further, the present invention encompasses present and future known equivalents to the known components referred to herein by way of illustration.

[0036] Prior to discussing the present invention, a brief discussion regarding the overall simulation and imaging process to be calibrated is provided. FIG. 1 illustrates an exemplary lithographic projection system 10. The major components are a light source 12, which may be a deep-ultraviolet excimer laser source, illumination optics which define the partial coherence (denoted as sigma) and which may include specific source shaping optics 14, 16a and 16b; a mask or reticle 18; and projection optics 16c that produce an image of the reticle pattern onto the wafer plane 22. An adjustable filter or aperture 20 at the pupil plane may restrict the range of beam angles that impinge on the wafer plane 22, where the largest possible angle defines the numerical aperture of the projection optics NA=sin(Θ_max).

[0037] In a lithography simulation system, these major system components can be described by separate functional modules, for example, as illustrated in FIG. 2. Referring to FIG. 2, the functional modules include the design layout module 26, which defines the target design; the mask layout module 28, which defines how the mask is laid out in polygons based on the target design; the mask model module 30, which models the physical properties of the pixilated and continuous-tone mask to be utilized during the simulation process; the optical model module 32, which defines the performance of the optical components of lithography system; and the resist model module 34, which defines the performance of the resist being utilized in the given process. As is known, the result of the simulation process produces, for example, predicted contours and CDs in the result module 36.

[0038] More specifically, it is noted that the properties of the illumination and projection optics are captured in the optical model module 32 that includes, but is not limited to, NA-sigma (σ) settings as well as any particular illumination source shape, where σ (or sigma) is outer radial extent of the illuminator. The optical properties of the photo-resist layer coated on a substrate—i.e. refractive index, film thickness, propagation and polarization effects—may also be captured as part of the optical model module 32, whereas the resist model module 34 describes the effects of chemical processes which occur during resist exposure, PEB and development, in order to predict, for example, contours of resist features formed on the substrate wafer. The mask model module 30 captures how the target design features are laid out in the reticle and may also include a representation of detailed physical properties of the mask, as described, for example, in U.S. Pat. No. 7,587,704. The objective of the simulation is to accurately predict, for example, edge placements and critical dimensions (CDs), which can then be compared against the target design. The target design is generally defined as the pre-OPC mask layout, and will be provided in a standardized digital file format such as GDSII or OASIS.

[0039] In general, the connection between the optical and the resist model is a simulated aerial image intensity within the resist layer, which arises from the projection of light onto the substrate, refraction at the resist interface and multiple reflections in the resist film stack. The light intensity distribution (aerial image intensity) is turned into a latent “resist image” by absorption of photons, which is further modified by diffusion processes and various loading effects. Efficient simulation methods that are fast enough for full-chip applications approximate the realistic 3-dimensional intensity distribution in the resist stack by a 2-dimensional aerial (and resist) image.

[0040] As should be therefore apparent from the above, the model formulation describes all of the known physics and chemistry of the overall process, and each of the model parameters preferably corresponds to a distinct physical or chemical effect. The model formulation thus sets an upper bound on how well the model can be used to simulate the overall lithography process. However, sometimes the model parameters may be inaccurate from measurement and reading errors, and there may be other imperfections in the system. With precise calibration of the model parameters, extremely accurate simulations can be done. In other words, the calibration of modem models is probably a larger factor in accuracy than the theoretical upper bound of accuracy.

[0041] There are various ways that the model parameters can be expressed. One efficient implementation of a lithography model is possible using the following formalism, where the image (here in scalar form, but which may be extended to include polarization vector effects) is expressed as a Fourier sum over signal amplitudes in the pupil plane. According to the standard Hopkins theory, the aerial image (AI) intensity may be defined by:

[0042] $\begin{array}{cc}\begin{array}{c}I\left(x\right)=\sum _k{A\left(k\right)\sum _{k^{\prime }}M\left(k^{\prime }-k\right)P\left(k^{\prime }\right)\exp \left(-j{k}^{\prime }x\right)}^2\\ =\sum _k{A\left(k\right)}^2\left\{\sum _{k^{\prime }}\sum _{k^{″}}M\left(k^{\prime }-k\right)P\left(k^{\prime }\right)M^{*}\left(k^{″}-k\right)P^{*}\left(k^{″}\right)\exp \left(-j\left(k^{\prime }-k^{″}\right)x\right)\right\}\\ =\sum _{k^{\prime }}\sum _{k^{″}}\left[\sum _k{A\left(k\right)}^{2}P\left(k+k^{\prime }\right)P^{*}\left(k+k^{″}\right)\right]\\ M\left(k^{\prime }\right)M^{*}\left(k^{″}\right)\exp \left(-j\left(k^{\prime }-k^{″}\right)x\right)\\ =\sum _{k^{\prime }}\sum _{k^{″}}TC{C}_{k^{\prime },k^{″}}M\left(k^{\prime }\right)M^{*}\left(k^{″}\right)\exp \left(-j\left(k^{\prime }-k^{″}\right)x\right)\end{array}& \left(\mathrm{Eq}.1\right)\end{array}$
where I (x) is the aerial image intensity at point x within the image plane (for notational simplicity, a two-dimensional coordinate represented by a single variable is utilized), k represents a point on the source plane, A (k) is the source amplitude from point k, k′ and k″ are points on the pupil plane, M is the Fourier transform of the mask image, P is the pupil function (M* and P* are the complex conjugates of M and P, respectively).

[0043] An important aspect of the foregoing derivation is the change of summation order (moving the sum over k inside) and indices (replacing k′ with k+k′ and replacing k″ with k+k″), which results in the separation of the Transmission Cross Coefficients (TCCs), defined by the term inside the square brackets in the third line in the equation and as shown in the fourth line, from the other terms. In other words:
TCC_k′,k″=Σ_kA(k)² P(k+k′) P*(k+k″).

[0044] These transmission cross coefficients are independent of the mask pattern and therefore can be pre-computed using knowledge of the optical elements or configuration only (e.g., NA and σ or the detailed illuminator profile). It is further noted that although in the given example (Eq. 1) is derived from a scalar imaging model for ease of explanation, those skilled in the art can also extend this formalism to a vector imaging model, where TE and TM polarized light components are summed separately.

[0045] For simplicity, note that the relationship between aerial image intensity and TCCs, i.e., (Eq. 1) can be expressed as a bilinear operator:
I(x)=M*TCC*M*

[0046] Furthermore, the approximate aerial image intensity I can be calculated by using only a limited number of dominant TCC terms, which can be determined by diagonalizing the TCC matrix and retaining the terms corresponding to its largest eigenvalues, i.e.,

[0047] $\begin{array}{cc}TC{C}_{k^{\prime },k^{″}}=\sum _{n=1}^N{\lambda }_n{\varphi }_n\left(k^{\prime }\right){\varphi }_n^{*}\left(k^{″}\right)& \left(\mathrm{Eq}.2\right)\end{array}$
where λ_n(n=1, . . . , N) denotes the N largest eigenvalues and φ_n(●) denotes the corresponding eigenvector of the TCC matrix. It is noted that (Eq. 2) is exact when all terms are retained in the eigenseries expansion, i.e., when N is equal to the rank of the TCC matrix. However, in actual applications, it is typical to truncate the series by selecting a smaller N to increase the speed of the computation process.

[0048] Thus, (Eq. 1) can be rewritten as:

[0049] $\begin{array}{cc}\begin{array}{c}I\left(x\right)=\sum _{k^{\prime }}\sum _{k^{″}}TC{C}_{k^{\prime },k^{″}}M\left(k^{\prime }\right)M\left(k^{″}\right)\exp \left(j\left(k^{\prime }-k^{″}\right)x\right)\\ =\sum _{k^{\prime }}\sum _{k^{″}}\sum _{n=1}^N{\lambda }_n{\varphi }_n\left(k^{\prime }\right){\varphi }_n^{*}\left(k^{″}\right)M\left(k^{\prime }\right)M^{*}\left(k^{″}\right)\exp \left(j\left(k^{\prime }-k^{″}\right)x\right)\\ =\sum _{n=1}^N{\lambda }_n\sum _{k^{\prime }}{\varphi }_n\left(k^{\prime }\right)M\left(k^{\prime }\right)\exp \left(j{k}^{\prime }x\right)\\ \sum _{k^{″}}{\varphi }_n^{*}\left(k^{″}\right)M^{*}\left(k^{″}\right)\exp \left(-j{k}^{″}x\right)\\ =\sum _{n=1}^N{\lambda }_n{{\Phi }_n\left(x\right)}^2\end{array}& \left(\mathrm{Eq}.3\right)\end{array}$
where

[0050] ${\Phi }_n\left(x\right)=\sum _{k^{″}}{\varphi }_n\left(k^{″}\right)M^{*}\left(k^{″}\right)\exp \left(j{k}^{″}x\right)$
and |●| denotes the magnitude of a complex number.

[0051] As should be apparent from the above, the lithography simulation model parameters are captured in the TCCs. Accordingly, calibrating model parameters in embodiments of the invention is achieved by obtaining highly accurate raw TCC terms (i.e. before diagonalization). However, the invention is not limited to this example embodiment.

Univariate Parameter-Sensitive Calibration Pattern Set Design

[0052] Note that a gauge comprises one or more calibration patterns. When there are more than one calibration pattern in one gauge, the multiple patterns are typically duplicates of the same pattern. This allows users to take several measurements and then perform averaging to reduce random measurement error.

[0053] In first embodiments of the invention, gauge design methods are provided to maximize sensitivity of a given metric against variations in a single parameter. When the sensitivity is maximized, the robustness against random measurement error is also maximized, which minimizes the number of measurements. To further illustrate this point, consider the following example. Assume that parameter P is calibrated from CD measurements, as in a typical calibration process. Suppose the true value of P is P₀, and the value of P is estimated from L gauges. Further assume using a brute-force approach to calibrate parameter P, i.e., there are a set S_P of possible values for P, denoted as
S_P={P₀, P₁, P₂, . . . }

[0054] For every element P′∈S_P, the resulting CD, denoted as CD_l(P′) can be simulated for the l-th calibration pattern (l=1, . . . , L), for example using the numerical modeling systems and methods as described in U.S. Pat. No. 7,003,758 and an article titled “Optimized Hardware and Software For Fast, Full Chip Simulation”, by Y. Cao et al., Proc. SPIE, Vol. 5754, 405 (2005),

[0055] Next the error between measured CDs and the simulated CDs is computed for every P′∈S_P, and the value P* is chosen that minimizes the error among all elements in S_P. Typically Mean Squared Error (MSE) is used to measure such error, in which case P* can then be formulated as

[0056] $\begin{array}{cc}{P}^{*}=\arg \underset{P^{\prime }\in S_p}{\min }\frac{1}{L}\sum _{l=1}^L{\left({\widetilde{\mathrm{CD}}}_l-{\mathrm{CD}}_l\left(P^{\prime }\right)\right)}^2& \left(\mathrm{Eq}.4\right)\end{array}$
where {tilde over (C)}{tilde over (D)}_l is the measured CD value of the l-th gauge (l=1, . . . , L), and the function arg min (or arg max) represents the argument of the minimum (or maximum), i.e., the value of the given argument for which the value of the given object function achieves its minimum (or maximum).

[0057] In the ideal world where there is no CD measurement error, and the computation of CD_l(P′) is carried out perfectly without any approximation or numerical error, such that
{tilde over (C)}{tilde over (D)}_l=CD_l(P₀)
for all L gauges, and P′=P₀ will always lead to the minimum error (0) between the measured CDs and simulated CDs from the true parameter value.

[0058] Unfortunately, in real applications, there is always some (random) deviation between the measured CDs and the simulated CDs. Then for the l-th gauge (l=1, . . . , L), the measured CD value is
{tilde over (C)}{tilde over (D)}_l=CD_l(P₀)+E_l
where E_l is the deviation for the l-th gauge.

[0059] For small parameter variations, it can be assumed that the relationship between CD and P is linear, i.e., for the l-th gauge,
CD_l(P)=CD_l(P₀)+K_l(P−P₀)
as disclosed in U.S. Patent Publ. No. 2009/0157360, commonly owned by the present assignee. Here the CD sensitivity (i.e. one example of a metric) against parameter P is defined as:

[0060] $K_l=\frac{\partial {\mathrm{CD}}_l\left(P\right)}{\partial P}{❘}_{P=P_0}$
which is a constant for each gauge when the variation of P is small. In this case, (Eq. 4) can be written as

[0061] $\begin{array}{c}{P}^{*}=\arg \underset{P^{\prime }\in S_p}{\min }\frac{1}{L}\sum _{l=1}^L{\left({\widetilde{\mathrm{CD}}}_l-{\mathrm{CD}}_l\left(P^{\prime }\right)\right)}^2\\ =\arg \underset{P^{\prime }\in S_p}{\min }\frac{1}{L}\sum _{l=1}^L{\left(K_l\left(P^{\prime }-P_0\right)-E_l\right)}^2\\ =\arg \underset{P^{\prime }\in S_p}{\min }\frac{1}{L}\left[{\left(P^{\prime }-P_0\right)}^2\sum _{l=1}^L{K}_l^2-2\left(P^{\prime }-P_0\right)\sum _{l=1}^L\left(K_l{E}_l\right)+\sum _{l=1}^L{E}_l^2\right]\end{array}$

[0062] Taking the derivative of the expression inside the square brackets above with respect to P′ and setting it to zero results in the following equation:

[0063] $2\left(P^{*}-P_0\right)\sum _{l=1}^L{K}_l^2-2\sum _{l=1}^L\left(K_l{E}_l\right)=0$
Therefore,

[0064] $\begin{array}{cc}{P}^{*}=P_0+\frac{\sum _{l=1}^L\left(K_l{E}_l\right)}{\sum _{l=1}^K{K}_l^2}.& \left(\mathrm{Eq}.5\right)\end{array}$

[0065] As can be noticed, because of the E_l term, the calibration value of the parameter is no longer its true value P₀. The calibration error is denoted as

[0066] ${\Delta }_P=\frac{\sum _{l=1}^L\left(K_l{E}_l\right)}{\sum _{l=1}^L{K}_l^2}.$

[0067] From (Eq. 5), it can be seen that the smaller Δ_P, the better the calibration result. If the absolute value of the CD sensitivity K_l can be increased, the calibration error can be reduced. For example, if the CD sensitivity for each gauge is increased to NK_l, then the calibration error is reduced to 1/N of the original. Or vice verse, if it is desired to maintain the same calibration precision, one can reduce the number of gauges and number of CD measurements, which results in lower cost. One can also do a simple probability analysis on (Eq. 5) where it is assumed that such deviation is an additive random variable with mean 0 and variance σ² and independent of each other. One can further assume that all the features have same CD sensitivity K_l, i.e., K_l=K, for l=1, . . . , L. Thus,

[0068] ${\Delta }_P=\frac{\sum _{l=1}^L{E}_l}{\mathrm{KL}}$
Then the variance of calibration error is

[0069] $E\left({\Delta }_P^2\right)=\frac{\sum _{l=1}^{L}E\left(E_l^2\right)}{K^2{L}^2}=\frac{{\sigma }^2}{K^{2}L}$

[0070] Again, increasing CD sensitivity of each pattern will reduce the variance of the calibration error.

[0071] Next, the concepts of Delta I and Delta TCC are introduced herein. More particularly, given a nominal condition, the nominal TCC, nominal aerial image intensity and nominal CD can be computed. To study the CD sensitivity, a small perturbation is added to a certain parameter, and then the CD change (or Delta CD) under such parameter perturbation is studied. Suppose the nominal TCC is TCC_NC, then the nominal aerial image intensity I_NC is
I_NC=M*TCC_NC*M^*
from (Eq. 1).

[0072] Then suppose the TCC with parameter p perturbed is TCC_p, then the perturbed aerial image intensity I_p is:
I_p=M*TCC_p*M^*.

[0073] In example embodiments described below, the focus is on aerial image CD change (i.e. CD sensitivity is the observed metric). Therefore, the focus can be on aerial image intensity change near the aerial image threshold point, as explained in more detail below. If a gauge has a large aerial image intensity change near the threshold point when the parameter is perturbed, it is very likely that it also has a large CD sensitivity against that parameter variation.

[0074] To further simplify the problem, the difference between the perturbed aerial image intensity (which corresponds to perturbed parameter values) and nominal aerial image intensity (which corresponds to nominal parameter values) is studied. This is called the Delta I or ΔI_p, i.e.,
ΔI_p=I_p−I_NC=M*(TCC_p−TCC_NC)*M^*

[0075] The term Delta TCC refers to the term in parentheses above, i.e.:
ΔTCC_p=TCC_p−TCC_NC

[0076] Delta TCC (i.e. ΔTCC_p in the equation above) is computed by taking the element-by-element matrix difference between TCC_p and TCC_NC. It should be noted that Delta TCC (i.e. ΔTCC_p), TCC_p and TCC_NC are all in “raw” matrix form, before any diagonalization.

[0077] Then the Delta I can be viewed as an aerial image intensity resulting from the Delta TCC and the original mask image, i.e., the bilinear operator on the mask image and the Delta TCC. In an embodiment, similar to (Eq. 2) and (Eq. 3), eigen-decomposition of the Delta TCC is performed to speed up the computation of Delta I. In the embodiment, the eigenvalues of the Delta TCC may be negative, so the eigenvalues with largest absolute values should be kept.

[0078] With this foundation, various embodiments for designing 1D maximum sensitive test gauges according to the invention will now be described. In a first example embodiment, the methodology includes finding a line/space pattern and associated assist features that maximize sensitivity for each parameter. A general description of this embodiment will be provided in connection with the flowchart in FIG. 4. However, further details of this and alternative embodiments will be provided following the description of FIG. 4.

[0079] First, in step S402, a nominal condition setting is identified, which is a set of parameter values for all relevant parameters. Next in step S404, a loop begins for each parameter. First, in step S406, a perturbation value for the parameter is applied, typically a small amount which yields a Delta TCC as described above. In some embodiments, the perturbation value is an amount that will cause the CD change of the most important pattern to be around 10 nm. Here, use is made of the fact that the relationship between CD change and small parameter variations is quite linear. Also, the CD sensitivity of a given pattern is defined as the difference between the nominal CD value and the perturbed CD value, i.e., the CD value when that particular parameter value is the sum of its nominal value and a perturbation. As set forth above, this can be inherently identified from the Delta I resulting from this Delta TCC.

[0080] So far, the actual design of the calibration patterns is not yet involved, as previous steps concentrated on the imaging. In step S408, a loop begins for each possible or allowed value of the parameters (such as line width or pitch of a calibration grating). In step S410, for the current line/space pattern, place an assisting feature to maximize the CD sensitivity against that parameter. Next, after looping over all possible main feature line width and pitch combinations (including the opposite tones) as described in more detail below, in step S412, the combination leading to the maximum (absolute) CD sensitivity is identified. This combination corresponds to a calibration pattern for a perturbed state. When used in combination with the calibration pattern corresponding to the nominal values of the parameters, the model can be calibrated.

[0081] After optimal calibration patterns with assisting features have been generated for all parameters, the set of resulting calibration patterns are output in step S414. This set includes a calibration pattern according to the nominal parameters to enable comparison. By using the same nominal values for all parameters, only a single calibration pattern can be used for respectively calibrating all parameters. In an embodiment, a plurality of calibration patterns is used to allow multiple measurements to be taken in order to improve accuracy of the calibration.

[0082] In this embodiment, it is assumed that the changes in aerial image intensity and thus the aerial image CD is representative of the changes in wafer CD. In terms of the model (e.g. FIG. 2), the resist model 34 is not used (as explained above, the connection between the optical model and the resist model is a simulated aerial image intensity within the resist layer.) Thus, one can identify the optimal assisting feature locations in step S410 purely based on TCCs. This assumption reduces the design complexity considerably, as will become more apparent from the descriptions below.

[0083] Further, the exhaustive search of the most sensitive calibration pattern sets over all possible main feature line width and pitch as set forth in step S408 may be too expensive. Accordingly, embodiments of the invention employ a few approximations to reduce the computational complexity.

[0084] First, for identifying the optimal main feature pitch, the present inventors recognize through simulations that it is not necessary to loop through all possible pitch values in step S408. The reason is, for each gauge, if only the central main feature is measured, other main features can also be viewed as assisting features to the central feature. Thus, it is theoretically optimal to start with an infinite pitch, i.e., an isolated main feature or a line space pattern with a very large pitch. This is also consistent with the observation that if the assisting feature placement methodology is optimal, then the more space to place assisting features, the higher the CD sensitivity, because more assisting features are contributing to the sensitivity. However, in real applications, it is often preferred that a gauge contains a few repetitive main features, so that measurements can be taken at different locations then averaged in order to reduce (random) measurement error. Thus it is desired to have a pitch ¼ or ⅕ of the gauge width (which is typically a few micrometers), so that the pitch is still quite large and four or five measurements can be taken for each gauge.

[0085] For the identification of most sensitive main feature line width, the present inventors recognize that it is possible to separate the process of identifying most sensitive CD from the process of adding assisting features without hurting the performance. More particularly, they note that it is reasonable to assume that the line width having the most sensitive CD without any assisting feature, when added with assisting features, is close to the line width having the most sensitive CD with assisting features. Accordingly, another approximation is to first search the most sensitive main feature CD without any assisting feature, and then add assisting features to the isolated main feature with that CD value.

[0086] The present inventors further recognize that it is also possible to pre-determine the most optimal line width for the main feature without an exhaustive search. Suppose the delta TCC matrix for a 2D image is ΔTCC(k_x1, k_y1, k_x2, k_y2) The Delta TCC matrix for a 1D Line Space pattern can be recomposed as ΔTCC(k_x1, k_x2)=ΔTCC(k_x1, 0, k_x2, 0), and the corresponding kernel in the image plane is ΔW(x, y). The indices k_x1, k_y1, k_x2, k_y2 are the two-dimensional indices corresponding to the frequency-domain indices k′, k″ as set forth above.

[0087] Then, by changing (Eq. 1) from pupil plane to image plane, or equivalently, changing from the frequency domain to the spatial domain using an inverse Fourier transform process, the Delta I for a Line Space mask pattern M(x) becomes
ΔI(x)=∫∫M(ξ₁)M^*(ξ₂)ΔW(x−ξ₁, x−ξ₂)dξ₁dξ₂

[0088] Assume that the line width in pattern M(x) is LW, i.e., M(x) is an isolated line with width of LW. Further, without loss of generality, assume a darkfield mask without attenuation (the results for other scenarios are very similar), then

[0089] $M\left(x\right)=\{\begin{array}{cc}1& \mathrm{if}x\le \mathrm{LW}/2\\ 0& \mathrm{otherwise}\end{array}$

[0090] The method above is only interested in the aerial image intensity near an aerial image threshold. This threshold corresponds to a minimum intensity to be received by the resist layer (photoactive layer) used in the lithographic process before it is activated. Here, for simplicity, it can be assumed that that point is close to the two mask image edge locations, x=±LW/2. Further, because of symmetry, one need only look at the single point x=LW/2, i.e., one can focus on ΔI (LW/2), which can be simplified as

[0091] $\begin{array}{c}\Delta I\left(\mathrm{LW}/2\right)={\int }_{{\xi }_1=-\mathrm{LW}/2}^{\mathrm{LW}/2}{\int }_{{\xi }_2=-\mathrm{LW}/2}^{\mathrm{LW}/2}\Delta W\left(\mathrm{LW}/2-{\xi }_1,\mathrm{LW}/2-{\xi }_2\right)d{\xi }_{1}d{\xi }_2\\ ={\int }_{{\xi }_1=0}^{\mathrm{LW}}{\int }_{{\xi }_2=0}^{\mathrm{LW}}\Delta W\left({\xi }_1,{\xi }_2\right)d{\xi }_{1}d{\xi }_2\end{array}$

[0092] This formula relates the CD sensitivity, as expressed by ΔI (LW/2), to the line width LW of the main feature, and thus enables the determination of the line width that maximizes the CD sensitivity. Accordingly, with the process for determining optimal line width and pitch being simplified as set forth above, steps S408 to S412 of the general methodology for calibration pattern design described in connection with FIG. 4 can be implemented by identifying the most sensitive line width and pitch for a given parameter (for both opposite tones) using the above equation, and, for the main features with identified line width and pitch, place assisting features to maximize the CD sensitivity against that parameter.

[0093] For this latter step, the present inventors have further developed advantageous methods for placing assisting features to maximize the CD sensitivity for a given main feature. In embodiments, this includes identifying an Assisting Feature Guidance Map (AFGM), which is similar to the SGM described in the PCT application published as WO 2008-151185. The following describes two alternative methods to compute AFGM.

[0094] The first method is called a single kernel approach. From (Eq. 3), one can also express the Delta I computation in the space domain as:
ΔI=L₁*(MF₁)²+L₂*(MF₂)²+ . . . +L_N*(MF_N)²
where: M is the mask image in the space domain (which is typically real); N is the number of eigenvalues of the Delta TCC; F₁ to F_N are the real-space filters corresponding to each TCC term (i.e., the real part of the inverse Fourier Transform of φ₁ to φ_N); L₁ to L_N are the corresponding eigenvalues of each Delta TCC term; “” means convolution, and “*” is the regular multiplication. One can assume without loss of generality that |L₁|≧|L₂| . . . ≧|L_N|.

[0095] In the single kernel approach, the emphasis is on the aerial image intensity from the kernel corresponding to the eigenvalues with the largest absolute value. Then by ignoring the scaling factor L₁:
ΔI≈(MF)²
where F=F(x,y) is a scalar field, and can be well approximated by F₁ in the “near-coherent” case, i.e., |L_n|/|L₁|<<1, for any n=2, 3, . . . , N.

[0096] For each field point x′ on the mask, this approach places a hypothetical point source δ(x−x′) as an assisting feature and studies the contribution from this point source to the change in aerial image intensity Delta I around the mask edges. If the contribution is positive, then it implies that the change (Delta AI or ΔI) in aerial image intensity will increase if the assisting features contain this point. This means that adding this point to the assist features for the calibration pattern corresponding to the perturbed parameter values contributes to the sensitivity of the set of calibration patterns to parameter changes. Thus, the assist features for the calibration pattern corresponding to the perturbed parameter values should comprise this point.

[0097] For each field point x′, the change (Delta I) in aerial image intensity with the point source is
ΔI_x′=((M+δ(x−x′))F)²

[0098] Note that the convolution operation is linear, thus the change related to the field point in the change (Delta I) related to different parameter values in aerial image intensity, caused by placing the point source is

[0099] $\begin{array}{c}\Delta I_{x^{\prime }}-\Delta I={\left(\left(M+\delta \left(x-x^{\prime }\right)\right)\otimes F\right)}^2-{\left(M\otimes F\right)}^2\\ {=2\left(\delta \left(x-x^{\prime }\right)\right)\otimes F)\left(M\otimes F\right)+\left(\delta \left(x-x^{\prime }\right)\right)\otimes F)}^2\\ {=2F\left(x-x^{\prime }\right)\left(M\otimes F\right)+\left(F\left(x-x^{\prime }\right)\right))}^2\end{array}$

[0100] Assuming a real mask, then the AFGM, which is the point source's contribution to all mask edge locations, is

[0101] $\begin{array}{c}AFGM\left(x^{\prime }\right)={\int }_x\frac{dM\left(x\right)}{dx}\left[\Delta I_{x^{\prime }}-\Delta I\right]dx\\ ={\int }_x\frac{dM\left(x\right)}{dx}\left[2F\left(x-x^{\prime }\right)\left(M\otimes F\right)+{\left(F\left(x-x^{\prime }\right)\right)}^2\right]dx\\ =2\left(\frac{dM\left(x\right)}{dx}\left(M\otimes F\right)\right)\otimes F\left(-x\right)+\frac{dM\left(x\right)}{dx}\otimes F^2\left(-x\right)\end{array}$

[0102] With Fourier transforms, one can replace convolutions in the space domain by multiplication in the frequency domain, such that:

[0103] $AFGM\left(x^{\prime }\right)=2IFFT\left\{FFT\left\{\frac{dM\left(x\right)}{dx}\left(M\otimes F\right)\right\}FFT\left\{F\left(-x\right)\right\}\right\}+IFFT\left\{FFT\left\{\frac{dM\left(x\right)}{dx}\right\}FFT\left\{F^2\left(-x\right)\right\}\right\}$
where FFT(.) is the Fourier Transform operation and IFFT(.) is the inverse Fourier Transform operation. One advantageous thing about the frequency domain operation is that FFT{F(−x) } and FFT{F²(−x)} are independent of the mask, thus as soon as the optical condition is fixed, they can be pre-computed.

[0104] A second embodiment to compute AFGM is called multi-kernel approach. In the multi-kernel approach, mask transmittance M(x) is separated into a pre-OPC component (M^T), an assisting feature (AF) component (M^A) and an OPC corrections component (M^C), i.e.:
M(x)=M^T(x)+M ^A(x)+M^C(x)
If
M^k(x)=M^T(x)+M^C(x)
represents the post-OPC layout transmittance, then by applying the inverse Fourier Transform (i.e., space domain representation) of (Eq. 1), the change in aerial image intensity (Delta I) is

[0105] $\begin{array}{c}\Delta I\left(x\right)=\int \left[M^K\left(x_1\right)+M^A\left(x_1\right)\right]\left[M^{K^{*}}\left(x_2\right)+M^{A*}\left(x_2\right)\right]\\ \Delta W\left(x-x_1,x-x_2\right)d{x}_{1}d{x}_2\\ =\int \left[M^K\left(x_1\right)M^{K^{*}}\left(x_2\right)+M^A\left(x_1\right)M^{K^{*}}\left(x_2\right)+M^K\left(x_1\right)M^{A*}\left(x_2\right)+M^A\left(x_1\right)M^{A*}\left(x_2\right)\right]\\ \Delta W\left(x-x_1,x-x_2\right)d{x}_{1}d{x}_2\\ =\Delta I^K\left(x\right)+\int \left[M^A\left(x_1\right)M^{K^{*}}\left(x_2\right)+M^K\left(x_1\right)M^{A^{*}}\left(x_2\right)\right]\\ \Delta W\left(x-x_1,x-x_2\right)d{x}_{1}d{x}_2\end{array}$
where ΔW(x,y) is the space domain representation of the Delta TCC and ΔI^T(x) is the Delta AI change in aerial image intensity (Delta I) without assisting features. In practice, the inventors note that the following term of the above equation can be ignored:
∫M^A(x₁)M^A*(x₂)ΔW(x−x₁, x−x₂)dx₁dx₂
Because M^A (associated with the AF component) is typically small compared to M^K.

[0106] Moreover, to derive the AFGM expression from the remaining terms, a unit source at x′ in the AF portion of the mask layout is assumed, i.e.,
M^A(x)=δ(x−x′).
This unit source at x′ contributes the following amount to the change in aerial image intensity (Delta I) at x:

[0107] $\begin{array}{c}\Delta I\left(x\right)-\Delta I^K\left(x\right)=\int \left[M^A\left(x_1\right)M^{K^{*}}\left(x_2\right)+M^K\left(x_1\right)M^{A^{*}}\left(x_2\right)\right]\\ \Delta W\left(x-x_1,x-x_2\right)d{x}_{1}d{x}_2\\ =\int \left[\delta \left(x_1-x^{\prime }\right)M^{K^{*}}\left(x_2\right)+M^K\left(x_1\right)\delta \left(x_2-x^{\prime }\right)\right]\\ \Delta W\left(x-x_1,x-x_2\right)d{x}_{1}d{x}_2\\ =\int M^{K^{*}}\left(x_2\right)\Delta W\left(x-x^{\prime },x-x_2\right)d{x}_2+\\ \int M^K\left(x_1\right)\Delta W\left(x-x_1,x-x^{\prime }\right)d{x}_1\end{array}$

[0108] The weighting of the vote from field point x to source point x′ is equal to the gradient of the pre-OPC image, such that

[0109] $\frac{d{M}^R\left(x\right)}{dx}=\frac{1}{2}\frac{d}{dx}\left[M^T\left(x\right)+M^{T^{*}}\left(x\right)\right]$

[0110] For AFGM, what needs to be determined is whether this point source, as an assisting feature, would enhance or weaken the change in aerial image intensity (Delta I) from main features only, at all locations near the aerial image intensity threshold. So at each location, the process multiplies the contribution from the point source by the change in aerial image intensity (Delta I) without any assisting feature. After summing this value over all aerial image intensity contour locations, a positive AFGM value implies that this point will enhance the CD sensitivity, and vice versa for a negative value. Assuming an OPC process is performed such that the aerial image intensity contour after OPC matches the pre-OPC edge locations, then one can sum the contributions from the point source over all points where the gradient of the pre-OPC mask image is nonzero. As a result, the AFGM value at x′ is equal to

[0111] $V\left(x^{\prime }\right)=\int {\frac{d{M}^R\left(x\right)}{dx}}^2\Delta I^K\left(x\right)\left[\Delta I\left(x\right)-\Delta I^K\left(x\right)\right]dx=\int {\frac{d{M}^R\left(x\right)}{dx}}^2\Delta I^K\left(x\right)\left[\int M^{K*}\left(x_2\right)\Delta W\left(x-x^{\prime },x-x_2\right)d{x}_2+\int M^K\left(x_1\right)\Delta W\left(x-x_1,x-x^{\prime }\right)d{x}_1\right]dx$

[0112] For simplicity, let

[0113] $R\left(x\right)={\frac{d{M}^R\left(x\right)}{dx}}^2\Delta I^K\left(x\right)$

[0114] Then with a change of variables, i.e., x=x′−ζ₁, x₂=x′−ζ₂, for the first integration inside the above brackets and x₁=x′−ζ₁, x=x′−ζ₂, for the second integration inside the above brackets, then

[0115] $V\left(x^{\prime }\right)=\int R\left(x\right)M^{K^{*}}\left(x_2\right)\Delta W\left(x-x^{\prime },x-x_2\right)dxd{x}_2+\int R\left(x\right)M^K\left(x_1\right)\Delta W\left(x-x_1,x-x^{\prime }\right)dxd{x}_1=\int R\left(x^{\prime }-{\xi }_1\right)M^{K^{*}}\left(x^{\prime }-{\xi }_2\right)\Delta W\left(-{\xi }_1,{\xi }_2-{\xi }_1\right)d{\xi }_{1}d{\xi }_2+\int R\left(x^{\prime }-{\xi }_2\right)M^K\left(x^{\prime }-{\xi }_1\right)\Delta W\left({\xi }_1-{\xi }_2,-{\xi }_2\right)d{\xi }_{1}d{\xi }_2$

[0116] Typically, the mask image M^K is real, thus
V(x′)=∫R(x′−ξ₂)M^K(x′−ξ₁)[ΔW(−ξ₂, ξ₁−ξ₂)+ΔW(ξ₁−ξ₂, −ξ₂)]dξ₁dξ₂   (Eq. 6)

[0117] The AFGM bilinear kernel (AFGK) can be related to the Delta TCC in the frequency domain:

[0118] $AFGK\left(k_1,k_2\right)=\int \left[\Delta W\left(-{\xi }_2,{\xi }_1-{\xi }_2\right)+\Delta W\left({\xi }_1-{\xi }_2,-{\xi }_2\right)\right]\exp \left(-i{k}_1{\xi }_1+i{k}_2{\xi }_2\right)d{\xi }_{1}d{\xi }_2=\int \Delta W\left({\xi }_1^{\prime },{\xi }_2^{\prime }\right)\exp \left(-i{k}_1\left({\xi }_2^{\prime }-{\xi }_1^{\prime }\right)-i{k}_2{\xi }_1^{\prime }\right)d{\xi }_1^{\prime }{d}_2^{\prime }+\int \Delta W\left({\xi }_1^{\prime },{\xi }_2^{\prime }\right)\exp \left(-i{k}_1\left({\xi }_1^{\prime }-{\xi }_2^{\prime }\right)-i{k}_2{\xi }_2^{\prime }\right)d{\xi }_1^{\prime }d{\xi }_2^{\prime }=\Delta \mathrm{TCC}\left(k_2-k_1,-k_1\right)+\Delta \mathrm{TCC}\left(k_1,k_1-k_2\right)$
Unlike TCCs, the Hermiticity of the AFGK is no longer guaranteed.

[0119] A practical difficulty is that if this formula is used directly, two raw TCCs appear simultaneously, which may be not feasible if the TCC is large (e.g., if each dimension of the TCC is 107 with float data type, then the total memory requirement exceeds 2G bytes). Therefore, it is desirable to make the computation “in-place.” To do so, the AFGK can be decomposed as
TCC₁(k₁,k₂)=ΔTCC(k₁,−k₂)
TCC₂(k₁,k₂)=TCC₁(k₁,k₂−k₁)=ΔTCC(k₁,k₁−k₂)
AFGK(k₁,k₂)=TCC₂(k₁,k₂)+TCC₂*(−k₁,−k₂)
where each step is in-place.

[0120] Another practical consideration is that TCCs are typically decomposed into convolution kernels, using an eigen-series expansion for computation speed and storage, because of their Hermiticity. Though AFGK is not necessarily Hermitian, one can apply Singular Value Decomposition (SVD) to it, i.e.,

[0121] $AFGK\left(k^{\prime },k^{″}\right)=\sum _{n=1}^N{\lambda }_n{\varphi }_n\left(k^{\prime }\right){\phi }_n^{*}(k^{″})$
where λ_n (n=1, . . . , N) denotes the N largest eigenvalues and φ_n(●) and φ_n* (●) denote the corresponding left and right eigenvector of the matrix, respectively. It is noted that (Eq. 6) is exact when all terms are retained in the SVD expansion, i.e., when N is equal to the rank of the AFGK matrix. However, in actual applications, it is typical to truncate the series by selecting a smaller N to increase the speed of the computation process.

[0122] Then, rewriting (Eq. 6) in frequency domain yields:

[0123] $V\left(x\right)={\int }_{{\xi }_1}{\int }_{{\xi }_2}R\left(x-{\xi }_2\right)M^K\left(x-{\xi }_1\right)\left[{\int }_{k^{\prime }}{\int }_{k^{″}}AFGK\left(k^{\prime },k^{″}\right)\exp \left(i\left(k^{\prime }{\xi }_1-k^{″}{\xi }_2\right)\right)d{k}^{\prime }d{k}^{″}\right]d{\xi }_{1}d{\xi }_2={\int }_{k^{\prime }}{\int }_{k^{″}}AFGK\left(k^{\prime },k^{″}\right)\left[{\int }_{{\xi }_2}R\left(x-{\xi }_2\right)\exp \left(-i{k}^{″}{\xi }_2\right)d{\xi }_2{\int }_{{\xi }_1}{M}^K\left(x-{\xi }_1\right)\exp \left(i{k}^{\prime }{\xi }_1\right)d{\xi }_1\right]d{k}^{\prime }d{k}^{″}={\int }_{k^{\prime }}{\int }_{k^{″}}AFGK\left(k^{\prime },k^{″}\right)\left(-k^{″}\right)K^{}\left(k^{\prime }\right)\exp \left(i\left(k^{\prime }-k^{″}\right)x\right)d{k}^{\prime }d{k}^{″}=\sum _{n=1}^N{\lambda }_n{\Phi }_n\left(x\right){\Gamma }_n^{*}\left(x\right)$
Where:
(k) and ^K(k) are the Fourier transforms of R(x) and M^K(x), respectively;
Φ_n(x)=∫φ_n(k)^K(k)exp(jkx)dk and
Γ_n(x)=∫φ_n(k)(k)exp(jkx)dk.

[0124] Note that one can first apply an OPC to the main feature, and then use this formula to generate AFGM, extract the assisting feature from AFGM, and apply another round of OPC to fix the print contour around the pre OPC edge locations.

[0125] To speed up the process, it is possible to skip the first OPC round and let M^K(x)=M^T(x) in the calculation of AFGM, in other words ignore the M^C(x) term in M^K(x)=M^T(x)+M^C(x), since M^C(x) is typically smaller compared to M^T(x).

[0126] After the AFGM (either single kernel or multi kernel AFGM) is computed, the assisting feature is extracted from this gray level map. Without any applicable constraints, one can simply place assisting features on every pixel with positive AFGM values. However, in reality, it may not be desirable to do so, since it may affect the manufacturability of the mask. For example, current technology does not allow very small assisting features.

[0127] For 1-Dimensional (1D) gauge design, there are three relevant Mask Rule Check (MRC) constraints that should be considered in placing assisting features according to embodiments of the invention: (1) there should be a minimum assisting feature width (W_min), i.e., for any assisting feature, its width should be no less than(W_min); (2) there should be a minimum spacing between assisting feature and main feature (S_main), i.e., the gap between any assisting feature and the main feature should be no less than (S_main); and (3) there should be a minimum spacing between any two neighboring assisting features (S_AF), i.e., for any two assisting features, the gap between them should be no less than (S_AF).

[0128] Next, there will be described a method of how to place assisting features to maximize the total AFGM value (thus maximize the CD sensitivity) under these MRC constraints according to embodiments of the invention. In some respects this method is similar to the SRAF rule generation based on SGM in U.S. patent application Ser. No. 11/757,805.

[0129] FIG. 5 is a diagram showing the coordinate system for generating assisting features using embodiments of AFGM according to the invention. Here both the main features 510 and assisting features (not shown) are assumed to have infinite length as compared to their width, i.e., they are all 1D patterns. The space between any two neighboring main features 510 is the specified space value between main features, which equals the difference between the pitch and main feature 510 line width. An AFGM is then generated for the main feature 510, with or without OPC.

[0130] As shown in this figure, a coordinate system is imposed on the main features 510, where the y-axis coincides with the boundary of an arbitrary main feature and the x-axis is perpendicular to the main feature 510. In the figure, at 502 x=0 and at 504 x=space. These points correspond to the boundaries of neighboring main features 510. For such one-dimensional patterns, the AFGM is also one dimensional. Thus only the AFGM value between any two neighboring main features 510 will be investigated, denoted as S(x)=AFGM (x,0), where x=[0, 1, . . . space], such as arbitrary point 506. The assisting feature extraction problem for these 1D features is then transformed into the problem of partitioning the interval [0, space] into n smaller intervals [x_1s, x_1e], [x_2s, x_2e], . . . [x_ns, x_ne], where 0≦x_1s<x_1e<x_2e<. . . <x_ns<x_ne≦space. Each interval represents an assisting feature, i.e., the i-th AF (1≦i≦n) can be described as x_is≦x≦x_ie.

[0131] Determining the optimal assisting feature placement is equivalent to maximizing the total AFGM value covered by all assisting features subject to MRC rules and possibly assisting feature printability constraints. Let S_i be the AFGM value covered by the i-th assisting feature (1≦i≦n), then the total AFGM value covered by all assisting features is

[0132] $\sum _{i=1}^n{S}_i=\sum _{i=1}^n\sum _{x=x_{\mathrm{is}}}^{x_{\mathrm{ie}}}S\left(x\right)$

[0133] There are five constraints on placing AFs in a layout: (1) minimum assisting feature width (W_min), i.e., for any i=(1, 2, . . . . n), x_ie−x_is≧W_min; (2) maximum assisting feature width (W_max), i.e., for any i=(1, 2, . . . . n), x_ie−x_is≦W_max; (For certain applications, there may be a finite constraint on the largest possible assisting feature width, for example, the assisting features should not print. If there is no such constraint, W_max can be considered equal to ∞; (3) minimum spacing between any assisting feature and any main feature (S_main), i.e., x_ls≧S_main and x_ne≦space−S_main; (4) minimum spacing between any two neighboring assisting features (S_AF), i.e., for any i=(1, 2, . . . . n), x_is−x_(i−1)e≧S_AF; and (5) For any i=(1, 2, . . . . , n), S_i≧0. (There is no need to place assisting features with negative AFGM value, even if its value is the largest possible).

[0134] Assuming the global optimal solution (partition) for [0,space] with constraints (W_min, W_max, S_main, S_AF) is Rule_opt={[x_1s, x_1e], [x_2s,x_2e], . . . . [x_ns,x_ne]}, then the i-th assisting feature (1≦i≦n) covers [x_is, x_ie]. What is more, for any i=(2, . . . , n), {[x_1s, x_1e], [x_2s,x_2e], . . . , [x_(i−1)s,x_(i−1)e]} is also the optimal partition for [0, x_is−S_AF] with the same constraints (otherwise, if there exists a better partition for [0, x_is−S_AF], then it can be combined with the i, i+1, . . . , n-th assisting feature placement in Rule_opt and land at a rule that is better than Rule_opt and still satisfies the constraints, which contradicts the optimality of Rule_opt).

[0135] Thus, the interval [0,space] is divided into smaller intervals and an algorithm is constructed based on dynamic programming. The summary of this algorithm follows, assuming space≧2S_main+W_min:

[0136]

INPUT: space, S(x) for x=[0,1,...,space], and constraints (Wmin, Wmax,
Smain, SAF)
INTERMEDIATE RESULTS:
NumAFArray[x] (x=[0,1,...,space−Smain]): an array which has a size of
space-Smain+1 and NumAFArray[x] stores the number of assisting features
of the optimal partition for [0,x]
AFAFGMArray[x] (x=[0,1,..., space−Smain]): an array which has a size of
space−Smain+1 and AFAFGMArray[x] stores the total AFGM covered by
assisting features of the optimal partition for [0,x]
AFLeftEndArray[x] (x=[0,1,..., space−Smain]): an array which has a size of
space−Smain+1 and AFLeftEndArray[x] stores the coordinate of the right
most assisting feature's left end of the optimal partition for [0,x]
(corresponds to the largest xis such that xie<x]
AFRightEndArray[x] (x=[0,1,..., space−Smain]): an array which has a size
of space−Smain+1 and AFRightEndArray[x]stores the coordinate of the
right most assisting feature's right end of the optimal partition for [0,x]
(corresponds to the largest xie such that xie<x
INITIALIZATION:
Set NumAFArray[x] and AFAFGMArray[x] to zero for all x=[0,1,...,
space−Smain]
AF COMPUTATION:
For i=Smain+Wmin to space−Smain, STEP=1 //For Constraint 3
	tempAFGMValue←AFAFGMArray[i−1]
	tempNumAF←NumAFArray[i−1]
	tempAFLeftEnd←AFLeftEndArray[i−1]
	tempAFRightEnd←AFRightEndArray[i−1]
	tempNewAFAFGM ⁢ ← ⁢ ∑ k = i - W min i ⁢ ⁢ S ⁡ ( k ) ⁢ // Candidate ⁢ ⁢ AF ' ⁢ s ⁢ ⁢ AFGM ⁢ ⁢ value
	for j=i−Wmin to max{i−Wmax, Smain}: STEP=−1
	//j: Candidate AF's left end.
	//The width of each AF is guaranteed to fall in [Wmin,Wmax]
		if(tempNewAFAFGM≧0) //For Constraint 5
			h←j−SAF
			if(h≧Smain+Wmin)
				PreviousAFGMValue←AFAFGMArray[h]
				PreviousNumAF← NumAFArray[h]
				//Optimal partition for [0, j−SAF]
			else
				PreviousAFGMValue←0
				PreviousNumAF← 0
			End
		if(tempNewAFAFGM+PreviousAFGMValue>tempAFGMValue)
				tempAFGMValue←
				tempNewAFAFGM+PreviousAFGMValue
				tempNumAF← PreviousNumAF+1
				tempAFLeftEnd←j
				tempAFRightEnd←i
			End
		End
		tempNewAFAFGM← tempNewAFAFGM+S(j−1)
	End
	AFAFGMArray[i] ←tempAFGMValue
	NumAFArray[i] ← tempNumAF
	AFLeftEndArray[i] ← tempAFLeftEnd
	AFRightEndArray[i] ← tempAFRightEnd //Update all
		intermediate results
End
OUTPUT: NumAFArray[space−SAF], AFLeftEndArray[x] (x=[0,1,...,
space−SAF]), and AFRightEndArray[x] (x=[0,1,..., space−SAF])

[0137] FIG. 6 shows an example of the 1D maximum sensitive gauges generated according to embodiments of the invention, where the patterns 602 are the periodic main features and the patterns 604 are assisting features placed based on AFGM. It should be noted that there can be, and preferably are, many more main features 602 than shown in FIG. 6.

[0138] FIG. 7 compares the CD sensitivities of line space patterns without any assisting feature 702 (in solid line) vs. those of the line space patterns with assisting features 704 (in dashed line) with respect to sigma variation. All of the main features' line width is 80 nm, while the x-axis shows the pitch in nm. The y-axis shows the AI CD variation in nm when sigma changes from 0.8 to 0.83. The assisting features improve the CD sensitivities significantly. Further, for small pitch, because of MRC, thus limited space for assisting feature placement, the CD sensitivities of the pure line space patterns are almost the same as those from gauges with assisting features. As pitch becomes considerably large (>1500 nm), then the CD sensitivities of the maximum sensitive gauges stabilize around the maximum value, which is consistent with the previous statement: it is preferable to choose a large pitch for the main feature.

[0139] The above describes how to identify the most sensitive main feature line width and how to use assisting features to strengthen the CD sensitivities from the main features, so that the resulting gauge has the most (absolute) CD sensitivity. The present inventors have noted that these methodologies are quite versatile in that, with a little change of sign, they can be used to design gauges for different purpose, such as gauges with maximum positive CD sensitivity, gauges with maximum negative CD sensitivity, and gauges with minimum CD sensitivity (i.e., most insensitive gauges).

[0140] Even though the above primarily discusses the design of 1D gauges, the AFGM can be computed for 2D main feature, such as contact, in the same way. Then it is possible to extract 2D assisting features from the 2D AFGM, subject to possible MRC constraints. For example, FIG. 8 shows an image with a contact 806 as the main feature and assisting features (black rings 802 and dots 804 in the surrounding area). The assisting features are designed to enhance the sensitivity (Contact's CD or energy) against Sigma variation. Those skilled in the art will recognize how to extend the 1D gauge design methodologies of the invention to designing 2D gauges based on the examples provided above.

[0141] It should be noted that these maximum sensitive gauge masks themselves can have all possible patterns designed using methodologies described above, covering a range of most possible optical setups and most possible parameter variations. However, for each specific lithography process, the actual SEM measurements do not need to be taken on all of them. The patterns that need to have SEM data taken can be analyzed after nominal parameters are known, e.g., NA, Sigma, film stack. Then the overall calibration flow can be as follows: print the Gauge Mask, determine the rough number of SEM measurements to be taken, then use a “maximum-sensitive-gauge-selection” algorithm to analyze the nominal model and select the set of gauge patterns to be measured by SEM, then take SEM measurements. The algorithm in this “selection” process can be based on methodologies described above. These selected patterns will represent the optimal set of SEM measurements.

[0142] In another embodiment, gauges are designed with different main feature line widths. This approach is called Spatial-Difference-Delta-TCC (SDDTCC) design. In this design approach, the method places two main features in each gauge and use the difference between the two CDs to calibrate the parameters. An advantage of this design is that it eliminates any measurement error in bias as long as these two main features are reasonably close. The two main feature line widths (denoted as LW1 and LW2) and the center-to-center distance between the main features (denoted as s) suffice to describe the main features. Again, for each parameter to be calibrated and each combination of (LW1, LW2, s), assisting features are added to improve the sensitivity of the delta CD (difference between the CDs of the two main features). Then the process loops over all possible (LW1, LW2, s) to identify the optimal combination. Further, one can apply the same approximations used in the previous embodiment to speed up the design. In particular, the larger the value of s, the better effect from assisting features. This is due to larger space to accommodate more assisting features and the main feature can also be viewed as assisting feature of each other. In this case, since s is quite large compared to LW1 and LW2, one can break down the SDDTCC design into two independent maximum sensitive gauge designs: the left half of the SDDTCC gauge is actually a gauge with maximum positive CD sensitivity; while the right half of the SDDTCC gauge is a gauge with maximum negative CD sensitivity.

[0143] It should be further noted that in some applications where intensity is known, e.g., certain mask model calibration, one can design intensity based gauges which maximize the intensity sensitivity against parameter variation. Here, the parameter variation can still be represented by Delta TCC. Similar to CD based calibration, increasing the intensity sensitivity also leads to better robustness against measurement error or less measurements. It is advantageous to design a mask which maximizes the intensity sensitivity at certain point. Without loss of generality, assume the point of interest is the origin point, i.e, x=0 in (Eq. 3). For mask image M(x), there should also be a constraint on the intensity norm to remove the scaling effect, for simplicity, assume that the norm satisfies

[0144] $M\left(x\right)=\sum _x{M\left(x\right)}^2\le 1.$

[0145] As is known, the eigenvectors φ_n (k) in (Eq. 3), which are resulted from eigenseries expansion of the Hermitian Delta TCC, are orthogonormal, i.e.,

[0146] ${\varphi }_n\left(k\right)=\sum _k{{\varphi }_n\left(k\right)}^2=1\mathrm{for}\mathrm{any}n$ ${\varphi }_{n_1}\left(k\right)·{\varphi }_{n_2}\left(k\right)=\sum _k{\varphi }_{n_1}\left(k\right){\varphi }_{n_2}^{*}\left(k\right)=0\mathrm{for}\mathrm{any}{n}_1\ne n_2$
Or, in other words:

[0147] $\sum _k{\varphi }_{n_1}\left(k\right){\varphi }_{n_2}^{*}\left(k\right)=\delta \left(n_1-n_2\right).$

[0148] Then a complete orthogonal coordinate system can be built in the space domain based on the eigenvectors, and the Fourier transform of any grayscale mask image can be represented as linear combination of all the eigenvectors, i.e.,

[0149] $\left(k\right)=\sum _n{c}_n{\varphi }_n\left(k\right).$

[0150] From the mask image intensity constraint and Parseval's Theorem,

[0151] $\sum _i{c_i}^2=\left(k\right)=M\left(x\right)\le 1.$

[0152] Without loss of generality, assume that the eigenvalues are ordered such that
|λ₁|≧|λ₂|≧|λ₃|≧. . .

[0153] Thus, the first eigenvector φ₁ (k) corresponds to the most significant eigenvalues of the Delta TCC. Then from (Eq. 3), the change (Delta I) in aerial image intensity at x=0 resulted from the mask image (k) is

[0154] $\begin{array}{c}\Delta I\left(0\right)=\sum _n{\lambda }_n{\sum _k{\varphi }_n\left(k\right)}^2\\ =\sum _n{\lambda }_n{\sum _k{\varphi }_n\left(k\right)\sum _m{c}_m{\varphi }_m^{*}\left(k\right)}^2\\ =\sum _n{\lambda }_n{\sum _m\sum _k{c}_m{\varphi }_n\left(k\right){\varphi }_m^{*}\left(k\right)}^2\\ =\sum _n{\lambda }_n{\sum _m{c}_m\delta \left(n-m\right)}^2\\ =\sum _n{\lambda }_n{c_n}^2\end{array}$

[0155] Since there exists the constraint:

[0156] $\sum _n{c_n}^2\le 1$
and the assumption that
|λ₁|≧|λ₂|≧|λ₃|≧. . . ,
then.

[0157] $\Delta I\left(0\right)=\sum _n{\lambda }_n{c_n}^2\le \sum _n{\lambda }_1{c_n}^2={\lambda }_1\sum _n{c_n}^2\le {\lambda }_1.$
Similarly,
ΔI(0)≧−|λ₁|.

[0158] Therefore, if a mask is designed which is the inverse Fourier transform of delta TCC's most significant kernel φ₁(k), i.e., c_n=δ(n−1), then it has the maximum (absolute) aerial image intensity sensitivity. Such mask pattern “resonates” with delta TCC.

Multivariate Calibration Pattern Design

[0159] The above describes methods to design gauges that are extremely sensitive (in CD, CD difference, or aerial image intensity, for example) to one specific optical parameter, which can be considered as a specific direction in the model parametric space.

[0160] In the above methodology, the process uses TCC's to compute the predicted CD and perform gauge design, and the parameter perturbations are expressed in Delta TCC's. So for each given parameter P, the methodology forms Delta TCC's by changing the single parameter P to P+ΔP while keeping the other parameters constant.

[0161] However, as mentioned earlier, for some parameters, the present inventors recognize that gauges designed using the methodologies above may exhibit quite similar behavior for several different parameter variations; therefore, it is extremely difficult to distinguish these parameters' variations using such a gauge set, and so the gauge set may be incomplete in calibrating these parameters separately (i.e. degeneracy among parameters). Using one simple example, suppose there are two parameters P₀ and P₁ and only one gauge, and the relationship between the predicted gauge CD and parameters P₀ and P₁ is CD=k+g(P₀+hP₁), where k, g, and h are arbitrary constants independent of P₀ and P₁ that will become more apparent from below. In a typical calibration, the process is to first measure the actual wafer CD, {tilde over (C)}{tilde over (D)}.

[0162] Then the process looks for the values for all parameters including P₀ and P₁ such that they minimize the error between the predicted CD and measured wafer CD, as defined in Eq. 4 above. For this example, any values of parameters P₀ and P₁ such that
P₀+hP₁=({tilde over (C)}{tilde over (D)}−k)/g
will lead to the minimum calibration error (i.e. predicted minus measured CD equals zero) for arbitrary values of g, h and k, where g is a scaling factor between CD and parameters. In other words, there is no way to determine unique values for P₀ and P₁; thus, the model built from this calibration process may only work well for this specific gauge, but it does not capture the behavior for all possible patterns. Therefore, it becomes apparent that it is critical to design gauges which maximize differential response to different parameter variations to achieve sufficient pattern coverage.

[0163] In general, the embodiments provided below provide several different methodologies to design gauges that respond vastly differently to different parameter variations in order to achieve complete pattern coverage. The focus of the discussion below is on CD sensitivities; however, those skilled in the art will recognize that all the embodiments described herein are versatile so that they can be applied to other measurements, such as CD difference sensitivities, or aerial image intensity sensitivities. The present inventors recognize that, in order to achieve best complete pattern coverage, the gauges should maximally and distinctly excite all the physics and chemistry parameters to be calibrated, making sure that the wafer data for the gauges can drive the model calibration to the optimal parameter values.

[0164] In a first embodiment, a set of gauges is created to maximize the pair-wise differential CD sensitivity to different parameter variation(s). The resulting gauges are referred to herein as Max-Delta-Parameter-Difference gauges (MDPDG). In particular, between each pair of the parameters, the invention creates pair-wise gauges that maximize the CD sensitivity with respect to the delta-parameter-difference. If there are N parameters, there are C_N²=N(N+1)/2 pairs or N(N+1) gauges. For each pair, this embodiment specifies that there are at least two gauges, represented as (ΔP₁, −ΔP₂) and (ΔP₁, ΔP₂), where ΔP₁ and ΔP₂ denote the perturbed amounts for parameters P₁ and P₂, respectively. More specifically, for (ΔP₁, −ΔP₂), this embodiment of the invention designs the gauge that maximizes the CD sensitivity when parameter P₁ changes from its nominal value P₁₀ to P₁₀+ΔP₁ and P₂ changes from its nominal value P₂₀ to P₂₀−ΔP₂, simultaneously.

[0165] In this embodiment, the gauge design method is preferably the same as that for gauges maximizing the CD sensitivity with respect to a single parameter as taught in the univariate calibration pattern design section. More particularly, the design of a gauge maximizing the CD sensitivity with respect to P₁ change as performed in the embodiments above relating to univariate calibration pattern design can be viewed as the design of MDPDG of (ΔP₁, 0)) according to the present embodiment (ΔP₁=P₁+ΔP, where ΔP is a small perturbation amount). The present embodiment extends the approach of the above embodiments to find a pair of MDPDG gauges (ΔP₁, ΔP₂) and (ΔP₁, −ΔP₂) that both maximize the CD sensitivity when parameter P₁ changes from its nominal value P₁₀ to P₁₀±ΔP₁ and P₂ changes from its nominal value P₂₀ to P₂₀±ΔP₂, simultaneously and respectively (ΔP₁ and ΔP₂ are small perturbation amounts, and may or may not be equal to each other).

[0166] For example, rather than designing a single gauge per parameter using a Delta TCC that captures only one parameter perturbation as in the embodiments above, the present embodiment designs a pair of gauges for each pair of parameters. The present embodiment therefore uses a respective pair of Delta TCCs that both capture two parameter perturbations as set forth above, namely one Delta TCC based on (ΔP₁, ΔP₂) and one Delta TCC based on (ΔP₁, −ΔP₂). These Delta TCCs are then used in the design gauge methodologies described in the above embodiments in place of the single-parameter Delta TCCs described therein.

[0167] This methodology directly maximizes the differential behavior to different parameter variations. Thus, unlike the gauges designed to maximize single parameters, MDPDGs should guarantee complete pattern coverage, given a sufficient number of gauges.

[0168] If the number of parameters to be calibrated is small, then the number of required MDPDGs is not large. For example, if there are N=5 optical parameters, the above method will require the design of N*(N−1)=20 gauges, which is easy to handle. However, if the number of parameters increases, the number of resulting gauges increases very fast and will be on the order of O(N²). In that case, in an embodiment manual intervention or artificial intelligence algorithms (e.g., Supporting Vector Machines) are used to separate the parameters into multiple less-correlated groups. Then, where inter-group degeneracy is weak, pair-wise MDPDGs are designed within each group.

[0169] Using the MDPDG design according to this embodiment for multivariate calibration pattern design the minimum required number of gauges is reduced to N*(N−1) gauges. In most practical situations this is already much better than “full-chip” pattern sets or other pattern sets that are representative of the pattern variety in the product. Additionally the method is independent of the variety in the actual products which are produced after a lithographic process is initiated based on the calibrated model. Theoretically speaking, assuming that there exists good linearity between CD sensitivity and parameter variation and no degeneracy in the transformation from model parametric space to CD-sensitivity space (this assumption is generally true from simulations and experiments. Also note that the degeneracy here is not the degeneracy from model parametric space to a very limited gauge set discussed previously, but the transformation from model parametric space to CD-sensitivity space, or all possible gauges patterns), and that wafer measurement error can be ignored, then there should exist N gauges that cover the N parameter variations sufficiently. However, MDPDG design requires at least N*(N−1) gauges, which suggests certain redundancy in the gauge design, and further in the amount of measurements required.

Multivariate Orthogonal Calibration Pattern Design

[0170] Additional embodiments will now be described according to a methodology to design gauges with sufficient pattern coverage. The number of resulting gauges in this approach is in the order of O(N). When there are only an order of O(N) gauges needed, less computation run time is required. Also, fewer gauges have to be measured physically, which reduces metrology costs and time constraints.

[0171] The inventors have noted that some parameters have strong degeneracy due to their close physical relationship, such as parameters related to focus and aerial image location. So many patterns exhibit very similar response to changes in these two parameters. Further, CD sensitivities for all parameters are typically correlated with aerial image intensity and/or resist image intensity slope at the contour, while the design goal of dose-sensitive gauges is to decrease the aerial image/resist image intensity slope at the contour. So dose-sensitive gauges are not only sensitive to dose variation, they are also quite sensitive to many other parameters and thus selected as gauges specific for these parameters. Therefore, it is extremely difficult to differentiate the behavior of one parameter from the other, and thus it is difficult to determine the precise values for these parameters in calibration. As a result, it would be advantageous to break such degeneracy in gauges' response to different parameter variations.

[0172] Accordingly, in the following design embodiments, one aspect is to determine N orthogonal directions in the N-dimensional model parametric space. Such gauges are thus referred to herein as orthogonal gauges. The first such embodiment depends on a large set of test calibration patterns representative of all aspects of product patterns which are to be produced based on the calibrated model. The second such embodiment is based on the change in TCC's (Delta TCC). Note that TCC's can be computed without knowledge of the pattern to be imaged. Hence, this leads to independence of the second embodiment from the requirement to follow a large set of test calibration patterns.

[0173] Before describing the details of gauge design according to this embodiment, a brief review of techniques in pattern coverage measurement is provided, which can help in understanding the present invention gauge design procedure. Pattern coverage can be measured by a certain cost function, for example the cost function described in, for example paragraph [0051] of U.S. Patent App. No. 61/140,812 (P-3365), commonly owned by the present assignee and incorporated herein by reference in its entirety, where the pattern coverage of a gauge set is measured by the remaining CD prediction uncertainty after applying this gauge set in calibration. In this embodiment, the invention preferably uses the same cost function definition in the gauge set design. The details of this cost function as well as those of the gauge design, to the extent helpful to teach those skilled in the art about the principles of this embodiment of the present invention, are described in more detail below.

[0174] Assuming there are a total number of G (G>>N) patterns, denote the sensitivity of the g-th gauge (g=1, . . . , G) to the model parameter P_n (n=1. . . N) by

[0175] $K_{\mathrm{ng}}=\frac{\partial {\mathrm{CD}}_g\left(P_n\right)}{\partial P_n}{❘}_{P_n=P_{n0}},$
where P_n0 is the nominal value for parameter P_n. Imagine a plot in the N-dimensional sensitivity space in which each point is a prospective test pattern, with coordinate K_ng for the g-th point on the n-th axis. Each axis corresponds to one parameter, and each point corresponds to the sensitivity of a particular pattern with respect to all N parameter changes. Similar to the approach for pattern selection in U.S. Patent App. No. 61/140,812 (P-3365), this gauge design can also be described as a maximization of an enclosed volume of designed gauges in sensitivity space. For a given target number L of gauges, this approach attempts to design the L patterns that results in the maximum enclosed volume in the sensitivity space. This approach prevents degeneracy because degeneracy will lead to small enclosed volumes. Below, the concrete procedures to find the patterns that achieve this maximization are explained.

[0176] A next step in this approach is to perform Principal Component Analysis (PCA, also known as Karhunen-Loève Transform (KLT), Hotelling Transform, or proper orthogonal decomposition (POD)) on the sensitivity matrix K_ng. In this mathematically rigorous approach, the principal components of the sensitivity matrix are identified. These are the directions in the N-dimensional model parametric space that capture the most variations in the sensitivities, and as a result of the PCA, these directions are orthogonal (perpendicular) to each other. This way, the maximum sensitivity requirement and the least degeneracy requirement can be satisfied at the same time.

[0177] An outline of part of this embodiment is as follows, whereas further details will be provided below. A next step in this approach is to find a Confidence region of model parameters. Then a cost function related to the confidence region is designed as the optimization target. One cost function that can be used is the trace of the inverse of the covariance matrix in the normal equation for the least-squares optimization (see descriptions below for the mathematical formulation). The present inventors recognize that minimization of this confidence region is equivalent to maximization of enclosed volume in sensitivity space. But some refinement will be needed to handle the interaction between the model parameters, since they have different units, and it is generally undesirable that the choice of units should result in some parameters being favored over the others in the gauge design.

[0178] After designing the cost function, a next step in this approach is to find a Confidence region of predictions of test pattern CD's. In order to deal with the problem of handling model parameters with different units, the algorithm calculates the model parameter uncertainties and simulates their impact on CD prediction uncertainty. The problem can be stated as follows: Design the set of L gauges that minimizes the sum-of-squared-errors uncertainty in the CD prediction for a target set of test patterns—which can be the set of the G training patterns themselves or another pattern set from an actual product layout—as estimated from the confidence region of the model parameters.

[0179] Certain aspects of a Mathematical Formulation of cost functions based on the confidence regions mentioned above according to an embodiment of the invention for designing gauges for lithographic model calibration will now be described in more detail.

[0180] The index p is further used to correspond to the actual physics and chemistry in the lithographical process, while m is used to correspond to nominal values of the parameters in the model description of the process. Furthermore, let P_p and P_m represent the physical parameter values and the nominal model parameter values, and use a common parameter index j to signify the conceptual correspondence between actual parameter values P_pj and nominal model parameter values P_mj.

[0181] Let CD denote the values for candidate test patterns from which the pattern coverage of the designed gauges for the actual model calibration is to be evaluated. More particularly, let CD _p denote the set of measurement values for these patterns and CD_m the set of nominal model predicted values for the same patterns. The task of model calibration will then be the determination of P_pj, given P_mj, CD_p, and CD_m as input.

[0182] In the context of model calibration (in fact, for many other applications, such as scanner matching, this approximation is still valid), an assumption can be made that the variations of physical and model parameters are within a small range during calibration. Accordingly, it can be assumed that a linear formulation is appropriate for the small value difference under consideration. Therefore, the simulated pattern sensitivity for a given pattern i to a parameter j can be described as

[0183] $K_{m,n,g}=\frac{\partial {\mathrm{CD}}_{m,g}}{\partial P_{m,n}}$

[0184] The pattern coverage measurement thus the gauge design will be based on this sensitivity matrix K_m. Assume that for the purpose of gauge design, the sensitivity matrix calculated from a reasonably accurate nominal model description (m) will be sufficient. as the reasonably accurate nominal model description brings the sensitivities in the linear range.

[0185] Starting from the linearized model calibration equation

[0186] ${\mathrm{CD}}_{p,g}-{\mathrm{CD}}_{m,g}=\sum _n{K}_{m,n,g}\left(P_{p,n}-P_{m,n}\right),$

[0187] This is usually an over-determined equation (i.e. the number of gauges is larger than the number of parameters) such that actual parameters P_pj can be fitted in the least-squares sense. Considering the fact that wafer measurements CD_p,i have inherent uncertainties (denoted σ_g), the equation may be rewritten in the matrix form Ax=b and solved using the normal equation formulation (A^TA)x=A^Tb where
A_ng=K_m,n,g/σ_g, x_n=P_p,n−P_m,n, b_g=(CD_p,g−CD_m,g)/σ_g

[0188] Matrix C=(A^TA)⁻¹ is closely related to the probable or standard uncertainties of the estimated parameters x, in that it is the covariance matrix of x (see, e.g., a book titled “Numerical Recipes in C: The Art of Scientific Computing”, by W. H. Press, Cambridge, U.K.: Cambridge Univ. Press, (1992)):
σ²(x_j)=C_jj, Cov(x_j,x_k)=C_jk

[0189] A confidence region for the parameters can be derived from the covariance matrix, for example as described in the book entitled “Numerical Recipes in C, The Art of Scientific Computing,” 2d ed., W. H. Press, S. A. Teukolsky et al., Chap. 15, Cambrige University Press, 1992. If it is chosen to minimize the ‘size’ of this confidence region (defined as the sum of squared principal axis dimensions of the confidence ellipsoid, or the sum of squared uncertainty level of individual parameters) at a certain confidence level, this cost function will in fact be equal to the trace of matrix C.

[0190] The cost function definition above contains an undesirable arbitrariness due to the artificial choice of parameter units. One alternative to avoid this is to project the estimated covariance matrix of the parameters to CD prediction errors over a large pattern space, and use the sum of squared error for all the pattern CD predictions as the cost function. In practice, a full-chip pattern set or a pre-determined set of representative patterns (i.e. a training pattern set) may be used that can emulate the full-chip pattern behavior. The point is that the pattern set (or set of patterns) should cover all possibly desired calibration patterns (and a pattern corresponding to nominal parameter values). If the ‘full-chip’ pattern set has the full sensitivity matrix of A , this updated cost function will take the form of

[0191] $\sum _{\mathrm{CD}\in \mathrm{full}-\mathrm{chip}\mathrm{pattern}\mathrm{set}}{\sigma }^2\left(CD\right)=\mathrm{Tr}\left[AC{A}^T\right]=\mathrm{Tr}\lfloor {A\left(A^{T}A\right)}^{-1}{A}^T\rfloor$
and the minimization of this function will be used to drive gauge design in a final implementation.

[0192] To achieve the minimization of the cost function described above, embodiments of the invention perform principal component analysis on the full sensitivity matrix A, and design gauges that maximize the CD sensitivity for the directions in the N-dimension sensitivity space along which most variations in A are distributed. The resulting gauges whose sensitivities are most aligned with these principal directions will be the most efficient ones in reducing the cost function.

[0193] Having provided the above details regarding the cost function to be used, a calibration gauge design methodology according to this embodiment of the invention will now be further described in connection with the flowchart provided in FIG. 9.

[0194] In an embodiment (FIG. 9), initial estimates of N model parameters are created (S902) by using, for example, any combination of machine setting data, on-tool metrology data and tool design data. As described above, let p describe the actual physics and chemistry in the lithographical process, while m is a nominal model description of the process. As further mentioned above, let P_p and P_m represent the physical parameters and the nominal model parameters, respectively, and use a common parameter index n to signify the conceptual correspondence between actual parameters P_pn and nominal model parameters P_mn.

[0195] In step S904, a large collection of patterns (e.g., ‘full-chip’ pattern set) that is representative of the pattern variety in the product on which the model is to be used is determined. For example, users can pool all the patterns that appear in a desired chip design into this collection.

[0196] In step S906, based on the initial model estimate, the sensitivities, K_ng, for all patterns g to each of the model parameters n=1, . . . , N that will be calibrated during the calibration process are calculated using the equation for K_ng set forth earlier in this application, wherein the predicted CD's in that equation are calculated using aerial image simulation techniques well known to those skilled in the art, and Delta TCCs corresponding to each parameter. In one example, each test pattern (g=1 . . . G, where G is the total number of test patterns) corresponds to a point, or equivalently, a vector, in the N-dimensional sensitivity space (as there are N model parameters). All these points can be viewed as cumulatively forming a sensitivity “cloud” for the full-chip pattern set.

[0197] In step S908, PCA is performed on the points in the sensitivity ‘cloud’ for the full-chip pattern set, using the matrix A=K/σ More particularly, the contribution to sensitivity variation from each of the principal components is calculated. The PCA of A can be determined using the eigen-decomposition of matrix A^TA, which yields the diagonal eigenvalues matrix D and the principal direction vectors (eigenvectors) V, i.e.,
A^TA=VDV^T, W^T=I
and the diagonal matrix is composed of the singular values, i.e.,
D=diag[λ₁²,λ₂², . . . , λ_N²],
where / is the identity matrix. Without loss of generality, it is assumed that |λ₁|≧|λ₂|≧ . . . |λ_N|.

[0198] In step S910, the PCA result of A^TA as expressed with the eigenvalues and eigenvectors as described above is viewed as an ellipsoidal surface (ES) embedded in the N-dimensional model parametric space.

[0199] In step S912, the principal direction corresponding to the largest singular value is chosen (λ₁) and this direction is defined as V_l, which is a direction in the N-dimensional model parametric space. Note that the eigenvector V_k can be written as V_k=[V_k1, V_k2, . . . , V_kN]^T.

[0200] In step S914, the index k, corresponding to the directions in the model parametric space to be calibrated, is initialized with 1. Moreover, the first current working sub-space (WSS) is created as the full N-dimensional model parametric space.

[0201] In step S916, one or more gauges are designed that maximize the sensitivity to the model parametric direction V_k. In this step, it is assumed that the CD sensitivity

[0202] $K_{m,n,g}=\frac{\partial {\mathrm{CD}}_{m,g}}{\partial P_{m,n}}$
(i.e., the g-th pattern's CD sensitivity against the n-th model parameter's variation) is computed using a finite difference method. For example, the CD change (denoted as ΔCD_m,n,g) of the g-th pattern when a small perturbation (denoted as ΔP_m,n) is added to the n-th model parameter is computed. Then K_m,n,g is computed as

[0203] $K_{m,n,g}=\frac{\Delta {\mathrm{CD}}_{m,n,g}}{\Delta P_{m,n}}.$

[0204] In order to design gauge(s) maximizing the sensitivity for direction V_k, first the delta TCC is computed for delta parameters (V_k1ΔP, V_k2ΔP, . . . , V_kNΔP). More particularly, this involves calculating the element-by-element difference between the perturbed TCC matrix when the n-th model parameter is perturbed by V_kjΔP for all j=1, 2, . . . , N and the nominal TCC matrix, where ΔP is a small perturbation value. For example, ΔP is set to a value such that CD sensitivity is around 10 nm. However, for small parameter variations (ΔP) the relationship with CD is linear. Then the methodology described in the above univariate parameter-sensitive calibration pattern design embodiments is used to design the gauge(s) maximizing the CD sensitivity for this delta TCC instead of the univariate (single-parameter) Delta TCC as described earlier. Note that, for each direction V_k, a few gauges can be designed with CD sensitivities slightly smaller than that of the most sensitive gauge. However, each of these different gauges will have a slightly different model parametric direction, which leads to more gauges and thus higher robustness against errors such as measurement errors. In the calibration pattern sets according to the embodiment there is optimal robustness against errors such as measurement errors.

[0205] In step S918, the actual sensitivity(ies) of the gauge(s) obtained in step S916, denoted as W_k are computed. This is done, for example, using the same equations as above, but with the actual gauge pattern(s) designed in S916. It should be noted that this sensitivity vector typically does not align with V_k perfectly.

[0206] In step S920, the current working sub-space WSS is reduced to a lower dimension by removing the sensitivity direction of W_k. In other words, the WSS is reduced to a subspace that is orthogonal to W_k.

[0207] In step S922, k is increased by 1. A new direction V_k is defined, which is the longest vector in the intersection between the ellipsoidal surface ES and the current working sub-space WSS.

[0208] In step S924, steps S916 through S922 are repeated for the rest of the pattern set until the dimensionality of the parametric space is fully exhausted.

[0209] By this point, the process has achieved coverage for every parameter, i.e., no parameter will be left unconstrained in the model calibration stage. Within the metrology time limitation, it would be desirable to include more gauges to enhance the coverage, e.g., improve the robustness against wafer measurement error. The total number of gauges can be increased by increasing the number of gauges designed in each step S916, thus increasing the subsequent directions. But the minimum number of gauges is N, i.e., the number of model parameters to be calibrated, which is consistent with the (theoretical) analysis given earlier.

[0210] The aforementioned methodology achieves the goal of designing O(N) gauges with complete pattern coverage. However, it may have certain drawback for some applications. In particular, PCA, and thus the parametric directions, rely on the identification of a large training pattern set, which may be full-chip patterns or a large collection of representative patterns as mentioned above. This may be difficult in real applications: first, it may be hard to determine the truly representative training patterns. Even if a full-chip pattern set is used, it may still become unmatched if the gauges are used for the calibration of the same lithography condition for a different design layout. Second, it may be expensive to compute and process the sensitivity matrix for a large pattern set. For these reasons, it is desired to devise a methodology which purely depends on the model and does not depend on the choice of initial pattern set.

[0211] Accordingly, another embodiment of a methodology for designing calibration gauges, based on identifying parametric directions from Delta TCCs, will now be described in connection with the flowchart in FIG. 10. In general, in the previous embodiment, optimization was based on the square root of the sum of square of CD sensitivities. In this embodiment, optimization is instead based on the Hilbert-Schmidt norm of the Delta TCC matrix, i.e., the square root of the sum of square of TCC matrix elements.

[0212] As in the previous embodiment, as shown in FIG. 10, in step S1002, initial estimates of N model parameters are created by using, for example, any combination of machine setting data, on-tool metrology data and tool design data. As described above, let p describe the actual physics and chemistry in the lithographical process, while m is a nominal model description of the process. As further described earlier, let P_p and P_m represent the physical parameters and the nominal model parameters, and use a common parameter index n to signify the conceptual correspondence between actual parameters P_pn and and nominal model parameters P_mn.

[0213] In step S1004, this embodiment first calculates the Delta TCCs for all parameter perturbations. This is done, for example, by an element by element difference between the nominal and perturbed TCC matrices. Denote Delta TCC_j as the Delta TCC resulting from a small perturbation (ΔP_n) on parameter P_n.

[0214] Next in step S1006, each Delta TCC matrix is rearranged into a vector. Howe this is done does not affect the outcome, however, all the Delta TCC matrices should be arranged in the same order. So the number of vectors is N, i.e., the number of model parameters. The Gram-Schmidt transform is then applied to these vectors, with the inner product defined to be consistent with the Hilbert-Schmidt norm (sum of product between corresponding matrix elements of the two Delta TCCs). In particular, the inner product between two Delta TCCs: ΔTCC₁=ΔTCC₁(k₁,k₂) and ΔTCC₂=ΔTCC₂(k₁,k₂) is:

[0215] $\langle \Delta TC{C}_1,\Delta TC{C}_2\rangle =\sum _{k_1,k_2}\Delta TC{C}_1^{*}\left(k_1,k_2\right)\Delta TC{C}_2\left(k_1,k_2\right),$
where * indicates complex conjugate.

[0216] The Hilbert-Schmidt norm of a vector can also be viewed as the square root of the inner product of the vector and itself, i.e., ∥ΔTCC∥=√{square root over (ΔTCC,ΔTCC)}. Further, the projection operator is defined by

[0217] ${\mathrm{proj}}_{\Delta {\mathrm{TCC}}_1}\Delta TC{C}_2=\frac{\langle \Delta TC{C}_1,\Delta TC{C}_2\rangle }{\langle \Delta TC{C}_1,\Delta TC{C}_1\rangle }\Delta TC{C}_1.$

[0218] In Gram-Schmidt transform, a new set of N orthogonal vectors (Z₁, Z₂, . . . , Z_N), is then generated. Each new vector has the same number of elements as Delta TCCs, thus the definitions of inner product, Hilbert-Schmidt norm, and projection operator for Delta TCCs also apply to the new vectors and Delta TCCs. Then (Z₁, Z₂, . . . , Z_N) are computed as:

[0219] $Z_1=\Delta TC{C}_1$ $Z_2=\Delta TC{C}_2-{\mathrm{proj}}_{z_1}\Delta TC{C}_2$ $⋮$ $Z_N=\Delta TC{C}_N-\sum _{n=1}^{N-1}{\mathrm{proj}}_{Z_n}\Delta TC{C}_N$

[0220] Finally, the Gram-Schmidt transform of Delta TCCs is defined as another set of N orthogonal vectors (Y₁,Y₂, . . . , Y_N) which are the normalization of (Z₁, Z₂, . . . , Z_N), i.e., Y_n=Z_n/∥Z_n∥, for any n=1, 2, . . , N. Note that since each of Z_n is a linear combination of Delta TCCs as a result of the orthogonalization process, each Y_n, which is Z_n divided by its length, is also a linear combination of Delta TCCs, i.e.,

[0221] $Y_n=\sum _{k=1}^N{h}_{\mathrm{nk}}\Delta TC{C}_k,$
where h_nk are constant coefficients which only depend on Delta TCCs. Then the parametric space is transformed using exactly the same coefficients. For example, a set of orthogonalized parametric vectors q_n, (n=1, . . . , N) is generated where q_n=(h_n1ΔP₁,h_n2ΔP₂, . . . , h_nNΔP_N). Each q_n vector represent a direction in the N-dimensional parametric space which leads to a Delta TCC with norm 1, and all the resulting Delta TCCs are orthogonal with each other, under the assumption that all Delta TCCs are linear with respect to small parametric variations.

[0222] In step S1008, this embodiment next defines a hypersphere HS in the q space with unit radius.

[0223] In step S1010, k is initialized to 1, and the current working sub-space WSS is initialized as the full N-dimensional model parametric q space, and V₁ to some arbitrary vector of unit length in the q space. The choice is arbitrary, now there is a hypersphere instead of an ellipsoidal surface.

[0224] Next in step S1012, the Delta TCC corresponding to V_k is computed. Similar to the embodiment described above in connection with step S916, gauge(s) are then designed for maximizing CD sensitivity for this Delta TCC, using the techniques of the above univariate parameter-sensitive gauge design embodiments. Several different gauges can be designed using this Delta TCC, further similar to the embodiment described above. Moreover, different gauges may be designed using a different V₁ in step S1010.

[0225] Further similar to the embodiment above with step S918, the actual sensitivity of the gauge is calculated in step S1014. In this case, the new gauge pattern must be projected onto the q space. The sensitivity vector, W_k, defined in the q space is then invoked. It may have components along directions other than V_k.

[0226] Similar to the embodiment above with step S920, in step S1016 the current working sub-space WSS is reduced to a lower dimension by removing the sensitivity direction of W_k. In other words, WSS is reduced to its subspace that is orthogonal to W_k.

[0227] In step S1018, k is increased by 1. Moreover, a new direction V_k is defined, which is an arbitrary vector in the intersection between the hypersphere HS and the current working sub-space WSS (which is a hypersphere of lower dimension).

[0228] In step S1020, steps S1012 to S1018 are repeated until the dimensionality of the q space is fully exhausted.

[0229] In an embodiment wherein more gauge patterns are allowed, the above steps are wrapped in a bigger loop. In each iteration of the bigger loop, there can be a different way for choosing a vector from a hypersphere. Of course, several gauges can also be generated for each direction, and then the number of subsequent directions can be increased. Again, the minimum number of gauges is N, i.e., the number of model parameters to be calibrated, which is consistent with the theoretical analysis.

[0230] There are many similarities between the embodiments described in FIGS. 9 and 10. However, there are two major differences:

[0231] 1. Note that the matrix g is not necessarily orthonormal. This is a major conceptual difference from the PCA-based procedure.

[0232] 2. The other difference from the procedure in FIG. 9 is the choice of the metric, the procedure in FIG. 10 is based on Delta TCCs rather than CD sensitivities from a chosen pattern set. This difference is independent of the previous difference, so one could also consider for example using the Gram-Schmidt procedure for CD sensitivities, or using PCA for Delta TCCs. Since the metric is based on Delta TCCs in the new procedure, thus it removes the pattern dependence. Moreover, the metric can also be based on aerial images or resist images, or 1D Delta TCCs as described in U.S. patent application No. 61/113,024 (P-3314), incorporated by reference herein in its entirety. Here, a 1D Delta TCC is defined as the 1D components of a Delta TCC.

[0233] In particular, for a 2D mask image M(k_x, k_y), the complete delta aerial image can be expressed as

[0234] $\Delta I\left(x,y\right)=\sum _{\left(\mathrm{kx}1,\mathrm{ky}1,\mathrm{kx}2,\mathrm{ky}2\right)}\Delta TCC\left(k_{x1},k_{y1},k_{x2},k_{y2}\right)M(’k_{x1},k_{y1})M^{*}\left(k_{x2},k_{y2}\right)\exp \left(-j\left(k_{x1}-k_{x2}\right)x-j\left(k_{x2}-k_{y2}\right)y\right)$