Three-state lithography model: an enhanced mathematical approach to predict resist characteristics in grayscale lithography processes

Abstract. Background: Physical modeling of grayscale lithography processes for the prediction of photoresist heights leads to complex mathematical algorithms. A promiment example is the numerical simulation of the photoresist shape after development through Dill’s equations. These grayscale lithography models exhibit accurate prediction quality but can not directly implemented into mask layout tools to simplify the layout procedure. Limited process windows, changes in the mask design, variations of the used materials or manufacturing tools lead to time-consuming and cost-intensive test procedures to adjust the photoresist model for sufficient results. Aim: The focus of this work is to enhance current grayscale lithography models for a straightforward method with the same precise prediction of remaining photoresist heights to simplify the mask layout process. Moreover, we aim for an uncomplicated optimization of the model to minimize the empirical analysis necessary for its use. Approach: Based on experimental results, we deploy a sectionally defined mathematical expression that includes the theory of Fraunhofer diffraction and illumination-dependent activation of the photo-sensitive component and its solubility in developer. Results: We produced pyramidal, spherical and chess field structures with exposure doses of 3000 and 15  ,  000  J  /  m2 on bare silicon substrates with 100-nm resolution and on silicon substrates with anti-reflective coatings, with accuracy as fine as 20 nm. Conclusion: The proposed three-state lithography model has been verified by experimental evaluation. It is able to operate in a wide process window and can be directly implemented in existing mask layout software. This model ensures a cost-efficient and precisely controlled production of three-dimensional topographies using grayscale lithography processes.


Introduction
Over the past three decades, grayscale lithography has steadily gained interest both in industrial and scientific community. The ability to create desired three-dimensional (3D) topographies on a surface through a cost-efficient method enabled grayscale lithography to improve the field of micro-optics. 1,2 The semiconductor industry exploits the advantages of grayscale lithography for the generation of special ion-implantation profiles, e.g., for the junction termination extension technique. 3 Investigations by Schneider et al. 4 revealed the possibilities of sharp 3D-ion implantation profiles on several simulations using semiconductor applications. However, recent research leverages grayscale lithography for production of blazed gratings in miniaturized *Address all correspondence to Bassem Badawi, bassem.badawi@emft.fraunhofer.de spectrometers, solid emitters for electrospray ionization applications, and liquid lenses for use in optogenetics. 5 The generation of arbitrary grayscales can be accomplished through various lithography technologies, i.e., e-beam, x-ray, or traditional UV exposure. This work utilizes UV exposure for structuring of photo-sensitive resists, as it is a common technique in the semiconductor industry. Due to complex mask preparation methods and their demanding manufacturing requirements, this technology becomes very time-consuming and cost-intensive for the production of 3D topographies. To achieve varying lateral height levels in a resist, the exposure dose must be modified laterally. The current resist technologies yield high contrasts, resulting in optimal imaging properties. This requires precise control of the energy exposure directed at the photosensitive resist. An effective principle to control the exposure intensity was published by Waits et al. 6 and utilizes the theory of Fraunhofer diffraction via sub-resolution patterns on a photomask. Prior to this, the mathematical background for this approach was presented in 1995 by Wagner et al. 7 As reported by Smith et al., 5 the usage of current simulation tools for grayscale lithography processes can lead to limited process windows in which only specific mask features for certain exposures yield acceptable results. Their work showed process windows for mask features with an areal density of 0.45 and 0.35 at SPR518 and SPR220 photoresist with possible exposure doses between 2000 and 3000 J∕m 2 and 2200 and 4000 J∕m 2 .
Standard lithography simulation tools apply Dill's equations, 8 which use the numerical solution of time-dependent differential equations to simulate the inhibitor concentration in thin, vertically stacked photoresist layers for lateral exposure intensity, calculated using the Fourier transformation of the mask and a given modulation transfer function. An empirically acquired development rate for specific inhibitor concentrations is used in conjunction with Dill's equations to calculate the resist solubility. The approach in this work simplifies this current mathematical process to simulate the photoresist solubility with an effortless method, which is accessible in mask layout software tools. Moreover, this mathematical concept enables an easy optimization procedure improving upon current simulation applications. This model has been given the name "three-state lithography model" due to the three different solubility phases that occur in the photoresist during exposure.

