Photolithographic process simulation including efficient result computation for multiple process variation values

[0001] This patent application is a continuation of U.S. Ser. No. 11/331,223, filed Jan. 11, 2006, now abandoned.

[0002] This patent specification relates to computer simulation of photolithographic processes as may be employed, for example, in sub-wavelength integrated circuit fabrication environments.

[0003] Computer simulation has become an indispensable tool in a wide variety of technical endeavors ranging from the design of aircraft, automobiles, and communications networks to the analysis of biological systems, socioeconomic trends, and plate tectonics. In the field of integrated circuit fabrication, computer simulation has become increasingly important as circuit line widths continue to shrink well below the wavelength of light. In particular, the optical projection of circuit patterns onto semiconductor wafers, during a process known as photolithography, becomes increasingly complicated to predict as pattern sizes shrink well below the wavelength of the light that is used to project the pattern. Historically, when circuit line widths were larger than the light wavelength, a desired circuit pattern was directly written to an optical mask, the mask was illuminated and projected toward the wafer, the circuit pattern was faithfully recorded in a layer of photoresist on the wafer, and, after chemical processing, the circuit pattern was faithfully realized in physical form on the wafer. However, for sub-wavelength circuit line widths, it becomes necessary to “correct” or pre-compensate the mask pattern in order for the desired circuit pattern to be properly recorded into the photoresist layer and/or for the proper physical form of the circuit pattern to appear on the wafer after chemical processing. Unfortunately, the required “corrections” or pre-compensations are themselves difficult to refine and, although there are some basic pre-compensation rules that a human designer can start with, the pre-compensation process is usually iterative (i.e., trial and error) and pattern-specific to the particular desired circuit pattern.

[0004] Because human refinement and physical trial-and-error quickly become prohibitively expensive, optical proximity correction (OPC) software tools have been developed that automate the process of pre-compensating a desired circuit pattern before it is physically written onto a mask. Starting with the known, desired circuit pattern, an initial mask design is generated using a set of pre-compensation rules. For the initial mask design, together with a set of process conditions for an actual photolithographic processing system (e.g., a set of illumination/projection conditions for a “stepper” and a set of chemical processing conditions), a simulation is performed that generates a simulated image of the pattern that would appear on the wafer after the photoresist was exposed and chemically processed. The simulated image is compared to the desired circuit pattern, and deviations from the desired circuit pattern are determined. The mask design is then modified based on the deviations, and the simulation is repeated for the modified mask design. Deviations from the desired circuit pattern are again determined, and so on, the mask design being iteratively modified until the simulated image agrees with the desired circuit pattern to within an acceptable tolerance. The accuracy of the simulated image is, of course, crucial in obtaining OPC-generated mask designs that lead to acceptable results in the actual production stepper machines and chemical processing systems of the actual integrated circuit fabrication environments.

[0005] Photolithographic process simulation generally comprises optical exposure simulation and chemical processing simulation. During optical exposure simulation, the operation of a photolithographic processing system (e.g., stepper) is simulated to compute an optical intensity pattern in the photoresist. The computed optical intensity can be purely two-dimensional, i.e., a function of two variables that treats the photoresist layer as a single plane. Alternatively, the optical exposure simulation can treat the photoresist layer as a volume, and compute two-dimensional optical intensities for a plurality of planar levels within the photoresist volume. In still another alternative, the optical exposure simulation can provide a three-dimensional optical intensity.

[0006] During chemical processing simulation, the effects of chemical processing (which may include, for example, treating the exposed/unexposed resist, washing away the exposed/unexposed resist, wet or dry etching, etc.) are simulated to compute a resultant intensity pattern in the resist and/or the wafer. In some cases the chemical processing simulation is carried out in a very simple manner, such as by simply thresholding a purely two-dimensional optical intensity pattern by a fixed, predetermined value to generate a resultant intensity that is a binary or continuous pattern. In other cases, the chemical processing simulation can be more complex by processing the optical intensity pattern in location-dependent, pattern-dependent, or depth-dependent manners, or in other ways. In still other cases, such as during the design or troubleshooting of the stepper machines themselves by a stepper machine manufacturer, the chemical processing simulation is omitted altogether (i.e., is a null step in the photolithographic process simulation), with the two-dimensional or three-dimensional optical intensity pattern itself being of primary interest to the user.

[0007] For purposes of clarity and not by way of limitation, the output of an optical exposure simulation is referenced herein as an optical intensity. As used herein, resultant intensity refers to the output of an overall photolithographic process simulation. In cases where the chemical processing simulation is a null step in the overall photolithographic process simulation, the resultant intensity corresponds to an optical intensity. Depending on the particular need, the resultant intensity can represent any of a variety of different images or physical quantities including, but not limited to, an optical intensity, an aerial image, a latent image in a resist film, a developed resist pattern, a final semiconductor pattern, and an etch depth profile within a final semiconductor.

[0008] It is crucial to accurately compute the optical intensity pattern in the photoresist and to accurately compute the chemical processing effects in order to enable the overall photolithographic process simulation result to be accurate. In addition to the OPC software context, accurate photolithographic process simulation is also useful in other contexts including, but not limited to, resolution enhancement techniques (RETs) additional to OPCs, in which the effectiveness of a particular RET strategy can be verified by simulation.

[0009] One or more issues arise in relation to the simulation of a photolithographic process, at least one of which is resolved by at least one of the preferred embodiments described herein. It is to be appreciated, however, that the scope of the preferred embodiments is also widely applicable to a variety of different optical exposure processes, with or without subsequent chemical processing steps, and is not limited to the particular environment of integrated circuit fabrication. One key issue relates to computation time. For example, as indicated in US2005/0015233A1, which is incorporated by reference herein, computation of a single optical intensity for a single resist plane and a single set of process conditions using conventional methods can take hours or even days. In practice, it is often desired to compute the optical intensity for many different values of one or more exposure system variations, and/or to compute the optical intensity for many different resist planes. It is further often desired to compute the resultant intensity for many different values of one or more exposure system or chemical processing system variations, and/or to compute the resultant intensity for many different planes within the photoresist layer or the processed wafer. This can drastically increase the already-substantial computation time in order to compute the desired overall set of results. Other issues arise as would be apparent to one skilled in the art upon reading the present disclosure.

[0010] A method, system, related computer program products, and related computer-readable numerical arrays for computer simulation of a photolithographic process are provided. In one preferred embodiment, a first resultant intensity for a first value of a process variation associated with the photolithographic process is computed by processing a mask transmittance function with a plurality of model kernels characterizing the photolithographic process. Subsequent resultant intensities for subsequent values of the process variation are computed in a manner that does not require additional processing of the mask transmittance function in conjunction with the model kernels, each subsequent resultant intensity being computed in a substantially shorter time than that needed to compute the first resultant intensity. Thus, where computation for a single process variation value requires a first time interval T, results for N different process variation values can be achieved within a second time interval that is substantially less than NT.

[0011] According to another preferred embodiment, a first resultant intensity for a first target depth is computed by processing a mask transmittance function with a plurality of model kernels characterizing the photolithographic process. Subsequent resultant intensities for subsequent target depths are computed in a manner that does not require additional processing of the mask transmittance function in conjunction with the model kernels, each subsequent resultant intensity being computed in a substantially shorter time than that needed to compute the first resultant intensity.

[0012] According to another preferred embodiment, a first resultant intensity for a first combination of values for at least two process variations is computed by processing a mask transmittance function with a plurality of model kernels characterizing the photolithographic process. Subsequent resultant intensities for subsequent combinations of values for the at least two process variations are computed in a manner that does not require additional processing of the mask transmittance function in conjunction with the model kernels, each subsequent resultant intensity being computed in a substantially shorter time than that needed to compute the first resultant intensity.

[0013] According to another preferred embodiment, a first resultant intensity for a first combination of values for a target depth and at least one process variation is computed by processing a mask transmittance function with a plurality of model kernels characterizing the photolithographic process. Subsequent resultant intensities for subsequent combinations of values for the target depth and the at least one process variation are computed in a manner that does not require additional processing of the mask transmittance function in conjunction with the model kernels, each subsequent resultant intensity being computed in a substantially shorter time than that needed to compute the first resultant intensity.

[0014] According to another preferred embodiment, a process model of a photolithographic process for use in a computer simulation thereof is computed, the photolithographic process having physical parameters associated therewith, the process model comprising model kernels. A transmission cross coefficient (TCC) is computed based on the physical parameters, the TCC characterizing the photolithographic processing system and incorporating a process variation and/or a target depth according to a series expansion having TCC basis spatial functions not dependent on that process variation and/or target depth and at least one TCC expansion coefficient dependent on that process variation and/or target depth. The model kernels are computed from the basis spatial functions.

[0015] According to another preferred embodiment, a process model of a photolithographic process for use in a computer simulation thereof is computed, the photolithographic process having physical parameters associated therewith, the process model comprising model kernels and model kernel weights. A mutual intensity function is computed based on the physical parameters, the mutual intensity function characterizing an illumination system that illuminates a mask. A plurality of coherent modes and a corresponding plurality of coherent mode weights associated with the mutual intensity function are computed. A point spread function based on the physical parameters is computed, the point spread function characterizing a projection system that projects the illuminated mask toward a target to expose the target. The point spread function may be further extended to include certain effects due to a chemical processing system that chemically processes the exposed target to generate a resultant intensity. The model kernels and model kernel weights are computed using the coherent modes, the coherent mode weights, and the point spread function. In other preferred embodiments, calibration methods are provided for calibrating process models to physical photolithographic processing systems.

[0041] FIG. 1A illustrates a simplified representation of a photolithographic processing system 102 for computer simulation according to one or more of the preferred embodiments described herein. Photolithographic processing system 102 comprises an optical exposure system 105 and a chemical processing system 114. Optical exposure system 105 comprises an illumination system 106 that illuminates a mask 108, such as a photolithographic mask, and a projection system 110 that projects the illuminated mask toward a target 112 to expose the target 112. The illumination system 106 includes an optical source 104. The optical exposure system 105 may represent a photolithographic stepper machine (“stepper”), and the target 112 can be a semiconductor wafer covered with photoresist (“resist”), although the scope of the present teachings is not so limited.

[0042] Chemical processing system 114 comprises devices that chemically process the target 112 after exposure, and also represents the photochemical reaction of resist molecules upon absorption of exposure photons. Subsequent to the exposure process, the chemical processing is often performed after the target 112 has been physically removed from the optical exposure system 105. For clarity of presentation and not by way of limitation, the chemical processing system 114 is conceptually illustrated by an enlarged arrow beneath the projection system 110 in FIG. 1A. Prior to chemical processing (i.e., above the enlarged arrow in FIG. 1A) the effects of the photolithographic system 102 on the target 112 can be characterized by an optical intensity pattern. As used herein, an optical intensity value can be a time-constant or time-varying instantaneous value and/or can be a time integral of the instantaneous value. Subsequent to chemical processing (i.e., below the enlarged arrow in FIG. 1A) the effects of the photolithographic system 102 on the target 112 can be characterized by a resultant intensity pattern corresponding to any of a variety of physical results including, but not limited to, an aerial image, a latent resist image, a developed resist pattern, a final semiconductor pattern, and an etch depth profile within a final semiconductor. If the chemical processing is a null step (i.e., if chemical processing system 114 is not present or turned off), the resultant intensity for the target 112 corresponds to the optical intensity for the target 112.

[0043] In general, the simulation process comprises receiving known physical parameters of the photolithographic processing system 102, generating mathematical models that represent its operation, and applying those mathematical models to proposed transmittance functions of the mask 108 to compute resultant intensities. For purposes of generality, the target 112 is illustrated in FIG. 1A as a volumetric object. However, it is to be appreciated that optical intensities and resultant intensities can be either two-dimensional profiles for single planes thereon or therein, or can be three-dimensional profiles for the entire volume or a subvolume therein.

[0044] Among the advantages of at least one of the preferred embodiments described herein, easier computation of I_TARGET(x_t,y_t;z_t;var) is provided, where z_t represents a target depth into the target 112, and where var represents one or more process variations. Process variations include optical exposure system variations for the optical exposure system 105 and chemical processing variations for the chemical processing system 114. Optical exposure system variations include illumination system variations for the illumination system 106 and projection system variations for the projection system 110. Examples of illumination system variations include, but are not limited to, intensity/position variations of multiple emitters, source chromatic variations, coherence variations, and source positioning variations. Examples of projection system variations include, but are not limited to, defocus variations, lens aberrations, immersion medium refractive index variations, immersion medium attenuation coefficient variations, stack film thickness variations, stack film material refractive index variations, and stack film material attenuation coefficient variations. A defocus variation can arise, for example, from mechanical errors in the vertical position of a wafer-supporting platform in a stepper machine relative to the optics of the projection system 110. Examples of chemical processing variations include, but are not limited to, resist development variations and chemical etching process variations. Although just the single term “var” is often used herein to denote process variations, it is to be appreciated that “var” can represent a plurality of process variations var1, var2, var3, etc.

[0045] As used herein, value of a process variation refers broadly to any quantitative descriptor, index, or indicator reflecting an amount, degree, class, type, magnitude, direction, nature, or other characteristic associated with the process variation. Although in many cases the value of a process variation will be a scalar number (e.g., a scalar distance metric for the case of defocus variations), in other cases the value of a process variation may be a vector or multidimensional expression, an index into a lookup table, a percentage, a color, a shape, or other indicator, including a null value if the particular process variation does not exist, has no effect, or is not to be considered. It is to be appreciated that a value of a process variation, as that term is used herein, might only be abstractly related to the physical underpinnings of the process variation itself, provided only that some kind of change in the process variation is associated with some kind of change in the value.

[0046] According to one of the computational advantages of at least one of the preferred embodiments described herein, the number of computationally expensive convolutions of the mask transmittance function with model kernels does not increase with the number of different value combinations for one or more process variations and/or target depths for which target intensity results are desired. Instead, according to at least one of the preferred embodiments herein, only a single set of convolutions of the mask transmittance function with model kernels is required for computing target intensity results for a number of different value combinations for the one or more process variations and/or target depths. Once this single set of convolutions is computed for a first value combination (for example, var=[defocus₁aberration₁] together with depth=depth₁), the same convolution results are then combined differently according to variation-dependent and/or depth-dependent weights for generating results for different value combination (for example, var=[defocus₂,aberration₂] together with depth=depth₂). Notably, combining the same convolution results with differing weights is substantially less computationally intensive than generating different convolution results for each different process variation value and/or target depth. Accordingly, when simulating a photolithographic process according to one or more of the preferred embodiments, substantial computation time is saved when resultant intensities are desired for several different value combinations for the one or more process variations and/or target depths.

[0047] For one or more of the preferred embodiments described herein, for an initial process variation value and/or an initial target depth, the number of model kernels with which the mask transmittance function is convolved can be larger than a number of optical kernels used in one or more prior art methods. Accordingly, for any particular mask transmittance function, if only a single aerial image is desired for a single process variation value for a single target depth, the computational advantages may not be readily apparent. However, for the more prevalent case in which results for several process variation values and/or target depths are desired, drastic computational advantages over the one or more prior art methods are realized.

