# Rigorous Model-Based Mask Data Preparation Algorithm Applied to Grayscale Lithography for the Patterning at the Micrometer Scale

Pierre Chevalier, Patrick Quemere, Sebastien Berard-Bergery, Jean-Baptist Henry, Charlotte Beylier, Jérôme Vaillant

<figure><img src="image_1.png" /></figure>

## To cite this version:

Pierre Chevalier, Patrick Quemere, Sebastien Berard-Bergery, Jean-Baptist Henry, Charlotte Beylier, et al. Rigorous Model-Based Mask Data Preparation Algorithm Applied to Grayscale Lithography for the Patterning at the Micrometer Scale. Journal of Microelectromechanical Systems, 2021, 30 (3), pp.442 - 455. 10.1109/jmems.2021.3067475. cea-04506422

HAL Id: cea-04506422

https://cea.hal.science/cea-04506422

Submitted on 15 Mar 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Rigorous Model-Based Mask Data Preparation Algorithm Applied to Grayscale Lithography for the Patterning at the Micrometer Scale

Pierre Chevalier<sup>10)</sup>, Patrick Quéméré, Sébastien Bérard-Bergery, Jean-Baptiste Henry, Charlotte Beylier, and Jérôme Vaillant

***Abstract***—Grayscale mask creation has for the most part been restricted to over-simplified optical and resist models usually based on a contrast curve approach. While this technique has proven to work for microstructures of large dimensions (ten to hundreds of micrometers), its capability has not been assessed for microstructures with smaller dimensions. In this paper, a rigorous lithographic model has been developed in Python to simulate the process of imaging, exposure and development of an i-line photoresist. Using this model, a mask data preparation algorithm capable of optimizing simultaneously both the size and position of the dots on a grayscale mask has been implemented. Experimental results after development of the photoresist confirm the capability of our mask data preparation algorithm to achieve microstructures with dimensions ranging between 1 to 3µm.
[2021-0018]

***Index Terms***—Grayscale lithography, microfabrication, data-preparation, optimization, Python.

## I. INTRODUCTION

GRAYSCALE lithography is an alternative lithography technique enabling the patterning of 3D microstructures in a photoresist material with application ranging from microfluidic [1], optics [2]–[5] and MEMS [6], [7]. Where conventional lithography seeks to produce purely vertical structures, grayscale lithography can generate slopes and smooth profiles with different heights with a single lithography and development step. Grayscale lithography consists in a combination between a low contrast photoresist and a spatially variable exposure dose that control the resist height after development. Such variable exposure dose is generated using a chrome on glass photo-mask designed with sub-resolutions dots arranged on a fixed grid. The dose modulation is obtained by adjusting the size of the chromium

dots that affects the mask diffraction. By using approximate models to compute the aerial image and the resist development this technique is capable of patterning vertically shaped microstructures with dimensions ranging between ten to a few hundreds of micrometers. The generation of a variable exposure dose is not restricted to grayscale lithography. Direct laser writing [8]–[10] or digital micromirror devices [11]–[15] are examples of maskless microfabrication techniques where the substrate surface is directly exposed without requiring an expensive and time-consuming photomask creation.

This paper proposes a rigorous approach to the modeling of grayscale lithography. A lithographic model compatible with i-line lithography has been developed. This model is used to create a new mask data preparation algorithm capable of optimizing both the size and the position of the dots for any resist target. The quality of this approach will be illustrated on experimental profiles of microstructures with vertical and lateral dimensions below 2µm.

The first part of the article introduces the working principles and limitations of existing grayscale lithography approaches. The second part thoroughly describes the optical and resist models that have been developed to simulate the resist profile for a given photomask and process conditions. Then the mask optimization algorithm principles are explained and illustrated for a microlens pattern. The last section presents experimental results for a variety of micro-structures obtained with a photomask designed with our mask data preparation algorithm. The final part concludes on the capability of our approach.

## II. BACKGROUND

The size, shape and disposition of the sub-resolutions dots on the mask impact its transmission, thus the light intensity received by the resist during exposure. Depending on the strategy used for the mask creation, grayscale publications can be split in three groups:

1. Pulse width modulation: the dots are set on regular grid on the mask. Their size can be changed but not their position.

2. Pulse density modulation: the dots size is fixed but their position can be changed.

3. A combination between the two.

For its simplicity of implementation, the first approach is predominantly represented in the literature. Figure 1 defines the

Manuscript received January 22, 2021; revised March 2, 2021; accepted March 7, 2021. Date of publication March 24, 2021; date of current version June 1, 2021. This work was supported by the STMicroelectronics and CEA Leti. Subject Editor P. M. Sarro. (Corresponding author: Pierre Chevalier.)

Pierre Chevalier and Charlotte Beylier are with the STMicroelectronics, Reticle Enhancement Technique, 38927 Crolles, France (e-mail: pierre.chevalier@st.com; charlotte.beylier@st.com).

Patrick Quéméré, Sébastien Bérard-Bergery, Jean-Baptiste Henry, and Jérôme Vaillant are with the University of Grenoble Alpes, CEA, Leti, F-38000 Grenoble, France (e-mail: patrick.quemere@cea.fr; sebastien.berardbergery@cea.fr; jean-baptiste.henry@cea.fr; jerome.vaillant@cea.fr).

Color versions of one or more figures in this article are available at https://doi.org/10.1109/JMEMS.2021.3067475.

Digital Object Identifier 10.1109/JMEMS.2021.3067475
<figure><img src="image_3.png" /><figcaption>Fig. 1. Example of a grayscale mask and the quantities of interest.</figcaption></figure>

quantities of interest for a pulse width modulation grayscale mask:

The density $D$, is defined as the area ratio between the chromium dot and the cell surface as follows:

$$D = \frac{S_{Cr}}{\text{pitch}^2} \qquad (1)$$

where $S_{Cr}$ is the dot surface and pitch corresponds to the cell size on Figure 1.

A chrome on glass photo-mask can be described by its transmission function $O(x, y)$. The chromium regions are opaque to the illumination thus their transmission is 0, whereas quartz regions are transparent and assumed to have a transmission of 1.

$$O(x, y) = \begin{cases} 1 & \text{if } (x, y) \text{ in quartz} \\ 0 & \text{if } (x, y) \text{ in chromium} \end{cases} \qquad (2)$$

The mask transmission is directly correlated to the density. If the density is unity, the transmission is zero (full chromium mask). For a density equals to zero, the transmission is unity (full quartz mask). For intermediate densities a model is required.

Originally, the transmission of a grayscale mask has been extrapolated from the transmission of an infinite line/space mask under coherent illumination. A line/space with a 1:1 ratio below the resolution limit always produce a constant light intensity of 0.25.<fnref n="1" /> By adjusting the ratio between the lines and spaces the average intensity can be adjusted. This observation lead Wagner and Quenzer [3], [4] to propose a similar approach for grayscale lithography where the light intensity is computed using the chromium density:

$$I = (1 - D)^2 \qquad (3)$$

We can note that for a 1:1 ratio grayscale mask the density is equal to $D = 0.5$ then $I = 0.25$. Equation (3) is used to compute the aerial image of the mask which is the light intensity impinging on the photoresist during exposure. There is a single intensity value computed for each cell on the mask, the transition between two intensity levels is abrupt.

