**Background:** Negative-tone development (NTD) photoresists are prone to shrinkage effects during lithographic processing. Along with deformation seen during the postexposure bake (PEB), there are additional effects during the development that cannot be fully explained by a conventional PEB shrinkage model alone.

**Aim:** Understand the impact of PEB shrinkage on the development rate. Develop a model that can help predict resist profiles after chemical development.

**Approach:** A PEB shrinkage model for NTD resists is introduced, which uses the thermal properties of the resist material to help simulate shrinkage. The deformed state of the resist is used as an input to the development rate equation to predict the final feature dimensions observed in experiments.

**Results:** The strain concentration within the resist bulk can have an influence on the stability of the resist during the development. The strain influences the development rate depending on the resist feature shape and contours.

**Conclusions:** The results from this study can help improve optical proximity correction (OPC) modeling performance and help better understand the deformation characteristics of NTD resist materials. The model also shows that the development shrinkage has an influence on the edge placement error.

## 1.

## Introduction

Traditional negative-tone development (NTD) resists are prone to shrinkage and deformation effects during various stages of the lithography process. Volume losses are clearly observed during the exposure and PEB for chemically amplified resists. In the past, various methods were employed to model photoresist shrinkage effects during PEB.^{1}^{–}^{3} These models used a correlation of the photochemical properties after the reaction diffusion step to the creation of voids in the resist. The resist CDs and dimensions were then estimated by a deformation model. The main benefit of these methods was a lower overall computation time compared to more rigorous methods. These approaches do not consider the thermomechanical interaction of the resist polymer with the heating step. The thermal properties of the resist provide information on its conductivity and heat capacity. From another perspective, the modeling technique introduced in this paper uses the thermal properties of the resist to help model the shrinkage effects during PEB. The physical material properties are then coupled together to calculate the overall stresses and strains in the resist bulk with the help of finite element modeling (FEM). The thermal strain is, therefore, the driving force behind deformation. The concentration of the deprotected sites is used to represent the mechanical properties of the resist bulk. This paper also compares the impact of various material properties on the deformation characteristics.

Modeling these undesirable shrinkage effects merely during the PEB stage does not give an accurate prediction of the final developed photoresist profile. In some cases, additional unexpected deformation effects of the resist are observed that cannot be predicted by the PEB shrinkage model alone. Some experiments show evidence of increased dissolution of the resist in certain areas during the development. It is, therefore, important and necessary to model the deformation after PEB considering supplementary mechanical aspects of the resist material.

Photoresist polymers have viscoelastic properties that result in plastic deformation after the initial elastic recovery. These materials are prone to irreversible viscous flow that causes deformation even after the application of load. Heating during the postexposure bake (PEB) is the main temperature load applied to the resist. The viscous flow of the polymer leads to broken bonds within the polymerized photoresist. The broken bonds together with the chemical developer penetration cause the resist to soften and develop away faster than expected. This leads to more resist material being washed away than expected and causes discrepancies between the measured and simulated dimensions of the final developed resist profile. This phenomenon can be considered as resist shrinkage during the development.

Some experimental observations (later discussed in Sec. 3.1) indicate that the shrinkage simulated after PEB does not capture all the deformation effects taking place during lithographic processing. In this paper, a new strain-based model is introduced, which helps in understanding the additional deformation effects in chemically amplified resists noticed during the development process. The development shrinkage model predicts the actual final resist profile much more closely compared to the standard deformation model. The developer liquid and resist polymer influence chemical and mechanical effects to produce unwanted artifacts in the resulting resist pattern. Figure 1 shows the various processing stages, in which shrinkage effects are prevalent. A slight material loss is possible during the exposure stage depending on the dose value. Dose and focus variations can lead to variations in resist sidewall angles. During PEB, there is an out-gassing of the volatile by-product of the acid-catalyzed deprotection reaction. The evaporation of this by-product leaves voids in the resist material and leads to losses in height, CDs, and volume. A more pronounced change in sidewall angles can also be observed. Our proposed simulation model can help throw some light on the further dimensional changes possible during development, which is a direct consequence to the PEB shrinkage. This is depicted in the last figure of the schematic and shows how the chemical development process could impact the lateral dimensions ($XY$ plane) of a given resist pattern.

Section 2 gives an insight into the PEB shrinkage model used together with the various material aspects of the photoresist important for simulating deformation. A finite element model is used as a basis for shrinkage modeling and is coupled with input coming from lithography process simulations to help determine resist dimensions after PEB. The model presented is also compared to results obtained from analytical estimates and an analysis of the various model parameters impacting deformation is performed. Several observations from experimental data have indicated at effects not entirely explained by the PEB shrinkage model alone. Section 3 throws some light on how the mechanical status of the resist after PEB could lead to changes in the development rate leading to a disparity in profile contours. A parametric study of the geometrical aspects of the profile exhibits important findings useful for the field of optical proximity correction (OPC). Finally, a summary of the crucial observations of this study is presented while providing an outlook for the future work.

## 2.

## Shrinkage during PEB

## 2.1.

### Mechanical and Thermal Properties of the Photoresist