[0048] Referring again to FIG. 1A, coordinates (x_s,y_s) are used for the plane of the source 104, coordinates (x_m,y_m) are used for the plane of the mask 108, and coordinates (x_t,y_t) or (x_t,y_t,z_t) are used for the target 112. FIG. 1B illustrates a conceptual view of the optical source 104, the mask 108, and the target 112. For clarity of presentation, image reversal and image reduction (e.g., 4:1) associated with the projection system 110 are not included in FIG. 1B and herein, as a person skilled in the art would be readily able to implement the preferred embodiments factoring in the appropriate image reversal/reduction in view of the present disclosure.

[0049] FIG. 1B illustrates the mask 108 and target 112 as subdivided for computational purposes into mask segments 119 and target segments 121, respectively. As known in the art, adjacent ones of mask segments 119 and target segments 121 may have slight overlaps according to an optical ambit of the projection system 110 and diffusion radius of the chemical processing system 114, respectively. Due to the spatially-limited optical ambit of the projection system 110, any particular one of the mask segments 119 (e.g., an exemplary mask segment 120) can be processed separately from the other mask segments 119 to compute a resultant intensity at a corresponding target segment 121 (e.g., an exemplary target segment 122 corresponding to the exemplary mask segment 120). By way of example, and not by way of limitation, for photolithographic simulation for chips of size (1 mm)² to (10 mm)² there may be between (100)² and (1000)² mask segments 119 and corresponding target segments 121, each having a size of about (10 μm)². In turn, each mask segment 119 may comprise M×M numerical values representing mask transmittances within that segment, with M being about 100-1000, with resultant intensity values being computed for each of M×M corresponding target segment locations. Another suitable representation for the mask 108 is binary or phased polygons specified by coordinates of vertices. Typical sampling grid spacing for both mask and target may be on the order of 10-100 nm. The optical ambit for such cases may be on the order of 1 μm, corresponding to between 10 sample points and 100 sample points along each dimension. It is to be appreciated, however, that such computational array sizes and sampling grids are highly dependent on the particular nature of the exposure process being simulated, the type and speed of computing hardware being used, and other factors, and therefore can vary by orders of magnitude without departing from the scope of the preferred embodiments.

[0050] The computations for the different mask segments 119 and their corresponding target segments 121 can be performed separately and in parallel on separate computing devices to reduce the overall simulation time required. Under an assumption of global shift-invariance of the optical exposure system 105, the different computing devices can be programmed with the same simulation algorithm. However, the scope of the preferred embodiments is not limited to the use of the same simulation algorithm for different mask segments. Moreover, the scope of the preferred embodiments can include scenarios in which different mask segments are of different size and/or shape. As illustrated in FIG. 1B, the exposure process can be conceptually broken down into the illumination of an exemplary mask segment 120 by the illumination system 106, and the projection of the illuminated mask segment 120 onto the corresponding target segment 122.

[0051] FIG. 2 illustrates simulating a photolithographic process according to a preferred embodiment. At step 202, a mutual intensity function is computed characterizing illumination of the mask 108 by the illumination system 106. At step 204, a plurality of coherent modes and weights associated with the mutual intensity function are computed. Steps 202 and 204 are discussed further infra in relation to FIG. 3.

[0052] At step 206, a point spread function is computed that characterizes the projection system and the chemical processing system (if the chemical processing is not a null step), the point spread function incorporating at least one process variation for at least one of the projection system and the chemical processing system. At step 208, a process model is computed using the coherent modes and the point spread function, the process model comprising model kernels. At step 210, the mask transmittance function is convolved with the model kernels to generate convolution results. At step 212, a resultant intensity in the target 112 is computed as a weighted sum of the squares of the convolution results. Steps 206-212 are discussed further infra in relation to FIGS. 4-5. For a scenario in which the chemical processing is a null process (e.g., if the user presumes that the chemical processing system is not present, or is turned off, or if the chemical processing system is simulated separately from the optical exposure system), the process model is an optical model, the model kernels are optical kernels, and the resultant intensity is an optical intensity as those terms are used herein.

[0053] FIG. 3 illustrates a portion of the photolithographic processing system 102 of FIG. 1A and expressions relating to describing photolithographic process simulation according to a preferred embodiment. In particular, FIG. 3 illustrates the illumination system 106 and mask 108 of FIG. 1A together with a plane 108⁻ immediately prior to the mask 108 and a plane 108⁺ immediately subsequent to the mask 108. Also shown in FIG. 3 is a mutual intensity function 301 representing incident electromagnetic radiation at plane 108⁻. Mutual intensity functions are discussed generally in M. Born and E. Wolf, Principles of Optics, 7^th ed., Cambridge University Press (1999) at Chapter 10. As known in the art, a mutual intensity function at a particular plane is a complex-valued function of four variables (the x and y locations of two points within the plane) and carries substantially “more information” about the electromagnetic radiation than a “regular” optical intensity function, which is a non-negative, real-valued function of a single point within the plane and which represents a special case of the mutual intensity function when its two points are on top of each other. Although it is usually only necessary to compute the “regular” optical intensity function at the target (related to the resultant intensity through chemical processing effects, if present), consideration of mutual intensities at intermediate locations is advantageous for properly computing that optical intensity in accordance with a preferred embodiment.

[0054] The mutual intensity function 301 can be computed at step 202 of FIG. 2, supra, using known methods based on the physical parameters of the illumination system 106 including the optical source 104. Nominal or “starter” values for the illumination system physical parameters and projection system physical parameters are usually provided with each stepper device. Notably, the illumination system physical parameters that are used to generate the mutual intensity function 301 can be generated, derived, assumed, etc., by any of a variety of different methods without departing from the scope of the preferred embodiments. For example, during a simulator calibration process described further with respect to FIG. 15B, infra, the illumination system physical parameters can themselves be the subject of an optimization process in which they are iteratively refined.

[0055] The mutual intensity function 301 is shift-invariant which, in conjunction with a spatially limited optical ambit, allows windowing to a smaller size for computational purposes to form a windowed mutual intensity function 302. The window size for the windowed mutual intensity function 302 can be, for example, about twice the size of the optical ambit on a side. Referring to an example scenario in which the optical ambit is 1 μm and in which the sampling grid is 50 nm, the windowed mutual intensity function 302 would be an M′×M′×M′×M′ array of complex values, where M′ is about 40.

[0056] Shown at FIG. 3 (Eq. {3c}) is a representation of the windowed mutual intensity function 302 in terms of its coherent modes 308 and coherent mode weights 306, as computed at step 204 of FIG. 2, supra. In one preferred embodiment, the coherent modes 308 are eigenvectors of a linear operator with the windowed mutual intensity function 302 as a Hilbert-Schmidt kernel, and can be generated therefrom using known, off-the-shelf eigenanalysis software. Referring to the above example scenario in which the windowed mutual intensity 302 is an M′×M′×M′×M′ array of complex values, each coherent mode 308 would be an M′×M′ array of complex values and each coherent mode weight 306 would be a scalar quantity. The expansion shown in FIG. 3 (Eq. {3c}) can be referred to as the coherent mode representation of the partially coherent illumination field. Although the number of coherent modes 308 and coherent mode weights 306 may be theoretically infinite, acceptable results can be achieved by limiting the expansion to a maximum number nmax, with an example value of nmax being about 10.

[0057] Also illustrated in FIG. 3 is a mutual intensity function 310 representing electromagnetic radiation at plane 108⁺. The mutual intensity function 310 is equal to the incident mutual intensity function 301 at plane 108⁻ modified by point-wise post-multiplication with the mask transmittance function 312 and point-wise pre-multiplication with the conjugate of the mask transmittance function 312 as indicated in FIG. 3 Eq. {3b}).

[0058] FIG. 4 illustrates a portion of the photolithographic processing system 102 of FIG. 1A and expressions relating to describing photolithographic process simulation according to a preferred embodiment. In particular, FIG. 4 illustrates the mask 108, the projection system 110, the chemical processing system 114, and the target 112 of FIG. 1A. FIG. 4 also conceptually illustrates at least one process variation 403 (“var”) representing projection system and/or chemical processing system variations. According to one preferred embodiment, the projection system 110 and chemical processing system 114 are characterized by a point spread function representing a response to a conceptual unit point emitter 402 at the mask plane 108, wherein the point spread function incorporates the at least one process variation. According to another preferred embodiment, the projection system 110 is characterized by a point spread function representing a response to the conceptual unit point emitter 402, wherein the point spread function incorporates a target depth for different target depths z_t. According to another preferred embodiment, the projection system 110 is characterized by a point spread function representing a response to the conceptual unit point emitter 402, wherein the point spread function incorporates both the at least one process variation and the target depth.

[0059] In accordance with the mask and target segmentations described supra in relation to FIG. 1B, for points having lateral locations (x_t,y_t) within the target segment 122 it is usually only necessary to consider the effects of conceptual unit point emitters at locations (x_m,y_m) within the mask segment 120, although the scope of the preferred embodiments is not so limited. Furthermore, it is only necessary to consider a small window of optical ambit with side equal to, say, twice the optical ambit, which has M′×M′ sampling points, with M′ being about 40, although the scope of the preferred embodiments is not so limited. Illustrated in FIG. 4 is a point spread function 406 that incorporates, in addition to the location arguments (x_m,y_m) and (x_t,y_t), the arguments z_t and var. Preferably, the point spread function 406 is characterized as being shift-invariant such that its value depends on the relative locations of the points (x_m,y_m) and (x_t,y_t) rather than their absolute locations, as illustrated at FIG. 4 (Eq. {4a}).

[0060] Referring again to step 206 of FIG. 2, the point spread function 406 can be computed using any of a variety of methods, whether analytical, numerical, or both, as a function of the spatial indices (x_m,y_m;x_t,y_t) and the depth and variation indices (z_t,var). There are several ways that this can be achieved, as would be apparent to one skilled in the art in view of the present disclosure. In one preferred embodiment, the point spread function 406 is characterized as being dependent on the process variations/target depth (var,z_t) according to a series expansion having basis spatial functions 410 not dependent on the process variations/target depth (var,z_t) and having at least one expansion coefficient 408 dependent on the process variations/target depth (var,z_t), as illustrated at FIG. 4 (Eq. {4b}). The particular basis spatial functions 410 of FIG. 4, which correspond to the PSF, may be called PSF basis spatial functions, and the expansion coefficients 408 may be called PSF expansion coefficients. The PSF basis spatial functions 410 are also termed partial response functions 410 hereinbelow, while the PSF expansion coefficients 408 are also termed partial response coefficients 408 hereinbelow. One specific formulation for achieving such characterization when there is a defocus variation z′ is described infra in relation to FIG. 8 (Eq. {8h}). Another formulation for achieving such characterization, which is for a more general case in which the process variations/target depth (var,z_t) can be modeled as affecting a model pupil function is described infra in relation to FIG. 17, although other methods may be used as well.

[0061] The terms basis spatial function and expansion coefficient are clarified here such that simple nomenclature adjustments will not cause departure from the scope of the preferred embodiments. Let a generic function f(x,y;var) be expandable in the form f(x,y;var)=g₁(var)h₁(x,y)+g₂(var)h₂(x,y)+g₃(var)h₃(x,y)+ . . . , where var represents one or more variation variables, and where (x,y) represent two or more spatial variables or indices into numerical arrays that depict spatially dependent values. Under nomenclature consistent with the description of FIG. 4 (Eq. {4b}), and as used herein, the functions h_k are termed basis spatial functions and the functions g_k shall be termed expansion coefficients.

[0062] FIG. 5 illustrates a portion of the photolithographic processing system 102 of FIG. 1A and expressions relating to describing photolithographic process simulation according to a preferred embodiment. In particular, FIG. 5 illustrates the plane 108⁺ immediately subsequent to the mask 108, as well as the projection system 110, the chemical processing system 114, and the target 112 of FIG. 1A. FIG. 5 also conceptually illustrates at least one process variation 403 (“var”) representing projection system and/or chemical processing system variations. FIG. 5 also illustrates equations based upon which the target intensity 502 may be computed according to one or more of the preferred embodiments.

[0063] Beginning with a first-principle quadruple-integral at FIG. 5 (Eq. {5a}) based upon the post-mask mutual intensity function 310 and the point response function 406, and incorporating the photolithographic processing system characterizations described supra in relation to FIGS. 1A-4, FIG. 5 illustrates how the target intensity 502 can be computed in one or more ways according to one or more of the preferred embodiments. In one preferred embodiment, with reference again to steps 208-212 of FIG. 2, a formulation shown at FIG. 5 (Eq. {5m}) can be used. For step 208, a process model comprising model kernels and model kernel weights is computed using the coherent modes 308 and the point spread function 406. For one preferred embodiment consistent with FIG. 5 (Eq. {5m}), model kernels 507 are generated by multiplying the point spread function 406 with each of the coherent modes 308 generated at step 204 of FIG. 2, and the model kernel weights are set equal to the respective coherent mode weights 306. For step 210, the mask transmittance function 312 is convolved with each model kernel 507 to generate convolution results 504. Finally, for step 212, the squared magnitudes of the convolution results 504 are generated and a weighted sum thereof is generated according to the coherent mode weights 306, which form the model kernel weights for this embodiment.

[0064] FIG. 6A illustrates simulating an exposure process according to a preferred embodiment. At step 602, a process model is computed, the process model comprising a plurality of model kernels and model kernel weights. At step 604, the mask transmittance function is convolved with each model kernel to generate a plurality of convolution results. At step 606, an initial resultant intensity for an initial process variation value is computed as a first combination of the convolution results, e.g., combined using an initial set of model kernel weight values. Finally, at step 608, resultant intensities are computed for subsequent process variation values as subsequent combinations of the convolution results, e.g., combined using subsequent sets of model kernel weight values that numerically depend on the degree of process variation. In one preferred embodiment, the process variations are projection system variations. According to another preferred embodiment, steps similar to those of FIG. 6A are carried out in a manner that includes different target depths, separately or in various combinations (var,z_t) with the different process variation values.

[0065] A model kernel at least partially characterizes the manner in which a photolithographic processing system generates a resultant intensity from a mask. For one preferred embodiment, model kernel refers to one of a set of model kernels that yield a resultant intensity when processed in conjunction with a mask transmittance function, the processing comprising convolving the mask transmittance function with each model kernel to generate a plurality of convolution results and coherently and/or incoherently combining the convolution results. Preferably, the combination of the convolution results is carried out according to one or more model kernel weights associated with the plurality of model kernels. The photolithographic processing system may be with or without a chemical processing system.

[0066] Examples of model kernels include imaging kernels and optical kernels. An imaging kernel at least partially characterizes the manner in which a photolithographic processing system generates a resultant intensity from a mask. For one preferred embodiment, imaging kernel refers to one of a set of imaging kernels that yield a resultant intensity when processed in conjunction with a mask transmittance function, the processing comprising convolving the mask transmittance function with each image kernel to generate a plurality of convolution results and computing a weighted sum of the squared magnitudes of the convolution results according to one or more imaging kernel weights associated with the imaging kernels.

