## Abstract

We have developed a three-dimensional simulation algorithm based on fast marching method that mimics the etching behavior of chalcogenide photoresists, especially for maskless interference lithography. This lithography exposure is characterized by continuous variation of the exposure intensity inside the photoresist, without step like variation. Furthermore, the chalcogenide photoresist has a “gray-scale” behavior, without definite threshold. The resulting etching process is very sensitive to exposure dose and etching time. The optimal relations between these parameters are determined both theoretically and experimentally. A very good agreement between calculation and experimental results is shown, opening the door to complex nanostructures engineering.

© 2007 Optical Society of America

## 1. Introduction

The exposure of thin amorphous chalcogenide layers to light of suitable wavelength and intensity can result in significant structural changes due to their very unique physicochemical properties [1–4]. One of the most interesting properties of amorphous chalcogenide thin layers is its chemical resistance to alkaline solvents that can be modified by light exposure [5, 6]. This phenomenon, together with near-atomic resolution capability and high transmission in the IR region enables them to be applied as a medium for the fabrication of nano-structured optical elements such as gratings [7, 8] and filters. Additionally, chalcogenide glasses possess a very high refractive index, ranging from 2.3 to 3.5 which make them attractive materials for fabrication of photonic crystals [9]. Recently our group has demonstrated double exposure two-dimensional photonic crystal and layer-by-layer three-dimensional (3D) photonic crystal (wood-pile) structures using maskless interference lithography [10–12]. The interference lithography method is characterized by long range order periodic pattern, very high resolution, and infinite depth of focus. Therefore it is very attractive for fabrication of dense gratings and 2D and 3D photonic crystal.

In interference lithography, the exposure intensity inside the photoresist varies continuously, without sharp steps as in mask lithography. On the other hand, the etching rate of most photosensitive chalcogenide does not present the definite threshold characteristics of commercial organic photoresists, but is rather near to linear for lower exposure energies and saturates at high energies. Therefore it presents a gray-scale behavior. The combination of gray-scale intensity patterning and gray-scale etching behavior of the glass makes interference lithography patterning in chalcogenide material difficult to predict and therefore to engineer.

Currently there is no report of quantitative study of the etching process in chalcogenide glasses. Although simulation tools developed for virtual linear photosensitivity photoresist [17] or commercial photoresist like SU-8 [13] can be adapted, the specific etching of chalcogenide glasses at the sub-micron level has to be addressed specifically. The objective of this study is to present a comprehensive approach to the etching of chalcogenide glasses for sub-micron patterning.

In the first part of this paper we present the development of a dedicated algorithm based on the fast marching method for the simulation of the etching process. In the second part we describe the experimental conditions for the quantitative study. In the third part we present a detailed comparison between simulation and experimental results. We then propose a simple predictive model that is tested against experiments.

## 2. Developing of the etching model

#### 2.1 Light interference pattern

The exposure model determines the distribution of absorbed energy within the resist. Four coplanar beams interfere: two primary beams coming from the laser and two reflections from the substrate. The primary sinusoidal interference pattern results from the contribution of two TE polarized coherent laser beams. This image is disturbed by absorption within the chalcogenide layer and by air and substrate reflections. The indexes of refraction and the coefficient of absorption involved are measured and, in the simulation, the complex Fresnel coefficient is calculated to determine the amplitude and phase of the reflected beams.

#### 2.2 Formulation of Etching Model

Simulation of the etching process is a typical problem of front propagation. In this work we use the fast marching method developed by Sethian and Osher [14–19]. Recently, implementation of the method to lithography in SU-8 was presented [20–21]. Instead of trying to determine the position of the solvent front at every time moment, we look for the time the front arrives at any given point. The problem then reduces to the solution of a partial differential equation as explained in the following.

In the case of the “wet etching” relevant here, the front propagates in its normal direction with a local etching time *t _{0}·γ(x,y,z). γ* is a parameter without units which determines the relative etching time (the etching time of the illuminated area divided by the etching time of the non-illuminated area) and