Several thermal and mechanical material properties need to be taken into consideration in order to model deformation and shrinkage effects using the finite-element method. Since the photoresist is basically a polymer, it has physical material properties such as the Young’s modulus $E$ (also known as elastic modulus), Poisson’s ratio $\nu $, density $\rho $, bulk modulus $K$, and shear modulus $G$. The Young’s modulus and Poisson’s ratio values are used to specify the physical material properties for linear isotropic elasticity. Typical values for photoresist materials in deep ultraviolet lithography are calculated from experiments and shown in Table 1. The bulk and shear modulii values are estimated using the Lamé parameters together with the Young’s modulus and Poisson’s ratio.^{4}

## Table 1

Typical material properties of NTD photoresists.

Property | Value |
---|---|

Young’s modulus $E$ | 3 GPa |

Poisson’s ratio $\nu $ | 0.4 |

Density | $1200\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{kg}/{\mathrm{m}}^{3}$ |

Thermal conductivity $\kappa $ | $0.2\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{W}/\mathrm{mK}$ |

Specific heat $c$ | $1.5\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{kJ}/\mathrm{kg}\text{\hspace{0.17em}}\mathrm{K}$ |

In addition to mechanical properties, the thermal properties are used to describe changes in the photoresist during the PEB process. Shrinkage effects as a result of heating/baking are seen due to the outgassing effect of the photoresist polymer. Consequently in this case, the PEB baking temperature and the thermal expansion co-efficient are decisive factors in influencing shrinkage. The heat conduction of a material is defined by the thermal conductivity $\kappa $. Materials with a higher $\kappa $ value are better conductors of heat and vice-versa. It is defined as the quantity of heat $Q$, transmitted in time $t$ through a thickness $L$, in a direction normal to a surface of area $A$, due to a temperature difference $\mathrm{\Delta}T$, under steady-state conditions and when the heat transfer is dependent only on the temperature gradient. Thermal conductivity equals to the heat flow rate $(\frac{Q}{t})$ times the distance $(L)$ divided by the cross sectional area $(A)$ and the temperature difference $(\mathrm{\Delta}T)$:

The specific heat capacity ($c$, also known as specific heat or SHC) of a substance is defined as heat capacity per unit mass. It is the amount of energy required to raise the temperature of 1 kg of a substance by 1 K. The unit for the specific heat capacity $c$ is $\frac{\mathrm{J}}{\mathrm{kg}\text{\hspace{0.17em}}\mathrm{K}}$. The strain generated due to temperature changes in a material is known as thermal strain ${\u03f5}^{\mathrm{th}}$. The secant or average thermal expansion coefficient $({\alpha}^{\mathrm{se}})$ is used to describe the general extension or contraction in length due to temperature $T$ changes based on a reference temperature ${T}_{\mathrm{ref}}$ and can be given as

To model shrinkage in our case, a negative value of ${\alpha}^{\mathrm{se}}$ is used. This parameter, therefore, has the ability to describe the amount of shrinkage in a material. Depending on the type of simulation (2D or 3D), ${\alpha}_{x}^{\mathrm{se}}$, ${\alpha}_{y}^{\mathrm{se}}$, and ${\alpha}_{z}^{\mathrm{se}}$ are specified accordingly. $T$ is the PEB temperature and ${T}_{\mathrm{ref}}$ is the room temperature. For isotropic materials, the volumetric thermal expansion coefficient is three times the linear coefficient. The ${\alpha}^{\mathrm{se}}(T)$ value is, therefore, multiplied by a factor of three to take into consideration shrinkage in all dimensions:

The material property values mentioned in Table 1 represent those for typical NTD resists and based on them corresponding values were chosen for the FEM simulations.

## 2.2.

### FEM Deformation Model

Figure 2 explains how the secant thermal expansion coefficient ${\alpha}^{\mathrm{se}}$ and the corresponding instantaneous value ${\alpha}^{\text{in}}$ are calculated. ${\alpha}^{\mathrm{se}}$ is used for the model since during the PEB process the baking temperature is gradually increased with respect to the room temperature till it reaches the desired value. This is shown in Fig. 2, where ${\alpha}^{\mathrm{se}}$ is the slope between the reference temperature (room temperature) and the PEB temperature. Based on Hooke’s law, stress is directly proportional to strain. The relation between the two can be represented by Eq. (4) with the help of a stiffness matrix formulation:

where $\mathbf{\sigma}$ is the stress vector, ${[\begin{array}{cccccc}{\sigma}_{x}& {\sigma}_{y}& {\sigma}_{z}& {\sigma}_{xy}& {\sigma}_{yz}& {\sigma}_{xz}\end{array}]}^{\mathrm{T}}$, $[\mathit{D}]$ is the stiffness matrix or the Young’s modulus, ${\mathbf{\epsilon}}^{\mathrm{el}}=\mathbf{\epsilon}-{\mathbf{\epsilon}}^{\mathrm{th}}$ is the elastic strain vector, $\mathbf{\epsilon}$ is the total strain vector = ${[\begin{array}{cccccc}{\epsilon}_{x}& {\epsilon}_{y}& {\epsilon}_{z}& {\epsilon}_{xy}& {\epsilon}_{yz}& {\epsilon}_{xz}\end{array}]}^{\mathrm{T}}$, and ${\mathbf{\epsilon}}^{\mathrm{th}}$ is the thermal strain vector.From Eq. (4), the strain $\mathbf{\epsilon}$ can be formulated as