[0067] For the particular case in which no chemical processing is simulated (e.g., if the user presumes that the chemical processing system is not present, or is turned off, or if the chemical processing system is simulated separately from the optical exposure system), imaging kernels operate as optical kernels. An optical kernel at least partially characterizes the manner in which an optical exposure system generates an optical intensity from a mask. For one preferred embodiment, optical kernel refers to one of a set of optical kernels that yield a resultant intensity when processed in conjunction with a mask transmittance function, the processing comprising convolving the mask transmittance function with each optical kernel to generate a plurality of convolution results and computing a weighted sum of the squared magnitudes of the convolution results according to one or more optical kernel weights associated with the optical kernels. The resultant intensity so generated using optical kernels is, in particular, an optical intensity.

[0068] Examples of model kernels further include partial kernels. A partial kernel at least partially characterizes the manner in which a photolithographic processing system generates a resultant intensity from a mask. For one preferred embodiment, partial kernel refers to one of a set of partial kernels that yield a resultant intensity when processed in conjunction with a mask transmittance function, the processing comprising convolving the mask transmittance function with each partial kernel to generate a plurality of convolution results and combining the convolution results using both coherent and incoherent weighted sums according to one or more partial kernel weights associated with the partial kernels. More particularly, the combining of the convolution results comprises computing a plurality of coherent weighted sums across related groupings of the convolution results, and computing an incoherent weighted sum of the squared magnitudes of the coherent weighted sums, each of the coherent weighted sums and incoherent weighted sums being weighted according to partial kernel weights. In one example, the convolution results are grouped for purposes of coherent weighted summation according to a coherent mode expansion index associated with expansion of a mutual intensity function into coherent modes.

[0069] The preferred embodiment of FIG. 5, supra, represents one example for which the model kernels are imaging kernels and the model kernel weights are imaging kernel weights. The preferred embodiment of FIG. 7, infra, represents one example for which the model kernels are partial kernels and the model kernel weights are partial kernel weights. The preferred embodiment of FIG. 12, infra, represents one example for which the model kernels are optical kernels and the model kernel weights are optical kernel weights. The preferred embodiment of FIG. 16A, infra, represents another example for which the model kernels are imaging kernels and the model kernel weights are imaging kernel weights.

[0070] FIG. 6B illustrates a computation time plot for illustrating one of the features and advantages of photolithographic process simulation according to a preferred embodiment. In particular, FIG. 6B illustrates plots of total computation time versus the number of process variation value/target depth combinations (var,z_t) for which resultant intensities are computed for a particular mask. Advantageously, once the resultant intensity is computed for an initial combination (var,z_t), the resultant intensity for subsequent additional combinations are readily and relatively quickly computed by re-using the convolution results associated with the initial combination (var,z_t), because the computational complexity for data convolutions is, generally speaking, substantially higher that the computational complexity for operations such as squaring, weighting, and addition for the sizes of data involved in photolithographic simulations such those presented herein. By way of example, where computation for an initial combination requires a time T, computation time for each subsequent combination will require a time Δ=T/F₁, where F₁ can be termed a subsequent computation time savings factor and is substantially greater than unity (e.g., between 5 and 10, for example, or even larger depending on the data sizes involved). Overall computation time for N different combinations using the preferred embodiment of FIG. 6A is therefore equal to T+(N−1)Δ, which is substantially less than NT, as indicated by the plots 654 and 650 in FIG. 6B.

[0071] Also shown in FIG. 6B is a plot 652 representing computation time for N different (var,z_t) combinations according to one or more kernel-based prior art methods. Generally speaking, it is not impossible that there might be fewer kernels used in the one or more prior art methods. In such cases, if they exist, there can be fewer initial convolutions necessary than for the preferred embodiment of FIG. 6A. For purposes of illustration, such cases are reflected in the plot 652 by an initial computation time P=T/F₂, with F₂ being greater than unity. Notably, however, the subsequent computation time savings factor F₁ for the preferred embodiment of FIG. 6A will usually be greater than and so the initial time deficit for the preferred embodiment of FIG. 6A will be compensated as N increases, as shown at the crossover locations of the plots 654 and 650. In the particular example of FIG. 6B, for which F₁ is assumed to be about 6 and F₂ is assumed to be about 3, the overall time savings is realized as N rises above about 6.

[0072] It is often the case that a photolithographic process simulation is repeated for different combinations of process variation values and/or different target depths for selected areas or positions on a chip rather than for the whole chip, the selected areas or positions being irregularly and sparsely distributed. It has been found that even further computational efficiencies can be realized using the preferred embodiment of FIG. 6A when resultant intensities need to be calculated at such irregularly and sparsely distributed target points. Generally stated, these computational efficiencies can have the effect of drastically reducing the slope of the plot 654 of FIG. 6B toward a horizontal orientation (i.e., drastically increase the value of the subsequent computation time savings factor F₁). This effect is discussed further infra in the present disclosure.

[0073] FIG. 7A illustrates a photolithographic processing system and expressions relating to describing photolithographic process simulation according to a preferred embodiment. More particularly, FIG. 7A illustrates a conceptual diagram of the photolithographic processing system of FIG. 1A in which the projection system 110 and chemical processing system 114 are in combined block form between the mask 108 and the target 112. Also illustrated are the illumination system 106 and the at least one process variation 403 (“var”) representing a projection system and/or chemical processing system variation.

[0074] Generally speaking, FIG. 7A illustrates an example of a computational characterization for a resultant intensity 702 that is consistent with the method of FIG. 6A, wherein the model kernels are partial kernels and the model kernel weights are partial kernel weights. It is to be appreciated, however, that other computational characterizations consistent with the method of FIG. 6A are within the scope of the preferred embodiments, such as those described further infra. Shown in FIG. 7 (Eq. {7a}) is a recast version of FIG. 5 (Eq. {5m}) with the point spread function 406 replaced by the series expansion expression of FIG. 4 (Eq. {4b}). Terms are rearranged in accordance with standard algebraic principles to yield the expression of FIG. 7 (Eq. {7b}).

[0075] With reference again to FIG. 6A in view of FIG. 7A, the process model referenced at step 602 can be generated by (i) computing a plurality of partial kernels 704 by multiplying coherent modes 308 (referenced supra with respect to FIG. 3) and partial response functions 410 (referenced supra with respect to FIG. 4), and (ii) forming the partial kernel weights as a union of coherent mode weights 306 (referenced supra with respect to FIG. 3) and partial response coefficients 408 (referenced supra with respect to FIG. 4). For step 604, convolution results 706 are generated by convolving the mask transmittance function 312 with each of the partial kernels 704, the convolution results 706 being advantageously applicable for computation of results for both initial and subsequent process variation values.

[0076] For step 606, resultant intensity results are computed for the initial process variation value by combining the convolution results 706, the combining comprising (i) forming a plurality of coherent weighted sums 708 using a first portion of the partial kernel weights corresponding to the partial response coefficients 408, and (ii) computing an incoherent weighted sum of the squared magnitudes of the coherent weighted sums 708 using a second portion of the partial kernel weights corresponding to the coherent mode weights 306. Stated another way, the convolution results are grouped by coherent mode expansion index “n” for purposes of computing the coherent weighted sums 708, and then the incoherent weighted sum of the squared magnitudes of the coherent weighted sums 708 is computed using the coherent mode weights 306.

[0077] For step 608, resultant intensity results are computed for subsequent process variation values by combining the convolution results 706 in a manner similar to that described for step 606, except that the partial response coefficients 408 that make up the second portion of the partial kernel weights are different because they are dependent on the process variation. Advantageously, as evidenced at least by FIG. 7 (Eq. {7b}), the resultant intensity can be computed for different process variation values using the same mask transmittance function convolution results 706 that were computed for the first process variation value. Importantly, it is the convolutions involving the mask transmittance function, i.e., the convolutions that yield the convolution results 706, that represent the bulk of the computational complexity needed for the simulation. Other computation steps, such as weighting, squaring, and summing of convolution results 706, are of generally minor computational complexity compared to such convolutions. Accordingly, the step 608 for each subsequent process variation value takes substantially less time to achieve than the step 606 for the initial process variation value, because the convolution results 706 have already been computed prior to step 608. Similar advantages apply in another preferred embodiment in which steps similar to those of FIG. 6A are carried out in a manner that includes different target depths, separately or in various combinations (var,z_t) with the different process variation values.

[0078] FIG. 7B illustrates a conceptual view of the formation of the partial kernels 704 referenced supra with respect to FIG. 7A. At the left side of FIG. 7B is an (Lmax+1)-member set of partial response functions 410 and an nmax-member set of coherent modes 308 as derived supra for a particular set of physical parameters of the photolithographic process. Referring again to the M′×M′ window of optical ambit scenario supra, the partial response functions 410 and coherent modes 308 would each be M′×M′ arrays of complex values. At the right side of FIG. 7B is an (Lmax+1)(nmax)-member set of partial kernels 704, any particular member PKRN_n,L comprising an point-wise product of the n^th coherent mode and the L^th partial response function. Importantly, in accordance with the features and advantages of at least one of the preferred embodiments, the partial kernels 704 do not depend on the process variation values, target depths, or mask transmittance function. Moreover, the convolution of a particular mask transmittance function with the partial kernels 704 only needs to be performed once for computing target intensities across many different process variation values and/or target depths.

[0079] FIG. 8 illustrates a portion of the photolithographic processing system 102 of FIG. 1A and expressions relating to describing photolithographic process simulation according to a preferred embodiment, for the purposes of illustrating one example of the manner in which a point spread function incorporating at least one process variation can be generated, for a particular example in which there is a defocus variation represented by z′. It is to be appreciated, however, that there are many ways for generating a point response function incorporating at least one process variation according to a preferred embodiment as would be readily apparent to a person skilled in the art in view of the present disclosure.

[0080] In particular, FIG. 8 illustrates the mask 108, the projection system 110, the chemical processing system 114, and the target 112 of FIG. 1A, with the projection system 110 being represented in “4F” format such that a model pupil aperture 802 is illustrated. For a variation-less projection system 110, the model pupil aperture 802 can be modeled by a pupil function 804 which is non-zero within a circular disc in the frequency domain, the pupil function 804 passing a maximum spatial frequency f_max related to a numerical aperture of the projection system 110. Under an assumption of shift-invariance, a point response function 806 simply corresponds to an inverse Fourier transform of the pupil function 804. For a defocus variation z′, however, there is a point response function 808 that is the inverse Fourier transform of the pupil function 804 as modified by a phase factor 810. The term η in FIG. 8 (Eq. {8b}) is a refractive index associated with material lying above the target 112 during the exposure process.

[0081] As illustrated in FIG. 8, the point response function 808 can be expressed as a sum of the responses through different annular portions 813 of the model pupil aperture 802 of the projection system 110 as modified according to the phase factor 810 corresponding to a process variation (defocus) of the projection system. Also, when each annular portion 813 is sufficiently narrow, the phase factor 810 may be approximated as a constant within each annular portion 813, as illustrated in FIG. 8 (Eq. {8d}), where f′_j is a constant magnitude of a spatial frequency in the j^th annular portion. Optionally, the phase factor 810 within each annular portion 813 may be further expanded into a truncated Taylor series with only a few terms. Alternatively, based on a Taylor series expansion of the phase factor 810 for the entire pupil aperture 802 (see FIG. 8 (Eq. {8e})), the point response function 808 can be characterized as being dependent on the defocus variation z′ according to a series expansion having partial response functions 814 not dependent on the defocus variation z′ and having at least one expansion coefficient 816 dependent on the defocus variation z′, as illustrated at FIG. 8 (Eqs. {8g} and {8h}). Although the number of partial response functions 814 and partial response coefficients 816 are theoretically infinite, acceptable results can be achieved by limiting the expansion to a maximum number Lmax, with an example value of Lmax being about 5.