The aerial image controls the remaining resist thickness after development. The contrast curve of the resist describes this effect from a purely vertical point of view. The contrast

<figure><img src="image_4.png" /><figcaption>Fig. 2. Theoretical contrast curves for three different resists with infinite, moderate and low contrast.</figcaption></figure>

curve is the theoretical or experimental relationship between the exposure dose<fnref n="2" /> and the remaining resist thickness after development. The contrast of the resist, written $\gamma$, is defined as the slope of the contrast curve. Figure 2 illustrates 3 different theoretical contrast curves for different resist behaviors.

If $\gamma = \infty$, the resulting contrast curve is stair-like which is highly desirable for conventional lithography because it creates a sharp transition between the non developed and the fully developed regions in the resist. For grayscale lithography, a low contrast resist is required because it offers a smooth transition between two consecutive resist heights. The ideal contrast curve for grayscale lithography would be a linear one. In reality, most contrast curves are approximated by a sigmoid function $f$ that defines the resist height after development $h$ as follows:

$$h = f(\text{dose}, \gamma) \qquad (4)$$

The contrast curve is at the center of many grayscale lithography publications for it is used as a transfer function between the resist target and the density on the mask. Mask creation with a pulse width modulation is achieved with the following steps:

1) The resist target is divided into square cells defined by the pitch in the $xy$ plane.

2) A representative height is computed for all cells (mean, median...).

3) This representative height is converted into a dose/light intensity using the contrast curve of the resist (4).

4) For all cells on the mask, the corresponding chromium dot is computed using (3).

Some approaches do not use (3) and skip step 3 by directly using the mask density [16], [17]. In such cases, a calibration mask with different chromium densities is used to obtain the transfer function between the density on the mask and the remaining resist height. This approach was improved by LeCompte with a linearization technique to limit the number of useful densities [18]. Indeed, two different densities may

<footnote n="1">¹Under coherent illumination, the resolution limit for an infinite line/space of ratio 1:1 is exactly $\lambda$, the wavelength of used light.</footnote>

<footnote n="2">²Exposure dose is simply the light intensity times the surfacic exposure energy in [J/m²].</footnote>
<table><caption>TABLE I
EXACT OPTICAL PARAMETERS OF THE PROJECTION
STEPPER USED IN THIS PAPER</caption><tr><td>Wavelength</td><td>365 nm</td><td>Coherence σ</td><td>0.7</td></tr><tr><td>Numerical aperture NA</td><td>0.57</td><td>Polarization</td><td>None</td></tr></table>

generate the same resist height. The linearization technique determines all densities that generates a unique resist height after development thus preventing redundant densities. This observation was confirmed by Waits [19]. More recently, Mosher introduced the double exposure grayscale lithography to extend the number of gray levels based on the formulation of equation (3) [7].

Yao is the first to propose an approach to account for the non linearity of the imaging process when using (3) [5]. By using a rigorous optical model the size and the position of the dots are adjusted to match a target aerial image. This aerial image is the scaled version of a microlens target between 0 and 1. The patterning of microlens with diameter of 40µm was reported by Yao.

Halftoning techniques for the mask creation have been investigated by Audran but this approach is also based on a contrast curve interpolation [20]. Experimental results do not prove the capacity of this approach to efficiently address microstructures with dimensions below 2µm.

So far, the different approaches suffer from serious limitations:

* Equation (3) does not account for the non linearity of the imaging process. This formula is valid for a coherent projection system but partially coherent illumination is used in reality. Besides, it supposes an ideal mask composed of an infinite array of dots with constant size.

* The horizontal aspect of development is not caught by the contrast curve which is essentially a vertical description of the development.

For large microstructures the aforementioned limitations do not limit the applicability of the pulse width modulation approach because the chromium dot size variation is small and over large distances on the mask. However, we will see that such approach cannot address microstructures with more challenging dimensions.

### III. MODELING OF GRAYSCALE LITHOGRAPHY

This section describes the mathematical framework used to simulate the lithographic process. The projection lithography stepper used in this article is an I-line Canon iZ+ with a circular geometry source. The main optical parameters of this equipment are presented in the Table I: Concerning the resist, a positive tone photoresist from JSR has been selected: the MFR530-R1. This photoresist is acryl based, and designed specially for grayscale applications. The experimental stack is depicted on Figure 3. It consists of a three layers structure. The substrate is a 300mm silicon wafer on top of which a first layer of JSR resist has been spin-coated and soft baked at 90°C for 90s. This layer is hard-baked at 160°C during 5 minutes.

<table><caption>Fig. 3. Experimental stack for the grayscale lithography.</caption><tr><td>JSR MFR530-R1</td></tr><tr><td>JSR MFR530-R1 (hardened)</td></tr><tr><td>Silicium</td></tr></table>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Main blocks of a projection lithography tool.</figcaption></figure>

This solidified layer is 1550nm thick. Another layer of JSR resist is spin-coated on top with a thickness of 1600nm and soft baked at 90°C for 90s. It is this top layer that will be patterned by grayscale lithography.

A TMAH solution concentrated at 1% is used for the resist development during 60 s. The first part of this section will describe the optical model used to compute the aerial image of an arbitrary grayscale mask. The second part details the modeling of the photoresist during exposure and development resulting in the final resist profile. All the equations are not only valid for grayscale lithography but also for conventional i-line lithography.

Both the optical and resist models have been implemented in Python.

#### A. Optical Model

The optical model computes the aerial image that corresponds to the light intensity received by the resist surface upon exposure. The aerial image is the result of complex interactions between an illumination, a photomask and a projection system. In this article all photo-masks are clear chrome on glass. In our case, we are considering diffracting elements on the mask that are in the order or larger than the wavelength of illumination $λ$. Consequently, the image formation will be described by the scalar diffraction theory and the Fourier Optics in the paraxial regime. A simplified view of the architecture of a projection lithography stepper is depicted on the Figure 4.

The optical elements are arranged in order to create a Köhler illumination on the mask [21]. This optical configuration is useful for several reasons:

1) Each point of the source generates a spectrum of plane waves on the mask. This ensures an uniform illumination of the mask from an angular point of view.

2) In the absence of mask the source is imaged by the condenser and collection lens at the entrance pupil of the projection lens. This image is called the effective source.

3) Upon illumination, the Fraunhofer diffraction pattern of the mask is projected at the entrance pupil of the
<figure><img src="image_6.png" /><figcaption>Fig. 5. Illustration of the Köhler illumination principle for 3 source points.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Angular acceptance of the projection pupil results in a frequency cutoff for the mask spectrum.</figcaption></figure>

projection lens due to the collection lens. Without any lens, this figure of diffraction would be projected at infinity.

The Figure 5 illustrates the light propagation for 3 sources points. The angle of incidence of the plane wave on the mask depends on the corresponding source point position regarding the optical axis. For a source point on the optical axis (green point) the illumination is parallel to the optical axis. However, as the source point moves away from the optical axis, the angle of incidence increases. The pupil plane of the projection lens is buried inside the lens in this representation.

Considering a Köhler illumination, the light intensity received at the top of the resist, $I(x, y)$, can be written as follows:

$$I(x, y) = \int_{-\infty}^{+\infty} \cdots \int_{-\infty}^{+\infty} \tilde{J}(f, g) \tilde{H}(f + f', g + g') \\ \times \tilde{H}^*(f + f'', g + g'') \tilde{O}(f', g') \tilde{O}(f'', g'') \\ \times e^{-2i\pi[(f'-f'')x+(g'-g'')y]} df dg df' dg' df'' dg'' \quad (5)$$

The derivation of this equation can be found in [22]–[24]. In its scalar form this equation is restricted to numerical aperture smaller than 0.5. However, we will see that the validity of this expression can be extended to higher numerical apertures using an obliquity factor. The tilde notation is used for functions expressed in frequency space expressed on the $f, g$ support. In the next paragraphs, all terms used in equation (5) are described.

$\tilde{O}(f, g)$ corresponds to the mask spectrum. It is obtained by Fourier Transform of the mask transmission function as follows:

$$\tilde{O}(f, g) = \mathcal{F}[O(x, y)](f, g) \quad (6)$$

In order to numerically compute the Fourier Transform the mask is expressed on a finite support. A pixel size in $x$ and $y$ direction is chosen in order to be representative of the features on the mask. The frequency support is directly computed from the spatial support. If the mask is composed of elements with existing analytical Fourier Transform (such as rectangles), the exact Fourier Transform can be computed without using Fast Fourier Transform libraries. The mask spectrum is the sum of all individual Fourier Transforms for each chromium dot. Analytical Fourier Transform are used in this paper.

The first term in the integral is $\tilde{J}(f, g)$ that represents the effective source illustrated on Figure 5. This function shouldn't be mistaken with the Fourier Transform of the source even if it is expressed in frequency space. The reason for this is that the figure of diffraction of the mask is imaged on the same plane as the effective source. Thus, it is more convenient to express both functions in the frequency space.

All the elements (points) of the effective source take part in the imaging formation. Each single point in the effective source corresponds to a spectrum of plane waves impinging on the photo-mask at a specific angle. In our case, the source is circular, its radius is defined by the coherence factor $\sigma$. The image of the source is also circular and corresponds to the radius of the real source filtered by the numerical aperture of the condenser lens that is equal to $\frac{\sigma NA}{\lambda}$. Thus, the effective source is written:

$$\tilde{J}(f, g) = \text{circ} \left( \frac{\sqrt{f^2 + g^2}}{\frac{\sigma NA}{\lambda}} \right) = \begin{cases} 1 & \text{if } \sqrt{f^2 + g^2} \le \frac{\sigma NA}{\lambda} \\ 0 & \text{otherwise} \end{cases} \quad (7)$$

The next term in equation (5) is the transfer function of the projection lens: $\tilde{H}(f, g)$. We have seen that the mask Fraunhofer diffraction pattern is projected on the pupil plane of the projection lens. This spectrum is recombined by the projection lens before being projected on the photoresist top surface. In the ideal case, the full spectrum is recombined and an exact reproduction of the mask is projected on the substrate. In reality, only a fraction of the mask spectrum can be recombined because of the finite dimension of the projection lens. The projection lens acts as a low pass filter for the mask spectrum. This behavior is described by the angle of acceptance of the projection lens, $U$, depicted on the Figure 6.

The angle of acceptance $U$ is related to the numerical aperture NA by the following formula:

$$NA = n \sin(U) \quad (8)$$

where $n$ is the optical index after the lens (air).

Small frequencies traveling close to the optical axis can go through the lens aperture unperturbed (purple light on Figure 6). Larger frequencies are cutoff by the pupil (red light on Figure 6). The frequency cutoff is governed by the projection lens cutoff frequency equals to $\frac{NA}{\lambda}$. The transfer
function is written as follows:

$$ \tilde{H}(f, g) = \text{circ} \left( \frac{\sqrt{f^2 + g^2}}{\frac{NA}{\lambda}} \right) = \begin{cases} 1 & \text{if } \sqrt{f^2 + g^2} \le \frac{NA}{\lambda} \\ 0 & \text{otherwise} \end{cases} \qquad (9) $$

The pupil function does not only describe the frequency cutoff but is also used to take into account other aspects of imaging such as aberrations.

Aberrations can be taken into account by adding a phase term to the pupil transfert function and an aberration function $W(f, g)$. The new transfer function is written as follows:

$$ \tilde{H}(f, g) = \begin{cases} e^{ikW(f,g)} & \text{if } \sqrt{f^2 + g^2} \le \frac{NA}{\lambda} \\ 0 & \text{otherwise} \end{cases} \qquad (10) $$

$W(f, g)$ can be defined using the Zernike polynomials [25]. For this article, we will only account for the defocus aberration that corresponds to a distance offset between the focus plane and the plane of best focus. In reality, this corresponds to a vertical shift of the substrate inside the stepper. For a defocus of length $D$, the aberration function is the following [26]:

$$ W(f, g) = D\sqrt{1 - (f^2 + g^2)^2} \qquad (11) $$

The defocus aberration affects the higher order of the mask spectrum.

Until this point we have been working in the paraxial approximation. There is no consensus on the domain of validity for this approximation. Paraxial approximation implies that rays travel close to the optical axis. The proximity of rays to the optical axis is directly related to the numerical apertures of the system. Paraxial approximation is assumed to give correct results for numerical aperture up to 0.5. However, larger numerical aperture require to switch to a vectorial description of light since polarization effects cannot be ignored anymore.

It is possible to extend the limits of paraxial approximation to higher numerical apertures using a correction term: the obliquity factor. The obliquity correction is defined at the pupil plane of the projection lens even if it is not an aberration. This factor has been introduced by Cole [27] and discussed by Gallatin [28], to extend the paraxial approximation to numerical apertures up to 0.6. The obliquity factor is noted $K$ and expressed as follows:

$$ K(f, g) = \left[ \frac{1 - \lambda^2 M^2 (f^2 + g^2)}{1 - \lambda^2 (f^2 + g^2)} \right]^{\frac{1}{4}} \qquad (12) $$

where $M$ is the stepper magnification (0.25 in our case). The higher the diffraction order, the higher the correction. The maximum correction is obtained at the boundary of the pupil plane, where $\sqrt{f^2 + g^2} = \frac{NA}{\lambda}$. This obliquity factor is added to the pupil transfer function as follows:

$$ \tilde{H}(f, g) = \begin{cases} K(f, g)e^{ikW(f,g)} & \text{if } \sqrt{f^2 + g^2} \le \frac{NA}{\lambda} \\ 0 & \text{otherwise} \end{cases} \qquad (13) $$

Equation (5) cannot be solved directly. Two implementations are possible to solve this equation: the Abbe and the Hopkins methods.

The Abbe's method [29] (also known as source point integration) consists in transforming the incoherent optical system in a sum of coherent optical systems. This transformation is possible since the source is not spatially coherent (it has only a temporal coherence because it is assumed to be monochromatic). The integration is carried on the effective source points. The Abbe's formulation of (5) is written:

$$ I(x, y) = \iint \sqrt{f^2 + g^2} \le \frac{NA}{\lambda} \tilde{J}(f, g) |E(f, g)|^2 df dg \qquad (14) $$

where $|E(f, g)|^2$ is the intensity produced by a source point of unit strength located at $(f, g)$:

$$ E(f, g) = \iint_{-\infty}^{+\infty} \tilde{H}(f + f', g + g') \tilde{O}(f', g') \times e^{-2i\pi[(f'x + g'y)]} df' dg' \qquad (15) $$

The aerial image is the sum of coherent images produced by all source points of the effective source.

The different numerical steps to compute the Abbe's method are given below:

1) Start by computing the Fourier Transform of the mask and the associated frequency support. Depending on the mask features one can use analytical formulas or a numeric Fast Fourier Transform.

2) Compute the effective source on the frequency support. In 2 dimensions, the effective source is a list of frequency pairs that respect the cutoff frequency of equation (7).

3) For a source point located in $(f, g)$ compute the shifted pupil function $H(f + f', g + g')$. This shifted pupil function is multiplied by the mask spectrum (that has not been shifted): $\tilde{H}(f + f', g + g')\tilde{O}(f', g')$.

4) Compute the numerical inverse Fourier transform of the product $\tilde{H}(f + f', g + g')\tilde{O}(f', g')$.

5) Compute the square of the absolute value of the inverse Fourier Transform.

6) Repeat steps 3 to 5 for all effective source points.

7) The aerial image is the sum of all coherent images contributions divided by the number of coherent images.

Hopkins method [24], [30] is more recent. It is based on the introduction of the transmission cross coefficients to separate the source contribution from the mask contribution. The Hopkins formulation is written as follows:

$$ I(x, y) = \iint \int \int \int \text{TCC}(f', g'; f'', g'') \tilde{O}(f', g') \tilde{O}^*(f'', g'') \times e^{-2i\pi[(f'-f'')x+(g'-g'')y]} df' dg' df'' dg'' \qquad (16) $$

where TCC are the transmission cross coefficients. These coefficients are used to described the optical system as a
<figure><img src="image_8.png" /><figcaption>Fig. 7. Transmission cross coefficient computation. (7a) No overlapping between source and pupils function results in a zero TCC coefficient while overlapping (7b) results in a non zero TCC coefficient.</figcaption></figure>

4 dimensional system. The TCC(f', g'; f'', g'') is written as:

$$TCC(f', g'; f'', g'') = \iint_{-\infty}^{+\infty} \tilde{J}(f, g) \tilde{H}(f + f', g + g') \times \tilde{H}^*(f + f'', g + g'') df dg \quad (17)$$

Calculating the TCC consists in exploring all possible positions for the two offset pupils $\tilde{H}(f + f', g + g')$ and $\tilde{H}(f + f'', g + g'')$ to find an overlap with the effective source function. This integration is illustrated graphically on the Figure 7.

Once a non zero TCC coefficient is obtained, it is multiplied by the value of the mask spectrum at the position (f', g') and the value of the conjugate mask spectrum at the position (f'', g''). The exponential term brings back the product to real space. The different numerical steps to compute the Hopkins method are given below:

1) Start by computing the Fourier Transform of the mask and the associated frequency support. Depending on the mask features one can use analytical formulas or a numeric Fast Fourier Transform.

2) Compute the TCC coefficients on the frequency support for all possible quadruplets (f', g', f'', g''). Note that the TCC coefficients can be complex in the presence of aberrations.

3) For all non zero TCC coefficients compute the product of equation (16).

4) The aerial image is obtained by summing of all non zero contribution.<fnref n="3" />

Both Abbe and Hopkins approaches are valid to compute the aerial image of a grayscale mask. Hopkins approach is the most commonly used optical model for optical proximity correction (OPC) because the TCC need to be calculated only once for a given illumination. For other applications such as source mask optimization, the Abbe's approach is preferred. In our case, the Hopkins method requires more complex numerical treatment to be more efficient than the Abbe's approach. Consequently, the Abbe's method has been selected.

## B. Exposure Model

During exposure, the aerial image computed with the optical model propagates into the photoresist, altering the chemical properties of the polymer. Exposure of an i-line photoresist has been described by a series of paper issued by Dill in the mid seventies [31]–[34]. A typical I-line photoresist is composed of the following elements:

* The Novolak polymer slightly soluble in the development solution.

* The photo-acid compound (PAC) that is fixed to the Novolak polymer and totally inhibits its dissolution.

* The solvent that binds the Novolak and PAC together.

During exposure, the PAC is converted into carboxylic acid increasing the novolak dissolubility. From an optical point of view the simulation of exposure consists in monitoring the evolution of the photoresist film absorption over time: $\alpha(x, y, z, t)$. Since the upper face of the resist is the first to be exposed, an absorption gradient is created inside the film: this effect is called bleaching.

The evolution of the film absorption is described by a 3 parameters model according to Dill's publications with the famous A, B and C parameters. The absorption value in the resist before exposure is equal to A + B. The absorption of a fully exposed resist is B. The C parameter corresponds to the PAC conversion rate. The absorption of the resist is directly related to the PAC concentration. The normalized PAC concentration is used: $m(x, y, z, t)$. It corresponds to the PAC concentration at a time t relative to the concentration at t = 0 (before exposure). $m(x, y, z, t)$ is also called the latent image. It is expressed in the [0, 1] interval. The absorption coefficient is written:

$$\alpha(x, y, z, t) = A \cdot m(x, y, z, t) + B \quad (18)$$

With this expression, before exposure, $m(x, y, z, t) = 1$ thus $\alpha(x, y, z, t) = A + B$. For a fully exposed resist $m(x, y, z, t = \infty) = 0$ thus $\alpha(x, y, z, t) = B$.

The first step in order to model exposure is to compute the 3D intensity inside the resist: $I(x, y, z)$. This intensity is called the bulk image. The simplest approach is to consider a vertical absorption of the aerial image with a Beer-Lambert model:

$$I(x, y, z; \alpha) = I(x, y, z = 0)e^{-\alpha z} \quad (19)$$

With this approach only one aerial image has to be computed. However it lacks precision because it does not account for the defocus of the aerial image inside the resist film. This defocus can be taken into account by using the scaled-defocus algorithm [26]. The aerial image at a position z inside the photoresist film is obtained by computing the aerial image with a scaled defocus:

$$I_{SD}(x, y, z; \alpha) = I\left(x, y, D - \frac{z}{n_{\text{resist}}}; \alpha\right) \quad (20)$$

where $n_{\text{resist}}$ is the resist optical index and D the defocus of the stepper. This formulation is more accurate, however it requires an important computational effort in comparison to the vertical absorption because multiple aerial images need to be computed.

<footnote n="3">³Depending on the convention, one need to normalize the result by the aerial image of a pure quartz mask.</footnote>
We decided to use the simple vertical absorption model instead of the scaled defocus algorithm to ensure fast computational time. We also neglected the reflection at the air/resist interface. The pixel size along the vertical direction has been set to 25nm resulting in 64 “resist slices”. It takes 64 times longer to compute the bulk image with the scaled defocus algorithm than with the simple vertical absorption.

The evolution of the latent image during exposure is related to the evolution of the bulk image by the Dill's equation [32]:

$$ \frac{\partial m(x, y, z, t)}{\partial t} = -C \cdot I(x, y, z, t) \cdot m(x, y, z, t) \quad (21) $$