## Eq. (5)

$$\mathbf{\epsilon}={\mathbf{\epsilon}}^{\mathrm{th}}+{[\mathit{D}]}^{-\mathbf{1}}\mathit{\sigma}.$$The total strain $\mathbf{\epsilon}$ is the sum of the elastic ${\mathbf{\epsilon}}^{\mathrm{el}}$ and thermal strain ${\mathbf{\epsilon}}^{\mathrm{th}}$ values:

## Eq. (6)

$${\mathbf{\epsilon}}^{\mathrm{th}}=\mathrm{\Delta}T{[\begin{array}{cccccc}{\alpha}_{x}^{\mathrm{se}}& {\alpha}_{y}^{\mathrm{se}}& {\alpha}_{z}^{\mathrm{se}}& 0& 0& 0\end{array}]}^{\mathrm{T}}.$$The flexibility or compliance matrix ${[\mathit{D}]}^{-1}$ is

## Eq. (7)

$${[\mathit{D}]}^{-1}=\left[\begin{array}{cccccc}1/{E}_{x}& -{\nu}_{xy}/{E}_{x}& -{\nu}_{xz}/{E}_{x}& 0& 0& 0\\ -{\nu}_{yx}/{E}_{y}& 1/{E}_{y}& -{\nu}_{yz}/{E}_{y}& 0& 0& 0\\ -{\nu}_{zx}/{E}_{z}& -{\nu}_{zy}/{E}_{z}& 1/{E}_{z}& 0& 0& 0\\ 0& 0& 0& 1/{G}_{xy}& 0& 0\\ 0& 0& 0& 0& 1/{G}_{yz}& 0\\ 0& 0& 0& 0& 0& 1/{G}_{xz}\end{array}\right],$$The thermal strain as formulated above is analogous to shrinkage strain since the deformation occurs due to changes in temperature/heat flux and is used directly as an input for the shrinkage model in the ANSYS simulation environment.^{5} Based on the concentration of protected/deprotected sites $P$ in the photoresist, the thermal strain ${\mathit{\epsilon}}^{\mathrm{th}}$ or shrinkage strain ${\mathit{\epsilon}}^{\text{shrink}}$ can be formulated as

## Eq. (8)

$${\mathit{\epsilon}}^{\text{shrink}}={\mathit{\epsilon}}^{\mathrm{th}}=3{\alpha}^{\mathrm{se}}\mathrm{\Delta}T=\eta (1-P).$$Here $\eta $ is the shrink parameter, which describes how strongly the resist material responds to shrinkage. From Eq. (8), the coefficient of thermal expansion can be obtained as

The coefficient of thermal expansion $({\alpha}^{\mathrm{se}})$ value is used as an input to model the shrinkage phenomenon using the FEM. ${\alpha}^{\mathrm{se}}$ takes negative values to help simulate the loss of resist volume during the PEB process.

Along with the other material parameters, the Young’s modulus value also varies in different regions of the photoresists. From experiments, it was observed that the exposed parts of thin-film resists are softer (i.e., having a lower Young’s modulus $E$) compared to the unexposed parts.^{6} This observation has been taken into consideration in the implementation of the FEM model. The stiffness ratio is, therefore, given as

## Eq. (10)

$${E}^{\prime}=\frac{\text{exposed}\text{\hspace{0.17em}}E}{\text{unexposed}\text{\hspace{0.17em}}E}.$$The stiffness ratio ${E}^{\prime}$ is a variable parameter in the FEM simulations and a value of 0.1 is used for simulating the test cases to best describe the experimental results. Atomic force microscopy (AFM) measurements, however, show that the stiffness ratio values can range from 0.1 to 1.

## 2.3.

### Analytical Approximation

The accuracy of the finite element model described in Sec. 2.2 can be assessed with the help of an analytical study of resist height shrinkage. This can be done by comparing the height loss of a fully deprotected resist between the analytical and FEM models. Periodic and fixed boundary conditions are set accordingly in the FEM model and it is assumed that the deformation during PEB is relatively small. The approximate height loss of the resist can be calculated by considering a cubic resist volume shrinking isometrically with fixed lateral boundary conditions. The effective height loss is, therefore, calculated with respect to no change in dimensions along the $X$ and $Y$ directions (substrate plane). This can be visualized as a three step process as shown in Fig. 3 with ${L}_{0}$ being the initial resist height and $H$ being the final height. The resulting change in height for a fully deprotected resist $(P=0)$ is dependent on the Poisson’s ratio $\nu $. This ratio is valid for small strains along a given co-ordinate axis and is generally in the range between 0 and 0.5 for the most materials:

## Eq. (11)