[0082] FIG. 9 illustrates a multi-emitter illumination system and expressions relating to describing photolithographic process simulation according to a preferred embodiment. In particular, FIG. 9 illustrates an illumination system 902 comprising an optical source 904 having multiple emitters R1, R2, . . . , RL, some or all of which may be off-axis, along with expressions characterizing simplified coherent mode computation according to a preferred embodiment. As shown in FIG. 9, the mutual intensity 906 of radiation incident to the mask plane 108 can be expressed as a sum of contributions 908 from each individual emitter R1, R2, . . . RL. The contribution 908 from each emitter can be expressed in terms of its coherent modes 910 and coherent mode weights 912, in accordance with a preferred embodiment. Under a first assumption of shift-invariance and a second assumption that the emitters are identical and translationally congruent, any particular k^th one of the of the emitter coherent modes 910 can be expressed as a phase-shifted version of a common coherent mode 914, i.e., the common coherent mode 914 multiplied by a phase factor 916 equal to exp[-iPHASE(x_sRK,y_sRk;x_m,y_m)]. As indicated by its arguments, the phase factor 916 is a function of mask point location (x_m,y_m) relative to the k^th emitter location (x_sRk,y_sRk. Computations involving the coherent modes are advantageously simplified by about a factor of L, the number of geometrically congruent emitters at the source 904.

[0083] FIG. 10 illustrates simulating an exposure process according to a preferred embodiment, the exposure process comprising illumination of a mask and projection of the illuminated mask toward a target by an optical exposure system. Notably, the steps of FIG. 10 represent one of several ways in which the method of FIG. 6A can be carried out in accordance with one or more of the preferred embodiments, for a case in which chemical processing is not simulated. At step 1002, a transmission cross-coefficient (TCC) is computed incorporating at least one exposure system variation in the optical exposure system according to a series expansion having basis spatial functions not dependent on the exposure system variation and at least one expansion coefficient dependent on the exposure system variation. Such basis spatial functions and expansion coefficients for a TCC may be termed TCC basis spatial functions and TCC expansion coefficients. At step 1004, optical kernels and associated optical kernel weights are computed from the TCC basis spatial functions. At step 1006, a mask transmittance function is convolved with each of the optical kernels to produce a set of convolution results. At step 1008, an optical intensity is computed for an initial exposure system variation value as a first combination of the convolution results. Finally, at step 1010, optical intensities are computed for subsequent exposure system variation values as subsequent combinations of the convolution results, e.g., combined using subsequent sets of optical kernel weight values that numerically depend on the process variation values. In one preferred embodiment, the process variations are optical variations of the projection system. According to another preferred embodiment, steps similar to those of FIG. 10 are carried out in a manner that includes different target depths, separately or in various combinations (var,z_t) with the different process variation values.

[0084] FIG. 11 illustrates simulating a photolithographic process according to a preferred embodiment, the photolithographic process comprising illumination of a mask, projection of the illuminated mask toward a target, and chemical processing of the exposed target by a chemical processing system. Notably, the steps of FIG. 11 represent one of several ways in which the method of FIG. 6A can be carried out in accordance with one or more of the preferred embodiments, for a case in which chemical processing simulation is included. At step 1102, a transmission cross-coefficient (TCC) is computed incorporating at least one process variation associated with the photolithographic process according to a series expansion having basis spatial functions not dependent on the process variation and at least one expansion coefficient dependent on the process variation. At step 1104, imaging kernels and associated imaging kernel weights are computed from the TCC basis spatial functions. At step 1106, a mask transmittance function is convolved with each of the imaging kernels to produce a set of convolution results. At step 1108, a resultant intensity is computed for an initial process variation value as a first combination of the convolution results. Finally, at step 1110, resultant intensities are computed for subsequent process variation values as subsequent combinations of the convolution results, e.g., combined using subsequent sets of imaging kernel weight values that numerically depend on the process variation values. In one preferred embodiment, the process variation is a projection system variation and/or a chemical processing system variation. According to another preferred embodiment, steps similar to those of FIG. 11 are carried out in a manner that includes different target depths, separately or in various combinations (var,z_t) with the different process variation values.

[0085] FIG. 12 illustrates portions of the photolithographic processing system 102 of FIG. 1A in block format, excluding the chemical processing system 114, along with expressions relating to describing optical exposure system simulation according to a preferred embodiment. Also shown in FIG. 12 are exposure system variations (“var”) in the form of projection system variations 403″. The optical exposure system is characterized by a TCC 1216 that, if used in a Hopkins formulation integral according to FIG. 12 (Eq. {12a} together with the mask transmittance function 312, would yield an output optical intensity 1218. Notably, according to a preferred embodiment, the integral of FIG. 12 (Eq. {12a}) is not itself required to be carried out. Although mask/target segmentations are not illustrated in FIG. 12, the computations described herein are preferably performed for each of a large number of mask segments and corresponding target segments similar to the mask segments 120 and target segments 122 of FIG. 1B, supra.

[0086] With reference to the steps of FIG. 10, supra, and FIG. 12 (Eq. {12b}) through FIG. 12 (Eq. {12h}), one particular method for accomplishing step 1002 is as follows. First, analytical and/or numerical expressions for the partial response functions 410 and partial response coefficients 408 are developed, wherein the partial response functions 410 are not dependent on the exposure system variations (var), and wherein the partial response coefficients 408 are dependent on the exposure system variations (var). These expressions could be readily developed by a person skilled in the art in view of the present disclosure, with one example for a defocus variation z′ being provided at FIG. 8 (Eq. {8h}), supra, and another example for a generalized variation provided at FIG. 17, infra. Second, analytical and/or numerical expressions for the windowed mutual intensity function 302 incident upon the mask 108 (see FIG. 3, supra) are generated. This can be achieved either using methods known in the art or, optionally, by the coherent mode method described in relation to FIG. 3, supra. Third, using the relationships shown at FIG. 12 (Eq. {12e}) and substituting counter variables as shown, a set of kmax TCC basis spatial functions 1219 not dependent on the optical variations (var) and TCC expansion coefficients 1220 dependent on the exposure system variations (var) are generated. The TCC 1216 can thereby be expressed as illustrated at FIG. 12 (Eq. {12g}) according to a series expansion having basis spatial functions 1219 not dependent on the exposure system variations and expansion coefficients 1220 dependent on the exposure system variations.

[0087] With reference to FIG. 12 (Eq. {12h}), step 1004 of FIG. 10 can be carried out by performing a Mercer expansion for each k^th TCC basis spatial function 1219 using known Mercer expansion methods. For each k^th TCC basis spatial function 1219, a set of Lmax optical kernels 1222 and coefficients 1224 are generated. Notably, the number Lmax of optical kernels 1222 and optical kernel weights 1224 can itself vary with k, i.e., there can be different numbers of terms in the Mercer expansions for different ones of the TCC basis spatial functions 1219.

[0088] With reference to FIG. 12 (Eq. {12j}), step 1006 is carried out by convolving the mask transmittance function 312 with each of the optical kernels 1222 to generate convolution results 1226. Step 1008 can then be carried out by computing a weighted sum of the squared magnitudes of the convolution results 1226 as weighted by a corresponding plurality of optical kernel weights 1228, the optical kernel weights 1228 being computed, in turn, from the TCC expansion coefficients 1220 and the coefficients 1224 associated with the Mercer expansion of the TCC basis spatial functions, for the initial projection system condition.

[0089] Notably, the convolution results 1226 are not dependent on the exposure system variations (var), whereas the optical kernel weights 1228 are dependent on the exposure system variations (var) by virtue of the TCC expansion coefficients 1220. Accordingly, for step 1010, target optical intensities can be computed for different values of the exposure system variation (var) simply by re-combining the same set of convolution results 1226 using the associated different optical kernel weights 1228. According to another preferred embodiment, steps similar to described with respect to FIGS. 10 and 12 are carried out in a manner that includes different target depths, separately or in various combinations (var,z_t) with the different process variation values.

[0090] FIG. 13 illustrates a conceptual view of the formation of the optical kernels 1222 referenced supra with respect to FIG. 12. At the left side of FIG. 13 is a (kmax+1)-member set of TCC basis spatial functions 1219, represented as cubes to denote that each is a function of four variables (x_m1,y_m1;x_m2,y_m2). For each k^th basis spatial function 1219, a set of optical kernels OPTKERN_k,1, OPTKERN_k,2, . . . OPTKERN_k,Lmax(k), which can be generated by Mercer expansion, for example, is illustrated. The number of optical kernels for each k^th TCC basis spatial function 1219 can vary with k, and therefore there are Lmax(k) elements for each k^th TCC basis spatial function 1219. Referring again to the example M′×M′ window of optical ambit scenario, the optical kernels OPTKERN_k,L would each be M′×M′ arrays of complex values. At the right side of FIG. 13 is the set of optical kernels 1222 reproduced as a group, which is precomputed for a particular optical exposure system. Importantly, in accordance with the features and advantages of at least one of the preferred embodiments, the optical kernels 1222 do not depend on the exposure system variations, and the convolution of a particular mask transmittance function with the optical kernels 1222 only needs to be performed once for computing resultant intensity values across many different exposure system variation values.

[0091] FIG. 14 illustrates a conceptual diagram of the use of a photolithographic process simulation engine 1402 according to a preferred embodiment. The photolithographic process simulation engine 1402 is used, for example, as part of a larger optical proximity correction (OPC) tool used by a semiconductor fabrication facility that is presented with a desired circuit pattern for a production run. The semiconductor foundry often possesses several photolithographic processing systems, each having their own particular variations and not processing mask patterns in precisely the same way. The characteristics for a typical photolithographic processing system are embodied in a set of model kernels and weights 1408 that have been developed or calibrated based on its particular physical parameters using one or more of the preferred embodiments described supra (and infra with respect to FIGS. 15A-15B). Based on the desired circuit pattern, the OPC tool has developed a particular mask pattern 1410, usually as part of a larger overall iterative process. As part of the larger overall iterative process, it is now desired to simulate how the current mask pattern 1410 will be treated by the particular photolithographic processing system in question, for the purposes of perfecting the mask pattern for producing the desired circuit pattern on the wafer using that photolithographic processing system.

[0092] The photolithographic process simulation engine 1402 comprises an optical exposure simulation engine 1404 and a chemical processing simulation engine 1406. Based on the model kernels and weights 1408 and the mask pattern 1410, the optical exposure simulation engine 1404 computes a target resist optical intensity pattern that would be generated by the particular stepper device in question, in accordance with one or more of the preferred embodiments described supra, for a mask pattern 1410. Based on the computed optical intensity pattern, the chemical processing simulation engine 1406 simulates the subsequent chemical processing of the exposed photoresist and the underlying semiconductor wafer so that a simulated image 1416 is generated.

[0093] Because a photolithographic processing system in question will inevitably incur process variations during a production run (for example, defocusing due to vertical placement errors of the wafer relative to the stepper device optics), and for other reasons, it may be desirable to compute simulated images for a variety of different process variations from an initial process condition. Also, to improve the accuracy of any particular simulated image, it may be desirable to compute resultant intensities for several different target depths. Advantageously, in accordance with a preferred embodiment, simulated images 1416 for a plurality of values for such variations 1414 (which includes both process variations and multiple resist depths) can be generated in an overall computation time that is not substantially greater than the computation time for a single variation value. For example, where computation for a single variation value requires a first time interval T, results for N different variation values can be achieved within a second time interval that is substantially less than NT.

[0094] FIG. 15A illustrates a conceptual diagram for the calibration of a photolithographic process simulation engine 1402, the photolithographic process simulation engine 1402 comprising an optical exposure simulation engine 1404 capable of simulating a photolithographic process for a plurality of different process and/or depth variations in an overall computation time that is not substantially greater than the computation time for a single process and/or depth condition in accordance with a preferred embodiment. As used herein, calibration refers to the generation of model kernels and weights customized for a photolithographic processing system for simulating the operation of that photolithographic processing system. For a photolithographic processing system, such calibration could be performed, for example, by the semiconductor fabrication facility upon receipt, modification, or repair of a stepper device, or at other times for other reasons as circumstances warrant.

[0095] Based on initial physical parameters 1502 for a subject physical photolithographic processing system 1509, as may be obtained from information sources including stepper device manufacturers, initial model kernels and weights are computed at block 1504 in accordance with one or more of the preferred embodiments. A mask pattern 1506 is provided to the simulation engine 1402, and is also physically written to a physical mask 1508. The physical mask 1508 is inserted into the subject photolithographic processing system 1509, and a physical result (actual image 1510) is generated using the photolithographic processing system. Meanwhile, a simulated image 1512 is generated by the simulation engine 1402 based on the test mask pattern 1506 and the initial model kernels and weights generated at block 1504. The actual image 1510 and simulated image 1512 are compared at block 1514. If the deviation between the actual image 1510 and the simulated image 1512 is within an acceptable tolerance, then the model kernels and weights 1518 are established for that photolithographic processing system for use in future simulations on actual proposed mask patterns for actual desired circuit designs. If the deviation is not acceptable, then at block 1516 a modification process is carried out that modifies the physical parameters and sends the modified physical parameters to block 1504 which, in turn, generates the modified model kernels and weights. The simulation engine 1402 then generates another simulated image based on the test mask pattern 1506 and the modified model kernels and weights, and the process is repeated until the simulated image 1512 is sufficiently similar to the actual image 1510.

[0096] FIG. 15B illustrates calibrating a first process model to a photolithographic processing system according to a preferred embodiment, the first process model corresponding to a first photolithographic process simulator and comprising a plurality of model kernels. An actual image (i.e., physical result) is generated by operation of the photolithographic processing system on a physical mask derived from a test mask pattern (blocks 1556, 1558, 1559, 1560). An initial parameter set 1552 at least partially characterizing the photolithographic processing system is received. An optimized parameter set is generated using the initial parameter set 1552, a second photolithographic process simulator 1590 different from the first photolithographic process simulator, the test mask pattern 1506, and the physical result 1560, using a comparison at block 1564 between computed resultant intensities 1562 and the physical result 1560 and making associated physical parameter modifications at block 1566. In one preferred embodiment, the second photolithographic process simulator 1590 uses an Abbe's formulation method for optical exposure simulation at block 1591. Finally, the plurality of model kernels along with associated model kernel weights is computed from the optimized physical parameters using, for example, one or more methods described hereinabove with respect to FIGS. 2-6A and FIGS. 7A-13 (blocks 1504 and 1570).

[0097] FIG. 16A illustrates the photolithographic processing system 102 of FIG. 1A in a block format and expressions relating to describing photolithographic process simulation according to a preferred embodiment, including process variations (var) 1614 associated with one or more of the illumination system, projection system, and chemical processing system. According to a preferred embodiment and as described further infra, a TCC 1616 that describes the photolithographic processing system and that also incorporates the process variations is characterized as a series expansion having basis spatial functions 1619 not dependent on the process variation and at least one expansion coefficient 1620 dependent on the process variations, even where it is not necessarily the projection system's point spread function that incorporates the process variations as with the preferred embodiment of FIG. 12, supra. Each TCC basis spatial function 1619 can then be Mercer-expanded to form imaging kernels 1622 and expansion coefficients 1624.

[0098] Thus, for a general process variation, a set of imaging kernels 1622 and imaging kernel weights 1628 are computed according to a preferred embodiment, the imaging kernel weights 1628 being computed, in turn, from the TCC expansion coefficients 1620 and the coefficients 1624 associated with the Mercer expansion of the TCC basis spatial functions, the imaging kernel weights 1628 being dependent on the process variations. For an initial process variation value, the resultant intensity 1618 can be computed by convolving the mask transmittance function 312 with each of the imaging kernels 1622 to generate a plurality of convolution results 1626, which are then combined according to the imaging kernel weights 1628 for the initial process variation value. For subsequent process variation values, resultant intensities 1618 can be computed simply by re-combining the same set of convolution results 1626 using the different values for imaging kernel weights 1628. According to another preferred embodiment, similar steps are carried out in a manner that includes different target depths, separately or in various combinations (var,z_t) with the different process variation values.

[0099] FIG. 16B illustrates a conceptual view of the formation of the imaging kernels 1622 referenced supra with respect to FIG. 16A. At the left side of FIG. 16B is a (kmax+1)-member set of TCC basis spatial functions 1619, represented as cubes to denote that each is a function of four variables (x_m,1,y_m1;x_m2,y_m2). For each k^th TCC basis spatial function 1619, a set of imaging kernels IMGKERN_l,1, IMGKERN_k,2, . . . IMGKERN_k,Lmax(k), which can be generated by Mercer expansion, for example, is illustrated. The number of imaging kernels for each k^th TCC basis spatial function 1619 can vary with k, and therefore there are Lmax(k) elements for each k^th TCC basis spatial function 1619. Referring again to the example M′×M′ window of optical ambit scenario, the imaging kernels IMGKERN_k,L would each be M′×M′ arrays of complex values. At the right side of FIG. 16B is the set of imaging kernels 1622 reproduced as a group, which is precomputed for a particular photolithographic processing system. Importantly, in accordance with the features and advantages of at least one of the preferred embodiments, the imaging kernels 1622 do not depend on the process variations, and the convolution of a particular mask transmittance function with the imaging kernels 1622 only needs to be performed once for computing resultant intensity values across many different process variation values and target depths.

[0100] FIG. 17 illustrates a portion of the photolithographic processing system 102 of FIG. 1A and expressions relating to describing photolithographic process simulation according to a preferred embodiment, for the purposes of illustrating another example, in addition to the example supra of FIG. 8, of the manner in which a point spread function incorporating at least one process variation can be provided. In contrast to the example supra of FIG. 8, FIG. 17 illustrates a generalized example in which the process variation can be expressed as a generalized frequency plane perturbation 1710 times the model pupil function 804. As shown in FIG. 17 (Eq. {17b}), the point spread function 1708 can be expressed as an inverse Fourier transform of the product of the frequency plane perturbation 1710 and the model pupil function 804. As illustrated in FIG. 17 (Eqs. {17c}-{17d}), the generalized frequency plane perturbation 1710 can be expressed as a summation of, for example, Zernike polynomials 1714 as weighted by expansion coefficients 1712, wherein the Zernike polynomials 1714 depend only on spatial frequencies, and wherein the expansion coefficients 1712 depend only on the process variations (var) and/or the target depth z_t. This results in the desired expression of FIG. 17 (Eq. {17g}) in which the point spread function 1708 is expressed as a series expansion having basis spatial functions 1716 that do not depend on the process variation or target depth and expansion coefficients 1712 dependent on the process variation and/or target depth. The basis spatial functions 1716 can be computed as the inverse Fourier transforms of products of the Zernike polynomials 1714 and the model pupil function 804. The expansion coefficients 1712 can be computed using the relationship shown in FIG. 17 (Eq. {17d}).

[0101] For one or more of the preferred embodiments described supra, features and advantages presented with respect to process variations can also be applied in the context of multiple target depths. For example, according to one preferred embodiment, simulation of a photolithographic process is provided in which computation time for computing the results for several different target depths is not significantly greater than the computation time for computing the results for a single target depth. Computation of model kernels and model kernel weights proceeds in accordance with one or more of the preferred embodiments supra. A mask transmittance function is convolved with each model kernel to generate a plurality of convolution results. A first resultant intensity for a first target depth is computed as a first combination of the convolution results, e.g., combined using a first set of model weight values. Resultant intensities are subsequently computed for different target depths as subsequent combinations of those convolution results, e.g., combined using subsequent sets of model weight values that numerically depend on the target depth.

[0102] According to another preferred embodiment, simulation of an exposure process is provided in which computation time for computing the results for several process variation values for each of a plurality of different process variations, as well as for several different resist planes, is not significantly greater than the computation time for computing the results for a single value for each of the process variations and a single target depth. Computation of model kernels and model kernel weights proceeds in accordance with one or more of the preferred embodiments supra. A mask transmittance function is convolved with each model kernel to generate a plurality of convolution results. A first resultant intensity for an initial value for each of the process conditions, as well as an initial target depth, is computed as a first combination of the convolution results, e.g., combined using a first set of model weight values. Resultant intensities are subsequently computed for subsequent values of each of the process variations and for different target depths as subsequent combinations of those convolution results, e.g., combined using subsequent sets of model weight values that numerically depend on the process variation values and the target depth.

[0103] Whereas many alterations and modifications of the present invention will no doubt become apparent to a person of ordinary skill in the art after having read the descriptions herein, it is to be understood that the particular preferred embodiments shown and described by way of illustration are in no way intended to be considered limiting. By way of example, references to convolution operations hereinabove also encompass mathematical operations including, but not limited to, transformation into the frequency domain, multiplication, and inverse transformation from the frequency domain, designed to achieve results similar to spatial domain convolution methods. Therefore, reference to the details of the preferred embodiments are not intended to limit their scope, which is limited only by the scope of the claims set forth below.

[0104] FIG. 20 illustrates a schematic diagram of a computer system capable of implementing the methods of the present invention. Methods of simulation in accordance with the present invention may be implemented in machine readable code (i.e. software or computer program product) and performed on a computer, such as, but not limited to, the type illustrated in FIG. 20, wherein a central processing unit (CPU) 2000 is connected by means of a bus, network or other communications means 2010 to storage devices such as disk(s) 2001, or memory 2002, as well as connected to various input/output devices, such as a display monitor 2004, a mouse 2003, a keyboard 2005, and other I/O devices 2006 capable of reading and writing data and instructions from machine readable media 2008 such as tape or compact disk (CD). In addition, there may connections via a network 2010 to other computers and devices, such as may exist on a network of such devices, e.g. the Internet 2009. Software to perform the simulation steps described herein may be implemented as a program product and stored on machine readable media 2008.

[0105] General Discussion

[0106] The feature size of integrated circuits (ICs) has been kept shrinking at a steady pace by technology advancements in photolithography and mask design, so to sustain the ongoing growth of circuit complexity and running speed of semiconductor chips. In the so-called deep-sub-micron regime, the features on wafers are so tiny (less than 100 nm) and so dense that it becomes increasingly difficult to create wafter patterns resembling the design patterns (created by IC design engineers/tools) with high fidelity. In a traditional IC design flow, design rule checking tools are used to check the design patterns against a set of design rules, and extraction tools are employed to estimate electric parameters of to-be-fabricated devices directly from the design patterns, the electric parameters are then used to perform timing analysis, signal integrity check, power and thermal analyses for the IC under design. In view of the inevitable difference between the design and wafer patterns, it makes more sense to employ a computer model simulating the manufacturing process, and to predict the wafer patterns using computer simulations, then to check the predicted wafer patterns against the design rules, and to extract electric parameters of devices based on the predicted wafer patterns.

[0107] Generally speaking, a currently available computer simulation process model can be considered as falling into one of two classes, as discussed in C. A. Mack, “30 years of lithography simulation,” Proc. SPIE, vol. 5754, pp. 1-12, 2005. One class is physically based, simulating the processes of optical imaging, resist development, and chemical etching in what can be called an ab initio (“from the beginning”) manner, or from first-principles. Such process models are fairly accurate but rather computation-intensive because they simulate the physics in a mechanical fashion, hence can often overload themselves with data manipulation and number crunching without the help of analytical or semi-analytical solutions. As a result, these process models are generally too slow to run full-chip simulations for most production applications, such as optical proximity correction (OPC) of masks and other resolution enhancement techniques (RETs), RET verification, etc. However, they are useful research tools to run simulations on a small chip area on the order of microns to tens of microns squared, to assist process development and to benchmark other simulation tools. The other class of process models is used in production tools, that is designed with high running speed and sacrificed simulation accuracy. A reason for such speed improvement and accuracy loss is that the models are not all ab initio, rather, they often are based upon or incorporate empirical assumptions about the processes of optical imaging, resist development, and chemical etching, that are not general principles governing the underlying physical and chemical processes. Apart from lower accuracy, the applicability of the process models for production tools is usually limited to a small region of the parameter space of process variations, because the models are calibrated using test data under a single, usually the best, process condition, and efforts are made to add and tweak empirical formulas and parameters until the best fit to the measured data is reached. A model built in this manner may reasonably predict the process behavior under the calibrated condition. However, the empirical assumptions often lack the generality to make reliable predictions when simulating process deviations away from the calibrated condition.

[0108] Consequently, existing OPC and RET verification products are often only capable of modelling a manufacturing process and simulating the wafer patterns under one (usually the best) process condition. For example, the simulations usually assume a process under the best focus of the optical exposure system (often referred as the optical stepper) and ideal intensity/position settings of multiple emitters. However, practical manufacturing inevitably involves process variations. Defocus, that is deviation of the wafer position from the best focal surface, and intensity/position variations of multiple emitters, are exemplary sources of process variations among others. Such process variations could degrade the quality of photolithographic pattern transfer, resulting in wafter patterns that fall out of certain design rules, which eventually lead to yield loss. Accordingly, there is a need for OPC and RET verification tools to be able to accurately model process variations and to sample over a large space of such variations by tuning continuously at least two process parameters, then to predict wafer patterns under variations of process parameters, eventually to calculate or optimize the so-called process windows. A hardware system has been available for full-chip process window qualifications, as discussed in M. J. Brodsky et al., “Process-window sensitive full-chip inspection for design-to-silicon optimization in the sub-wavelength era,” Proc. SPIE, vol. 5756, pp. 51-60, 2005. A software counterpart for a similar purpose could be equally useful, more cost-effective, more convenient to use, while offering more and readily expandable functionalities with faster turnaround. It would be desirable to provide a numerical algorithm that provides both accuracy and speed simultaneously to run full-chip simulations over a large space of process variations. It would be further desirable to provide a numerical algorithm that supports continuous sampling of at least two process variations. One approach that has actually been tried in several commercial products is to employ multiple empirical models calibrated separately at discrete process conditions, as discussed in J. Schacht et al., “Calibration of OPC models for multiple focus conditions,” Proc. SPIE, vol. 5377, pp. 691-702, 2004; J. A. Torres et al., “Process window modeling using compact models,” Proc. SPIE, vol. 5567, pp. 638-648, 2004; J. Word et al., “Lithography yield enhancement through optical rule checking,” Proc. SPIE, vol. 5645, pp. 142-153, 2005; N. Lafferty et al., “Physically-based compact models for fast lithography simulation,” Proc. SPIE, vol. 5754, pp. 537-542, 2005; A. Borjon et al., “High accuracy 65 nm OPC verification: full process window model vs. critical failure ORC,” Proc. SPIE, vol. 5754, pp. 1190-1201, 2005; M. Cote et al., “Process centering OPC using design intent to improve yield,” Proc. SPIE, vol. 5756, pp. 331-339, 2005; and P. Martin et al., “Exploring new high speed, mask aware RET verification flows,” Proc. SPIE, vol. 5853, pp. 114-123, 2005. However, a collection of discretely calibrated models may not support continuous tuning of process parameters because it becomes questionable when extrapolating empirical models into conditions other than where they are calibrated. Furthermore, multiple models generally require significantly more effort to calibrate, and may induce higher computational complexity to applications using them.

[0109] For the known existing software approaches identified above, one disadvantage is that, if N is the number of discrete sampling points of process variations, for example the number of different defocus values, or different values of lens abberations, and combinations of values of multiple process variations, then those known existing software approaches need to generate and use N sets of process models, each process model comprising a set of model kernels, which collectively characterize the photolithographic process from a mask to a resultant image of the mask (e.g., an aerial image, or a planar image at a certain depth in resist, or a three-dimensional (3D) distribution of optical intensity in resist, or a latent image in resist, or a set of developed resist patterns, or a set of final silicon patterns). For each discrete sampling point of process variation, each model kernel may be convolved with a mask transmittance function and the convolution results may be squared and combined in a weighted manner to yield the resultant image of the mask. The resultant image of the mask may be used to make predictions about the developed resist patterns or the final silicon patterns. Such software approaches generally suffer from stiff increases of computational complexity when the number N of discrete sampling points of process variations becomes large. The computational complexity of a computational task may be measured by the total computation time of all computers used for the task, so that the complexity measure is substantially independent of the number of computers used, so long as all computers are substantially equivalent in computational power. When the same or similar computational task is repeated under different values of process variations, excluding the computation time needed for initial preparations, such as reading files and preparing mask data, let T(1) denote the computational complexity for simulating one sampling point of process variations, and T(N), N≧2 be the total computational complexity for simulating N sampling points of process variations, then the following formula holds true for such software approaches,
T(N)=T(1)×N, N≧2.  (1)
That is, the computational complexity is directionally proportional to the number of sampling points of process variations, when the same or similar computational task is repeated. By repeating the same or similar computational task under process variations, it is meant that the computational task is repeated with the same or similar specifications including but not limited to input data and desired output data, except for changing the values of a set of process variations. For an example scenario, the distribution of resultant intensity within a chosen area on a wafer substrate may be calculated repeatedly under different values of defocus, where the fixed output data may be discretely sampled intensity values on a fixed mesh of target points over the chosen area on wafer, whereas the fixed input data may include the mask layout, especially the portion that corresponds to the chosen area on wafer, parameters related to chemical processing, and the configurations of the stepper machine such as the illumination setting and the numerical aperture of the projection lens, except for the distance from the last projection lens to the wafer substrate, namely the defocus value, which changes from one repetition to another. A known existing software tool tackles such simulation task for each value of defocus separately and independently, thus incurs a complexity cost that is directly proportional to the number of repetitions, that is, the number of different values of defocus. Such high computational complexity imposes a substantial limitation to such software approaches in simulating photolithographic processes under variations.

[0110] Another scenario may be that different mask layouts or different portions of the same mask layout are input for different sampling points of process variations, namely, repetitions of a computational task with different values of a set of selected process variations, while other (non-selected) parameters of stepper configurations and chemical processing remain the same. The output data may no longer be fixed under such scenario. Rather, each repetition may output resultant intensities of different target areas, corresponding to different portions of a mask layout or simply different mask layouts. Furthermore, the resultant intensity may be irregularly and sparsely chosen by a user or application software, and the irregular distribution of chosen target points may vary from one repetition, namely, one sampling point of process variations, to another. Consequently, the total computation time may vary a lot from one to another sampling point of process variations. However, it is often the case that the total computation time divided by the total number of target points, called the average cost of computational complexity per target point, or cost per point (CPP) in short, remains approximately the same from one sampling point of process variations to another. Therefore, the CPP may serve as another suitable measure for the efficiency or computational complexity of an algorithm or software for photolithographic simulations under process variations, especially with irregularly and sparsely distributed target points. The CPP is generally dependent on the method for computing convolutions of model kernels and a mask transmittance function. One method uses fast Fourier transform (FFT) and computes the convolutions as simple function multiplications in the spatial frequency domain, which is suitable for calculating resultant intensities at regularly and densely sampled target points, and achieves a rather low CPP by virtue of dividing the total computation time by a large number of target points. However, when the target points are irregularly and less densely distributed, such FFT-based method incurs an increase in CPP, because convolution results may have to be evaluated at regularly and densely distributed target points anyway, even though the calculation of resultant intensities is required for a much less number of target points that are irregularly and sparsely distributed. Another method is polygon-based, suitable for the case of irregularly distributed target points with low density, which has amplitude distributions pre-computed and tablized for a predetermined finite set of “primitive” mask shapes. During a photolithographic simulation with a mask layout, for each individual target point, such polygon-based method firstly lists the polygons of the mask layout that falls into a so-called “window of optical ambit” centered at the target point, and decomposes the polygons into primitive mask shapes, then obtains amplitude values for the target point with respect to the primitive mask shapes by looking up the pre-computed tables, and coherently sums the amplitude values associated with each model kernel, finally squares the coherent sums and sums the squares together to get a resultant intensity. The CPP with such polygon-based method is independent of the density and distribution of target points, rather, the CPP is proportional to the average number of primitive mask shapes within a window of optical ambit. For a known existing software of photolithographic simulation under process variations, let CPP be the cost per point regardless of which method is used for convolutions, and let Σ be the average number of target points per repetition, namely, per sampling point of process variations, then the total computation time T(N) for the known existing software to simulate N≧2 sampling point of process variations is given by,
T(N)=CPP×Σ×N, N≧2.  (2)
Once again, T(N) is roughly proportional to the number N of sampling points of process variations, which represents a high computational cost for known existing softwares to simulate the photolithographic effects of process variations.

[0111] Provided in accordance with one or more of the preferred embodiments is a fast ab initio simulation engine that calculates the 3D image in a resist, and simulates the processes of resist development and chemical etching, all from first-principles. The simulation engine is highly accurate yet sufficiently fast to perform full-chip simulations within standard turn-around-time (on the order of hours) of production tools. Incorporated into one or more of the preferred embodiments is a finding that many computational problems in optical imaging, resist development, and chemical etching are actually amenable to analytical or semi-analytical solutions, which help to greatly reduce the amount of necessary numerical calculations. In particular, the defocus behavior and the 3D image distribution in the resist of a partially coherent imaging system may be represented semi-analytically in terms of series expansions, where a typical term comprises one function of spatial coordinate(s) that does not depend on the defocus distance and/or depth into the resist film and another function that depends on the defocus distance and/or depth into the resist film. The functions of spatial coordinate(s) need only to be numerically calculated once, then the 3D image distribution under variable defocus may be more readily computed by series summation. The provided algorithm evaluates 3D image under arbitrary defocus (as well as abberations and other process variations) without incurring a computational complexity much higher than a conventional calculation of two-dimensional (2D) image under a fixed defocus. Also incorporated into one or more of the preferred embodiments is a finding that 3D image calculation is the major contributor of computational cost for a process simulation. Once a 3D latent image is available, the subsequent ab initio simulations of resist development and chemical etching induce a generally insignificant increase of computational cost. Also provided in accordance with one or more of the preferred embodiments is an efficient and common tree structure for mask data as well as image fields and wafer patterns, which increases data locality and tightly correlates numerical operations with their relevant data, and reduces data dependency among computations on different chip areas. Such data structure is well suited for distributed processing. Also provided in accordance with one or more of the preferred embodiments is efficient data processing by the combination of FFT and tree codes. Resultant amplitudes are calculated once by convolving mask data with pre-computed model kernels using FFT, and results are stored on regular grid points on a tree structure of image points on wafer. Such grid points storing image data are on a higher hierarchical level of the tree structure, where points on the lower hierarchical levels including the leaf nodes have no image information yet. In a latter stage, when image values are needed for a selected set of points on lower hierarchical levels of the tree, image data are pushed down the hierarchy using tree codes.

[0112] Provided in accordance with one or more of the preferred embodiments are algorithms that enable fast, first-principle, and full-chip (F³) simulations of photolithographic processes. It is believed that F³ engine will benefit the following applications as well as others: a) yield-oriented silicon fidelity checking (SiFiC), which computes the common process window of a chip or cell by sampling over a space of variable process parameters and checking predicted wafer patterns against a set of design rules; when a probability distribution of process variations is available, such yield-oriented SiFiC may then estimate a yield probability by integrating the probability density over the calculated process window; b) yield-oriented optical proximity correction, which uses a yield-oriented SiFiC-like routine to estimate the process window and yield probability, and selects mask corrections to optimize the process window and yield; c) chip-specific optimization of stepper machines and photolithographic processes, which runs a yield-oriented SiFiC-like routine over the whole chip of a specific chip design, repeatedly under changing stepper configurations and process parameters, until the best configurations and parameters are identified which optimize the predicted yield or process windows; d) virtual device and chip testing based on ab initio simulations, which predicts wafer patterns and extracts electric parameters for devices, then uses a circuit simulator to test the functionality and performance of the IC under design; e) benchmark process models embedded in physical design tools, which are fast simulators of manufacturing processes to guide design optimizations for design for manufacturing (DFM) and design for yield (DFY).