Experimental
For the experimental procedure, we utilize a postive AZ ECI 3027 photoresist for all test processes and the pyramidal, spherical, and chess field sample structures. The target shape designs of the pyramidal, chess field, and spherical structures are shown in Figs. 1-3. Operating on a Suess Gamma spin coating and developing system, resist thicknesses of 3 μm were coated on 200-mm silicon substrates. The UV exposures were performed using a Canon FPA-3000 i4 Stepper with an i-line (λ ¼ 365 nm) illumination system, a numerical aperture NA ¼ 0.62, and an object to image magnification ratio M ¼ 5∶1. No post-exposure bakes were executed before developing in AZ 726 MIF. To flatten the resist irregularities on the pyramidal resist structures, a subsequent hard-bake at 115°C was performed for 120 s. Thickness measurements of the contrast curve were carried out using a KLA-Tencor UV1280. The resultant pyramidal and chess field structures were investigated by KLA Tencor P-17 Stylus Profiler, while the spherical structures are inspected using a Digital Instruments Dimension 5000 Atomic Force Microscope (AFM) system.
To handle the high amount of data necessary for mask design, a combination of MATLAB and Tanner L-Edit scripts were used for the automatic generation of mask features and their positioning on the mask based on the three-state photoresist model. For this purpose, a test reticle was produced containing squared fields, with side lengths of 2.7 mm, consisting of quadratic shapes with a pitch of 1.5 μm and fill factors from 0.04 to 1. For sufficient validation of the process window of the proposed model, similar test layouts were fabricated using exposure doses of 3000 and 15;000 J∕m 2 . For additional verification of the model, all sample structures were created in both lightfield and darkfield ambients to prove the development uniformity under different exposure conditions.