$$\frac{H}{{L}_{0}}=\frac{{L}_{1}}{{L}_{0}}-\nu (1-\frac{{L}_{1}}{{L}_{0}})-(1+\nu )\left\{\frac{1-\nu}{2}\right(\frac{{L}_{1}}{{L}_{0}})-\sqrt{\frac{{(1-\nu )}^{2}}{4}{\left(\frac{{L}_{1}}{{L}_{0}}\right)}^{2}+\nu (\frac{{L}_{1}}{{L}_{0}}-1)\frac{{L}_{1}}{{L}_{0}}}\}.$$Equation (11) was derived and can be found in Appendix A. Figure 4 shows the influence of the shrink factor $\eta $ on the final resist height. The approximate nature of the derived Eq. (11) is the reason behind the small deviations in the final resist heights between the FEM model and derived formula. These deviations increase with larger Poisson’s ratios and volumetric shrink factors. In contrast to the analytically approximated formula, the FEM model can be applied for all valid $\nu $ and $\eta $ values and is used to obtain simulation results shown in this paper.

## 2.4.

### Simulation Results for Real Profiles

A number of 2D test cases were analyzed in order to verify and better understand the shrinkage effects resulting from a variation in the input parameters. Information regarding the 2D test cases can be summarized as follows.

• Dense feature: 50-nm line with 100-nm pitch.

• Semidense feature: 100–nm line with 350-nm pitch.

• Isolated feature: 250-nm line with 1200-nm pitch.

• PEB temperature is 130°C and $\eta $ is $-0.5$.

• CD cutline band is between 60% and 80% of resist thickness. This band is used for extracting CD values along the resist thickness. CD values are calculated at specific locations along the resist thickness and then averaged to get a quantitative value. For the deformed profile, the cutline band is in the same range but with respect to the new deformed thickness.

Deformed protection group concentration profiles are shown in Fig. 5. They show the overall deformation of the exposed regions. Larger exposed parts on the mask in the case of bright field masks leads to a much greater volume loss of the resist compared to dark field masks. The profiles were simulated using the model explained in Sec. 2.2.

A number of material parameters such as the Young’s modulus $E$, Poisson’s ratio, and shrink factor $\eta $ were varied to investigate their impact on the overall deformation of the photoresist. For the 2D test cases, the shrink factor $\eta $ has the largest impact on the final resist dimensions and CD values. The results of these parameter variations can be seen from Fig. 6 for dense, semidense, and isolated lines and spaces features.

A small change in CD is observed due to a variation of the Poisson’s ratio. Since the shrink factor is directly proportional to the amount of shrinkage as described in Eq. (8), it has maximum effect on both the CD and height. The $\mathrm{\Delta}\mathrm{CD}$ values are larger for the isolated case than the dense one, whereas the change in height $\mathrm{\Delta}H$ is more pronounced for the dense feature. This is observed while varying all the parameters. Based on the results obtained for the 2D case as seen in Fig. 6, it is clear that a change in the Young’s modulus value has almost no impact on the CD and height for a given feature density. This is an important observation, especially when we start looking at resist materials used in extreme ultraviolet lithography (EUVL). The Young’s modulus values $E$ in EUVL can be $<10\%$ of DUV resist materials.^{7} This would mean that there would be comparable volume losses as seen here even in EUVL.

Another crucial finding is that the change in CD is directly proportional to the feature linewidth, i.e., larger $\mathrm{\Delta}\mathrm{CD}$ values for a larger linewidth. On the other hand, the change in height is inversely proportional to the linewidth as seen in Figs. 6(d)–(f). This means that for the larger linewidth region with a higher stiffness a lower amount of longitudinal deformation compared to lateral deformation is observed. Looking at Figs. 6(b) and (e), we observed that as the Poisson’s ratio of the resist material increases, there is a linear increase in the height loss compared to a small decrease in $\mathrm{\Delta}\mathrm{CD}$ values. This is mainly due to the boundary conditions, whereas the resist is free to shrink at the top rather than in the lateral $XY$ plane. The results shown above are in general agreement with the experimental data shown in Refs. 1 and 2.

## 3.

## Impact of Shrinkage on the Development Rate

## 3.1.

### Influence of Strain

The mechanical properties of the resist depend on the molecular weight, chemical composition, and how the optical parameters influence the polymer material. The formation of free volume during PEB has an impact on the physical and chemical stability of the photoresist during the development. So far, the model described in Sec. 2 was capable of predicting resist profiles after PEB shrinkage well^{8}^{,}^{9} in combination with the kinetic rate equation (Mack model)^{10} for the development. We, however, soon came across a case where this approach could not provide a satisfactory estimate of the final developed resist profile. This is evident from the simulated 3D resist profile shown in Fig. 7(b). The wriggling contours along the spatial sidewalls shown in Fig. 8 are not captured by the deformation model. The simulated profile actually shows an outward ridge or indentation in the spatial direction in the regions close to the contact holes. This noncompliance to the experimental data prompted us into investigating the development model more closely. The deformation phenomenon being a mechanical phenomenon requires additional considerations to model its effect on the chemical development step.

The photoresist is elastic in nature during PEB since the baking temperature is below the glass transition temperature of the polymer ${T}_{\mathrm{PEB}}<{T}_{G}$.^{3}^{,}^{11} Also the glass transition temperature is directly proportional to the Young’s modulus value.^{11} This change in volume induces localized stresses and strains within the deformed resist material. The behavior of the chemical developer with these localized regions can have an influence on the contour and shape of the final developed resist. The deformation behavior before bond breakage will depend on the bond strength, temperature, and entanglement density, i.e., regions that are weakly bonded will undergo larger strains. The Young’s modulus or stiffness of the resist depends on the protection group concentration based on our PEB deformation model.