[0113] Mathematical Theory of Optical Imaging

[0114] Optical Setup of a Stepper

[0115] A stepper machine has a rather complicated optical system, which however, may be modelled as comprising of a spatially extensive illumination source placed on the focal plane of a condenser lens, a mask being illuminated by the condensed light, and a projection lens forming an M times reduced image of the mask onto a wafer, as shown by the left picture of FIG. 18. Many systems are with the reduction factor M=4. Even though a real projection system may have tens of lens pieces, it may be well represented by a two-lens setup, with one lens having a focal length F and the other MF, M≧1. The mask is placed at the input focal plane of the MF lens, the output focal plane of which coincides with the input focal plane of the F lens, where the Fourier transform of the mask transmittance function is formed. The semiconductor wafer is placed at the output focal plane of the F lens to receive the image of the mask. Some advanced designs have the last projection lens and the wafer immersed in a liquid, such as water, with an index of refraction η₀>1. No matter how the imaging system may be designed, it is usually the numerical aperture N.A.=η₀ sin θ_max of the last projection lens that limits the resolution of the system, where θ_max is the largest angle made by an optical ray with respect to the z-axis, which comes off the last projection lens and propagates toward the semiconductor wafer. The z-axis is the optical axis of the lens system, which is normal to the mask and wafer planes. The projection system from the mask down may be modelled by an equivalent 4F and 1:1 system as depicted by the right picture of FIG. 18, where the mask is reduced by M times from its actual size, a pupil may be placed on the Fourier plane to limit the numerical aperture, and the whole system, including the mask, is immersed in the liquid with index η₀. We stress that such completely immersed, 4F and 1:1 projection system is not the actual configuration of any real machine, rather it serves as a good model for real-world stepper systems, as far as the physics of optical imaging is concerned.