*t*is the etching time of a unit length non-illuminated sample. Applying γ on the illumination pattern

_{0}*I*provides the relative local etching time

_{(x,y,z)}*γ*. The arrival time (

_{(x,y,z)}*T*) is then governed by the Eikonal equation:

The value of *t _{0}* is measured experimentally. The numerical scheme appropriate for the solution of this equation must be able to treat the possible cases when a sharp corner develops in the front, the front separates into a few pieces or different parts of the front merge. We use both the first-order and the second-order schemes. The first-order one is defined by the following finite-difference equation:

where ${D}_{\mathrm{ijk}}^{+x}T=\frac{{T}_{i+1,jk}-{T}_{\mathrm{ijk}}}{{x}_{i+1,jk}-{x}_{\mathrm{ijk}}}$ and ${D}_{\mathrm{ijk}}^{-x}T=\frac{{T}_{\mathrm{ijk}}-{T}_{i-1,jk}}{{x}_{\mathrm{ijk}}-{x}_{i-1,jk}}$ are the forward and the backward partial x-derivatives respectively; *D ^{±y}_{ijk}*

*T*and

*D*

^{±z}_{ijk}*T*are similarly defined for

*j*and

*k*indices respectively. The key point here is the correct choice of the right and the left derivatives. The second-order scheme is built on the same principles.

In the fast marching method Eq. (2) is solved for every *T _{ijk}* assuming that the neighboring grid values are given. An examination of the equation shows that

*T*depends only on the neighboring values which are smaller than it. This means that the numerical information flows only in the direction of the time growth. The points with smallest T are the initial points where

_{ijk}*T*=0. In the beginning the photoresist forms a homogeneous slab while the solvent is on top of it. The method then starts from the boundary, builds a narrow zone of trial points in its vicinity and propagates it in the direction of growing T. The first step of the algorithm is to label the boundary points as

*known*. The second step is to compute the values in the points neighboring to the

*known*ones. The values in these points are considered as

*trial*because they still can change in the following steps. The next step is to find the trial point with the smallest T. This point may be then frozen (considered

*known*) because its value can not change due to the updates in the neighboring points. Next if necessary, the algorithm updates the narrow zone of

*trial*points around the new

*known*point. It is important to notice that at this stage the values in all trial points surrounding the new

*known*one must be computed or re-computed since the new path may be faster. From now on the algorithm proceeds with the same scenario until all the grid points become

*known*. As soon as the arrival time T(x,y,z) is computed, the position of the front at any given moment of time is easily found through the isosurfaces of T(x,y,z).

Figure 1 presents a movie showing the time-development simulation of the etching process for a sinusoidal illumination of a chalcogenide layer. The simulation grid resolution along the substrate normal direction is 9.375 nm. The simulation grid resolution along the grating periodicity direction is 21.25 nm. The all volume size is 1.7×1.7×1.7 µm^{3}

The reflected beams are responsible for the modulations on the grating walls. The modulation pitch size is 100 nm (half wavelength inside the photoresist).

## 3. Experimental setup and photosensitivity measurement

#### 3.1 Samples preparation and photosensitivity measurement

For the fabrication of thin chalcogenide films (750 nm thickness) we thermally evaporated glassy 5(As_{2}S_{3}):1(As_{2}Se_{3}) in a high-vacuum evaporation chamber at a pressure of 5×10^{-6} mbar, and with a rate, measured with a quartz microbalance, of approximately 2–3 nms^{-1}. The absorption coefficient and index of refraction were measured using an Ellipsometer (V.A.S.E Series of J. A. Wollam Co.). At the writing wavelength 532 nm the index of refraction of 5(As_{2}S_{3}):1(As_{2}Se_{3}) composition is 2.79 and the absorption is approximately 10,000cm^{-1}.