The chemical composition of the resist with regards to its bond strength, density, and molecular weight plays a pivotal role in determining the physical and mechanical properties of the resist.^{12} Negative-tone crosslinking resists like SU8 have lower ultimate strains as the crosslink density increases.^{13} Also regions in the resist bulk that exhibit higher strains tend to have weaker molecular bonding. Dissolution rates at lower exposure areas have an impact on the defect characteristics of the resist profile.^{14} The thermodynamic characteristics of the photoresist material influence the dissolution rate due to the fact that the chemical developer depends on how fast the solvent molecule penetrates the polymer network.^{15} The solubility characteristics of a particular polymer depends on its chemical structure and it tends to dissolve in solvents with similar solubility parameters and polarities. The residual solvent in the photoresist after PEB has an impact on the dissolution rate of the polymer.^{16} A rapid baking of PMMA polymer gives rise to “extra free volume” and strain in the polymer leading to accelerated dissolution rates.^{17} Moreover, the dissolution rate of the photoresist resin depends on the extent of hydrogen bonding of the hydroxyl sites. Hydrogen bonding is known to inhibit dissolution in diazonaphthoquinone (DNQ)-novolac type resist polymers.^{18} Due to tensile deformation, hydrogen bond breaking occurs at higher strain levels and temperatures^{19} and in turn aids in polymer dissolution. Although some of these arguments were derived for PMMA and DNQ-type photoresists, they indicate a certain sensitivity of dissolution rates of photopolymers to mechanical strain. We consider the Cartesian components $i,j=\mathrm{1,2},3$ of the strain tensor and it is represented as

## Eq. (12)

$${\u03f5}_{ij}=\frac{1}{2}(\frac{\partial {u}_{i}}{\partial {x}_{j}}+\frac{\partial {u}_{j}}{\partial {x}_{i}}).$$For a three-dimensional case, the strain tensor can be represented as

## Eq. (13)

$$\left[\begin{array}{ccc}{\u03f5}_{xx}& {\u03f5}_{xy}& {\u03f5}_{xz}\\ {\u03f5}_{yx}& {\u03f5}_{yy}& {\u03f5}_{yz}\\ {\u03f5}_{zx}& {\u03f5}_{zy}& {\u03f5}_{zz}\end{array}\right]=\left[\begin{array}{ccc}\frac{\partial {u}_{x}}{\partial x}& \frac{1}{2}(\frac{\partial {u}_{x}}{\partial y}+\frac{\partial {u}_{y}}{\partial x})& \frac{1}{2}(\frac{\partial {u}_{x}}{\partial z}+\frac{\partial {u}_{z}}{\partial x})\\ \frac{1}{2}(\frac{\partial {u}_{y}}{\partial x}+\frac{\partial {u}_{x}}{\partial y})& \frac{\partial {u}_{y}}{\partial y}& \frac{1}{2}(\frac{\partial {u}_{y}}{\partial z}+\frac{\partial {u}_{z}}{\partial y})\\ \frac{1}{2}(\frac{\partial {u}_{z}}{\partial x}+\frac{\partial {u}_{x}}{\partial z})& \frac{1}{2}(\frac{\partial {u}_{z}}{\partial y}+\frac{\partial {u}_{y}}{\partial z})& \frac{\partial {u}_{z}}{\partial z}\end{array}\right].$$Here $u$ represents the displacement along the principal co-ordinate $X$, $Y$, or $Z$ axes. ${\u03f5}_{xx}$, ${\u03f5}_{yy}$, and ${\u03f5}_{zz}$ are the normal strains shown on the trace of the strain tensor and the off-diagonal components ${\u03f5}_{ij}$ represent the shear strains along the $ij$ planes. The strain tensor is symmetric: ${\u03f5}_{ij}={\u03f5}_{ji}$. Consider a material cube having sides $a$, $b$, $c$ with $a=b=c$. The relative change in volume at any given point in the resist material can be given by

## Eq. (14)

$$\frac{\mathrm{\Delta}V}{V}=\frac{(a+\mathrm{\Delta}a)(b+\mathrm{\Delta}b)(c+\mathrm{\Delta}c)-abc}{abc}=(1+{\u03f5}_{xx})(1+{\u03f5}_{yy})(1+{\u03f5}_{zz})-1\phantom{\rule{0ex}{0ex}}\approx {\u03f5}_{xx}+{\u03f5}_{yy}+{\u03f5}_{zz}.$$Normal strains induce a change in volume while shear strains induce a change in shape. During PEB, there is a significant change in volume. The volume loss causes strains within the resist bulk, which can be represented by the expression in Eq. (14), which approximates to the trace of the strain tensor. Figure 9 shows the individual strain components of the strain tensor for the test case considered in this section. Strains generated are calculated after implementing the PEB shrinkage model introduced in Sec. 2.2.

## 3.2.

### Impact of Strain on the Development Rate