[0116] The condenser transforms the light of each source point indexed by (x,y)=ρ(cos α, sin α), with ρ=(x²+y²)^1/2 and α∈[0,2π), into a coherent plane wave √{square root over (I(x,y))}exp(i2πf′·r), with the spatial frequency f′^def(f′_z,f′_y,f′_z)=λ⁻¹(cos θ,sin θ cos α,sin θ sin α), where I(x,y) denotes the intensity of the source point, r^def(x,y,z) is the spatial coordinate, λ is the light wavelength in vacuum, θ=arctan(ρ/F_c) is the angle between the z-axis and the propagation direction of the plane wave, F_c is the focal length of the condenser lens, and

[0117] $f^{\prime }·r\stackrel{\mathrm{def}}{=}\left(f_x^{\prime }x+{f_y}^{\prime }y+f_z^{\prime }z\right)$
is the scalar product between f′ and r. Because of the one-to-one correspondence between the location of a source point and the spatial frequency of the plane wave that it generates after the conderser lens, it is possible (and industrial custom) to index the source points by the spatial frequency coordinate. It is further possible to place the frequency-indexed source pattern on the Fourier plane of the projection system as shown in the picture on the right side of FIG. 18, after scaling up (f′_z,f′_y) of the source “coordinate” to (f_z,f_y)=(Mf′_z,Mf′_y), because the mask has been sized down by M times. Then it makes sense to talk about the relative positions and sizes of the source and pupil aperture of the projection system, all in a bounded region on the Fourier plane. The bound in the frequency domain reflects the bandlimiting nature of the system. In particular,

[0118] ${\sigma }_{\max }\stackrel{\mathrm{def}}{=}\max \left\{{\lambda \left(f_x^2+f_y^2\right)}^{1/2}\right\}$
corresponds to the outmost source point from the optical axis, which generates the most off-axis illuminating beam, making an angle sin⁻¹(σ_max/M) with respect to the optical axis.

[0119] General Theory of Partially Coherent Imaging

[0120] There are two mathematical formulations commonly used for partially coherent imaging. One is Abbe's formulation, the other is Hopkins' formulation, as discussed generally in M. D. Smith and C. A. Mack, “Methods for benchmarking photolithography simulators,” Proc. SPIE, vol. 5040, pp. 57-68, 2003; and R. J. Socha, “Resolution enhancement techniques,” in B. G. Eynon, Jr. and B. Wu, Photomask Fabrication Technology, McGraw-Hill, New York, 2005, at Chapter 7. In Abbe's formulation, a spatially extensive illumination source is divided into non-overlapping pixels that are naturally incoherent, where each pixel is sufficiently small to be treated as a single coherent point source. Through the condenser lens, each pixel as a point source generates a parallel beam, or called a plane wave, to illuminate the mask, with the incident angle of the beam determined by the position of the source pixel. Such plane wave multiplied by the mask transmittance function generates an input amplitude field to the projection system, which then forms an output amplitude distribution on the target, whose magnitude squared becomes a contributory resultant intensity distribution for the mask image due to the source pixel. Since different source pixels are mutually incoherent, all contributory intensity distributions due to individual source pixels are summed together to yield the total intensity distribution for the mask image on the target. By contrast, Hopkins' formulation does not divide the illumination source into pixels, instead it uses a Hilbert-Schmidt kernel function called the transmission cross coefficient (TCC) to incorporate the effects of both the illumination and the projection systems together. See, generally, H. H. Hopkins, “On the diffraction theory of optical images,” Proc. Roy. Soc. Lond. A, vol. 217, pp. 408-432, 1953; and M. Born and E. Wolf, Principles of Optics, 7th ed., Cambridge University Press, 1999, at Chapter 10.

[0121] To establish a general framework with Hopkins' formulation, we start from the mutual intensity of two points r₁ and r₂ on a surface S₁, which is illuminated by light coming from a surface S₀ of incoherent source points,

[0122] $\begin{array}{cc}{J}_1\left(r_1,r_2\right)\stackrel{\mathrm{def}}{=}\langle E_1^{*}\left(r_1,t\right)E_1\left(r_2,t\right)\rangle ={\int }_{S_0}{K}_{01}^{*}\left(s,r_1\right)K_{01}\left(s,r_2\right)I_0\left(s\right)ds,& \left(3\right)\end{array}$
where E₁(r,t) is the light electric field at a point r=(x,y,z), denotes time and statistical average, I₀(s) is the intensity per unit area of the source at a point s∈S₀, K₀₁(s,r) is the optical propagation amplitude from the source point s to the target point r. See, generally, as referenced supra, H. H. Hopkins; and M. Born and E. Wolf. The symbol “*” denotes the complex conjugate for a scaler quantity, and more generally the Hermitian conjugate for a vector, a matrix, or an operator. The following describes the propagation of mutual intensity from one surface S₁ to another surface S₂,
J₂(q₁,q₂)=∫_s1∫_s1K*₁₂(r₁,q₁)J₁(r₁,r₂)K₁₂(r₂,q₂)dr₁dr₂,  (4)
where q₁, q₂∈S₂, and K₁₂(r,q) is the optical propagation amplitude from r∈S₁ to q∈S₂.

[0123] Now let S₁ be the plane where a mask is located, and S₂ be the image surface (that is the wafer surface for lithography). Equation (3) may be used to calculate the mutual intensity immediately before the mask, which is denoted as J₁⁻(r₁,r₂), r₁,r₂∈S₁. The mask is assumed to multiply the input light field E₁(r,t), r∈S₁ by a complex transmittance function u(r), r∈S₁. So the mutual intensity immediately after the mask is given by,

[0124] $\begin{array}{cc}\begin{array}{c}{J}_1^{+}\left(r_1,r_2\right)=\langle u^{*}\left(r_1\right)E_1^{*}\left(r_1,t\right)E_1\left(r_2,t\right)u\left(r_2\right)\rangle \\ =u^{*}\left(r_1\right)u\left(r_2\right){\int }_{S_0}{K}_{01}^{*}\left(s,r_1\right)K_{01}\left(s,r_2\right)I_0\left(s\right)ds\\ =u^{*}\left(r_1\right)J_1^{-}\left(r_1,r_2\right)u\left(r_2\right).\end{array}& \left(5\right)\end{array}$
According to equation (4), the mutual intensity for two points q₁,q₂∈S₂ would be,
J₂(q₁,q₂)=∫_s1∫_s1u*(r₁)K*₁₂(r₁,q₁)J₁⁻(r₁,r₂)K₁₂(r₂,q₂)u(r₂)dr₁dr₂.  (6)
Instead of the general mutual intensity, we are actually interested in the resultant intensity distribution I₂(q)=J₂(q,q)=E*₂(q,t)E₂(q,t), for which the formula is slightly simpler,
I₂(q)=∫_s1∫_s1u*(r₁)K*₁₂(r₁,q)J₁⁻(r₁,r₂)K₁₂(r₂,q)u(r₂)dr₁dr₂.  (7)
Imaging systems under our consideration are usually made shift-invariant in space, namely,
J₁⁻(r₁,r₂)≡J₁⁻(r₁−r,r₂−r), ∀r₁,r₂r∈S₁,  (8)
K₁₂(r,q)≡K₁₂(r−q|S₁,0), ∀r∈S₁, ∀q∈S₂,  (9)
where q|S₁ is the point on S₁ that is optically conjugate to q, namely, q|S₁ is the source point on S₁ whose geometrical optics image on S₂ is q. For plane-to-plane and 1:1 imaging, q|S₁^def(x,y,z₁)∈S₁, ∀q=(x,y,z₂)∈S₂. With such shift invariance, equation (7) may be rewritten as,

[0125] $\begin{array}{cc}{I}_2\left(q\right)={\int }_{S_1}{\int }_{S_1}{u}^{*}\left(r_1\right)K_{12}^{*}\left(r_1-q|S_1,0\right)J_1^{-}\left(r_1-q|S_1,r_2-q|S_I\right)\times \hspace{1em}{K}_{12}\left(r_2-q|S_1,0\right)u\left(r_2\right)d{r}_{1}d{r}_2.& \left(10\right)\end{array}$
With substitutions r₁→r₁+q|S₁, r₂→r₂+q|S₁, equation (10) may be further rearranged into,

[0126] $\begin{array}{cc}\begin{array}{c}{I}_2\left(q\right)={\int }_{S_1}{\int }_{S_1}{u}^{*}\left(r_1+q|S_1\right)K_{12}^{*}\left(r_1,0\right)J_1^{-}\left(r_1,r_2\right)K_{12}\left(r_2,0\right)u(r_2+\\ q|S_1)d{r}_{1}d{r}_2\\ ={\int }_{S_1}{\int }_{S_1}{u}^{*}\left(r_1+q|S_1\right)T\left(r_1,r_2\right)u\left(r_2+q|S_1\right)d{r}_{1}d{r}_2\\ ={\int }_{S_1}{\int }_{S_1}{u}_{q|S_1}^{*}\left(r_1\right)T\left(r_1,r_2\right)u_{q|S_1}\left(r_2\right)d{r}_{1}d{r}_2,\end{array}& \left(11\right)\end{array}$
with T(r₁,r₂) being the so-called TCC and represented as,

[0127] $\begin{array}{cc}\begin{array}{c}T\left(r_1,r_2\right)\stackrel{\mathrm{def}}{=}{K}_{12}^{*}\left(r_1,0\right)J_1^{-}\left(r_1,r_2\right)K_{12}\left(r_2,0\right)\\ ={\int }_{S_0}{K}_{12}^{*}\left(r_1,0\right)K_{01}^{*}\left(s,r_1\right)K_{01}\left(s,r_2\right)K_{12}\left(r_2,0\right)I_0\left(s\right)ds,\end{array}& \left(12\right)\end{array}$
and u_q(r) being the transmittance function of the same mask with shifted coordinate,

[0128] $\begin{array}{cc}{u}_q\left(r\right)\stackrel{\mathrm{def}}{=}u\left(r+q\right),\forall r,q\in S_1.& \left(13\right)\end{array}$

[0129] To simplify mathematical notations, we shall adopt the scale-normalized (Seidel) coordinate system and consider an equivalent 1:1 imaging system. We shall index points on any plane that is perpendicular to the z-axis by the same set P={(x,y):x∈R, y∈R}. P is actually a projective plane, in which an element q= (x,y,0) is the equivalent class of all 3D points (x,y,z) with arbitrary z. For convenience in dealing with mathematical convergence, P is further assumed to be bounded, that is,

[0130] $P\stackrel{\mathrm{def}}{=}\left\{\left(x,y\right):x\le L,y\le L\right\},$
where L>0 is sufficiently large so that the imaging system is not affected by the finiteness of the plane. A typical element q∈P is herein referred to as a planar spatial coordinate. Both surfaces S₁, S₂, and more that may be introduced later, are equivalent co-sets to P, and shall be identified with P as a convention. Mathematical formulas may be simplified by such collapse of parallel surfaces into one projective plane and indexing all points by the same projective planar coordinate. However, when the actual z coordinate is needed, the information may be obtained from the context. In particular for our propagation amplitudes K_dji(·,·), i,j=0,1, . . . , the subscripts i and j imply the propagation of light from one surface at z=z_i to another at z=z_j. Under said convention, equation (11) may be rewritten as,
I₂(q)=∫_P∫_Pu*_q(r₁)T(r₁,r₂)u_q(r₂)dr₁dr₂, ∀q∈P,  (14)
and the TCC may be calculated as,
T(r₁,r₂)=K*₁₂(0,−r₁)J₁⁻(r₁,r₂)K₁₂(0,−r₂),  (15)
where the propagation amplitude K₁₂(0,r), r∈P is actually the impulse response, also known as the point-spread function (PSF), of the optical projection system. It is to be appreciated that in one or more of the preferred embodiments, the PSF could be rendered (i.e., generalized) to incorporate some effects of photochemistry of resist molecules and post-exposure chemical processing and thereby characterizes the response of both the projection and chemical processing systems. By way of example and not by way of limitation, the PSF may be associated with the formation of any of the following resultant images: an aerial image, an optical intensity distribution in resist, a latent image in resist, a developed resist pattern, a final silicon pattern, and other responses. The PSF may be 2D, 3D, four- and even higher-dimensional where a dependence on time, defocus and/or other parameters is to be considered. For clarity of explanation in the present example, the PSF is considered as to form an optical intensity distribution (i.e., the chemical processing step is a null step) at one depth in the resist layer, hence to be 2D.