Modeling
Reliable height predictions require precise knowledge of the contrast function for a given resist and developer combination. Moreover, this function or contrast curve includes specific process  parameters, such as the optical absorption and reflection properties of the substrate, and as applied soft-bake and developing process parameters. To compensate for the asymmetric course of the contrast curve, the exposure dose of the measured heights were corrected under consideration of the resist's inner light absorption by the Lambert-Beer law using the boundary condition of a constant absorption coefficient α during exposure. Subsequently, the corrected data set was fitted to a modified sigmoid function (henceforth designated s-function), which matches best to the contrast curve data of our measurments and directly correlates to the production of the soluble photoresist component (carboxylic acid) during exposure. This correction determines a threshold energy dose needed to start the chemical transformation of the photo-sensitive component.
As a result, the mathematical expression for the measured contrast curve can be written as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 6 0 4 where a and b are the parameters extracted from the corrected data, whereas D 0 represents the initial resist thickness. As presented in the work of Wagner et al., 7 rectangular sub-resolution mask patterns with a width x and length y separated through a constant pitch p can be applied to alter the intensity distribution in the areal image. Out of these parameters, a fill factor GS for mask features is expressed as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 4 8 8 Considering the Rayleigh resolution limit, we have set p to 300 nm for the used i-line projection system. The physical basis of intensity variations for specific mask features is a limited entrance pupil of the projection system, which acts as an optical low-pass filter. The diffraction image plane in a Köhler illumination system is located in the entrance pupil position. 9 Considering ray optics, the object-sided opening angle θ can be given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 3 8 6 where NA is the numerical aperture and M represents the magnification ratio. The mathematical treatment of the intensity distribution is derived from the theory of Fraunhofer diffraction 10 and describes the propagation of an incident plane wave on a lattice in the far field. For a twodimensional quadratic aperture, the intensity I d of the diffraction pattern is proportional to E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 2 9 6 where α and β represent the angular divergence of the diffracted light in x and y direction and λ is the illumination wavelength. For a squared-entrance pupil geometry, the two-dimensional integral of Eq. (4) over the angles of divergence determines the intensity in the aerial image plane and is expressed as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 1 9 7 where θ∕2 acts as spatial cut-off frequency of the optical low-pass filter. The impact of altering intensities on the contrast curve in Eq. (1) can be taken into account by k F (henceforth designated Fourier factor), which includes the contribution of accessible intensity propagating through the entrance pupil. The Fourier factor is expressed using Eq. (2) as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 6 ; 1 0 6 where the expression SiðÞ represents the sine integral, k x ¼ πx sinðαÞ∕λ and k y ¼ πy sinðβÞ∕λ.
A closer examination on Eq. (6) reveals that the transmitted light is not merely proportional to the squared fill factor, but also dependent on the transmitted wavelength, system entrance pupil, and the pattern size. The part of illumination passing the entrance pupil of the optical system applied in this work is shown in Fig. 4 and is dependent on the size of the mask features. To ensure the Fourier factor is dependent exclusively on the fill factor, mask features are limited to quadratic shapes. This limitation provides the additional benefit of lowering the calculation efforts needed for the model. In this case, k F becomes one-dimensional. The course of k F for quadratic mask features compared to a simple quadratic dependency is shown in Fig. 5.
As a consequence of lower intensities for narrower mask patterns, the threshold for chemical transformation in the resist shifts to higher energy doses. This shifted threshold E shift is dependent on GS and is added to the initial threshold energy dose E in Eq. (1) with respect to the expression E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 6 ; 3 3 1 Fig. 4 Contribution of diffraction to k F in dependence on the standardized mask pattern dimensions x ∕p and y ∕p. where c and d represent characteristic parameters for each type of photoresist. These parameters are empirically determined by exposing a test mask to a variable exposure dose. During exposure, the optical properties of the positive photoresist change due to the chemical formation of carboxylic acid. The chemical transformation, which primarily impacts the extinction coefficient, causes the photoresist to appear transparent to UV radiation. This effect, known as bleaching, is strongly dependent on the chemical compounds of the photoresist and affects the slope of the contrast curve individually for any mask pattern fill factor. To a good approximation, the soluble photoresist component near the resist-developer surface is modeled by utilizing Eqs. (1) and (8) without taking Lambert-Beer into account. The solubility of deeperlying material within the photoresist follows Eq. (7). To create a continuous transition between the shallow and deep photoresist solubilities, an intermediate phase is mathematically expressed through a spline-interpolation function and is defined individually, dependent on the fill factor. Therefore, the three-state photoresist model follows the sectionally defined expression E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 6 ; 5 8 0 (9) where the boundary points P 1 and P 2 represent the intersections between the phases. The symbols e, f, g, and h operate as fit parameters for the spline function. To decrease calculation complexity, the implicit expressions in Eqs. (1) and (7) are expressed in their explicit forms in Eq. (9), dependent on the photoresist thickness T.