The impact of the strains induced within the photoresist bulk on the photoresist development rate will be explained in this section. So far, we used the Mack model^{10} to simulate a developed resist profile. The only modification made was to the protection group concentration input, which was inverted to account for the effects of a negative tone developer.^{20} This approach, however, cannot provide results seen from experiments. On investigating the causes for these discrepancies further, it was observed that the volumetric strain had a much higher value in the regions close to the contact holes as shown in Fig. 10(a).

Standard negative-tone development (NTD) can be modeled by substituting the normalized protection group concentration $M$ as $(1-M)$ and inserting it into the kinetic rate equation. ^{10}^{,}^{20} ${r}_{\mathrm{max}}$ is the maximum dissolution rate for a fully exposed resist, whereas ${r}_{\mathrm{min}}$ is the minimum dissolution rate for an unexposed resist ($M=1$). $n$ is called the dissolution selectivity parameter and it varies depending on the photoresist contrast. ${r}_{\mathrm{min}}$ and ${r}_{\mathrm{max}}$ differ for positive-tone development (PTD) and NTD.^{21} This is shown in Fig. 11 with higher ${r}_{\mathrm{min}}$, lower ${r}_{\mathrm{max}}$, and $n$ values for negative-tone developers compared to positive-tone developers. Using this kinetic rate equation directly resulted in a resist profile as seen in Fig. 7. To describe the impact of strain on the development rate, we propose the following extension to the Mack development rate model:

## Eq. (15)

$$r={r}_{\mathrm{max}}\frac{(a+1){(1-M)}^{n}}{a+{(1-M)}^{n}}+{r}_{\mathrm{min}}+{(1-{M}_{\mathrm{th}})}^{n}\xb7\u03f5\xb7{r}_{\text{strain}}.$$The parameter values chosen in Table 2 are in accordance to the experimental observations of Tarutani et al.,^{21} where he observed that the ${r}_{\mathrm{max}}$ values for NTD resists were smaller compared to PTD resists, whereas the ${r}_{\mathrm{min}}$ values were larger. The additional strain component ${(1-{M}_{\mathrm{th}})}^{n}\xb7\u03f5\xb7{r}_{\text{strain}}$ was chosen to help model the effect of strain on the overall development rate. ${r}_{\text{strain}}$ represents the rate at which the strained regions are developed away and can be used as an adjustable model parameter.

## Table 2

Parameter values for the modified development rate equation.

Parameter | Value |
---|---|

$n$ | 10 |

${r}_{\mathrm{min}}$ | $0.01\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{nm}/\mathrm{s}$ |

${r}_{\mathrm{max}}$ | $50\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{nm}/\mathrm{s}$ |

${M}_{\mathrm{th}}$ | 0.5 |

${r}_{\text{strain}}$ | $2\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{nm}/\mathrm{s}$ |

Figure 12 shows the final resist profile after implementing the modified development rate equation. It demonstrates the improvement in the contour prediction compared to the standard rate equation. Going back to the previous contour estimation from Fig. 7(b), we see the clear improvement over the old development rate model. The wriggling and indentations of the spatial regions are captured well as can be seen in the SEM images of Fig. 13.

The steepness of the inflection curve is governed by the parameter $n$ and has a lower value for modeling negative-tone developer action.^{20} Due to the lower $n$ value, the other parameters in the modified rate equation namely ${r}_{\text{strain}}$ and $\u03f5$ have a greater influence on the dissolution rate of the photoresist as shown in Fig. 14. The dissolution rate basically increases in regions that are not developed away and that have a lower protection group concentration. Also an increase in $n$ reduces the impact on the dissolution rates of these two parameters. Since we are inverting the protection group concentration for NTD, these slightly higher development rates influence the final resist shape. This approach has helped us in modeling the impact of strain on the overall chemical development process.

## 3.3.

### Parametric Study

Developed profiles simulated using the modified rate equation can show significant differences with respect to the ones simulated without considering the impact of shrinkage effects on the development rate. This is particularly the case with numerous 3D resist profiles similar to the one shown in Fig. 7. A combination of optical, chemical, and mechanical effects influences the final shape and contours of the developed resist profiles. Optical proximity effects are known to be dependent on the development parameters,^{22} and in Sec. 3, we have shown that mechanical effects can impact the shape and contours of the resist profile significantly as well.

A 3D resist profile is generated using the mask layout shown in Fig. 15 with varying space feature lengths in the horizontal direction. Figure 16 exhibits the simulated impact of the size of the space on simulated concentrations of the deprotection group, strain, and photoresist profiles without/with the additional term in the development rate [Eq. (15)]. As the space feature length increases, there is an increase in the width of the horizontal space as shown in Figs. 16(e) and 16(f). The volumetric strain plots from Figs. 16(c) and 16(d) exhibit an increased strain response in certain areas. The development rate [Eq. (15)] when used with the strain concentration input results in an increase in vertical space CD as seen in Figs. 16(g) and 16(h). This change in CD is mainly attributed to the strain distribution within the resist bulk and also optical proximity effects. Denser features are shown be particularly vulnerable to these additional effects during the development. Looking again closely to the strain distribution plots in Figs. 16(c) and 16(d) shows us that as the size of the horizontal space feature decreases and the distance between the vertical trench increases, the strain concentration between the two decreases. When the two features are close to each other like in Fig. 16(c), the strains around them tend to build-up/overlap and in turn affect the development rate in those regions. The increased rate of development causes an inward indentation in the sidewalls along the vertical spatial feature shown in Fig. 16(g).