[0131] Formulations Using Functional Analysis

[0132] Using the language of functional analysis, both the amplitude of light on a surface and the transmittance function of a mask may be regarded as elements in the Hilbert space L₂(P), that is the set of square Lebesgue-integrable functions on P with an inner product ·∥· represented as,

[0133] $\begin{array}{cc}\langle u||v\rangle \stackrel{\mathrm{def}}{=}{\int }_P{u}^{*}\left(r\right)v\left(r\right)dr,\forall u,v\in L^2\left(P\right),& \left(16\right)\end{array}$
where the so-called “bra” and “ket” notation, as widely used in quantum mechanics, has been adopted. Namely, the ket symbol |u denotes a vector u in the Hilbert space, whereas the bra symbol u| denotes the vector conjugate to u in the dual space. We will also be concerned with linear operators, or called linear maps, from L²(P) to itself. For any linear operator O:L²(P)→L²(P), there is an adjoint operator, denoted by O*, that acts on the dual Hilbert space, in the bra and ket notation,

[0134] $\begin{array}{cc}|u\rangle \stackrel{\mathrm{def}}{=}u\rangle ,\langle u\stackrel{\mathrm{def}}{=}{*}^{}\langle u\stackrel{\mathrm{def}}{=}\langle \left(u\right),\forall u\in L^2\left(P\right).& \left(17\right)\end{array}$
Regarded as Hilbert-Schmidt kernels, the TCC T(r₁,r₂) and propagation amplitudes K_ij(r,q), i,j=0,1, . . . may be associated with linear operators, called the TCC operator and propagators, denoted by T and K_ij, i,j=0,1, . . . respectively, as in the following,

[0135] $\begin{array}{cc}𝒯❘u\rangle \left(r_1\right)\stackrel{\mathrm{def}}{=}{\int }_{P}T\left(r_1,r_2\right)u\left(r_2\right)d{r}_2,& \left(18\right)\\ {𝒦}_{\mathrm{ij}}u\rangle \left(q\right)\stackrel{\mathrm{def}}{=}{\int }_{P}u\left(r\right)K_{\mathrm{ij}}\left(r,q\right)dr,& \left(19\right)\end{array}$
∀u∈L²(P), ∀i,j=0,1, . . . . The TCC operator T is Hermitian, and the propagators K_ij are unitary operators, namely, T*=T, K*_ij=K_ij⁻¹, where K_ij⁻¹ denotes the inverse operator of K_ij. In terms of the associated Hilbert-Schmidt kernels, we have the following,
T(r₁,r₂)=T*(r₂,r₁), ∀r₁,r₂∈P,  (20)
K_ij⁻¹(q,r)=K*_ij(r,q), ∀r,q∈P.  (21)
Moreover, the TCC operator T is positive definite, as equation (14) is a positive definite quadratic form of elements u_q∈L²(P), or in the abstract notation,
I₂(q)=u_q|T|u_q=u_q∥Tu_q=Tu_q∥u_q≧0, ∀q∈P.  (22)
Being compact, Hermitian, and positive definite, the TCC operator T has positive eigenvalues λ₁≧λ₂λ≧ . . . >0 with corresponding eigenvectors φ_n∈L²(P), n=1,2, . . . . Moreover, T may be represented as a sum of projection operators,

[0136] $\begin{array}{cc}=\sum _{n\ge 1}{\lambda }_n{\varphi }_n\rangle \langle {\varphi }_n,& \left(23\right)\end{array}$
and equivalently, there is the Mercer expansion for the Hilbert-Schmidt kernel,

[0137] $\begin{array}{cc}T\left(r_1,r_2\right)=\langle r_1r_2\rangle =\sum _{n\ge 1}{\lambda }_n{\varphi }_n\left(r_1\right){\varphi }_n^{*}\left(r_2\right),& \left(24\right)\end{array}$
with which, equation (22) becomes,

[0138] $\begin{array}{cc}{I}_2\left(q\right)=\sum _{n\ge 1}{\lambda }_n\langle u_q||{\varphi }_n\rangle \langle {\varphi }_n||u_q\rangle =\sum _{n\ge 1}{\lambda }_n{\langle {\varphi }_n||u_q\rangle }^2\ge 0,\forall q\in P.& \left(25\right)\end{array}$
This is the so-called “decomposition of a partially coherent imaging system into an Incoherent Sum Of Coherent Systems” (ISOCS), which may be used to speedup the calculation of resultant intensities with the eigenvectors {φ_n(r)}_n≧1 pre-computed. Imaging kernels refer herein to mathematical functions that characterize a photolithographic process, or a general optical imaging process, where each imaging kernel may be convolved (or correlated using inner product) with a mask transmittance function to yield an amplitude distribution as a convolution result, then the convolution result may be squared to get a contributory resultant intensity distribution, finally the contributory intensities of different imaging kernels may be summed together to generate the intensity distribution for the resultant image of the mask. It follows that the functions {φ_n(r)}_n≧1 are imaging kernels of the partially coherent imaging system. Moreover, the imaging kernels {φ_n(r)}_n≧1 are mutually orthogonal, and they form an orthogonal basis for the functional space L²(P). The problem of computing resultant intensity reduces to calculating convolutions of a mask transmittance function with pre-computed imaging kernels. The convolutions may be efficiently computed as simple function multiplications in the spatial frequency domain, with the help of FFT to transform between the two domains of real space and spatial frequency. Working in the spatial frequency domain will be discussed in details below. Other approaches of resultant intensity calculation rely on decomposing mask polygons into primitive geometric shapes, such as corners, rectangles, or triangles, pre-computing and storing convolution results, which may be called pre-images, of the primitive shapes with the imaging kernels, and calculating the amplitudes of a mask image by summing up pre-images of all primitive shapes constituting the mask. Then the resultant intensity is simply the sum of the amplitudes squared. In such approaches, imaging kernels may not be directly stored numerically in physical devices, rather the pre-images of primitive shapes could be stored and transferred between physical devices of computer memory. However, the imaging kernels may be recovered by differentiating the pre-images, or calculating the image of a tiny pixel on mask, for example a 1 nm×1 nm square, using the pre-images.

[0139] It may be noted that the inner product representation,
φ_n∥u_q=∫_Pφ*_n(r)u(r+q)dr,  (26)
for a typical n≧1, is not a convolution in the strict traditional sense, rather it is usually called the correlation function of φ_n(r) and u(r). However, by a change of variable r→−r, the inner product representation of equation (26) may be interpreted as a convolution between φ*_n(−r) and u(r) strictly in the conventional sense,
φ_n∥u_q=∫_Pφ*_n(−r)u(q−r)dr.  (27)
Such change of variable and convolutional interpretation are applicable throughout this description to all inner products between imaging kernels (as well as other model kernels) and a coordinate-shifted mask transmittance function u_q(r). As a consequence, in this description we shall not differentiate between the conventional convolution and correlation function of two functions. Rather we will use the two representations interchangeably and refer to them commonly as convolution.

[0140] It is also possible to reformulate the mathematical equations to manifest convolution structures explicitly in the conventional sense. To that end, one may start from equation (7) and work with the scale-normalized coordinate system,

[0141] $\begin{array}{cc}\begin{array}{c}{I}_2\left(q\right)={\int }_P^{}{\int }_P^{}{u}^{*}\left(r_1\right)K_{12}^{*}\left(r_1,q\right)J_1^{-}\left(r_1,r_2\right)K_{12}\left(r_2,q\right)u\left(r_2\right)d{r}_{1}d{r}_2\\ ={\int }_P^{}{\int }_P^{}{u}^{*}\left(r_1\right)K_{12}^{*}\left(0,q-r_1\right)J_1^{-}\left(r_1,r_2\right)K_{12}(0,q-\\ r_2)u\left(r_2\right)d{r}_{1}d{r}_2\\ ={\int }_P^{}{\int }_P^{}{u}^{*}\left(r_1\right)K_{12}^{*}\left(0,q-r_1\right)J_1^{-}\left(q-r_2,q-r_1\right)K_{12}(0,q-\\ r_2)u\left(r_2\right)d{r}_{1}d{r}_2\\ ={\int }_P^{}{\int }_P^{}{u}^{*}\left(r_1\right)T_{-}\left(q-r_1,q-r_2\right)u\left(r_2\right)d{r}_{1}d{r}_2,\end{array}& \left(28\right)\end{array}$
where T_(r₁,r₂) is a reformulated TCC,

[0142] $\begin{array}{cc}{T}_{-}\left(r_1,r_2\right)\stackrel{\mathrm{def}}{=}T\left(-r_1,-r_2\right)=K_{12}^{*}\left(0,r_1\right)J_1^{-}\left(r_2,r_1\right)K_{12}\left(0,r_2\right).& \left(29\right)\end{array}$
Using reformulated functions,
φ_n(r)=φ*_n(−r), ∀n≧1,  (30)
T_(r₁,r₂) may be Mercer-expanded as,

[0143] $\begin{array}{cc}{T}_{-}\left(r_1,r_2\right)=T\left(-r_1,-r_2\right)=\sum _{n\ge 1}{\lambda }_n{\varphi }_n\left(-r_1\right){\varphi }_n^{*}\left(-r_2\right)=\sum _{n\ge 1}{\lambda }_n{\Phi }_n^{*}\left(r_1\right){\Phi }_n\left(r_2\right),& \left(31\right)\end{array}$
which may be plugged into equation (28) to yield,

[0144] $\begin{array}{cc}{I}_2\left(q\right)=\sum _{n\ge 1}{\lambda }_n{\int }_P{\int }_P{u}^{*}\left(r_1\right){\Phi }_n^{*}\left(q-r_1\right){\Phi }_n\left(q-r_2\right)u\left(r_2\right)d{r}_{1}d{r}_2=\sum _{n\ge 1}{\lambda }_n{u\left(q\right){\mathrm{★\Phi }}_n\left(q\right)}^2,\mathrm{with}& \left(32\right)\\ u\left(q\right){\mathrm{★\Phi }}_n\left(q\right)={\Phi }_n\left(q\right)★u\left(q\right)\stackrel{\mathrm{def}}{=}{\int }_{P}u\left(r\right){\Phi }_n\left(q-r\right)dr={\int }_P{\Phi }_n\left(r\right)u\left(q-r\right)dr,{\forall }_n\ge 1,& \left(33\right)\end{array}$
being convolutions represented strictly in the conventional sense.

[0145] Formulations in the Spatial Frequency Domain

[0146] Sometimes it may be convenient to work in the frequency domain. One advantage of working in the frequency domain is that the plane waves are eigenvectors of the spatial shift-invariant propagators, so that the action of a propagator on a plane wave is simply a scaling of amplitude and a shift of phase. A spatial frequency (f_x,f_y,f_z) may be associated with a plane wave exp[i2π(f_xx+f_yy+f_zz)], with (x,y,z) coordinating the 3D real space. If x, y, and z are measured in units of λ, the wavelength of the source light in vacuum, then (f_x,f_y,f_z) lies on a spherical surface with radius η, that is, f_x²+f_y²+f_z²=η², where η is the refractive index of the medium transmitting the wave. In a spherical coordinate system, the spatial frequency may be represented as (f_x,f_y,f_z)=η(cos θ, sin θ cos φ, sin θ sin φ), θ∈[0,π], φ∈[0,2π). When restricted to a plane with z=z₀, the plane wave becomes e^i2πfxz0|f with |f being the ket notation for a wave function,

[0147] $\begin{array}{cc}\langle r❘f\rangle \stackrel{\mathrm{def}}{=}\exp \left(\mathrm{i2\pi }f·q\right)=\exp \left[i2\pi \left(f_{x}x+f_{y}y\right)\right],& \left(34\right)\end{array}$
where

[0148] $r=\left(x,y\right)\in P,f\stackrel{\mathrm{def}}{=}\left(f_x,f_y\right)\in F,\mathrm{with}F\stackrel{\mathrm{def}}{=}\left\{\left(f_x,f_y\right):f_x^2+f_y^2\le f_{\max }^2\right\}$
being the space of bounded 2D spatial frequencies, where

[0149] $f_{\max }\stackrel{\mathrm{def}}{=}N.A.+{\sigma }_{\max }.$
It is obvious that |f∈L²(P) because of the finiteness of P in our convention. It should be noted that there is never an actual 2D wave corresponds to any |f, f∈F, which is just a slice of a 3D wave on a plane with constant z. Light propagates with the phase velocity c/η in a medium with index η, where c is the speed of light in vacuum. With this understanding, however, we shall often use |f with f=(f_x,f_y) to denote either a 3D wave exp[i2π(f_xx+f_yy+f_zz)] or a slice of it on a plane with constant z, so long as the context provides information about which surface z=constant is concerned, or which medium with what refractive index is transmitting the wave. Such convention is backed up by the fact that the x- and y-directional spatial frequency f_x and f_y are conserved by reflection/refraction at interfaces perpendicular to the z-axis between media with different refractive indices, due to or manifested by Snell's law of reflection/refraction.

[0150] We shall use frequently the following pair of Fourier and inverse Fourier transforms and their nature extensions to higher-dimensional spaces,

[0151] $\begin{array}{cc}U\left(f\right)=ℱ\left[u\left(x\right)\right]\stackrel{\mathrm{def}}{=}\int u\left(x\right)\exp \left(-i2\pi \mathrm{fx}\right)dx,& \left(35\right)\\ u\left(x\right)={ℱ}^{-1}\left[U\left(f\right)\right]\stackrel{\mathrm{def}}{=}\int U\left(f\right)\exp \left(i2\pi \mathrm{fx}\right)df.& \left(36\right)\end{array}$
By applying the operator identity ∫_F|ff|df=1 twice, equation (22) may be rewritten as,

[0152] $\begin{array}{cc}{I}_2\left(q\right)\begin{array}{c}=\langle u_qu_q\rangle \\ ={\int }_F{\int }_F\langle u_q|f_1\rangle \langle f_1f_2\rangle \langle f_2|u_q\rangle d{f}_{1}d{f}_2\\ ={\int }_F{\int }_F{\bigcup }_q^{*}\left(f_1\right)T\left(f_1,f_2\right){\bigcup }_q\left(f_2\right)d{f}_{1}d{f}_2,\end{array}& \left(37\right)\end{array}$
where,