Each set contained 17 samples that were thermally evaporated together, for maximum uniformity and for reproducibility analysis. In principle, chalcogenide etching properties depend not only on the total illumination energy, but also on the illumination power [4]. In order to avoid the use of additional parameter in this study, we kept the illumination time fixed at 200 sec, and we determined γ for different values of the illumination intensity. To measure γ, seven samples thermally evaporated together were used. Six of them were illuminated for 200 sec. The illumination intensity was varied and the total dose was in the 3–30 J/cm^{2} range. The etching time was measured for all the samples. Then, the ratio between the etching time of illuminated and non-illuminated samples was calculated. In Fig. 2 an experimental curve of γ values versus exposure dose is shown. This curve is slightly non-linear at low intensity and saturates at high intensity. In the simulation we use an analytic function that fits the experimental results. The best fit is given by

With *a=0.416, b=0.00716, n=1.92*. The exposure dose units are *J/cm ^{2}*.

#### 3.2 Writing setup

The set-up used to write the gratings is based on the interference of two CW laser beams (frequency doubled Nd:YAG, COHERENT VERDI-5). The incoming beam is split by a polarizing beam splitter into two equal intensity beams that form the arms of the interferometer with a sample mounted at the point where the two beams interfere. A λ/2 plate is added to one of the arms to set TE polarization in both beams. The angle at which the two beams overlap determines the spacing of the fringe pattern and sets the grating pitch Λ=*λ*(2·sin(*θ* 2)), where λ=0.532 µm is the laser wavelength and *θ* is the angle between the arms. Each beam is spatially filtered and then collimated to a diameter of 20 mm (FWHM). The sample size is 8×8 mm^{2}, smaller than the beam size, exposed only to the central flat part of the beam. The laser output power is limited to 900 mWatt (450 mWatt per beam) to prevent damage to the irises of the spatial filters.

The 0.532 µm writing wavelength of the laser is in the absorption tail of the 5(As_{2}S_{3}):1(As_{2}Se_{3}) chalcogenide material in order to enable effective photo-structural process as well as deep penetration of the laser light into the photoresist film. For example, 1 µm thick layer absorbs 60% of the light.

## 4. Experimental results and comparison with the simulation

The aim of this study is to experimentally validate the simulation algorithm, so that it can be used as a predictive tool for the fabrication of complex structures. Therefore, we have chosen two criteria for the comparison. The first criterion for comparison is the geometry of the structures. We compared the filling factor (the fraction of resist that remains after developing) and the shape (parallel walls or trapezium shape). The second criterion is the period of time during which the etched structure is stable inside the etching agent. For appropriate selection of the illumination intensity and material composition, the model predicts a stable period of time (dwell time) during which the structure remains almost unchanged during etching. After that, the walls get thinned until the structure breaks down.

For the comparative study, several sets of samples are illuminated and etched for different illumination energy doses and etching times. Then SEM cross sections of the etched samples are performed and analyzed. Finally we simulate the etching process for the same parameters using the function shown on Fig. 6. The following table summarizes the samples used in this work:

A sensitive factor for the comparison of the simulation with the experiment is the unavoidable set-up mechanical vibrations. These vibrations smooth the actual light intensity pattern and create an effective background that reduces the maximum possible contrast. We model these mechanical vibrations by convoluting the laser intensity profile in the sample with a Gaussian kernel of standard deviation of 30 nm.

#### 4.1 Shape and etching time comparison

In a first step, we focused on fully etched samples (Group II). Experimental and simulation results for one such sample are displayed in Fig. 3. We measured the wall thickness in the upper and lower surfaces, and compared the results with simulations. For a set of samples with total laser energy ranging from 12 J/cm^{2} to 22 J/cm^{2} the mismatch in the upper surface wall thickness was less then 10%. Comparing the lower surface walls thickness yields the same agreement.

To study the time evolution of the etching process we used non-completely etched samples (Group I). For each experimental condition, we run a simulation and recorded the time at which the simulated structure best fit the experimental SEM cross-section (Fig. 4). This time defined a simulated etching time that could be directly compared with the experimental etching time. The variation between these values was less then 10%.