The disparity between the simulated CD values with respect to the values without considering shrinkage effects is shown in Fig. 17. The $\mathrm{\Delta}\mathrm{CD}$ value represents the difference between the CD values of the model considering various shrinkage effects and a simple model (baseline) not including any of those effects. There is a significant difference in CDs between the two development rate models. The blue line represents CD values considering the standard development rate model and the green one represents the CDs considering the additional effects of strain on the development. The $\mathrm{\Delta}\mathrm{CD}$ values calculated for the pure PEB shrinkage model are similar to the values obtained by the baseline model but fluctuate quite a bit. This is mainly due to the OPC of the mask shown in Fig. 15. A combination of mechanical and OPC effects leads to a significant difference in the CDs obtained after additional consideration of the development shrinkage model. It is, therefore, crucial to also consider the impact of strain on the overall development rate to help predict resist profiles more accurately.

## 4.

## Conclusions

Shrinkage and deformation of NTD photoresists during the PEB step have significant effects on the feature CD, height, and volume. Changes in volume during this step can be as high as 25% with CD and height losses being in the range of several nanometers. The finite-element method can be seen as a robust modeling tool to predict the deformation behavior of photoresists. The main material parameter influencing shrinkage is the Poisson’s ratio. This together with other simulation parameters such as the PEB temperature and shrink factor $\eta $ has an influence on the resist volume, sidewall angles, CD, and height. For thin film resists, the stiffness ratio of the exposed regions is lower than the unexposed regions. A quantitative assessment of the PEB shrinkage model for the 2D test cases shows that changes in CDs and height depends on the feature pitch. Isolated features tend to show larger CD variations compared to denser features. The shrink factor parameter adjustment allows for flexibility in modeling the shrinkage and deformation effectively.

It was demonstrated that merely modeling shrinkage effects during PEB is not enough and there could be other additional causes to pattern deformation. An approach involving the investigation of the significance of the volumetric strain on the development step was formulated. Shrinkage and deformation seen during PEB can have an impact on the final development rate. A modified kinetic rate equation is used to model additional shrinkage effects during the development. The amount of volume losses during the development depends mainly on the pattern size, density, and complexity. The strain concentrations vary based on these aspects and impact the overall development rate accordingly. A combination of mechanical and optical proximity effects is shown to contribute to the amount of volume losses in the resist polymer. This can lead to localized regions with lower CDs and undercut profiles causing pattern instability/collapse during further processing steps.

Future studies could involve the investigation of the mechanical stability of various photoresist materials to negative-tone chemical developers. A subsequent analysis of these resist losses on pattern collapse behavior could help in understanding the true impact on lithographic performance.

## 5.

## Appendix A: Derivation of Height Loss Approximation

The derivation of the height loss approximation formula for a fully deprotected resist starts with the basic understanding of how the lateral strain $\frac{\mathrm{\Delta}d}{d}$ and longitudinal strain $\frac{\mathrm{\Delta}L}{L}$ are related to each other:

$\nu $ is the ratio of lateral strain to longitudinal strain and can be represented by the equation shown above. The relation of the volumetric shrink factor with respect to the resist height loss can be deduced, which is also used in Eq. (8):

Since $\eta $ is negative in our case:

To derive the effective height of a fully deprotected resist, we first assume an isotropic shrinkage of the cube resulting in side length ${L}_{1}$ along all three co-ordinate axes. A stretch along the $X$ direction $(\mathrm{\Delta}{X}_{2})$ yields a contraction along $Y$ $(\mathrm{\Delta}{Y}_{2})$ and $Z$ $(\mathrm{\Delta}H)$ that depends on $\nu $ given by

## Eq. (18)

$${L}_{0}=L1+\mathrm{\Delta}{X}_{2}-\mathrm{\Delta}{Y}_{3}\phantom{\rule{0ex}{0ex}}{L}_{0}=L1+\mathrm{\Delta}{X}_{2}-\nu \mathrm{\Delta}{Y}_{3}\left(\frac{{L}_{1}+\mathrm{\Delta}{X}_{2}}{{L}_{1}-\nu \mathrm{\Delta}{X}_{2}}\right).$$Similarly, a stretch along $Y$ yields a contraction along the height and the $X$ direction given by

## Eq. (19)

$${L}_{0}={L}_{1}-\nu \mathrm{\Delta}{X}_{2}+\mathrm{\Delta}{Y}_{3}\phantom{\rule{0ex}{0ex}}{Y}_{3}={L}_{0}-{L}_{1}+\nu \mathrm{\Delta}{X}_{2}.$$Substituting the value of $\mathrm{\Delta}{Y}_{3}$ in Eq. (18) results in the following equation:

## Eq. (20)