[0153] $\begin{array}{cc}\begin{array}{c}T\left(f_1,f_2\right)\stackrel{\mathrm{def}}{=}\langle f_1f_2\rangle ={\int }_P{\int }_P\langle f_1|r_1\rangle \langle r_1r_2\rangle \langle r_2|f_2\rangle d{r}_{1}d{r}_2\\ \begin{array}{c}={\int }_P{\int }_P{e}^{-i2\pi f_1·r_1}T\left(r_1,r_2\right)e^{i2\pi f_2·r_2}d{r}_{1}d{r}_2\\ ={\int }_P{\int }_P{e}^{-i2\pi f_1·r_1}{K}_{12}^{*}\left(0,-r_1\right)J_1^{-}\left(0,r_2-r_1\right)K_{12}\left(0,-r_2\right)\\ e^{\mathrm{i2}\pi f_2·r_2}d{r}_{1}d{r}_2\\ ={\int }_P{\int }_P{e}^{i2\pi f_1·r_1}{K}_{12}^{*}\left(0,r_1\right)J_1^{-}\left(0,r_1-r_2\right)\\ K_{12}\left(0,r_2\right)e^{-\mathrm{i2}\pi f_2·r_2}d{r}_{1}d{r}_2\\ ={\int }_F{\int }_P{\int }_P{e}^{\mathrm{i2}\pi f_1·r_1}{K}_{12}^{*}\left(0,r_1\right)J\left(f\right)e^{i2\pi f·\left(r_1-r_2\right)}{K}_{12}\left(0,r_2\right)\\ e^{-i2\pi f_2·r_2}d{r}_{1}d{r}_{2}df\\ ={\int }_F{K}_{12}^{*}\left(f+f_1\right)J\left(f\right)K_{12}\left(f+f_2\right)df,\end{array}\end{array}& \left(38\right)\end{array}$
and the Fourier transforms are represented as,

[0154] $\begin{array}{cc}{U}_q\left(f\right)\stackrel{\mathrm{def}}{=}ℱ\left[u_q\left(r\right)\right]={\int }_P\langle f❘r\rangle \langle r❘u_q\rangle dr={\int }_P{u}_q\left(r\right)e^{-i2\pi f·r}dr,& \left(39\right)\\ K_{12}\left(f\right)\stackrel{\mathrm{def}}{=}ℱ\left[K_{12}\left(0,r\right)\right]={\int }_P{K}_{12}\left(0,r\right)e^{-i2\pi f·r}dr,& \left(40\right)\\ J\left(f\right)\stackrel{\mathrm{def}}{=}ℱ\left[J_1^{-}\left(0,r\right)\right]={\int }_P{J}_1^{-}\left(0,r\right)e^{-i2\pi f·r}dr.& \left(41\right)\end{array}$
The inverse Fourier transforms are of course,
u_q(r)=F⁻¹[U_q(f)]=∫_FU_q(f)(r)e^i2πf·rdf,  (42)
K₁₂(0,r)=F⁻¹[K₁₂(f)]=∫_FK₁₂(f)e^i2πf·rdr,  (43)
J₁⁻(0,r)=F⁻¹[J(f)]=∫_FJ(f)e^i2πf·rdf.  (44)
K₁₂(f) is called the (frequency) transfer function of the projection system. As mentioned before, there is a one-to-one relationship between the coordinate s∈S₀ of a source point and the spatial frequency f∈F of the parallel beam generated by it. When σ_max/M is small, the relationship may be approximated as linear, namely, s≈fF_c/M, with which equation (3) becomes,
J₁⁻(r₁,r₂)=J₁⁻(0,r₂−r₁)=constant×∫_FI₀(fF_c/M)e^{i2πf·(r2−r1)}df,  (45)
and evidently,
J(f)=constant×I₀(fF_c/M).  (46)
The function T(f₁,f₂) may be regarded as the Hilbert-Schmidt kernel of the linear operator FTF⁻¹:L²(F)→L²(F), which is still Hermitian and positive definite, and has the same eigenvalues λ₁≧λ₂≧ . . . >0, with the corresponding eigenvectors,

[0155] $\begin{array}{cc}{\psi }_n\left(f\right)\stackrel{\mathrm{def}}{=}ℱ\left[{\varphi }_n\left(r\right)\right]={\int }_P{\varphi }_n\left(r\right)e^{-i2\pi f·r}dr,n=1,2,\dots ,& \left(47\right)\end{array}$
since F maps an orthonormal basis onto an orthonormal basis. Indeed, we have the following,

[0156] $\begin{array}{cc}T\left(f_1,f_2\right)=\langle f_1f_2\rangle =\sum _{n\ge 1}{\lambda }_n\langle f_1|{\varphi }_n\rangle \langle {\varphi }_n|f_2\rangle =\sum _{n\ge 1}{\lambda }_n{\psi }_n\left(f_1\right){\psi }_n^{*}\left(f_2\right),& \left(48\right)\end{array}$
and for the resultant intensity distribution,

[0157] $\begin{array}{cc}{I}_2\left(q\right)=\sum _{n\ge 1}{n}_{}{\langle {\varphi }_nu_q\rangle }^2=\sum _{n\ge 1}{\lambda }_n|{\langle {\mathrm{ℱ\varphi }}_nℱu_q\rangle }^2=\sum _{n\ge 1}{\lambda }_n{\langle {\psi }_n||U_q\rangle }^2=\sum _{n\ge 1}{\lambda }_n{{\int }_F{\psi }_n^{*}\left(f\right)U_0\left(f\right)e^{i2\pi f·q}df}^2=\sum _{n\ge 1}{\lambda }_n{{-1}^{}\left[{\psi }_n^{*}\left(f\right)U_0\left(f\right)\right]}^2,& \left(49\right)\end{array}$
where we have used,
U_q(f)=∫_Pu(r+q)e^i2πf·rdr=e^i2πf·qU₀(f).  (50)

[0158] For numerical simulations, J(f) is calculated in the frequency domain according to equation (46), then the mutual intensity J₁⁻(r₁,r₂)=J₁⁻(0,r₂−r₁) in the spatial domain may be obtained through an inverse FFT. It is also convenient to describe a projection system by its transfer function K₁₂(f) for several reasons: a) the bandlimited system is well characterized by a pupil function within a region bounded by N.A.; b) abberation effects may be readily included in the pupil function; c) wave propagation between two planes perpendicular to the z-axis is simply an f_z-dependent phase shift, with the (f_x,f_y) components of spatial frequency conserved so long as the media are uniform in the x and y directions. Then the PSF K₁₂(0,r) may be obtained through an inverse FFT. However, it may not be efficient to calculate the TCC in the frequency domain directly using equation (38), because equation (38) involves 2D numerical integration. It is more efficient to compute T(r₁,r₂) using equation (15). The TCC operator may then be eigen-analyzed to get the eigenvalues λ₁≧λ₂≧ . . . >0 and eigenvectors (imaging kernels) {φ_n(r)}_n≧1. When imaging kernels in the frequency domain {ψ_n(f)}_n≧1 are needed, they may be readily obtained by applying FFT to {φ_n(r)}_n≧1 as in equation (47). As noted above, one way to compute the resultant intensity distribution I₂(q) is using equation (49) and FFT: firstly computing the mask spectrum U₀(f), then multiplying U₀(f) with the frequency representations of the imaging kernels and performing inverse Fourier transform, finally squaring the amplitudes and summing the squares together.

[0159] A Different Approach—Decomposing J₁⁻(r₁,r₂) into Coherent Modes

[0160] To calculate the resultant intensity at any target point q=(q₁,q₂)∈P, it is sufficient to consider J₁⁻(r₁,r₂) and J₁⁺(r₁,r₂) in a 2a₀×2a₀ square region,

[0161] $\begin{array}{cc}{A}_q\stackrel{\mathrm{def}}{=}\left\{\left(x,y\right)\in P:q_1-a_0<x<q_1+a_0,q_2-a_0<y<q_2+a_0\right\},& \left(51\right)\end{array}$
where a₀ is a predetermined constant, called the “diffraction radius” that characterizes the critical distance between two points on a mask (usually scaled), beyond which the optical proximity effect between the two points may be neglected. The constant a₀ has been called by various names besides diffraction radius, including “optical ambit” or simply “ambit”, and “halo radius” or simply “halo”. In this description, the two terms “diffraction radius” and “optical ambit” are used, and in a synonymous and interchangeable manner. Because of shift-invariance, it is sufficient to consider J₁⁻(r₁,r₂) restricted in the region A₀^defA_q=(0,0), namely,

[0162] $\begin{array}{cc}{J}_1^{-}\left(r_1,r_2\right){|}_{A_0}=\{\begin{array}{cc}{J}_1^{-}\left(r_1,r_2\right)& r_1,r_2\in A_0,\\ 0,& \mathrm{otherwise},\end{array}& \left(52\right)\end{array}$
and using the following to calculate J₁⁺(r₁,r₂) with r₁,r₂∈A_q,
J₁⁺(r₁+q,r₂+q)=u*_q(r₁)J₁⁻(r₁,r₂)u_q(r₂), ∀q∈P, r₁,r₂∈A₀.  (53)
More generally, J₁⁻(r₁,r₂) may be multiplied by a smoother window function to yield a smoother version of windowed mutual intensity. In particular, a gradually varying and non-negative (or equivalently, non-positive) window function W(r) may be employed, which traces approximately the envelope of K₁₂(0,r) and satisfies,

[0163] $\begin{array}{cc}\{\begin{array}{cc}W\left(r\right)>0,& \forall r\in A_0,\\ W\left(r\right)=0,& \mathrm{otherwise}.\end{array}& \left(54\right)\end{array}$
For example, W(r) may be a truncated Gaussian function

[0164] $\begin{array}{cc}W\left(r\right)=\{\begin{array}{cc}\exp (-{\beta }_1{r}^2/a_0^2,& \forall r\in A_0,\\ 0,& \mathrm{otherwise},\end{array}& \left(55\right)\end{array}$
or a raised cosine window function represented as,

[0165] $\begin{array}{cc}W\left(r\right)=\{\begin{array}{cc}1,& 0\le r<\left(1-{\beta }_2\right)b_0,\\ \frac{1}{2}+\frac{1}{2}\cos \pi \left[\frac{r-\left(1-{\beta }_2\right)b_0}{2{\beta }_2{b}_0}\right],& \left(1-{\beta }_2\right)b_0\le r<\left(1+{\beta }_2\right)b_0,\\ 0,& r\ge \left(1+{\beta }_2\right)b_0,\end{array}& \left(56\right)\end{array}$
where a₀ is the previously mentioned diffraction radius, b₀=a₀/(1+β₂), β₁>0 and 0≦β₂≦1 are constant parameters controlling the shape and width of the windows. Another example of W(r) may be a gradually varying function which coincides with the local maximum points of |K₁₂(0,r)| inside A₀ and vanishes outside A₀. Namely, if r′₀,r′₁,r′₂, . . . are the locations inside A₀ at which |K₁₂(0,r)| reaches local maximum, then W(r) may be a gradually varying function satisfying W(r′_n)=|K₁₂(0,r′_n)| for all n=0,1,2, . . . . In any case, once a window function W(r) is chosen, a mutual intensity restricted by W(r) may be calculated as,

[0166] $\begin{array}{cc}{J}_1^{-}\left(r_1,r_2\right){❘}_W\stackrel{\mathrm{def}}{=}W\left(-r_1\right)J_1^{-}\left(r_1,r_2\right)W\left(-r_2\right),& \left(57\right)\end{array}$
and a PSF restricted by W(r) may be represented as,

[0167] $\begin{array}{cc}{K}_{12}\left(0,r\right){|}_w\stackrel{\mathrm{def}}{=}\{\begin{array}{cc}{K}_{12}\left(0,r\right)/W\left(r\right),& \forall r\in \left(1-{\beta }_3\right)A_0,\\ 0,& \mathrm{otherwise},\end{array}& \left(58\right)\end{array}$
where β₃∈[0,1) is a small number to clip a small portion of K₁₂(0,r) close to the edge of A₀, as merely one way to prevent potential amplitude runaway of K₁₂(0,r)|_W around the edge. Equations (57) and (58) may be plugged into equation (15) as well to calculate a TCC restricted by W(r),

[0168] $\begin{array}{cc}T\left(r_1,r_2\right)\approx T\left(r_1,r_2\right){❘}_W\stackrel{\mathrm{def}}{=}{K}_{12}^{*}\left(0,-r_1\right){❘}_W{J}_1^{-}\left(r_1,r_2\right){❘}_W{K}_{12}\left(0,-r_2\right){❘}_W.& \left(59\right)\end{array}$
The approximation T(r₁,r₂)≈T(r₁,r₂)|_W may be justified by the fact that |K₁₂(0,r)| becomes negligibly small outside A₀. Note that J₁⁻(r₁,r₂)|_A0 in equation (52) may be regarded as a special case of J₁⁻(r₁,r₂)|_W in equation (57) with a hard-cut window function,

[0169] $\begin{array}{cc}W\left(r\right)\stackrel{\mathrm{def}}{=}\{\begin{array}{cc}1,& r\in A_0,\\ 0,& \mathrm{otherwise}.\end{array}& \left(60\right)\end{array}$
Regarded as a Hilbert-Schmidt kernel, a windowed mutual intensity J₁⁻(r₁,r₂)|_W may also be associated to a linear operator J₁⁻:L²(P)→L²(P) as in the following,
r₁J₁⁻|r₂=J₁⁻(r₁,r₂)|_W, ∀r₁,r₂∈P.  (61)
It follows that J₁⁻ is Hermitian because J₁⁻(r₂,r₁)|_W=[J₁⁻(r₁,r₂)|_W]* and positive definite as a result of equation (3). As a consequence, J₁⁻ has eigenvalues μ₁≧μ₂≧ . . . >0, and corresponding eigenvectors ε₁(r),ε₂(r), . . . , which decomposes the Hilbert-Schmidt kernel J₁⁻(r₁,r₂)|_W as,

[0170] $\begin{array}{cc}{J}_1^{-}\left(r_1,r_2\right){|}_w=\sum _{n\ge 1}{u}_n{\xi }_n\left(r_1\right){\xi }_n^{*}\left(r_2\right),\forall r_1,r_2\in P.& \left(62\right)\end{array}$
Such Mercer's expansion for J₁⁻(r₁,r₂)|_W is called the coherent-mode representation of the partially coherent illumination field, where different modes are statistically independent (i.e., incoherent), while the spatial amplitude distribution is perfectly coherent within each mode. Furthermore, any linear (possibly lossy) propagation of the illumination field may still be represented by an incoherent combination of the coherent modes linearly propagated. As far as the second-order coherence theory is concerned, the partially coherent illumination field is physically indistinguishable from a group of lasers each oscillating with one of the coherent modes, and with the relative phase (as a time series) between any two lasers uniformed distributed in [0,2π). The coherent-mode representation of J₁⁻(r₁,r₂)|_W may be plugged into equation (59) to decompose the TCC,

[0171] $\begin{array}{cc}T\left(r_1,r_2\right)\approx T\left(r_1,r_2\right){|}_w=\sum _{n\ge 1}{\mu }_n{\phi }_n\left(r_1\right){\phi }_n^{*}\left(r_2\right),\forall r_1,r_2\in P,& \left(63\right)\\ {\phi }_n\left(r\right)\stackrel{\mathrm{def}}{=}{\left[K_{12}\left(0,-r\right){|}_w\right]}^{*}{\xi }_n\left(r\right),\forall n\ge 1,\forall r\in P.& \left(64\right)\end{array}$
Even though {ψ_x(r)}_n≧1 may no longer be mutually orthogonal, nor they may be eigenvectors of the TCC operator T with {μ_n}_n≧1 as eigenvalues, the kernel decomposition of equation (63) could nonetheless be applied to equation (14) to yield,