$$ \frac{\partial I(x, y, z, t)}{\partial z} = -[A \cdot m(x, y, z, t) + B] I(x, y, z, t) \quad (22) $$

It is actually possible to compute an analytical solution to solve this set of equations. However, such solution requires to neglect the bleaching of the resist ($B=0$) and the absence of defocus of the aerial image. Interested reader can find the analytical solution in [35], [36].

The partial derivative must be solved iteratively taking into account the kinetic aspects of the exposure. The exposure time is divided into small time steps $\Delta t$. Latent and bulk images are updated one after the other assuming the other one does not change during the integration time $\Delta t$. With such hypothesis it is possible to integrate (21) and (22) as follows:

$$ m(x, y, z, t + \Delta t) = m(x, y, z, t)e^{-C \cdot I(x, y, z, t) \Delta t} \quad (23) $$

$$ I(x, y, z, t + \Delta t) = I(x, y, z; \alpha = A \cdot m(x, y, z, t) + B) \quad (24) $$

Of course the light intensity kinetic is faster than the chemical kinetic of the latent image. The bulk image at a time $t + \Delta t$ does not depend on the bulk image at an earlier time $t$ but only on the latent image at this time. On the other side, the latent image at $t + \Delta t$ requires both latent and bulk image at $t$.

The numerical steps to model exposure are the following:

1) Instantiate exposure time to zero: $t \to 0$

2) Compute the initial bulk image, $I(x, y, z, t = 0)$ with (19) or (20).

3) Update the latent image with (23).

4) Update the bulk image with (24).

5) Increment the exposure time: $t \to t + \Delta t$

6) Repeat steps 3 to 5 until $t = t_{\text{expo}}$

At the end of the exposure simulation, the final bulk and latent images are obtained: $I(x, y, z, t = t_{\text{expo}})$ and $m(x, y, z, t = t_{\text{expo}})$. The next step is to compute the resist development.

### C. Development Model

Depending on the latent image at the end of exposure, $m(x, y, z, t = t_{\text{expo}})$, the resist is selectively removed during the development step resulting in the final profile.

The simplest approach to obtain the resist profile is to use a threshold development model. With this approach the resist profile corresponds to the envelope of the final latent image at a specified value called threshold. The threshold is written $m_{th}$ and is between 0 and 1. This is the fastest technique to compute a resist profile after exposure. However

this method does not account for all the complexity behind resist development.

Sethian [37] proposed to model the action of the developer by the propagation of a front perpendicular to the resist surface. The front position at any time in the resist is governed by the eikonal equation:

$$ |\nabla T(x, y, z)| R(x, y, z) = 1 \quad (25) $$

where $R(x, y, z)$ is the front speed in nm/s also called the development rate. $T(x, y, z)$ is the time required for the front to reach the point $(x, y, z)$. Equation (25) can be rewritten as follows:

$$ |\nabla T(x, y, z)| = \frac{1}{R(x, y, z)} \quad (26) $$

or:

$$ \sqrt{\left(\frac{\partial T}{\partial x}\right)^2 + \left(\frac{\partial T}{\partial y}\right)^2 + \left(\frac{\partial T}{\partial z}\right)^2} = \frac{1}{R(x, y, z)} \quad (27) $$

The surface $T_{dev}(x, y, z, t = t_{dev})$ gives the resist profile at the end of development. Many different development rate functions exist but we only consider here the 4 parameters Mack model [38] function (resist supplier provided with the corresponding model parameters) written as follows:

$$ R(x, y, z) = R_{max} \frac{(a+1)(1-m(x, y, z))^n}{a+(1-m(x, y, z))^n} + R_{min} \quad (28) $$

where $R_{max}$ is the maximum development rate, $n$ the developer selectivity and $a$ a constant that depends on the threshold concentration $m_{th}$ (value of $m$ at the inflection point of the development rate function) and the developer selectivity $n$:

$$ a = \frac{n+1}{n-1}(1-m_{th})^n \quad (29) $$

The resolution of the eikonal equation is not straightforward. It is possible to neglect the partial derivative along $x$ and $y$ direction. The development is assumed to be purely vertical. The resolution of (26) is simply an integral along $z$ direction.

In order to account for all three directions, [39] developed an algorithm called the fast marching method. This technique is already implemented in open source libraries such as the sci-kit fast marching method. This solution has been used for the model. Interested reader may find additional information about mathematical implementation of the fast marching method in [37].

### D. Model Calibration

The development model parameters are given by the resist supplier. However, in order to match our development process they have been calibrated using a grayscale calibration mask. Model calibration has already been discussed in [40]. The calibrated parameters are given in Table II.

## IV. GRAYSCALE MASK OPTIMIZATION

The lithography model described in the last section can be used to simulate the resist profile for a grayscale mask, illumination and process conditions. In this section, we discuss the model inversion in order to find the optimum grayscale
<table><caption>TABLE II
CALIBRATED PARAMETERS FOR THE MACK DEVELOPMENT MODEL</caption><thead><tr><th>Parameters</th><th>r<sub>max</sub></th><th>r<sub>min</sub></th><th>m<sub>th</sub></th><th>n</th></tr></thead><tbody><tr><td>Values</td><td>51.03</td><td>3.03</td><td>0.0022</td><td>3.11</td></tr></tbody></table>

mask for a target resist profile. Starting from a neutral mask, both the size and the position of the chromium dots on the mask will be optimized. This optimum mask should generate a resist profile as close as possible to the profile target.

This inversion problem can be computed by iteratively solving the direct problem. This approach requires a numeric optimizer that can efficiently explore the parameters space (size and positions). The lithographic model computation times are of paramount importance because the optimizer routine will make many calls to the model during the optimization.

The total mask optimization is a sequence of two optimizations. The idea is to divide the mask optimization in two steps. Firstly by optimizing the mask with a target aerial image (raw but fast optimization) and switch to a resist profile optimization (slow but accurate optimization). The first optimization only requires the optical model to work while the second optimization requires the complete model including the development model.

The prerequisite for such optimization to work is to be able to find an aerial image that produces a resist profile as close as possible to the resist profile target. Indeed, the idea of the first optimization is to save time for the second optimization by returning a mask that is already close to the final result. Determination of the aerial image is described in the next section.

### A. Resist Profile Target to Aerial Image Conversion Routine

This routine requires to measure or simulate the contrast curve of the photoresist. The target resist profile is interpolated with the contrast curve to produce a first aerial image. Such aerial image could be used directly for the optimization but it is possible to improve it.

This aerial image is iteratively corrected by looping on the exposure and development model. If the altitude of a point in the simulated profile is too high, the corresponding point in the aerial image is increased and vice versa. The correction function is a sigmoid function in order to apply a smooth correction (the maximum correction does not exceed 1% of the pixel each time). This correction loop stops if the maximum number of iterations is reached (500), or if the aerial image stops being corrected (smallest difference with the target profile).

The Figure 8 presents the results for a microlens of radius 1400nm with a period of 1800nm. The resulting aerial image is the black dotted lines. The resist target is the black line and the simulated resist profile with the resulting aerial image is the red curve.

The simulated profile is not exactly matching the target resist profile. However, it is close enough to use the aerial image and start the mask optimization. One can notice that there are two important discontinuities in the aerial image for