## 5. Discussion

The goal of the previous section was to establish the validity of the simulation tool and to link it with the experimental results. In order to build a predictive tool, we conducted a series of simulations in which we recorded the time tE at which the solvent reaches the substrate through the etching process and the time tB at which the structure beaks down, for different illumination doses.

We then plotted tE and tB as a function of the dose (Fig. 6). These two graphs define a region of parameters, shaded in yellow on Fig. 6, for which the structure is stable. In this region, the dwell time t_{B} - t_{E} is positive and characterizes the structure stability: larger dwell time meaning better tolerance to parameter fluctuations. It can be seen that at about 30 J/cm^{2}, the dwell time is maximal. In order to be less sensitive to parameters fluctuations, the etching time should be chosen in the vicinity of (t_{E}+t_{B})/2.

Group I experimental samples (uncompleted etching) are plotted in this graph using square red dots. These dots are outside the yellow shaded region. Group II experimental samples (complete etching) are plotted using green dots and are located within the shaded region. Since group I and group II samples have been used for the simulation calibration of the mechanical vibrations, the good agreement between the simulation and the experiment is expected.

In order to test the simulation tool predictive power, we have used the third group of samples (group III). These samples have been prepared during previous studies and were not influenced by the present one. They corresponded to exposure dose and etching times that led to successful fabrications (not successful samples had been thrown away). We therefore plotted the parameters used for the samples preparation on the same graph (round white dots). They all fell within the shaded region.

## 6. Conclusion

There is an increased interest in chalcogenide materials as high index, high resolution photoresists for the fabrication of complex optical nanostructures such as two or three-dimensional photonic crystals. Naïve approaches based on intensity threshold for the prediction of the final shape of the nanostructure miss the dramatic effects many process parameters have on the final structure geometry and etching stability. We have shown that there is a delicate interplay between these parameters and that proper structure design should take these effects into account.

In order to do so, we have developed a simulation tool based on a fast marching method and compared the simulation with experimental results. We found an excellent agreement between both, provided that we introduced empirically mechanical vibrations of about 30 nm amplitude. In addition, we defined a predictive chart based on the simulation that allows fast determination of the optimal illumination and etching parameters for obtaining a stable process.

This chart has been developed for the case of simple sinusoidal illumination using continuous wave laser. However, the presented method is general and can be adapted to different chalcogenide composition and different types of physical mechanism for the light absorption, such as single or multi photon absorption, continuous or pulsed illumination. In addition, it can be easily adapted to more complicated interference structures such as two or three-dimensional photonic crystals.

## Acknowledgment

The authors express their gratitude to Prof. V. Lyubin and Dr. M. Klebanov from Physics Department, Ben-Gurion University, Beer-Sheva, Israel for providing the samples and for fruitful discussions.

## References and links

**1. **K. Shimakawa, A. Kolobov, and S. R. Elliott, “Photoinduced effects and metastability in amorphous semiconductors and insulators,” Adv. Phys. **44**, 475 (1995). [CrossRef]

**2. **K. TanakaA.V. Kolobov, ed.,*Photo-Induced Metastability in Amorphous Semiconductors* (Wiley, Weinheim, 2003) pp 69. [CrossRef]

**3. **A. Ozols and K. Shvarts, “Photosensitivity of amorphous semiconductor As-S and As-Se films under CW, nanosecond and picosecond laser irradiation,” Cryst. Latt. Def. and Amorph. Mat. **17**, 235–239 (1987).

**4. **G. Rosenblum, B. G. Sfez, Z. Kotler, V. Lyubin, and M. Klebanov, “Nonlinear optical effects in chalcogenide photoresists,” Appl. Phys. Lett. **75**, 3249 (1999). [CrossRef]

**5. **V. M. Lyubin, A. M. Sedikh, N. N. Smirnova, and V. P. Shilo, Microelectronica18, 523 (1989).