Verification and Discussion
In this work, a dose range from 600 to 3000 J∕m 2 discloses the contrast curve shape. The measured resist thickness T dependent on the exposure dose E, the data after the correction and the s-functions are shown in Fig. 6  the intermediate phase reaches its maximum depth, before the intermediate phase shrinks for higher fill factors. This physical phenomenon is called reciprocity failure, and was investigated in detail by Sheats 11 in 1988. Accordingly, the chemical meta-stable product absorbance during the carboxylic acid production is shown to be responsible for this effect. To validate the simulated contrast curves, Fig. 7 shows measurement data produced with the test mask out of Sec. 2.1 for fill factors GS ¼ 0.28, 0.32, and 0.40, which follow the simulated contrast curve courses. Accordingly, the model enables the accurate prediction of remaining resist heights in a wide range of low to high exposure doses. Figure 8 shows two instances of the photoresist pyramidal structure at different positions (site I and site II) on a wafer compared to the target design topography. Furthermore, Fig. 8 also illustrates the topography after a hard-bake procedure for each energy dose and in both lightfield and darkfield ambients.
The formation of a standing wave within the photoresist results in a minimum height step of ≈100 nm, which corresponds to the constructive interference locations within the standing wave, as seen on the faces of the pyramidal structures. The topographies show discrepancies between the designed and resultant topographies of up to 200 nm. In darkfield ambient, both exposure doses yield similar profiles. Thickness discrepancies between lightfield and darkfield ambients  of 40 to 60 nm were observed for 15;000 J∕m 2 , while the lower dose of 3000 J∕m 2 led to variations of 100 to 200 nm. The higher instability at lower exposure doses against ambient is caused by the decreased exposure time, resulting in less irregularly formed carboxylic acid. Additionally, this effect led to better roughness properties for higher exposure doses. At an exposure dose of 3000 J∕m 2 , the produced pyramid slope exhibited inaccurate angles and sensitivity to ambient conditions. Conversely, an exposure dose of 15;000 J∕m 2 yielded more desirable results. In particular, the cavities produced at 3000 J∕m 2 showed inconsistent depths. The added hard-bake process led to smooth surfaces of the pyramidal structure. Due to a missing height step at a photoresist thickness of ≈700 nm, the topographies at exposure doses of 15;000 J∕m 2 still show irregularities. To achieve a smooth surface without hard-bake post processing the silicon surface is coated with a BARC to avoid the formation of a standing wave within the photoresist. Due to this substrate modification, the contrast curve function changes, as shown in Fig. 6. To limit errors from the mask design for a contrast curve on bare silicon, the target sphere structure in Fig. 3 is fabricated only with an exposure dose of 3000 J∕m 2 . An AFM line measurement over the sphere center is shown in Fig. 9. Figure 10 shows a 30 μm × 30 μm map over the sphere.
Although the mask design is not based on the substrate surface characteristics, the resultant resist sphere shows a satisfactory level of consistency compared to the target sphere design. The base of the sphere at ≈200 nm could not be printed into the resist due to the significantly steeper contrast curve course for thinner photoresist heights with BARC coating.
The manufactured chess field structures are shown in Figs. 11(a)-11(c). Here the simulated step designs for BARC/silicon sufaces are adjusted in consideration of the changed contrast curve in Fig. 6.  As seen previously, height variation under 100 nm are not resolvable on a bare silicon substrate [ Fig. 11(a)] due to standing wave effects within the photoresist. The BARC/silicon substrate optimization reveals a vertical resolution limit of ≈(20 to 40 nm) for 15;000 J∕m 2 [ Fig. 11(b)], if structure edges are neglected. This leads to a continuous height variation and not into a stepwise vertical change of ∼8 nm between the fields. The continuous slope is caused by the chemical resolution limit of the photoresist, which is usually applied for lateral critical dimensions around 350 nm. Moreover, higher exposure intensities [ Fig. 11(c)] create rougher surfaces, which further increases the resolution limit.

Conclusion
The mathematical procedure proposed in this paper describes the change of the soluble photoresist component dependent on mask features and fixed illumination settings for the production of arbitrary 3D topographies. The capabilities of the model were demonstrated through the creation of pyramidal, spherical and chess field photoresist structures with a prediction accuracy down to 20 to 40 nm. An additional parameter optimization process using the three-state lithography model will lead to an increasing prediction quality. The mathematical approach was applied in a wide process window between 3000 and 15;000 J∕m 2 , dependent on given optical system parameters and mask features with fill factors between 0.17 and 0.81. Process variations or changes in fabrication tools and materials are easily accounted for in the procedure and do not require large empirical tests to examine the change of the development rate for specific inhibitor concentrations, as presented by Dill et al. 8 Moreover, this straightforward method enables an adequate prediction of remaining photoresist after development without numerical techniques used in commercial simulation software, such as PROLITH, 12 to solve time dependent differential equations. The simplified mathematical description of the photoresist solubility allows for easy integration within standard mask layout tools. Thus, the mask layout for grayscale applications is feasible without additional simulation tools. The three-state lithography model enables a cost-efficient method for the production of 3D topographical structures in micro-electromechanical systems (MEMS), micro-opto-electro-mechanical systems (MOEMS) and semiconductor applications with a high degree of accuracy. Due to substantial simplifications in the mathematical description of the photoresist and light physics, the proposed model is unable to simulate diffusion processes or interference effects. Therefore, the three-state lithography model is targeted for the use in grayscale lithography applications and is not developed to replace the available, sophisticated numerical simulation tools. Fig. 11 Resulting example step topographies of the target chess field in Fig. 2 for an exposure of 3000 J∕m 2 in (a) on bare silicon substrate and on BARC/silicon substrate for an exposure of 15000 J∕m 2 in (b) and 3000 J∕m 2 in (c). First two steps in (b) and (c) on structure edges are neglected.