<figure><img src="image_9.png" /><figcaption>Fig. 8. Aerial image obtained for a microlens resist profile. Simulated resist profile with the aerial image and a simple vertical absorption (19).</figcaption></figure>

$x \pm 1400nm$. These coordinates correspond to the edge of the microlens where the spatial derivative is not defined. It is already understandable that such behavior will be difficult to reproduce with a grayscale mask.

### B. Mask Optimization

The optimization of the mask is also implemented in Python. The scipy least squares function is used for it is suited to local optimization problems. Besides, it offers the possibility to fix boundaries to the variables being optimized. The initial guess for the optimizer is written $x_0$, it is a list composed of the size and position for all polygons of the mask.

$$x_0 = [s_1, x_1, y_1, s_2, x_2, y_2, \dots] \quad (30)$$

where $s_i$ is the size, $x_i$ the x position and $y_i$ the y position of the $i^{th}$ dot on the mask. The starting point of the optimization is a grayscale mask made of regularly spaced chromium dots along x and y directions. The initial dot size is set to half the pitch.

For each dot on the mask, both the size and the position have to be optimized. Thus, the number of parameters to consider is equal to three times the number of dots. This number can grow very large for complex masks and results in extended computational time. It is possible to relieve the computational effort by taking into account the symmetries of the resist target.

Indeed, if the targeted resist profile presents a certain degree of symmetry in the xy plane, this symmetry can also be found at the mask level. In other words there are some dots on the mask that are identical from the target resist profile perspective. For a centro-symmetric shape, such as a microlens, only an eighth of the total dots is required to generate the complete mask. This principle is illustrated on the Figure 9 which illustrates the mask symmetries for the microlens targeted resist profile.

Only the dots overlapping the blue area are required to construct the complete mask. For illustration purposes, the dots circle in red represent a 1/8 symmetry whereas the dots circled in green represent a 1/4 symmetry.
<figure><img src="image_10.png" /><figcaption>Fig. 9. Mask symmetries for the microlens target. The blue area corresponds to 1/8 of the total mask.</figcaption></figure>

During optimization the Mask Design Rules (MRC) have to be respected. The most important rules to respect are the following:

1) The minimum CD that corresponds to the smallest feature that the maskshop can write.

2) The minimum space that corresponds to the smallest space between two features that the maskshop can write.

Respecting the first rule is straightforward. One only needs to provide this value as the lower boundary for the size parameter when setting the optimizer. However, keeping the space rule respected is more complex and a specific solution has been implemented. During the optimization process, every time an intermediate mask is generated, a buffer is computed for all dots on the mask. This buffer corresponds to the points located at a distance smaller or equal to half the minimum space distance for the specific dot.

Buffer are used to control the respect of the minimum space condition as illustrated on the Figure 10.

During optimization the total area of the mask buffers is computed as well as the area of the union of all buffers. The ratio written $r_{MRC}$ is computed:

$$r_{MRC} = \frac{\sum_n A_{b_n}}{\bigcup_n A_{b_n}} \qquad (31)$$

where $A_{b_n}$ is the area of the $n^{th}$ buffer.

If there is no overlap this ratio is equal to one: all dots respect the minimum space condition. When two or more buffers overlap that means that some dots on the mask do not respect the minimum space condition: the coefficient becomes greater than one. The $r_{MRC}$ ratio is used to modify the cost function during optimization.

The Table III summarizes the cost function $C_i$ at the iteration $i$, for both aerial image and resist profile optimization. $I_{tgt}$ and $P_{tgt}$ are the aerial image and the targeted resist profile respectively. $I_i$ and $P_i$ are the simulated aerial image and resist profile at iteration $i$ during optimization respectively. If the minimum space is not respected, the cost function is equal to the cost function of the previous iteration, $C_{i-1}$ times the $r_{MRC}$ ratio. By doing this we “inform” the optimizer

<figure><img src="image_11.png" /><figcaption>Fig. 10. Buffer creation for minimum space verification during mask optimization.</figcaption></figure>