$${L}_{0}={L}_{1}+\mathrm{\Delta}{X}_{2}-\nu ({L}_{0}-{L}_{1}+\nu \mathrm{\Delta}{X}_{2})\left(\frac{{L}_{1}+\mathrm{\Delta}{X}_{2}}{{L}_{1}-\nu \mathrm{\Delta}{X}_{2}}\right)\phantom{\rule{0ex}{0ex}}\mathrm{\Delta}{X}_{2}^{2}-{L}_{1}\mathrm{\Delta}{X}_{2}\left(\frac{1-\nu}{\nu}\right)-{L}_{1}({L}_{1}-{L}_{0})\frac{(1-\nu )}{(1+\nu )\nu}=0.$$The equation for the final resist is given by

## Eq. (21)

$$H={L}_{1}-\nu \mathrm{\Delta}{X}_{2}-\nu \mathrm{\Delta}{Y}_{3}={L}_{1}-\nu \mathrm{\Delta}{X}_{2}-\nu ({L}_{0}-{L}_{1}+\nu \mathrm{\Delta}{X}_{2})\phantom{\rule{0ex}{0ex}}={L}_{1}-\nu ({L}_{0}-{L}_{1})-\nu \mathrm{\Delta}{X}_{2}(1+\nu ).$$On solving Eq. (20) for $\mathrm{\Delta}{X}_{2}$ and inserting it in the above equation, it gives us:

## Eq. (22)

$$\mathrm{\Delta}{X}_{2}=\frac{1-\nu}{2\nu}{L}_{1}-\sqrt{\frac{{(1-\nu )}^{2}}{4{\nu}^{2}}{L}_{1}^{2}+({L}_{1}-{L}_{0})\frac{{L}_{1}}{\nu}}.$$Substituting the value of $\mathrm{\Delta}{X}_{2}$ into Eq. (21) yields the equation for the final resist height with respect to the initial height:

## Eq. (23)

$$\frac{H}{{L}_{0}}=\frac{{L}_{1}}{{L}_{0}}-\nu (1-\frac{{L}_{1}}{{L}_{0}})-(1+\nu )\left\{\frac{1-\nu}{2}\right(\frac{{L}_{1}}{{L}_{0}})-\sqrt{\frac{{(1-\nu )}^{2}}{4}{\left(\frac{{L}_{1}}{{L}_{0}}\right)}^{2}+\nu (\frac{{L}_{1}}{{L}_{0}}-1)\frac{{L}_{1}}{{L}_{0}}}\}.$$Expressing the above equation in terms of the volumetric shrink factor $\nu $ from Eq. (17) is

## Eq. (24)

$$\frac{H}{{L}_{0}}=(1-\frac{\eta}{3})-\nu \left(\frac{\eta}{3}\right)-(1+\nu )\left\{\frac{1-\nu}{2}\right(1-\frac{\eta}{3})-\sqrt{\frac{{(1-\nu )}^{2}}{4}{(1-\frac{\eta}{3})}^{2}+\nu [{(1-\frac{\eta}{3})}^{2}-(1-\frac{\eta}{3}\left)\right]}\}.$$Results have been plotted for relatively low values of $\eta $, because Eq. (24) holds true for small deformations only. Larger $\eta $ values (for modeling larger deformation during PEB) can cause the value inside the square root to turn negative.

## Acknowledgments

This project work was partly funded by Synopsys GmbH (Germany). We would like to thank and gratefully acknowledge the experimental input from Uzodinma Okoroanyanwu (Enx Labs, USA) and the valuable SEM image data from our colleagues at IMEC (Belgium).

## References

## Biography

**Sean D’Silva** received his BEng degree in mechanical engineering from the University of Mumbai and his MSc degree in computational engineering from the University of Erlangen-Nuremberg. He is a doctoral candidate at Fraunhofer IISB/University of Erlangen-Nuremberg. His field of research includes the study of photoresist material properties, modeling of photoresist deformation, and the implementation of machine learning techniques for predicting stochastic effects.

**Thomas Mülders** received his PhD in physics in 1999 from the Technical University Aachen, Germany. From 1999 to 2001, he worked at the CNRS Orléans, France, focusing on the development of simulation methods for many particle systems. Since 2001, he has been working in the field of lithography modeling, first in the Optics Simulation Group at Carl Zeiss, later from 2003 to 2008 at Infineon Technologies. In 2008, he became a member of the Silicon Engineering Group of Synopsys.

**Hans-Jürgen Stock** received his PhD in physics in 1998 from the Universität Bielefeld, Germany. From 1997 to 2004, he worked at Varetis AG in the area of speech automatization for DA call center. Since 2004, he has been working in the field of lithography simulations. Initially, he was in the Solid C/E Team of Sigma-C, afterward as a R&D manager in the Silicon Engineering Group of Synopsys.

**Andreas Erdmann** is the head of the Fraunhofer IISB Computational Lithography and Optics Group and teaches as “Privatdozent” at the University of Erlangen. He has more than 20 years of experience in optical and EUV lithography. He chaired SPIE conferences on “optical microlithography” and “optical design” and is an organizer of the International Fraunhofer Lithography Simulation Workshop. He contributed to the development of several advanced lithography simulators including the development and research lithography simulator Dr. LiTHO. He is a fellow of SPIE.