- Email: [email protected]
Modelling of relief of phase reflection diffraction grating
Applied Mathematical Modelling 27 (2003) 1035–1049 www.elsevier.com/locate/apm
Modelling of relief of phase reflection diffraction grating Mindaugas Margelevicius
a,*
, Viktoras Grigali unas b, Henrikas Pranevicius a, Juozas Margelevicius c
a
Department of Business Informatics, Kaunas University of Technology, Studentu st. 50, LT-3031/3032 Kaunas, Lithuania Institute of Physical Electronics, Kaunas University of Technology, Savanoriu st. 271, LT-3009 Kaunas, Lithuania c Department of Graphic Communication Engineering, Kaunas University of Technology, Studentu st. 50, LT-3031 Kaunas, Lithuania
b
Received 7 January 2002; received in revised form 16 May 2003; accepted 10 June 2003
Abstract This article describes theoretical and experimental results on phase reflection diffraction gratings. Based on Fourier optics the mathematical formulation describing diffraction dispersion of light from a relief grating of the trapezium profile is derived. We propose a method that lets us estimate the gratingÕs geometric parameters in a versatile modelling system. We have designed an original programme that estimates diffraction intensities and calculates diffraction efficiency. The estimated intensities are used to reconstruct the gratingÕs geometrical properties using our mathematical model. The precision of the method is evaluated as the deviation of obtained results from microscopy data. 2003 Elsevier Inc. All rights reserved. Keywords: Diffraction grating; Diffractive optical element; Fourier optics; Levenberg–Marquardt
1. Introduction Recently, there has been a great deal of interest in microstructures [1] which exhibit periodic dielectric properties. In such materials, strong diffraction effects can inhibit the propagation of electromagnetic waves of certain frequencies. With appropriate three-dimensional symmetry and
*
Corresponding author. E-mail address: [email protected] (M. Margelevicius).
0307-904X/$ - see front matter 2003 Elsevier Inc. All rights reserved. doi:10.1016/S0307-904X(03)00137-9
1036
M. Margelevicius et al. / Appl. Math. Modelling 27 (2003) 1035–1049
sufficiently large spatial modulation of the dielectric, such systems can manifest a full photonic band gap; that is, they possess a spectral range over which propagation of the electromagnetic waves is inhibited regardless of the propagation direction within the periodic medium. It is the photonic band-gap crystals that are currently in the focus of our attention [2–4]. The advance of fabrication methods has given rise to new challenging problems. The smallest details need to be checked with, e.g., scanning electron microscopy (SEM) or atomic force microscopy (AFM), both of which require expensive equipment and a preprocessing of the samples. An alternative is to use non-destructive non-imaging methods, such as ellipsometry or optical scatterometry. These methods belong to the class of inverse scattering techniques: the collected information from the scattered radiation, e.g., visible light, is processed and properties of the scatterer are deduced. Scatterometry only measures the intensity of the scattered light and does not require expensive equipment. However, in this case a mathematical reconstruction of the microstructuresÕ geometry becomes necessary. There exist different mathematical models of dielectric microstructures. Methods for geometrical optics, although inaccurate, are often used as first approximations. This is because raytracing methods are computationally less demanding compared to algorithms used in wave optics. As diffractive optical elements (DOEs) get more complicated, the application of wave diffraction theory becomes necessary. Wave diffraction theory and scalar Fourier optics are becoming increasingly more important since the current fabrication methods can already produce structures with smallest features of the order of tens of nanometers. To solve wave diffraction theory equations in practice, however, an approximation is necessary. In one of them, distant wavefield may be described as a scalar if smallest features in the DOE are several times larger than one wavelength of light [5]. Determination of the geometrical properties of the microstructure from the intensities of the scattered light belongs to the class of inverse problems. Such problems are often more difficult to solve than the corresponding direct problems. However, the special data training and mining can be feasible in inverse modelling as this is shown in [5,6], where geometrical properties of gratings are estimated by means of the neural network. Another approach presented in [7,8] invokes Monte Carlo simulations for inelastic scattering of secondary electrons from the surface of a diffraction grating. The method is very accurate but it needs SEM equipment to produce experimental data and also a database of precompiled Monte Carlo scans. Unlike solution of inverse problems, solution of direct problems often conforms to analytical derivation of the equations [9]. Modified integral method [10] for relief gratings is accurate, and as it is shown in the reference, the method can be perfectly applicable for prediction of diffraction efficiency. However, sometimes it can be impossible to apply this method for inverse modelling. In this work, relief reflection DOEs are studied. A profile of its smallest topological features varies in form from a rectangular to trapezium with the oblique sides. In the article, we discuss a wave optics theory shortly. With a reference to Fraunhofer diffraction theory, we present the derived mathematical expressions describing the energy diffracted from a diffraction grating of the trapezium profile. We sketch the experimental setup and environment to be used in the investigations described in the article. We compare the results obtained by our method with SEM data.
M. Margelevicius et al. / Appl. Math. Modelling 27 (2003) 1035–1049
1037
2. Angular spectrum of plane waves In this section, we discuss Fraunhofer diffraction shortly. We describe the Fourier spectrum of diffraction field. It is shown that the Fourier spectrum can be expressed in terms of the reflection coefficient of infinite homogeneous periodic structure. We also show that diffraction intensity depends on the Fourier spectrum. We consider several reflection models and derive the expressions for intensities of a particular diffraction order. The derived expressions are used in a mathematical module we developed to find the geometrical parameters of diffraction gratings. On the basis of an electromagnetic theory starting with the Maxwell differential equations, it is possible to derive scalar approximations of microrelief diffraction gratings [11], while these can be developed into the methods of modern Fourier optics. Meanwhile kinematic approaches and dynamic diffraction theory offers its own accurate determination of light dispersion and distribution with the boundary conditions [12]. Modelling and solving of diffraction gratings by electromagnetic diffraction theory become a little simpler if we use finite calculus methods, methods of finite differences, Galerkin methods or methods of finite elements [13]. We will pass to the far-field formulation and proceed to the approximation of electromagnetic theory, namely Fraunhofer approximation. The research in Ref. [14] shows that results on the basis of vector and scalar diffraction theory are identical, if the diffraction fields are analysed at distances that are large compared to the dimensions of the DOEs of a grating. In the object and image plane, we consider small patches surrounding the point we analyse and associate them with the Fourier expansion, a superposition of plane waves [15]. Time coordinates are omitted in distribution we considered; the detailed formulation of space–time field distribution may be found in reference [16]. In the prediction of the surface morphology of microstructures, the field intensity values may be used in neural networks model [5,6]. The concept of intensity is tightly related to space-frequency distribution. A widely employed space-frequency representation is the Wigner distribution [17,18]. Let us consider the diffraction of the coherent light. Let us also suppose the light is being reflected from the homogeneous diffraction grating. Say that the beam of the plane waves falls obliquely to the grating, and the plane of the fall is arbitrarily oriented with respect to the strips of diffraction grating. The picture of such a scheme is shown in Fig. 1, where incidence of the beam is defined by angles h and /. The angular spectrum of the diffracted light is described by the expression [19]:
z
ηn
h
Θ
a
b
φ φn
x
Θn y Fig. 1. Geometrical reflection scheme of a diffraction grating of the trapezium profile.
1038
M. Margelevicius et al. / Appl. Math. Modelling 27 (2003) 1035–1049
F0 ðkx ; ky Þ ¼
Z Z
1 2
ð2pÞ
1
u0 ðx; y; 0Þ exp½iðkx x þ ky yÞ dx dy:
ð1Þ
1
The complex amplitude of the wave in the plane z ¼ 0 is described by the function u0 ðx; y; 0Þ. In fact, in the plane z ¼ 0 the integrand of (1) expresses the complex amplitude of the plane harmonious wave with the components of the wave vector kx , ky and qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð2Þ kz ¼ k 2 kx2 ky2 : The complex spectrum F0 ðkx ; ky Þ depends upon the distribution of direction (time-term multiplier expðixtÞ is omitted). Variables kz ¼ k cos /, ky ¼ k cos h, kz ¼ k cos g, where cos /, cos h, cos g are the direction cosines of the wave-front normal and k ¼ 2p=k is the wave number. The initial distortion of the angular spectrum depends upon the range of kx and ky , i.e. upon the width of the angular spectrum F0 ðkx ; ky Þ in the plane z ¼ 0. Let us introduce function Rðx; yÞ that characterizes reflection of the grating. Every particular grating is characterized by its own Rðx; yÞ. The function Rðx; yÞ is also supposed to be dependent upon the incidence angle of the electromagnetic waves. If we define the field of the falling wave in the plane z ¼ 0 as uf ðx; yÞ, then the field behind the grating is defined by the formula [19] u0 ðx; y; 0Þ ¼ uf ðx; yÞRðx; yÞ;
ð3Þ
and the angular spectrum behind the grating is equal to the convolution of the angular spectrum of the falling field and the spectrum of the reflection coefficient: Z Z 1 Ff ðn; gÞUðkx n; ky gÞ dn dg: ð4Þ F0 ðkx ; ky Þ ¼ Ff ðkx ; ky Þ Uðkx ; ky Þ ¼ 1
If the plane wave of unit amplitude is falling obliquely upon the grating, where the lurch is defined by the angles h and /, then its angular spectrum is described [19,20] Ff ðkz ; ky Þ ¼ dðkz k cos /Þdðky k cos hÞ;
ð5Þ
where dðxÞ is the Dirac delta function. As the plane wave reflects from the screen, its spectrum becomes wider. The Fourier spectrum of the reflection coefficient of infinite homogeneous periodic structure will gain the expression 1 X 2pn an d k x ð6Þ dðky Þ; Uðkx ; ky Þ ¼ d n¼1 where d is the period of the grating and an is the coefficient of the Fourier series of the reflection function of the grating. In the considered case the grating is periodic along one axis only, thus the function of the reflectance RðxÞ depends upon one spatial parameter. Consequently, we express coefficient an as Z 1 d=2 an ¼ RðxÞ exp½ið2pnx=dÞ dx: ð7Þ d d=2
M. Margelevicius et al. / Appl. Math. Modelling 27 (2003) 1035–1049
1039
Substituting Eqs. (5) and (6) into Eq. (4) and performing the integration, the angular spectrum of the field reflected from the grating is obtained [20]: 1 X 2pn an d kx k cos / ð8Þ dðky k cos hÞ: F0 ðkx ; ky Þ ¼ d n¼1 Comparing expression (8) with expression (5) we can see that the direction cosines of the diffracted waves might be related as sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ffi nk nk : ð9Þ cos /n ¼ cos / ; cos hn ¼ cos h; cos gn ¼ 1 cos2 h cos / d d The last relative formula cos gn is found from relation (2) for diffraction wave vector kn . From Eq. (9) it is seen that, for any diffraction order n, the equality kn j ¼ cos h holds true, where jð0; 1; 0Þ is the unit vector in the y-axis direction. Thus, all the diffraction orders draw up on the surface of the cone whose axis coincide with the ordinate of the analysed coordinate system, and the angle of the apex is equal to 2h. If we watch the diffraction view in the plane, perpendicular to the y-axis (i.e. to the strips of the grating) and at a distance l from the origin of the co-ordinates (where the beam falls into the grating), then all the diffraction orders will draw up in a semicircle of the radius r ¼ l tan h. In the plane y ¼ l, the coordinates xln and zln of the intensity of diffraction order n are expressed: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ffi, nk nk cos h: ð10Þ cos h; zln ¼ l sin2 h cos / xln ¼ l cos / d d The distributed quantity of the orders is found from Eq. (10) by applying the condition that the expression under the square root sign must be positive. Let us calculate the intensity of the diffracted light at the special case when / ¼ 90. In the case of such a model of illumination, the quantity RðxÞ might be defined quite easy [21]. The analysed grating is a phase structure formed from a non-transparent plate, where a required trapezium microrelief is formed after the etching process [2]. The reflection coefficients of the upper and lower surfaces of the microrelief of the periodic structure are RV eiuV and RA eiuA . The reflection coefficient of the oblique part of width a of the profile (Fig. 1) in an ordinary case is equal to zero (absolutely absorbent sides of the inclination). The phase difference of the waves reflected from the upper and lower surfaces of the microrelief is calculated D¼
4ph cos h ðuV uA Þ: k
ð11Þ
Then
8 < RV ; RðxÞ ¼ 0; : RA ; eiD
0 < x < B; B < x < ðB þ AÞ and ð1 AÞ < x < 1; ðB þ AÞ < x < ð1 AÞ;
ð12Þ
where B ¼ b=d, A ¼ a=d and b, a are the length of the upper base of the profile of the structure and the length of the projection of the inclination (Fig. 1), respectively; h is the height of the profile of the relief, k is the wavelength of the irradiated light.
1040
M. Margelevicius et al. / Appl. Math. Modelling 27 (2003) 1035–1049
The angular distribution of the intensity of the diffracted light is found from the angular 2 spectrum equation [22] according to the expression Iðkx ; ky Þ ¼ jF0 ðkx ; ky Þj . Then the intensity of 2 the nth diffraction order of the diffracted light In ¼ jan j . Defining the initial model of the reflection of the phase diffraction grating, according to Fraunhofer approach we can deduce a formula of the intensity of the diffracted light for the diffraction order of n: 1 ½R2V sin2 pnB þ R2A sin2 pnð2A þ BÞ 2RV RA sin pnB sin pnð2A þ BÞ cos D: ð13Þ In ¼ 2 ðpnÞ According to Eq. (10) we can state that, when the angle h of incidence decreases, the distance Dxln ðhÞ between the adjacent orders contracts in the plane y ¼ l. We will use this feature while estimating the diffraction orders whose dispersion of distribution is so big that the digital camera registering the image of diffraction does not record enough of them. The above defined reflection model (12) will introduce errors that will depend upon the length of projection on the x-axis of the profile side-walls (parameter a: see Fig. 1). Let us consider such a reflection model RðxÞ which would calculate on that the side-walls of the profile will change the phase of the incident light; let us also suppose that the amplitude of the reflection coefficient along the axes x and y are distributed evenly (RV ¼ RA ¼ RS ¼ R, where RS denotes the reflection amplitude of the side-walls of the profile). Let us also consider / ¼ 90, 8 R; > h i 0 6 x < B; > > > R exp i 4ph ðxBÞ ; B 6 x < B þ A; < k 4ph A RðxÞ ¼ ð14Þ B þ A 6 x < 1 A; R exp i > > h k ; i > > : R exp i 4ph ð1xÞ ; 1 A 6 x < 1: k A For two- and three-dimensional grating models, dependencies of derivatives of the amplitude and phase components on spatial orientation of an incidence beam are described and depicted in Ref. [23]. Using the above formulas we can deduce the exact expression of intensity of the order n for the given reflection model: ( R2 1 1 sin 2pnB cosð2H pnBÞ sin pnð2A þ BÞ þ sin 2ðH þ pnAÞ In ¼ 2 2 2ð1 þ H =pnAÞ ðpnÞ 2 1 sinðH pnAÞ cosðH pnð2B þ AÞÞ 1 H =pnA 1 þ ½cos 2pnB 1 sinð2H pnBÞ sin pnð2A þ BÞ 2 1 ½1 cos 2ðH þ pnAÞ þ 2ð1 þ H =pnAÞ 2 ) 1 ; ð15Þ sinðH pnAÞ sinðH pnð2B þ AÞÞ 1 H =pnA where H ¼ 2ph=k.
M. Margelevicius et al. / Appl. Math. Modelling 27 (2003) 1035–1049
1041
Fig. 2. The equipment used in the experiments; DSIC––digital still image camera.
3. Experimental setup Formulas (13) and (15) are essential for modelling process and are used to compose the main part of the mathematical module we use. The values of intensities in the formulas must be known, and a diffraction image must be analysed to estimate the diffraction intensities to get the values. To illustrate how the system works, we present the experimental setup. The experiment equipment (Fig. 2) consists of a He–Ne laser (k ¼ 632:8 nm), a digital camera registering the image, a computer, and a diffraction grating which is the specimen of the investigation. The experiment is performed as follows: (a) the laser irradiates a grating, (b) laser light diffracted from the grating is captured by the digital camera, (c) the image from the camera is transferred to the computer for further treatment of the data and production of the results. For data analysis we have created programme DIFFRAIN that analyses the image of diffraction and estimates the values of diffraction intensities. These values are as an input data for the mathematical module of the modelling system. This module is a computer programme that calculates/modulates the geometry of the diffraction grating using the estimated values of the diffraction intensities in its system of equations. To evaluate the accuracy of the computer-aided estimation of the diffraction intensities, in addition the intensities are estimated with a sensitive photo-sensor that gives indication in the changes of voltage.
4. The numerical results of estimation and calculation In this section, we discuss an inverse analysis shortly and experimentally examine accuracy of the proposed method. For each specimen we compare SEM data to the obtained modelling results when the diffraction intensities are estimated (1) with a photo-sensor and (2) by a computer programme. We assess reliability of the computer programme by comparing results to those
1042
M. Margelevicius et al. / Appl. Math. Modelling 27 (2003) 1035–1049
obtained when intensities were being estimated by the photo-sensor. We calculate and plot the deviations of the modelling results from the data observed by SEM. 4.1. The inverse analysis Let us introduce functions In ðA; B; H ; RÞ that define the right-hand sides of equations (13) or (15). With the help of In we can describe the ratio as Wmn ðA; B; HÞ ¼ Im ðA; B; H ; RÞ=In ðA; B; H ; RÞ. Here m and n are the values of the diffraction orders for which the ratio of intensities to be calculated. In a non-linear optimization problem with constraints where the basis is formed from the functions of Wmn ðA; B; H Þ, there are three unknowns. Therefore the ratio functions of 2/1, 3/1, 3/2 orders were included in the system of equations that composes the optimization problem. The highest diffraction orders are included for the following reasons: • most of the energy is often concentrated in those orders (lower relative error), • smaller amount of solutions is obtained in the data domain [22], • registration for the higher-order intensities is simpler. Solving the global optimization problem (Fig. 3) with constraints is based on the Levenberg– Marquardt (quasi-Newton) method [24,25]. We define the least squares error function as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi EðA; B; H Þ ¼ ðW21 ðA; B; H Þ z21 Þ2 þ ðW31 ðA; B; H Þ z31 Þ2 þ ðW32 ðA; B; H Þ z32 Þ2 ; ð16Þ where zmn ¼ Im =In is the ratio of the values of intensities measured with the equipment (programme or photo-sensor) at the mth and nth diffraction orders. The plot of this function is presented in Fig. 4. For the depiction of this distribution, the fourth sample is analysed (see Table 1) where the ratios of intensities zmn were estimated by the programme. The graphics is made when the quantity H was constant and it was chosen the height of the grating profile to be h ¼ 2:164 lm. Such a value was obtained as the height component of the solution (see Table 3). 4.2. The inverse problem of grating profile Six different diffraction gratings with different geometrical parameters were used for the experiments. In this research, three relative parameters of the grating were investigated A ¼ a=d, B ¼ b=d, H ¼ 2ph=k (Fig. 1), where d is the grating period. The values of the geometrical parameters A, B, H , d of the gratings investigated experimentally are listed in Table 1. The geometrical parameters of the gratings listed in the table were measured using SEM. The mean square deviations of the inverse problem solutions were computed with respect to the SEM data. 4.2.1. The registration of intensities with a photo-sensor The photo-sensor provides us with the intensities of the diffracted beams directly, there is no need for the analysis of diffraction image. However, measuring intensities with a photo-sensor is a
M. Margelevicius et al. / Appl. Math. Modelling 27 (2003) 1035–1049
1043
Y
Registration of a diffraction image
Computing the least squares error
Diffraction image Scan next diffraction image
Analysis of the diffraction image
The solution is precise
N
Y
Estimation of the diffraction intensities
Comparison of the modelling results to SEM data
N
Values of the intensities
Evaluation of ratio of the required intensities
Mathematical model based on Fourier optics
Solving the Global Optimization Problem
Ratios of the intensities
0.56
0.00
0.19
0.28
0.37
0.74
0.56
B
0.93
Fig. 3. The structural scheme of an analysis process of light diffraction.
0.00 A
0.84
Fig. 4. The distribution of the error function EðA; B; H ¼ Hc Þ when the height of the profile is h ¼ 2:164 lm. The intensity values are calculated for the sample no. 4.
1044
M. Margelevicius et al. / Appl. Math. Modelling 27 (2003) 1035–1049
Table 1 The parameters of the diffraction gratings used in the experiments Sample no.
Period, d, lm
SEM A
B
h, lm
1 2 3 4 5 6
5.5 5.5 4 4 2 2
0.01 0.045 0.01 0.01 0.059 0.01
0.5 0.218 0.5507 0.667 0.5 0.375
0.5 1.846 0.5 2.077 0.5 1.692
The first column represents the sample number, the second represents the grating period, and the parameters of the diffraction gratings measured by SEM follows next.
time- and effort- consuming process, and also requires special hardware as well as a stationary registration stand. The photo-sensor was used in the experiments to provide the mathematical module with the values of intensities. The obtained modelling results when the intensities were registered with the photo-sensor are listed in Table 2. The incidence angle h of the laser beam is indicated for every sample. The error of the gradient method (GME) is the error function EðA; B; HÞ calculated according to (16); the result error (RE) is the deflection (norm) of the solution in the data space from the original point with the coordinates ðA0 ; B0 ; H 0 Þ. If for the same sample we assign the obtained solution to the vector a with its components ðA; B; hÞ and SEM data to the vector b with its components ðA0 ; B0 ; h0 Þ, then the result error is calculated as e ¼ ka bk. All the solutions obtained using the formula (15) for calculation of the diffraction intensities of order n are listed in Table 2. For the first sample, when h ¼ 87, the least squares error E is equal to 0.195. This point with the coordinates (A ¼ 0:01, B ¼ 0:5, h ¼ 0:616) is close to the values measured by the SEM, and e ¼ 0:116. For the case of incidence angle of h ¼ 70, the result is worse for the same sample. However, E ¼ 0 for this case. The result error is not very small; and this could be explained by the oblique incidence of the beam, where according to the reflection model (14) angle h is ignored, while the angle of the abscissa remains / ¼ 90. Table 2 The results and A, B, h solutions for each sample, when the angle of incidence of the laser beam is perpendicular or oblique S. no.
h (deg)
Solution of the inverse diffraction problem A
B
h, lm
1 1 2 3 4 4 5 6 6
87 70 87 70 87 70 70 87 70
0.01 0.055 0.045 0.128 0.021 0.032 0.099 0.115 0.041
0.5 0.467 0.218 0.596 0.631 0.554 0.505 0.385 0.38
0.616 0.342 1.859 0.585 1.901 2.208 0.45 1.667 1.73
GME, E
RE, e
0.195 0 0 0 0 0 0 0 0
0.116 0.168 0.023 0.152 0.179 0.174 0.049 0.109 0.05
M. Margelevicius et al. / Appl. Math. Modelling 27 (2003) 1035–1049
1045
The values of error are minimal for the second and fifth samples. For the first of them h ¼ 87 and for the second h ¼ 70. For the third sample, the obtained meanings of B and h are close to the real measurements. However A is rather large: the distance between two adjacent protuberant segments of the profile is estimated to be close to the width of the sloping part of the trapezium profile and equal to 0.148. For all the experiments where the result is distant from the real meanings, the value of coordinate A influences the error e for the most part. 4.2.2. The estimation of intensities by the programme DIFFRAIN The programme DIFFRAIN analyses the recorded diffraction image and provides the mathematical module with the values of intensities. The solutions of the inverse diffraction problem when the intensities were estimated by the programme are given in Table 3. For the first sample, the least squares error E is larger than zero and the result error decreases to its smallest value e ¼ 0:116, which is close to the next smallest result e ¼ 0:119. The second solution is exact and the found values are not too different from the real ones: the greatest deviation is for the A component that makes up to 6.5% of the grating period. Using the scatterometry method, we obtain a ¼ 0:4125 lm, while the actual width of the sloping part of the profile is a ¼ 0:055 lm (a near rectangular shape). In this case the microrelief of the grating might be not of an exact shape what causes a different dispersion of the diffraction energy. For the second sample, the component B ¼ 0:266 of the third result is slightly diverged from the true value, the parameters A and h almost match the original meanings. The component B ¼ 0:219 of the second result matches the meaning of the parameter measured by the SEM, and the Table 3 The solutions of A, B, h for each sample, when the angle of incidence of the laser beam is perpendicular or oblique; diffraction intensities were evaluated by the computer programme DIFFRAIN S. no. 1 1 2 2 2 3 3 3 4 4 5 5 5 6 6 6* 6 6
h (deg) 90 90 90 90 90 90 90 90 90 90 90 90 90 90 90 56.25 56.25 56.25
Solution of the inverse diffraction problem A
B
h, lm
0.007 0.075 0.065 0.078 0.047 0.076 0.064 0.049 0.032 0.007 0.087 0.088 0.061 0.011 0.038 0.047 0.026 0.05
0.493 0.532 0.239 0.219 0.266 0.601 0.608 0.542 0.643 0.667 0.541 0.5 0.5 0.375 0.367 0.389 0.413 0.387
0.615 0.595 1.814 1.804 1.832 0.55 0.559 0.603 2.164 2.203 0.5 0.55 0.577 1.586 1.593 1.597 1.694 1.68
GME, E
RE, e
0.354 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0.116 0.119 0.043 0.053 0.05 0.097 0.099 0.11 0.093 0.126 0.05 0.058 0.077 0.106 0.103 0.103 0.042 0.044
1046
M. Margelevicius et al. / Appl. Math. Modelling 27 (2003) 1035–1049
Fig. 5. SEM photograph of the microrelief of the third sample; mark size 10 lm.
Fig. 6. The diffraction image of the third sample. Three highest diffraction orders are indicated with numbers, the central one is the intensity of the zeroth order.
parameters A and h are more diverged compared to the third result. For the second sample, the values of result error e are close to each other. The enlarged view of the third diffraction grating is shown in Fig. 5, and the recorded diffraction image to be computer-aided analysed is seen in Fig. 6. Three diffraction intensities of the highest order are enumerated in the image. This number of diffraction orders is sufficient for us to solve the inverse diffraction problem with three unknowns. The values of A and B obtained for the third result of the third grating are close to the real parameters but the height of the grating h are evaluated with a larger error. Meanwhile the values of the first result are such that the height is approximately equal to real data, and the components A and B are more diverged. However, the values of error e are close to each other. For the sixth sample, in the case of perpendicular illumination, the dispersion of the diffracted energy was greater and the diffraction intensities could be seen only up to the second order in the registered diffraction image. For this sample including only two diffraction orders, the calculations produced the solutions summarized in Table 3. At the incidence angle of h ¼ 56150 of the monochromatic laser beam, a distribution of the diffraction intensities was more condensed (10). The result denoted with 6 was obtained according to the formula (15) including more diffraction orders. The error e is not smaller than that received in the calculations including two values of the diffraction intensities. This might be explained by the fact that the assumed increasing precision while including three diffraction orders was diminished by the phenomenon of oblique radiation
M. Margelevicius et al. / Appl. Math. Modelling 27 (2003) 1035–1049 35
30 25
A
20 15
h
deviation,%
deviation,%
35 B
10 5
30 25
A
20
h
B
15 10 5
0
(a)
1047
0
1
2
3 4 Sample no.
5
6
(b)
1
2
3 4 Sample no.
5
6
Fig. 7. The percentage of deviation of the solutions from the real parameters for each sample, when diffraction intensities are estimated (a) with the photo-sensor, and (b) by the programme DIFFRAIN.
(the grating reflection model implicates its own errors), thus the result became similar to the previous. The last two results of the sixth sample were obtained by modelling the grating according to the reflection model (12), where the intensity values are expressed by Eq. (13). The values of error e of these last solutions are found to be the smallest ones. Let us conclude this chapter by comparing the results obtained by both trends (Fig. 7). Let us express deviation of the solutions from the real measurements in percentage. We express the deflections of the components A and B in percent of the length of the grating period. We also express the deflection of the height component h in percent of its real value. The calculated deflections from the observed data are depicted in Fig. 7. The figure represents the deviations for those solutions that are highlighted in Tables 2 and 3. For the first sample, in both cases (different estimation of intensities) agreement with the real data is observed with the largest error. In the case of photo-sensor estimation (Fig. 7(a)), the best solutions are for the second (d ¼ 5:5 lm) and the sixth samples (d ¼ 2 lm); and in the case of the programme estimation (Fig. 7(b)), they are for the second (d ¼ 5:5 lm), fourth (d ¼ 4 lm), and sixth samples (d ¼ 2 lm). For the relative quantity A the greatest deviation is 11.8% and 6.6% for the first and second direction of investigation respectively. For the relative quantity B the greatest deviation is 4.53% and 5.03% for the first and second direction respectively; and for the height quantity h the greatest deviation is 31.6% and 19% for the first (a) and second (b) direction respectively. The smallest deviations of the solved quantities are placed with pairs of values for the cases (a) and (b) respectively: DA ¼ f0%; 0:3%g, DB ¼ f0%; 0%g, Dh ¼ f0:7%; 0%g.
5. Conclusions We propose a method of estimation of geometrical parameters of relief reflection diffraction gratings by means of which the measuring of grating becomes non-destructive and faster process. The method is based on the modelling of relief of diffraction gratings and on solving an inverse scattering problem. The mathematical equations are derived only for the diffraction gratings of trapezium (and rectangular) profile. Consequently, the inverse modelling is valid and can be applied only for the gratings of that type. Despite such constraint on diffraction gratingsÕ shape,
1048
M. Margelevicius et al. / Appl. Math. Modelling 27 (2003) 1035–1049
we are able to obtain the explicit mathematical formulation for inverse modelling that can be easily solved by global optimization methods. The advantage of the proposed method is a possibility to automate the measurement process of optical control providing the communication between the computer programme for diffracted energy estimation, mathematical module, and global optimization module. We used a He–Ne laser stand in the experimental setup. The gratings were exposed to the laser beam, wavelength of which was 632.8 nm. It was determined that the modelling results correspond to SEM data with small error––average width deviation of the protuberant strips of the gratings with respect to grating period is 2.6%. When estimating the diffraction intensities by the computer programme DIFFRAIN, we obtained the modelling results that are not worse than those obtained by means of a photo-sensor.
References [1] D.M. Mittleman, J.F. Bertone, P. Jiang, et al., Optical properties of planar colloidal crystals: dynamical diffraction and the scalar wave approximation, J. Phys. Chem 111 (1) (1999) 345–354. . Meskinis, M. Margelevicius, I. Mikulskas, R. Tomasi [2] V. Grigali unas, V. Kopustinskas, S unas, Replication technology for photonic band gap crystals, Opt. Mater. 17 (1–2) (2001) 15–18. . Meskinis, M. Margelevicius, R. Tomasi [3] V. Kopustinskas, V. Grigali unas, D. Jucius, S unas, a-C:H Films for photonic structures fabrication, Proc. SPIE 4318 (2001) 295–300. . Meskinis, M. Margelevicius, R. Tomasi [4] V. Kopustinskas, V. Grigali unas, S unas, Formation of photonic band gap crystal by SF6 reactive ion beam etching, Lithuanian J. Phys. 40 (1–3) (2000) 140–143. [5] I. Kallioniemi, J. Saarinen, E. Oja, Optical scatterometry of subwavelength diffraction gratings: neural network approach, Appl. Opt. 37 (1998) 5830–5835. [6] I. Kallioniemi, J. Saarinen, E. Oja, Extraction of opto-geometrical properties of gratings from scattering data by means of the neural network, in: Diffractive Optics Õ97 12, EOS Topical Meetings Digest Series, 1997, pp. 202–203. [7] M. Davidson, N. Sullivan, Monte Carlo simulation for CD SEM calibration and algorithm development, Proc. SPIE 2439 (1995) 334–344. [8] M. Davidson, A. Vladar, Inverse scattering approach to SEM linewidth measurements, Proc. SPIE 3677 (1999) 640–649. [9] D. Content, Diffraction grating groove analysis used to predict efficiency and scatter performance, Proc. SPIE 3778 (1999) 19–30. [10] L.I. Goray, Modified integral method for weak convergence problems of light scattering on relief grating, Proc. SPIE 4291 (2001) 1–12. [11] J. Turunen, Diffraction theory of microrelief gratings, in: H.P. Herzig (Ed.), Micro-Optics: Elements, Systems and Application, Taylor & Francis, London, 1997, pp. 31–52. [12] B.A.
M. Margelevicius et al. / Appl. Math. Modelling 27 (2003) 1035–1049
1049
[20] B.B. Bokrod, K.K. Uepacbmob, B.B. Ragaeb, et al., BpM epeybe papM epO B 'keM eyT O B BC c yxeT O M peakmyO uO gpoabkz T padkeybz gO lbaparwbO yyO q rapT bye, Mbrpo'kerT poybra 13 (1) (1984) 64–73. [21] C.T..U. Typredbx, OgT bra lbaparwboyymI X 'keveyT O d b cbcnev, Maibyocnpoeybe, Leningrad (St. Peterburg), 1986. [22] M. Margelevicius, V. Grigali unas, H. Pranevicius, Modeling of lamellar transmission diffraction grating by solving the reverse diffraction problem, Inform. Technol. Contr. 18 (1) (2001) 52–59. [23] L. E'qn, Teopbz paccezvbz ya gpodokoxymI X peienrax b cenrax, in: A.H. Rokmouopob, C.G. Hobbrob (Eds.), XbckeymI e M eT O lmI T eO pbb lbaparwbb, Mbp, Moscow, 1982, pp. 113–142. [24] A. Tarantola, Inverse Problem Theory: Methods for Data Fitting and Model Parameter Estimation, Elsevier, Amsterdam, 1987. [25] D.M. Trujillo, H.R. Busby, Practical Inverse Analysis in Engineering, CRC Press, Boca Raton, FL, 1997.
Modelling of relief of phase reflection diffraction grating Mindaugas Margelevicius
a,*
, Viktoras Grigali unas b, Henrikas Pranevicius a, Juozas Margelevicius c
a
Department of Business Informatics, Kaunas University of Technology, Studentu st. 50, LT-3031/3032 Kaunas, Lithuania Institute of Physical Electronics, Kaunas University of Technology, Savanoriu st. 271, LT-3009 Kaunas, Lithuania c Department of Graphic Communication Engineering, Kaunas University of Technology, Studentu st. 50, LT-3031 Kaunas, Lithuania
b
Received 7 January 2002; received in revised form 16 May 2003; accepted 10 June 2003
Abstract This article describes theoretical and experimental results on phase reflection diffraction gratings. Based on Fourier optics the mathematical formulation describing diffraction dispersion of light from a relief grating of the trapezium profile is derived. We propose a method that lets us estimate the gratingÕs geometric parameters in a versatile modelling system. We have designed an original programme that estimates diffraction intensities and calculates diffraction efficiency. The estimated intensities are used to reconstruct the gratingÕs geometrical properties using our mathematical model. The precision of the method is evaluated as the deviation of obtained results from microscopy data. 2003 Elsevier Inc. All rights reserved. Keywords: Diffraction grating; Diffractive optical element; Fourier optics; Levenberg–Marquardt
1. Introduction Recently, there has been a great deal of interest in microstructures [1] which exhibit periodic dielectric properties. In such materials, strong diffraction effects can inhibit the propagation of electromagnetic waves of certain frequencies. With appropriate three-dimensional symmetry and
*
Corresponding author. E-mail address: [email protected] (M. Margelevicius).
0307-904X/$ - see front matter 2003 Elsevier Inc. All rights reserved. doi:10.1016/S0307-904X(03)00137-9
1036
M. Margelevicius et al. / Appl. Math. Modelling 27 (2003) 1035–1049
sufficiently large spatial modulation of the dielectric, such systems can manifest a full photonic band gap; that is, they possess a spectral range over which propagation of the electromagnetic waves is inhibited regardless of the propagation direction within the periodic medium. It is the photonic band-gap crystals that are currently in the focus of our attention [2–4]. The advance of fabrication methods has given rise to new challenging problems. The smallest details need to be checked with, e.g., scanning electron microscopy (SEM) or atomic force microscopy (AFM), both of which require expensive equipment and a preprocessing of the samples. An alternative is to use non-destructive non-imaging methods, such as ellipsometry or optical scatterometry. These methods belong to the class of inverse scattering techniques: the collected information from the scattered radiation, e.g., visible light, is processed and properties of the scatterer are deduced. Scatterometry only measures the intensity of the scattered light and does not require expensive equipment. However, in this case a mathematical reconstruction of the microstructuresÕ geometry becomes necessary. There exist different mathematical models of dielectric microstructures. Methods for geometrical optics, although inaccurate, are often used as first approximations. This is because raytracing methods are computationally less demanding compared to algorithms used in wave optics. As diffractive optical elements (DOEs) get more complicated, the application of wave diffraction theory becomes necessary. Wave diffraction theory and scalar Fourier optics are becoming increasingly more important since the current fabrication methods can already produce structures with smallest features of the order of tens of nanometers. To solve wave diffraction theory equations in practice, however, an approximation is necessary. In one of them, distant wavefield may be described as a scalar if smallest features in the DOE are several times larger than one wavelength of light [5]. Determination of the geometrical properties of the microstructure from the intensities of the scattered light belongs to the class of inverse problems. Such problems are often more difficult to solve than the corresponding direct problems. However, the special data training and mining can be feasible in inverse modelling as this is shown in [5,6], where geometrical properties of gratings are estimated by means of the neural network. Another approach presented in [7,8] invokes Monte Carlo simulations for inelastic scattering of secondary electrons from the surface of a diffraction grating. The method is very accurate but it needs SEM equipment to produce experimental data and also a database of precompiled Monte Carlo scans. Unlike solution of inverse problems, solution of direct problems often conforms to analytical derivation of the equations [9]. Modified integral method [10] for relief gratings is accurate, and as it is shown in the reference, the method can be perfectly applicable for prediction of diffraction efficiency. However, sometimes it can be impossible to apply this method for inverse modelling. In this work, relief reflection DOEs are studied. A profile of its smallest topological features varies in form from a rectangular to trapezium with the oblique sides. In the article, we discuss a wave optics theory shortly. With a reference to Fraunhofer diffraction theory, we present the derived mathematical expressions describing the energy diffracted from a diffraction grating of the trapezium profile. We sketch the experimental setup and environment to be used in the investigations described in the article. We compare the results obtained by our method with SEM data.
M. Margelevicius et al. / Appl. Math. Modelling 27 (2003) 1035–1049
1037
2. Angular spectrum of plane waves In this section, we discuss Fraunhofer diffraction shortly. We describe the Fourier spectrum of diffraction field. It is shown that the Fourier spectrum can be expressed in terms of the reflection coefficient of infinite homogeneous periodic structure. We also show that diffraction intensity depends on the Fourier spectrum. We consider several reflection models and derive the expressions for intensities of a particular diffraction order. The derived expressions are used in a mathematical module we developed to find the geometrical parameters of diffraction gratings. On the basis of an electromagnetic theory starting with the Maxwell differential equations, it is possible to derive scalar approximations of microrelief diffraction gratings [11], while these can be developed into the methods of modern Fourier optics. Meanwhile kinematic approaches and dynamic diffraction theory offers its own accurate determination of light dispersion and distribution with the boundary conditions [12]. Modelling and solving of diffraction gratings by electromagnetic diffraction theory become a little simpler if we use finite calculus methods, methods of finite differences, Galerkin methods or methods of finite elements [13]. We will pass to the far-field formulation and proceed to the approximation of electromagnetic theory, namely Fraunhofer approximation. The research in Ref. [14] shows that results on the basis of vector and scalar diffraction theory are identical, if the diffraction fields are analysed at distances that are large compared to the dimensions of the DOEs of a grating. In the object and image plane, we consider small patches surrounding the point we analyse and associate them with the Fourier expansion, a superposition of plane waves [15]. Time coordinates are omitted in distribution we considered; the detailed formulation of space–time field distribution may be found in reference [16]. In the prediction of the surface morphology of microstructures, the field intensity values may be used in neural networks model [5,6]. The concept of intensity is tightly related to space-frequency distribution. A widely employed space-frequency representation is the Wigner distribution [17,18]. Let us consider the diffraction of the coherent light. Let us also suppose the light is being reflected from the homogeneous diffraction grating. Say that the beam of the plane waves falls obliquely to the grating, and the plane of the fall is arbitrarily oriented with respect to the strips of diffraction grating. The picture of such a scheme is shown in Fig. 1, where incidence of the beam is defined by angles h and /. The angular spectrum of the diffracted light is described by the expression [19]:
z
ηn
h
Θ
a
b
φ φn
x
Θn y Fig. 1. Geometrical reflection scheme of a diffraction grating of the trapezium profile.
1038
M. Margelevicius et al. / Appl. Math. Modelling 27 (2003) 1035–1049
F0 ðkx ; ky Þ ¼
Z Z
1 2
ð2pÞ
1
u0 ðx; y; 0Þ exp½iðkx x þ ky yÞ dx dy:
ð1Þ
1
The complex amplitude of the wave in the plane z ¼ 0 is described by the function u0 ðx; y; 0Þ. In fact, in the plane z ¼ 0 the integrand of (1) expresses the complex amplitude of the plane harmonious wave with the components of the wave vector kx , ky and qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð2Þ kz ¼ k 2 kx2 ky2 : The complex spectrum F0 ðkx ; ky Þ depends upon the distribution of direction (time-term multiplier expðixtÞ is omitted). Variables kz ¼ k cos /, ky ¼ k cos h, kz ¼ k cos g, where cos /, cos h, cos g are the direction cosines of the wave-front normal and k ¼ 2p=k is the wave number. The initial distortion of the angular spectrum depends upon the range of kx and ky , i.e. upon the width of the angular spectrum F0 ðkx ; ky Þ in the plane z ¼ 0. Let us introduce function Rðx; yÞ that characterizes reflection of the grating. Every particular grating is characterized by its own Rðx; yÞ. The function Rðx; yÞ is also supposed to be dependent upon the incidence angle of the electromagnetic waves. If we define the field of the falling wave in the plane z ¼ 0 as uf ðx; yÞ, then the field behind the grating is defined by the formula [19] u0 ðx; y; 0Þ ¼ uf ðx; yÞRðx; yÞ;
ð3Þ
and the angular spectrum behind the grating is equal to the convolution of the angular spectrum of the falling field and the spectrum of the reflection coefficient: Z Z 1 Ff ðn; gÞUðkx n; ky gÞ dn dg: ð4Þ F0 ðkx ; ky Þ ¼ Ff ðkx ; ky Þ Uðkx ; ky Þ ¼ 1
If the plane wave of unit amplitude is falling obliquely upon the grating, where the lurch is defined by the angles h and /, then its angular spectrum is described [19,20] Ff ðkz ; ky Þ ¼ dðkz k cos /Þdðky k cos hÞ;
ð5Þ
where dðxÞ is the Dirac delta function. As the plane wave reflects from the screen, its spectrum becomes wider. The Fourier spectrum of the reflection coefficient of infinite homogeneous periodic structure will gain the expression 1 X 2pn an d k x ð6Þ dðky Þ; Uðkx ; ky Þ ¼ d n¼1 where d is the period of the grating and an is the coefficient of the Fourier series of the reflection function of the grating. In the considered case the grating is periodic along one axis only, thus the function of the reflectance RðxÞ depends upon one spatial parameter. Consequently, we express coefficient an as Z 1 d=2 an ¼ RðxÞ exp½ið2pnx=dÞ dx: ð7Þ d d=2
M. Margelevicius et al. / Appl. Math. Modelling 27 (2003) 1035–1049
1039
Substituting Eqs. (5) and (6) into Eq. (4) and performing the integration, the angular spectrum of the field reflected from the grating is obtained [20]: 1 X 2pn an d kx k cos / ð8Þ dðky k cos hÞ: F0 ðkx ; ky Þ ¼ d n¼1 Comparing expression (8) with expression (5) we can see that the direction cosines of the diffracted waves might be related as sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ffi nk nk : ð9Þ cos /n ¼ cos / ; cos hn ¼ cos h; cos gn ¼ 1 cos2 h cos / d d The last relative formula cos gn is found from relation (2) for diffraction wave vector kn . From Eq. (9) it is seen that, for any diffraction order n, the equality kn j ¼ cos h holds true, where jð0; 1; 0Þ is the unit vector in the y-axis direction. Thus, all the diffraction orders draw up on the surface of the cone whose axis coincide with the ordinate of the analysed coordinate system, and the angle of the apex is equal to 2h. If we watch the diffraction view in the plane, perpendicular to the y-axis (i.e. to the strips of the grating) and at a distance l from the origin of the co-ordinates (where the beam falls into the grating), then all the diffraction orders will draw up in a semicircle of the radius r ¼ l tan h. In the plane y ¼ l, the coordinates xln and zln of the intensity of diffraction order n are expressed: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ffi, nk nk cos h: ð10Þ cos h; zln ¼ l sin2 h cos / xln ¼ l cos / d d The distributed quantity of the orders is found from Eq. (10) by applying the condition that the expression under the square root sign must be positive. Let us calculate the intensity of the diffracted light at the special case when / ¼ 90. In the case of such a model of illumination, the quantity RðxÞ might be defined quite easy [21]. The analysed grating is a phase structure formed from a non-transparent plate, where a required trapezium microrelief is formed after the etching process [2]. The reflection coefficients of the upper and lower surfaces of the microrelief of the periodic structure are RV eiuV and RA eiuA . The reflection coefficient of the oblique part of width a of the profile (Fig. 1) in an ordinary case is equal to zero (absolutely absorbent sides of the inclination). The phase difference of the waves reflected from the upper and lower surfaces of the microrelief is calculated D¼
4ph cos h ðuV uA Þ: k
ð11Þ
Then
8 < RV ; RðxÞ ¼ 0; : RA ; eiD
0 < x < B; B < x < ðB þ AÞ and ð1 AÞ < x < 1; ðB þ AÞ < x < ð1 AÞ;
ð12Þ
where B ¼ b=d, A ¼ a=d and b, a are the length of the upper base of the profile of the structure and the length of the projection of the inclination (Fig. 1), respectively; h is the height of the profile of the relief, k is the wavelength of the irradiated light.
1040
M. Margelevicius et al. / Appl. Math. Modelling 27 (2003) 1035–1049
The angular distribution of the intensity of the diffracted light is found from the angular 2 spectrum equation [22] according to the expression Iðkx ; ky Þ ¼ jF0 ðkx ; ky Þj . Then the intensity of 2 the nth diffraction order of the diffracted light In ¼ jan j . Defining the initial model of the reflection of the phase diffraction grating, according to Fraunhofer approach we can deduce a formula of the intensity of the diffracted light for the diffraction order of n: 1 ½R2V sin2 pnB þ R2A sin2 pnð2A þ BÞ 2RV RA sin pnB sin pnð2A þ BÞ cos D: ð13Þ In ¼ 2 ðpnÞ According to Eq. (10) we can state that, when the angle h of incidence decreases, the distance Dxln ðhÞ between the adjacent orders contracts in the plane y ¼ l. We will use this feature while estimating the diffraction orders whose dispersion of distribution is so big that the digital camera registering the image of diffraction does not record enough of them. The above defined reflection model (12) will introduce errors that will depend upon the length of projection on the x-axis of the profile side-walls (parameter a: see Fig. 1). Let us consider such a reflection model RðxÞ which would calculate on that the side-walls of the profile will change the phase of the incident light; let us also suppose that the amplitude of the reflection coefficient along the axes x and y are distributed evenly (RV ¼ RA ¼ RS ¼ R, where RS denotes the reflection amplitude of the side-walls of the profile). Let us also consider / ¼ 90, 8 R; > h i 0 6 x < B; > > > R exp i 4ph ðxBÞ ; B 6 x < B þ A; < k 4ph A RðxÞ ¼ ð14Þ B þ A 6 x < 1 A; R exp i > > h k ; i > > : R exp i 4ph ð1xÞ ; 1 A 6 x < 1: k A For two- and three-dimensional grating models, dependencies of derivatives of the amplitude and phase components on spatial orientation of an incidence beam are described and depicted in Ref. [23]. Using the above formulas we can deduce the exact expression of intensity of the order n for the given reflection model: ( R2 1 1 sin 2pnB cosð2H pnBÞ sin pnð2A þ BÞ þ sin 2ðH þ pnAÞ In ¼ 2 2 2ð1 þ H =pnAÞ ðpnÞ 2 1 sinðH pnAÞ cosðH pnð2B þ AÞÞ 1 H =pnA 1 þ ½cos 2pnB 1 sinð2H pnBÞ sin pnð2A þ BÞ 2 1 ½1 cos 2ðH þ pnAÞ þ 2ð1 þ H =pnAÞ 2 ) 1 ; ð15Þ sinðH pnAÞ sinðH pnð2B þ AÞÞ 1 H =pnA where H ¼ 2ph=k.
M. Margelevicius et al. / Appl. Math. Modelling 27 (2003) 1035–1049
1041
Fig. 2. The equipment used in the experiments; DSIC––digital still image camera.
3. Experimental setup Formulas (13) and (15) are essential for modelling process and are used to compose the main part of the mathematical module we use. The values of intensities in the formulas must be known, and a diffraction image must be analysed to estimate the diffraction intensities to get the values. To illustrate how the system works, we present the experimental setup. The experiment equipment (Fig. 2) consists of a He–Ne laser (k ¼ 632:8 nm), a digital camera registering the image, a computer, and a diffraction grating which is the specimen of the investigation. The experiment is performed as follows: (a) the laser irradiates a grating, (b) laser light diffracted from the grating is captured by the digital camera, (c) the image from the camera is transferred to the computer for further treatment of the data and production of the results. For data analysis we have created programme DIFFRAIN that analyses the image of diffraction and estimates the values of diffraction intensities. These values are as an input data for the mathematical module of the modelling system. This module is a computer programme that calculates/modulates the geometry of the diffraction grating using the estimated values of the diffraction intensities in its system of equations. To evaluate the accuracy of the computer-aided estimation of the diffraction intensities, in addition the intensities are estimated with a sensitive photo-sensor that gives indication in the changes of voltage.
4. The numerical results of estimation and calculation In this section, we discuss an inverse analysis shortly and experimentally examine accuracy of the proposed method. For each specimen we compare SEM data to the obtained modelling results when the diffraction intensities are estimated (1) with a photo-sensor and (2) by a computer programme. We assess reliability of the computer programme by comparing results to those
1042
M. Margelevicius et al. / Appl. Math. Modelling 27 (2003) 1035–1049
obtained when intensities were being estimated by the photo-sensor. We calculate and plot the deviations of the modelling results from the data observed by SEM. 4.1. The inverse analysis Let us introduce functions In ðA; B; H ; RÞ that define the right-hand sides of equations (13) or (15). With the help of In we can describe the ratio as Wmn ðA; B; HÞ ¼ Im ðA; B; H ; RÞ=In ðA; B; H ; RÞ. Here m and n are the values of the diffraction orders for which the ratio of intensities to be calculated. In a non-linear optimization problem with constraints where the basis is formed from the functions of Wmn ðA; B; H Þ, there are three unknowns. Therefore the ratio functions of 2/1, 3/1, 3/2 orders were included in the system of equations that composes the optimization problem. The highest diffraction orders are included for the following reasons: • most of the energy is often concentrated in those orders (lower relative error), • smaller amount of solutions is obtained in the data domain [22], • registration for the higher-order intensities is simpler. Solving the global optimization problem (Fig. 3) with constraints is based on the Levenberg– Marquardt (quasi-Newton) method [24,25]. We define the least squares error function as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi EðA; B; H Þ ¼ ðW21 ðA; B; H Þ z21 Þ2 þ ðW31 ðA; B; H Þ z31 Þ2 þ ðW32 ðA; B; H Þ z32 Þ2 ; ð16Þ where zmn ¼ Im =In is the ratio of the values of intensities measured with the equipment (programme or photo-sensor) at the mth and nth diffraction orders. The plot of this function is presented in Fig. 4. For the depiction of this distribution, the fourth sample is analysed (see Table 1) where the ratios of intensities zmn were estimated by the programme. The graphics is made when the quantity H was constant and it was chosen the height of the grating profile to be h ¼ 2:164 lm. Such a value was obtained as the height component of the solution (see Table 3). 4.2. The inverse problem of grating profile Six different diffraction gratings with different geometrical parameters were used for the experiments. In this research, three relative parameters of the grating were investigated A ¼ a=d, B ¼ b=d, H ¼ 2ph=k (Fig. 1), where d is the grating period. The values of the geometrical parameters A, B, H , d of the gratings investigated experimentally are listed in Table 1. The geometrical parameters of the gratings listed in the table were measured using SEM. The mean square deviations of the inverse problem solutions were computed with respect to the SEM data. 4.2.1. The registration of intensities with a photo-sensor The photo-sensor provides us with the intensities of the diffracted beams directly, there is no need for the analysis of diffraction image. However, measuring intensities with a photo-sensor is a
M. Margelevicius et al. / Appl. Math. Modelling 27 (2003) 1035–1049
1043
Y
Registration of a diffraction image
Computing the least squares error
Diffraction image Scan next diffraction image
Analysis of the diffraction image
The solution is precise
N
Y
Estimation of the diffraction intensities
Comparison of the modelling results to SEM data
N
Values of the intensities
Evaluation of ratio of the required intensities
Mathematical model based on Fourier optics
Solving the Global Optimization Problem
Ratios of the intensities
0.56
0.00
0.19
0.28
0.37
0.74
0.56
B
0.93
Fig. 3. The structural scheme of an analysis process of light diffraction.
0.00 A
0.84
Fig. 4. The distribution of the error function EðA; B; H ¼ Hc Þ when the height of the profile is h ¼ 2:164 lm. The intensity values are calculated for the sample no. 4.
1044
M. Margelevicius et al. / Appl. Math. Modelling 27 (2003) 1035–1049
Table 1 The parameters of the diffraction gratings used in the experiments Sample no.
Period, d, lm
SEM A
B
h, lm
1 2 3 4 5 6
5.5 5.5 4 4 2 2
0.01 0.045 0.01 0.01 0.059 0.01
0.5 0.218 0.5507 0.667 0.5 0.375
0.5 1.846 0.5 2.077 0.5 1.692
The first column represents the sample number, the second represents the grating period, and the parameters of the diffraction gratings measured by SEM follows next.
time- and effort- consuming process, and also requires special hardware as well as a stationary registration stand. The photo-sensor was used in the experiments to provide the mathematical module with the values of intensities. The obtained modelling results when the intensities were registered with the photo-sensor are listed in Table 2. The incidence angle h of the laser beam is indicated for every sample. The error of the gradient method (GME) is the error function EðA; B; HÞ calculated according to (16); the result error (RE) is the deflection (norm) of the solution in the data space from the original point with the coordinates ðA0 ; B0 ; H 0 Þ. If for the same sample we assign the obtained solution to the vector a with its components ðA; B; hÞ and SEM data to the vector b with its components ðA0 ; B0 ; h0 Þ, then the result error is calculated as e ¼ ka bk. All the solutions obtained using the formula (15) for calculation of the diffraction intensities of order n are listed in Table 2. For the first sample, when h ¼ 87, the least squares error E is equal to 0.195. This point with the coordinates (A ¼ 0:01, B ¼ 0:5, h ¼ 0:616) is close to the values measured by the SEM, and e ¼ 0:116. For the case of incidence angle of h ¼ 70, the result is worse for the same sample. However, E ¼ 0 for this case. The result error is not very small; and this could be explained by the oblique incidence of the beam, where according to the reflection model (14) angle h is ignored, while the angle of the abscissa remains / ¼ 90. Table 2 The results and A, B, h solutions for each sample, when the angle of incidence of the laser beam is perpendicular or oblique S. no.
h (deg)
Solution of the inverse diffraction problem A
B
h, lm
1 1 2 3 4 4 5 6 6
87 70 87 70 87 70 70 87 70
0.01 0.055 0.045 0.128 0.021 0.032 0.099 0.115 0.041
0.5 0.467 0.218 0.596 0.631 0.554 0.505 0.385 0.38
0.616 0.342 1.859 0.585 1.901 2.208 0.45 1.667 1.73
GME, E
RE, e
0.195 0 0 0 0 0 0 0 0
0.116 0.168 0.023 0.152 0.179 0.174 0.049 0.109 0.05
M. Margelevicius et al. / Appl. Math. Modelling 27 (2003) 1035–1049
1045
The values of error are minimal for the second and fifth samples. For the first of them h ¼ 87 and for the second h ¼ 70. For the third sample, the obtained meanings of B and h are close to the real measurements. However A is rather large: the distance between two adjacent protuberant segments of the profile is estimated to be close to the width of the sloping part of the trapezium profile and equal to 0.148. For all the experiments where the result is distant from the real meanings, the value of coordinate A influences the error e for the most part. 4.2.2. The estimation of intensities by the programme DIFFRAIN The programme DIFFRAIN analyses the recorded diffraction image and provides the mathematical module with the values of intensities. The solutions of the inverse diffraction problem when the intensities were estimated by the programme are given in Table 3. For the first sample, the least squares error E is larger than zero and the result error decreases to its smallest value e ¼ 0:116, which is close to the next smallest result e ¼ 0:119. The second solution is exact and the found values are not too different from the real ones: the greatest deviation is for the A component that makes up to 6.5% of the grating period. Using the scatterometry method, we obtain a ¼ 0:4125 lm, while the actual width of the sloping part of the profile is a ¼ 0:055 lm (a near rectangular shape). In this case the microrelief of the grating might be not of an exact shape what causes a different dispersion of the diffraction energy. For the second sample, the component B ¼ 0:266 of the third result is slightly diverged from the true value, the parameters A and h almost match the original meanings. The component B ¼ 0:219 of the second result matches the meaning of the parameter measured by the SEM, and the Table 3 The solutions of A, B, h for each sample, when the angle of incidence of the laser beam is perpendicular or oblique; diffraction intensities were evaluated by the computer programme DIFFRAIN S. no. 1 1 2 2 2 3 3 3 4 4 5 5 5 6 6 6* 6 6
h (deg) 90 90 90 90 90 90 90 90 90 90 90 90 90 90 90 56.25 56.25 56.25
Solution of the inverse diffraction problem A
B
h, lm
0.007 0.075 0.065 0.078 0.047 0.076 0.064 0.049 0.032 0.007 0.087 0.088 0.061 0.011 0.038 0.047 0.026 0.05
0.493 0.532 0.239 0.219 0.266 0.601 0.608 0.542 0.643 0.667 0.541 0.5 0.5 0.375 0.367 0.389 0.413 0.387
0.615 0.595 1.814 1.804 1.832 0.55 0.559 0.603 2.164 2.203 0.5 0.55 0.577 1.586 1.593 1.597 1.694 1.68
GME, E
RE, e
0.354 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0.116 0.119 0.043 0.053 0.05 0.097 0.099 0.11 0.093 0.126 0.05 0.058 0.077 0.106 0.103 0.103 0.042 0.044
1046
M. Margelevicius et al. / Appl. Math. Modelling 27 (2003) 1035–1049
Fig. 5. SEM photograph of the microrelief of the third sample; mark size 10 lm.
Fig. 6. The diffraction image of the third sample. Three highest diffraction orders are indicated with numbers, the central one is the intensity of the zeroth order.
parameters A and h are more diverged compared to the third result. For the second sample, the values of result error e are close to each other. The enlarged view of the third diffraction grating is shown in Fig. 5, and the recorded diffraction image to be computer-aided analysed is seen in Fig. 6. Three diffraction intensities of the highest order are enumerated in the image. This number of diffraction orders is sufficient for us to solve the inverse diffraction problem with three unknowns. The values of A and B obtained for the third result of the third grating are close to the real parameters but the height of the grating h are evaluated with a larger error. Meanwhile the values of the first result are such that the height is approximately equal to real data, and the components A and B are more diverged. However, the values of error e are close to each other. For the sixth sample, in the case of perpendicular illumination, the dispersion of the diffracted energy was greater and the diffraction intensities could be seen only up to the second order in the registered diffraction image. For this sample including only two diffraction orders, the calculations produced the solutions summarized in Table 3. At the incidence angle of h ¼ 56150 of the monochromatic laser beam, a distribution of the diffraction intensities was more condensed (10). The result denoted with 6 was obtained according to the formula (15) including more diffraction orders. The error e is not smaller than that received in the calculations including two values of the diffraction intensities. This might be explained by the fact that the assumed increasing precision while including three diffraction orders was diminished by the phenomenon of oblique radiation
M. Margelevicius et al. / Appl. Math. Modelling 27 (2003) 1035–1049 35
30 25
A
20 15
h
deviation,%
deviation,%
35 B
10 5
30 25
A
20
h
B
15 10 5
0
(a)
1047
0
1
2
3 4 Sample no.
5
6
(b)
1
2
3 4 Sample no.
5
6
Fig. 7. The percentage of deviation of the solutions from the real parameters for each sample, when diffraction intensities are estimated (a) with the photo-sensor, and (b) by the programme DIFFRAIN.
(the grating reflection model implicates its own errors), thus the result became similar to the previous. The last two results of the sixth sample were obtained by modelling the grating according to the reflection model (12), where the intensity values are expressed by Eq. (13). The values of error e of these last solutions are found to be the smallest ones. Let us conclude this chapter by comparing the results obtained by both trends (Fig. 7). Let us express deviation of the solutions from the real measurements in percentage. We express the deflections of the components A and B in percent of the length of the grating period. We also express the deflection of the height component h in percent of its real value. The calculated deflections from the observed data are depicted in Fig. 7. The figure represents the deviations for those solutions that are highlighted in Tables 2 and 3. For the first sample, in both cases (different estimation of intensities) agreement with the real data is observed with the largest error. In the case of photo-sensor estimation (Fig. 7(a)), the best solutions are for the second (d ¼ 5:5 lm) and the sixth samples (d ¼ 2 lm); and in the case of the programme estimation (Fig. 7(b)), they are for the second (d ¼ 5:5 lm), fourth (d ¼ 4 lm), and sixth samples (d ¼ 2 lm). For the relative quantity A the greatest deviation is 11.8% and 6.6% for the first and second direction of investigation respectively. For the relative quantity B the greatest deviation is 4.53% and 5.03% for the first and second direction respectively; and for the height quantity h the greatest deviation is 31.6% and 19% for the first (a) and second (b) direction respectively. The smallest deviations of the solved quantities are placed with pairs of values for the cases (a) and (b) respectively: DA ¼ f0%; 0:3%g, DB ¼ f0%; 0%g, Dh ¼ f0:7%; 0%g.
5. Conclusions We propose a method of estimation of geometrical parameters of relief reflection diffraction gratings by means of which the measuring of grating becomes non-destructive and faster process. The method is based on the modelling of relief of diffraction gratings and on solving an inverse scattering problem. The mathematical equations are derived only for the diffraction gratings of trapezium (and rectangular) profile. Consequently, the inverse modelling is valid and can be applied only for the gratings of that type. Despite such constraint on diffraction gratingsÕ shape,
1048
M. Margelevicius et al. / Appl. Math. Modelling 27 (2003) 1035–1049
we are able to obtain the explicit mathematical formulation for inverse modelling that can be easily solved by global optimization methods. The advantage of the proposed method is a possibility to automate the measurement process of optical control providing the communication between the computer programme for diffracted energy estimation, mathematical module, and global optimization module. We used a He–Ne laser stand in the experimental setup. The gratings were exposed to the laser beam, wavelength of which was 632.8 nm. It was determined that the modelling results correspond to SEM data with small error––average width deviation of the protuberant strips of the gratings with respect to grating period is 2.6%. When estimating the diffraction intensities by the computer programme DIFFRAIN, we obtained the modelling results that are not worse than those obtained by means of a photo-sensor.
References [1] D.M. Mittleman, J.F. Bertone, P. Jiang, et al., Optical properties of planar colloidal crystals: dynamical diffraction and the scalar wave approximation, J. Phys. Chem 111 (1) (1999) 345–354. . Meskinis, M. Margelevicius, I. Mikulskas, R. Tomasi [2] V. Grigali unas, V. Kopustinskas, S unas, Replication technology for photonic band gap crystals, Opt. Mater. 17 (1–2) (2001) 15–18. . Meskinis, M. Margelevicius, R. Tomasi [3] V. Kopustinskas, V. Grigali unas, D. Jucius, S unas, a-C:H Films for photonic structures fabrication, Proc. SPIE 4318 (2001) 295–300. . Meskinis, M. Margelevicius, R. Tomasi [4] V. Kopustinskas, V. Grigali unas, S unas, Formation of photonic band gap crystal by SF6 reactive ion beam etching, Lithuanian J. Phys. 40 (1–3) (2000) 140–143. [5] I. Kallioniemi, J. Saarinen, E. Oja, Optical scatterometry of subwavelength diffraction gratings: neural network approach, Appl. Opt. 37 (1998) 5830–5835. [6] I. Kallioniemi, J. Saarinen, E. Oja, Extraction of opto-geometrical properties of gratings from scattering data by means of the neural network, in: Diffractive Optics Õ97 12, EOS Topical Meetings Digest Series, 1997, pp. 202–203. [7] M. Davidson, N. Sullivan, Monte Carlo simulation for CD SEM calibration and algorithm development, Proc. SPIE 2439 (1995) 334–344. [8] M. Davidson, A. Vladar, Inverse scattering approach to SEM linewidth measurements, Proc. SPIE 3677 (1999) 640–649. [9] D. Content, Diffraction grating groove analysis used to predict efficiency and scatter performance, Proc. SPIE 3778 (1999) 19–30. [10] L.I. Goray, Modified integral method for weak convergence problems of light scattering on relief grating, Proc. SPIE 4291 (2001) 1–12. [11] J. Turunen, Diffraction theory of microrelief gratings, in: H.P. Herzig (Ed.), Micro-Optics: Elements, Systems and Application, Taylor & Francis, London, 1997, pp. 31–52. [12] B.A.
M. Margelevicius et al. / Appl. Math. Modelling 27 (2003) 1035–1049
1049
[20] B.B. Bokrod, K.K. Uepacbmob, B.B. Ragaeb, et al., BpM epeybe papM epO B 'keM eyT O B BC c yxeT O M peakmyO uO gpoabkz T padkeybz gO lbaparwbO yyO q rapT bye, Mbrpo'kerT poybra 13 (1) (1984) 64–73. [21] C.T.