<table><caption>TABLE III
COST FUNCTIONS FOR THE AERIAL IMAGE AND RESIST PROFILE OPTIMIZATION</caption><thead><tr><th>Optimization</th><th>Cost function</th></tr></thead><tbody><tr><td>Aerial image</td><td>C<sub>i</sub> = {√(I<sub>tgt</sub> - I<sub>i</sub>)<sup>2</sup> if r<sub>MRC</sub> = 1<br>C<sub>i-1</sub> · r<sub>MRC</sub> if r<sub>MRC</sub> > 1</td></tr><tr><td>Resist profile</td><td>C<sub>i</sub> = {√(P<sub>tgt</sub> - P<sub>i</sub>)<sup>2</sup> if r<sub>MRC</sub> = 1<br>C<sub>i-1</sub> · r<sub>MRC</sub> if r<sub>MRC</sub> > 1</td></tr></tbody></table>

not to explore this parameters space direction since it is not compliant with the design rules.

The Figure 11 presents the optimised mask resulting from the aerial optimization of the microlens of Figure 8.

This mask is used as starting point for the resist profile optimization. The Figure 12 presents the mask resulting from this final optimization for the microlens.

Interestingly, the dots on the mask naturally arranged in a circular shape coherent with the microlens target. The optimization quality can be assessed by simulating the resist profile obtained with the mask on Figure 12. The Figure 13 presents horizontal and diagonal cutlines for the target and simulated microlens.

The simulated resist profile is close to the target. The final cost function is equal to 14.96 nm. We can note the presence
<figure><img src="image_12.png" /><figcaption>Fig. 11. Optimized mask after aerial image optimization for the microlens target.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Optimized mask after resist profile optimization for the microlens target.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Target and simulated cutlines for the microlens obtained with the mask of Figure 12.</figcaption></figure>

of oscillations at the resist surface and the difficulty to obtain smooth slopes on the sides of the microlens.

## V. EXPERIMENTAL RESULTS

Using the mask data-preparation algorithm a photo-mask has been designed with hundreds of different pixelizations to

<figure><img src="image_15.png" /><figcaption>Fig. 14. SEM image of the microlens array.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 15. AFM scan for the microlens pattern.</figcaption></figure>

address a wide variety of microstructures. In this section we present some experimental results obtained with this mask on the experimental stack described on Figure 3. The exposure dose has been set to 1500 J/m² and no defocus is applied. For all patterns, a SEM image and an AFM profile extracted from the central chip will be provided. The SEM is a G3 SEMVision from Applied Material and the AFM is an Insight 3D from Bruker.

Experimental cutlines of the horizontal and diagonal axis are extracted from the AFM profiles. Each cutline corresponds to the mean envelope extracted from 5 different chips on the wafer: the central chip and its four direct neighbors (left, right, up and down). The envelope is an indication on the profile reproducibility.

### A. Microlens

The Figure 14 is a SEM image of the microlens array obtained with the mask of Figure 12 showing the presence of standing waves:

The Figure 15 is an AFM scan of a microlens.

The scanning direction is along the y axis. We can notice that the microlens is bridged with its four direct neighbors. Horizontal and diagonal cutlines for the microlens pattern are extracted and plotted on the Figure 16.
<figure><img src="image_17.png" /><figcaption>Fig. 16. Horizontal (top) and diagonal (bottom) cutline for the microlens pattern.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 17. Grayscale mask for the pyramid target.</figcaption></figure>

The cutlines confirm that the microlens profiles are bridged with the direct neighbors at an altitude approximately equals to 250 nm. This effect has been underestimated by the lithographic model. However, it is not visible for the diagonal cutline because there is more space between the neighbors for the resist to be removed.

## B. Pyramid

A pyramid pattern has been integrated with a square base of 2800 nm, a height of 1400 nm and a period of 2800 nm. The mask resulting from the optimization is illustrated on the Figure 17.

<figure><img src="image_19.png" /><figcaption>Fig. 18. SEM image of the pyramid array.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 19. AFM scan for the pyramid pattern.</figcaption></figure>

The Figure 18 is a SEM image of the pyramid array obtained with the mask of Figure 17 showing the presence of standing waves.

The Figure 19 is an AFM scan of a pyramid.

Horizontal and diagonal cutlines for the pyramid pattern are extracted and plotted on the Figure 20.

The junction with the eight neighbors is higher than expected in all directions. Over dosing could be a solution but the overall profile decrease which is not desirable. All the resist has not been properly removed.

## C. Inverted Microlens

An inverted microlens has been integrated with a radius of 1400 nm and period of 3200 nm. It basically corresponds to the extrusion of the microlens pattern. The mask resulting from the optimization is illustrated on the Figure 21.

The Figure 22 is an AFM scan of an inverted microlens.

Horizontal and diagonal cutlines for the inverted microlens pattern are extracted and plotted on the Figure 23.

Experimental profiles for this pattern are in very good agreement with the model and the target for the diagonal cutline. Figure 23 can be compared with Figure 16. It would seem that it is easier to pattern an inverted microlens rather than a microlens itself.
<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 20. Horizontal (top) and diagonal (bottom) cutline for the microlens pattern.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 21. Grayscale mask for the inverted microlens.</figcaption></figure>

### D. Dual Microlens Configuration

Finally, a hybrid configuration with two different microlens has been integrated. It is composed of a small microlens with radius 700 nm and a large microlens with a radius of 1400 nm.

This pattern has been selected because it represents a configuration where two microstructures with different heights are patterned at the same time in a dense configuration. The mask resulting from the optimization is illustrated on the Figure 25.

The Figure 26 is a SEM image of the dual microlens array obtained with the mask of Figure 25.

The Figure 27 is an AFM scan of a dual microlens array.

<figure><img src="image_24.png" /><figcaption>Fig. 22. AFM scan of an inverted microlens pattern.</figcaption></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /><figcaption>Fig. 23. Horizontal (top) and diagonal (bottom) cutlines for the inverted microlens pattern.</figcaption></figure>

Horizontal and diagonal cutlines for the dual microlens configuration are extracted and plotted on the Figure 28.

For this challenging structure it is not possible to achieve a zero-gap between two adjacent microlens. This experimental cutlines are in good agreement with the model but highlight its limitations at the same time.

### E. Discussion

We have confirmed the capacity of grayscale lithography to pattern microstructures of challenging dimensions using a rigorous lithography model and an new mask optimization approach. Multiple conclusions can be drawn from the experimental results presented here:

1) The presence of standing waves is totally normal because there is no anti-reflective layer in the stack. The smooth-
<figure><img src="image_27.png" /><figcaption>Fig. 24. Profile target for the dual microlens configuration.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Fig. 25. Grayscale mask for the dual microlens configuration.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Fig. 26. SEM image of the dual microlens array.</figcaption></figure>

ing of the standing waves could be achieved during the heat treatment which is an essential step for freezing the resist profile after development.

<figure><img src="image_30.png" /><figcaption>Fig. 27. AFM scan of an hybrid microlens pattern with two microlens geometries.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>Fig. 28. Horizontal (top) and diagonal (bottom) cutline for the dual microlens pattern.</figcaption></figure>

2) The junction with direct neighbors has not been anticipated by the model. The simplest reason for that is the fact that defocus of the aerial image during its propagation in the resist thickness is not taken into account when simulating the exposure with the simple vertical absorption model (19). The bulk image intensity at the bottom of the resist is overestimated which results in the resist being developed. One should consider using the scaled defocus algorithm (20) to compute the bulk image in order to improve the model predictability. However, switching to the scaled defocus will increase the computational effort and thus the optimization time.

3) Interestingly, it is easier to pattern an extrusion of a microlens rather than a microlens itself. This observation lead us to think to that using a negative tone photoresist could improve some experimental results.

## VI. CONCLUSION

In this paper we have demonstrated the capacity of grayscale lithography for the patterning of complex photoresist microstructures with different heights and dimensions ranging between 1 to 3µm. A new mask data-preparation algorithm has been developed in order to compute the most suitable grayscale mask for any arbitrary microstructure target. These
achievements have been made possible by introducing new considerations for the mask creation. In contrast to existing techniques relying on simplified optical and resist models, our approach is based on rigorous lithographic simulations. Where classical approaches only optimize the chromium dot sizes on a fixed grid, we have proposed an optimization technique that handles both the size and the position of the chromium dots. There is still room for improvements as the experimental profiles show a deviation with the model in the lowest parts of the microstructures suggesting an important effect of defocus that should be accounted for to improve the results.

## ACKNOWLEDGMENT

The authors would like to thank Rémi Coquand for his time and implication in the AFM metrology campaigns. They also are grateful to Xavier Buch and Jérôme Serrand from JSR Micro for providing the materials and resist parameters.

## REFERENCES

[1] V. Nock and R. J. Blaikie, "Fabrication of optical grayscale masks for tapered microfluidic devices," *Microelectron. Eng.*, vol. 85, nos. 5-6, pp. 1077-1082, May 2008.

[2] B. Morgan, C. M. Waits, J. Krizmanic, and R. Ghodssi, "Development of a deep silicon phase Fresnel lens using gray-scale lithography and deep reactive ion etching," *J. Microelectromech. Syst.*, vol. 13, no. 1, pp. 113-120, Feb. 2004.

[3] B. Wagner, H. J. Quenzer, W. Henke, W. Hoppe, and W. Pilz, "Microfabrication of complex surface topographies using grey-tone lithography," *Sens. Actuators A, Phys.*, vol. 46, nos. 1-3, pp. 89-94, Jan. 1995.

[4] H. J. Quenzer, W. Henke, W. Hoppe, W. Pilz, K. Reimer, and B. Wagner, "Fabrication of relief-topographic surfaces with a one-step UV-lithographic process," *Microsyst. Technol.*, vol. 1, no. 4, pp. 196-201, Sep. 1995.

[5] J. Yao et al., "Coding gray-tone mask for refractive microlens fabrication," *Microelectron. Eng.*, vol. 53, nos. 1-4, pp. 531-534, Jun. 2000.

[6] C. Waits, B. Morgan, M. Kastantin, and R. Ghodssi, "Microfabrication of 3D silicon MEMS structures using gray-scale lithography and deep reactive ion etching," *Sens. Actuators A, Phys.*, vol. 119, no. 1, pp. 245-253, Mar. 2005.

[7] L. Mosher, C. M. Waits, B. Morgan, and R. Ghodssi, "Double-exposure grayscale photolithography," *J. Microelectromech. Syst.*, vol. 18, no. 2, pp. 308-315, Apr. 2009.