**6. **A. Arsh, M. Klebanov, V. Lyubin, L. Shapiro, A. Feigel, M. Veinger, and B. Sfez, “Glassy *m*As_{2}S_{3}·nAs_{2}Se_{3} photoresist films for interference laser lithography,” Opt. Mater. **26**, 301–304 (2004). [CrossRef]

**7. **M. Vlcek, P. J. S. Ewen, and T. Wagner, “High efficiency diffraction gratings in As-S layers,” J. Non-Cryst. Solids **227**–**230**, 743 (1998). [CrossRef]

**8. **A. V. Stronski, M. Vlcek, A. Sklenar, P. E. Shepeljavi, S. A. Kostyukevich, and T. Wagner, “Application of As40S60-xSex layers for high-efficiency grating production,” J. Non-Cryst. Solids **266**–**269**, 973 (2000). [CrossRef]

**9. **S. Wong, M. Deubel, F. Pérez-Willard, S. John, G. A. Ozin, M. Wegener, and G. von Freymann, “Direct Laser Writing of Three- Dimensional Photonic Crystals with a Complete Photonic Bandgap in Chalcogenide Glasses,” Adv. Mater. **18**, 265–269 (2006). [CrossRef]

**10. **A. Feigel, M. Veinger, B. Sfez, A. Arsh, M. Klebanov, and V. Lyubin, “Two dimensional photonic band gap pattering in thin chalcogenide glassy films,” Thin Solid Films **488**, 185–188 (2005). [CrossRef]

**11. **A. Feigel, Z. Kotler, B. Sfez, A. Arsh, M. Klebanov, and V. Lyubin, “Chalcogenide glass-based three-dimensional photonic crystals,” Appl. Phys. Lett. **77**, 3221 (2000). [CrossRef]

**12. **A. Feigel, M. Veinger, B. Sfez, A. Arsh, M. Klebanov, and V. Lyubin, “Three-dimensional simple cubic woodpile photonic crystals made from chalcogenide glasses,” Appl. Phys. Lett. **83**, 4480 (2003). [CrossRef]

**13. **R. C. Rumpf and E. G. Johnson, “Fully three-dimensional modeling of the fabrication and behavior of photonic crystals formed by holographic lithography,” J. Opt. Soc. Am. A **21**, 1703–1713 (2004). [CrossRef]

**14. **J. A. Sethian, “A fast marching level set method for monotonically advancing fronts,” Proc. Nat. Acad. Sci. **93**, 1591–1595 (1996). [CrossRef] [PubMed]

**15. **J. A. Sethian, *Level Set Methods and Fast Marching Methods*, (Cambridge Univ. Press, 2nd ed., 1999).

**16. **S. Osher and R. Fedkiw, *Level Set Methods and Dynamic Implicit Surfaces* (Springer, 2003).

**17. **J. A. Sethian and D. Adalsteinsson, “Overview of level set methods for etching, deposition, and lithography development,” IEEE Transactions on Semiconductor Devices **10**, 167–184 (1997). [CrossRef]

**18. **D. Adalsteinsson and J. A. Sethian, “A level set approach to a unified model for etching, deposition, and lithography. I. Algorithms and two-dimensional simulations,” J. Comp. Phys. **120**, 128–144 (1995). [CrossRef]

**19. **D. Adalsteinsson and J. A. Sethian, “A level set approach to a unified model for etching, deposition, and lithography II: three-dimensional simulations [integrated circuits],” J. Comp. Phys. **122**, 348–366 (1995). [CrossRef]

**20. **R. C. Rumpf, P. Srinivasan, and E. G. Johnson, “Modeling the fabrication of nano-optical structures,” Proc. SPIE **6110**, 611004 (2006). [CrossRef]

**21. **R. C. Rumpf and E. G. Johnson, “Comprehensive modeling of near-field nano-patterning,” Opt. Express **13**, 7198 (2005). [CrossRef] [PubMed]