[8] F. Lima, I. Khazi, U. Mescheder, A. C. Tungal, and U. Muthiah, "Fabrication of 3D microstructures using grayscale lithography," *Adv. Opt. Technol.*, vol. 8, nos. 3-4, pp. 181-193, Jun. 2019.

[9] I. Khazi, U. Muthiah, and U. Mescheder, "3D free forms in c-Si via grayscale lithography and RIE," *Microelectron. Eng.*, vol. 193, pp. 34-40, Jun. 2018.

[10] J. Loomis, D. Ratnayake, C. McKenna, and K. M. Walsh, "Grayscale lithography—Automated mask generation for complex three-dimensional topography," *J. Micro/Nanolithography, MEMS, MOEMS*, vol. 15, no. 1, Art. no. 013511, 2016.

[11] Q. Zheng, J. Zhou, Q. Chen, L. Lei, K. Wen, and Y. Hu, "Rapid prototyping of a dammann grating in DMD-based maskless lithography," *IEEE Photon. J.*, vol. 11, no. 6, pp. 1-10, Dec. 2019.

[12] Q. Deng, Y. Yang, H. Gao, Y. Zhou, Y. He, and S. Hu, "Fabrication of micro-optics elements with arbitrary surface profiles based on one-step maskless grayscale lithography," *Micromachines*, vol. 8, no. 10, p. 314, Oct. 2017.

[13] N. Luo and Z. Zhang, "Fabrication of a curved microlens array using double gray-scale digital maskless lithography," *J. Micromech. Microeng.*, vol. 27, no. 3, Mar. 2017, Art. no. 035015.

[14] X. Ma et al., "Experimental study of numerical optimization for 3-D microstructuring using DMD-based grayscale lithography," *J. Microelectromech. Syst.*, vol. 24, no. 6, pp. 1856-1867, Dec. 2015.

[15] K. Totsu, K. Fujishiro, S. Tanaka, and M. Esashi, "Fabrication of three-dimensional microstructure using maskless gray-scale lithography," *Sens. Actuators A, Phys.*, vols. 130-131, pp. 387-392, Aug. 2006.

[16] Y. Oppliger, P. Sixt, J. M. Stauffer, J. M. Mayor, P. Regnault, and G. Voirin, "One-step 3D shaping using a gray-tone mask for optical and microelectronic applications," *Microelectron. Eng.*, vol. 23, nos. 1-4, pp. 449-454, Jan. 1994.

[17] D. R. Purdy, "Fabrication of complex micro-optic components using photo-sculpting through halftone transmission masks," *Pure Appl. Opt., J. Eur. Opt. Soc. A*, vol. 3, no. 2, pp. 167-175, Mar. 1994.

[18] M. LeCompte, X. Gao, and D. W. Prather, "Photoresist characterization and linearization procedure for the gray-scale fabrication of diffractive optical elements," *Appl. Opt.*, vol. 40, no. 32, pp. 5921-5927, 2001.

[19] C. M. Waits, A. Modafe, and R. Ghodssi, "Investigation of gray-scale technology for large area 3D silicon MEMS structures," *J. Micromech. Microeng.*, vol. 13, no. 2, p. 170, 2002.

[20] S. Audran, "Etude des propriétés et des mécanismes de mise en forme des résines photolithographiques pour une application capteur d'images CMOS avancés," Ph.D. dissertation, Université Louis Pasteur, Strasbourg, France, 2007.

[21] A. Köhler, "Ein neues beleuchtungsverfahren für mikrophotographische zwecke," *Zeitschrift Für Wissenschaftliche Mikroskopie Und Für Mikroskopische Technik*, vol. 10, no. 4, pp. 433-440, 1893.

[22] A. K. Wong, *Resolution Enhancement Techniques in Optical Lithography*. Bellingham, WA, USA: SPIE, Mar. 2001.

[23] A.-K. Wong, *Optical Imaging in Projection Microlithography*. Bellingham, WA, USA: SPIE, 2005.

[24] J. W. Goodman, *Statistical Optics 2E*. Hoboken, NJ, USA: Wiley, 2015.

[25] M. Born et al., *Principles of Optics*. Cambridge, U.K.: Cambridge Univ. Press, 1999.

[26] C. Mack, *Fundamental Principles of Optical Lithography: The Science of Microfabrication*. Hoboken, NJ, USA: Wiley, 2008.

[27] D. C. Cole, "Extending scalar aerial image calculations to higher numerical apertures," *J. Vac. Sci. Technol. B, Microelectron. Nanometer Struct.*, vol. 10, no. 6, p. 3037, Nov. 1992.

[28] G. M. Gallatin, "High-numerical-aperture scalar imaging," *Appl. Opt.*, vol. 40, no. 28, p. 4958, Oct. 2001.

[29] E. Abbe, "Beiträge zur theorie des mikroskops und der mikroskopischen wahrnehmung," *Archiv Für mikroskopische Anatomie*, vol. 9, no. 1, pp. 413-418, 1873.

[30] H. H. Hopkins, "On the diffraction theory of optical images," *Proc. Roy. Soc. A*, vol. 217, no. 1130, pp. 408-432, May 1953.

[31] F. H. Dill, A. R. Neureuther, J. A. Tuttle, and E. J. Walker, "Modeling projection printing of positive photoresists," *IEEE Trans. Electron Devices*, vol. 22, no. 7, pp. 456-464, Jul. 1975.

[32] F. H. Dill, "Optical lithography," *IEEE Trans. Electron Devices*, vol. ED-22, no. 7, pp. 440-444, Jul. 1975.

[33] K. L. Konnerth and F. H. Dill, "In-situ measurement of dielectric thickness during etching or developing processes," *IEEE Trans. Electron Devices*, vol. ED-22, no. 7, pp. 452-456, Jul. 1975.

[34] T. Dillon, M. Zablocki, J. Murakowski, and D. Prather, "Processing and modeling optimization for grayscale lithography," *Proc. SPIE*, vol. 6923, Apr. 2008, Art. no. 6923B.

[35] C. E. Herrick, "Solution of the partial differential equations describing photodecomposition in a light-absorbing matrix having light-absorbing photoproducts," *IBM J. Res. Develop.*, vol. 10, no. 1, pp. 2-5, Jan. 1966.

[36] S. V. Babu and E. Barouch, "Exact solution of Dill's model equations for positive photoresist kinetics," *IEEE Electron Device Lett.*, vol. ED-7, no. 4, pp. 252-253, Apr. 1986.

[37] J. A. Sethian, "Fast-marching level-set methods for three-dimensional photolithography development," in *Optical Microlithography IX*, G. E. Fuller, Ed. Bellingham, WA, USA: SPIE, Jun. 1996.

[38] C. A. Mack, *Modeling Solvent Effects in Optical Lithography*. Austin, TX, USA: Univ. Texas at Austin, 1998.

[39] J. A. Sethian, "Level set methods: Evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science," *Cambridge Monographs on Applied and Computational Mathematics*. Cambridge, U.K.: Cambridge Univ. Press, 1996.

[40] P. Chevalier et al., "Beyond contrast curve approach: A grayscale model applied to sub-5µm patterns," *Proc. SPIE*, vol. 10958, Mar. 2019, Art. no. 109581E.