Effective and flexible modeling approach to investigate various 3D Talbot carpets from a spatial finite mask

We present an effective modeling approach for a fast calculation of the Talbot carpet from an initially 2-dimensional mask pattern. The introduced numerical algorithm is based on a modiﬁed angular-spectrum method, in which it is possible to consider the border effects of the Talbot region from a mask with a ﬁnite aperture. The Bluestein’s fast Fourier transform (FFT) algorithm is applied to speed up the calculation. This approach allows as well to decouple the sampling points in the real space and the spatial frequency domain so that both parameters can be chosen independently. As a result an extended three-dimensional Talbot-carpet can be calculated with a minimized number of numerical steps and computation time, but still with high accuracy. The algorithm was applied to various 2-dimensional mask patterns and illumination setups. The inﬂuence of speciﬁc mask patterns to the resulting ﬁeld intensity distribution is discussed.


INTRODUCTION
The self-imaging properties of the Talbot effect was discovered in the early nineteenth century [1] and theoretically interpreted for the first time by Rayleigh [2].The Talbot effect occurs when a two-dimensional periodic mask is illuminated with a monochromatic light source.Due to diffraction effects, a series of self-images of the initial periodic structure occurs in the adjacent region behind the mask in multiples of a welldefined distance.
In optics, the self-imaging properties of the Talbot effect have attracted significant technological interest, as it allows the replication of an even high-frequent periodic structure without the need of any additional optical elements.
The unique characteristics of the Talbot effect opens up a broad spectrum of different applications.For example, Talbot interferometry allows to study the surface profile of transparent objects [3], the slope contours of bent plates [4] and also was used for collimation testing [5].Additionally, the selfimaging properties of the Talbot effect are also suitable for sub-wavelength focusing of light [6] and offer new solutions for spectroscopy [7,8].
An essential application field concerns lithography.In comparison to alternative laser based lithographic techniques, the Talbot lithography allows to combine advantages of both, se-rial and parallel processes, so that also micro-optical structures can be realized in a single exposure step [9].The utilization of two-photon effects in Talbot lithography allows the manufacturing of well-defined three-dimensional (3D) nanostructures [10].The application of Talbot Lithography in a mask aligner system with an adapted illumination system allows the manufacturing of small pitched structures and sub-micron resolution also in an industrial fabrication environment [11,12].With the introduction of displacement techniques, Talbot lithography offers the possibility for high-resolution patterning of large areas [13].
Based on the wide variety of the application spectrum of the Talbot effect it is necessary to get an accurate and fast prediction of the expected three-dimensional Talbot carpet.In the past, typical Talbot carpets are established along the x-and z-direction, which occurs from a one dimensional transfer function of the mask.Proceedings from these simulation results is the basic characteristic of the Talbot effect discussed like the fractional characteristic or typical carpet patterns behind amplitude or phase masks.The next stage of the developing the simulation algorithm was the calculation of the intensity distribution behind a twodimensional mask in a special distance z.But there were only x -y intensity layers presented not a complete three dimensional pattern [14]- [16].The challenge is the dealing with the high number of sampling points to gauge the transfer function of the mask and to sample the resulting x -y layers, which can be overlay to generate three dimensional profiles.
The here introduced algorithm offers the simulation of three dimensional intensity profiles of light behind a mask with a finite aperture.The base of the algorithm is the angular spectrum method.We combine the angular spectrum method with the Bluestein's FFT algorithm and show that the algorithm is able to calculate typical Talbot intensity patterns with less numerical effort.That means in detail that the sampling rate can be chosen independently in the spatial and frequency domain and also along the x and y axis as it required.In addition, it is possible to generate the intensity profile at the edge of the mask.With this feature it is possible to observe an adequate intensity profile for real lithography applications.

THE CHIRP-Z-TRANSFORMATION FOR MODELING FINITE TALBOT CARPETS
In the basic Talbot setup a two-dimensional mask is illuminated by an incoming monochromatic plane wave with a wavelength λ, which causes a complex field distribution in the adjacent region behind the mask.The mask possesses a finite aperture oriented parallel to the xy-plane in a Cartesian coordinate system, and the light wave propagates in the direction of the positive z-axis.
The complex field distribution which is generated directly behind the mask (z = 0) is represented by the scalar field U(x 0 , y 0 ; 0).In order to find the field distribution in any plane z behind the mask, a free-space propagation has to be performed, which is done in Fourier space.Figure 1 shows the schematic representation in both the real space and the spatial frequency domain.
We use the classical angular spectrum approach in combination with a chirp-Z-transformation [17] to be able to take finite boundaries of the mask aperture into account.Following the basic angular spectrum approach (cf.[18]), the initial field U(x 0 ,y 0 ;0) is, first, transferred into the spatial frequency domain by the two-dimensional Fourier transform: Here f x and f y represent the two axes in the spatial frequency domain.Secondly, with the background of the Helmholtz equation, the Fourier components F( f x , f y ; z) at z > 0 are accessible by the initial Fourier components at z = 0 where α and β are the direction cosines of the propagating waves, which read: Eqs.
(2) define the angular spectrum at any transverse plane perpendicular to the propagation direction z: The resulting field distribution in the real space at a position z > 0 can be calculated by superposing the propagated components, which is a transformation from the spatial frequency domain back to the real space: To find the extended three-dimensional Talbot carpet generated by arbitrary mask structures, numerical simulations have to be performed, in which the Fourier integral is commonly approximated by a discrete fast Fourier transform (FFT).Consequently, in the mask plane (z=0) and for each field position z > 0, all coordinate axes, both in real space (x, y) and in the spatial frequency domain ( f x , f y ), are specified by appropriate sampling points N x and N y .
In the modification of the standard free-space propagation, which we have solved using the angular spectrum method, two essential aspects are added in the presented calculation procedure: first, Bluestein's FFT algorithm, and, second, a fast convolution procedure.
The Bluestein's FFT algorithm comprises, in principle, a simple binomial equation for the x-and y-direction in the real and in the spatial frequency domain: The direct insertion of both binomial Eqs. ( 5) and ( 6) into Eq.( 4) yields: Due to the fact that now the differences x 0 − f x and y 0 − f y occur, the Fourier transformation can be interpreted in terms of a convolution applied to the two-dimensional space as follows where the functions g 1 and g 2 for the Bluestein algorithm in forward transformation are defined from Eq. ( 7): Because of the fact that g 2 is a product of either dimension, the 2D-convolution is separable, and, instead, a 1D-convolution can be applied successively to each of the real space coordinates of Eq. ( 2).In consequence, not a two-dimensional matrix convolution but a more convenient vector convolution has to be calculated, which reduces the numerical effort significantly.Rewriting Eq. ( 7) with the convolution from Eqs. ( 8)-( 10) results in the chirp-Z-tranformation of the initial field distribution: Based on this last formulation and in analogy to the preceding general procedure, the complex field distribution U(x, y; z) at a definite plane z > 0 is accessible by calculating the corresponding propagated Fourier components and applying the convolution approach for the reverse Fourier transformation.The calculation speed is significantly increased as compared to the 2D-Fourier transform without convolution.
This approach allows avoiding the correlated contributions of the real space and the spatial frequency domain, so that sampling points of the coordinate system for either domain can be chosen independently, and, in particular, the choice is not restricted to integer multiples of 2, as it is the case for the FFT.
For each of the coordinate pairs in the real space and in the spatial frequency domain, (x, y),(x 0 , y 0 ) and f x , f y , resp., we define specific variables M f , M 0 and M with: With this approach and the inverse transform we obtain the complex field distribution at z > 0: For the numerical calculation the initial complex field distribution U(x, y; z) at z = 0 has to be sampled.If a plane, monochromatic wave is incident perpendicular to the mask, the field distribution is assumed to be equal to the mask structure.With the chirp-Z-transformation outlined above, the sampling point N and range of the spatial domain can be adapted to the specific region of interest within the initial field and frequency domain.That means, it is possible to increase or reduce the number of the necessary sampling points with respect to the periodic structure of the initial field.
From this approach follows an optimum number of sampling points for each initial distribution, and, in particular, the sampling points in the frequency domain can be adapted to the specific modes appearing in the mask structure and effects of finite mask aperture can be taken into account without increasing the sampling points in the spatial domain.As a consequence, the computation time will be further reduced, and still the required accuracy of the Talbot pattern can be achieved.
An essential numerical aspect concerns the computation of the convolution.Here, an effective approach is based on the convolution theorem, which indicates that the convolution of two functions is identical to the inverse Fourier transform of the product of the Fourier transform of these two functions.In the discrete case of a numerical calculation, the convolution theorem can then be written as: Here g 1 and g 2 are the original functions in the real space.F is the discrete Fourier transform and F −1 is the discrete inverse Fourier transform, resp.In the discrete one-dimensional sum notation the Fourier transformation and the inverse transformation are: Here N represents the number of the sampling points, k is indicating the variable of the real space and n is the variable of the wave number or frequency domain, respectively.
The application of the discrete convolution theorem allows us to calculate the two-dimensional problem of the complex field distribution at z > 0 with a product of three subsequent FFT terms: 13004-3 J. Maaß, et al.The resulting field distribution U(x, y; z) we obtain again by an inverse Fourier transformation: U (x, y; z) = exp (iπ • M) Eq. ( 18) is the basis for our numerical simulation of Talbot carpets generated by finite mask structures.

SIMULATION OF 3D TALBOT CARPETS FROM FINITE MASKS
With the numerical algorithm derived in the previous section, we will demonstrate in the following the capability of our modelling approach.Thereby we will show that it is capable of simulating complex three-dimensional Talbot intensity profiles, taking into account boundary effects from the finite mask when calculating the Talbot carpet.We demonstrate the simulation with unequal sampling points along the x-and ydirection both in space and in the spatial frequency domain.Additionally, it is possible to enlarge the three dimensional intensity profile to get the essential information.
First of all we prove the correctness of the above presented algorithm by applying it to the example of referee [19].The ap-propriate mask has a square amplitude transfer function with a period of 280 µm and is illuminated with a plane wavefront of 632 nm.The classical Talbot distance is calculated with the equation z T = 2p 2 λ to be z T = 248 mm, for this example.The Figure 2a shows this mask pattern with squares of 70 x 70 µm.
The other three intensity profiles of Figures 2(b-d The FFT extends the mask aperture generally to infinity along both axes.However, the chirp-Z-transformation allows a defined finite aperture size.For the FFT calculation from  Z-transformation (Figure 2(c)), there is no significant difference in the alignment and shape of the intensity spots.Slight variations at the edge of the squared maxima become visible.The difference between the Figures 2(c,d) is considerable.With increasing distance to the centre, the spot maximum becomes more blurred.
In a second example, we applied the chirp-Z-transformation to calculate a complete three-dimensional (3D) intensity distribution in the adjacent region of a two-dimensional structured amplitude mask.The mask has a sinusoidal amplitude profile with a period of 20 µm in each transverse direction and is illuminated by a plane wave with a wavelength of λ=800 nm.For the calculation, the structure is covered with a grid combined of 1001 sampling points for both axes along the aperture of 1000 µm x 1000 µm see Figure 3(a).
Since a sinusoidal amplitude profile generates only the first and zeroth order, it is not necessary to regard the entire frequency spectrum.Accordingly, it is sufficient to take a frequency range of ±0.05 in the direction cosines with only 200 sampling points into account.In that case, the computation time can be reduced significantly compared to the classical FFT algorithm.The Figure 3(b) illustrates a homogenous periodic intensity distribution.The maxima energy is in the multiple of the Talbot distance.Between these maxima exists subtleties which can be influences the structure profile.To observe the 3Dintensity distribution in detail it is the profile of Figure 3(b) enlarged in Figure 4.An intensity maximum is visible in every half of the Talbot distance.These maxima are shifted by half of the period in the transverse direction similar to a body-centred cubic pattern.In addition, the volume structure posses interconnection between the half of the Talbot distance.In this case, the profile can be used for generating volume gratings in a lithographic process [20].Such structure features cannot be distinguished in a general view, like that of the Figure 3.For a practical application of a spatial finite mask in Talbot lithography, it is essential to precisely know the intensity distribution generated.In calculating this profile, only 110 sampling points are necessary.Hence, there is less numerical effort necessary for the simulation.With the algorithm it is possible to allocate more sampling points for a smaller observation range.
In addition to the ability to zoom into an intensity structure, the algorithm allows to perform simulations with independent sampling rates along the x-and y-axis in the spatial domain and f x -and f y -axis in the frequency domain.This feature increases the flexibility of the algorithm, additionally.It is particularly suitable for micro-optical components, which change their periodicity in only one direction.
A grating with a period of 20 µm possessing a binary amplitude transfer function is displayed in Figure 5.The mask is built with 650 sampling points along the x-axis and 250 points along the y-axis with a size of 400 µm x 160 µm.Therefore the mask has a restricted aperture.Such kinds of masks are used to structure, for example, optical fibers or their end surface [21].The small aperture of the mask strongly influences the resulting free-space intensity distribution.
The resulting frequency domain is also modulated in one direction only.Consequently, there are less numerical steps in the frequency domain along the y-direction necessary.This means for the calculation that the frequency space is scanned with 800 points along the x-and with 200 points along the yaxis.
Figure 6 shows the resulting intensity profile behind the red rectangle of Figure 5.This 3D-profile reaches from ±25 µm up to ±80 µm with 100 and 200 sampling points, resp.For this simulation it is useful to use different sampling rates in the spatial domain of the mask plane, the frequency domain and the resulting three dimensional spatial intensity profile.These fixed parameters have been selected to make the boundary effects of the spatial 3D-intensity distributions visible, which particularly occur behind a finite mask.The volume structure generated by the local intensity becomes smaller with larger distance along the z axis.These facts have to be considered in the setup of the experimental illumination assembly.

CONCLUSION
The chirp-Z-transformation provides an effective numerical tool to simulate the three-dimensional intensity distribution behind a lithography mask of finite aperture.The freedom of choosing the sampling rates in the transverse directions of both the spatial and frequency domain provides the ability to deal with very flexible mask geometries.This is particularly useful for the prediction of lithographic structures.The numerical algorithm presented here has specially been developed for the Talbot lithography and the generation of periodic microstructures.In addition, it is also useful for all other applications of generating Talbot interference patterns for illumination.
In the field of lithography or structured illumination it is essential to observe the three-dimensional range of homogeneity of the intensity profile to generate high-quality illumination conditions.Moreover, the algorithm has the potential to be easily extended to simulate the light propagation through an optical material, such as a photoresist, e.g., in which the refraction and absorption of light and the solubility of the resist must be taken into account.With this extension our algorithm can predict micro-or nanostructures in a polymer.

FIG. 1
FIG. 1 Coordinate systems and variables for both the real space and the spatial frequency domain.These systems are related to each other by Fourier transformation.

FIG. 2
FIG. 2 Mask z = 0 (a), FFT calculation by z = z T /2 (b) and the chirp-Z-transformation simulation with huge aperture (c) and small aperture (d) in the half of the Talbot distance.
) are simulated for the half of the Talbot distance.The Figure 2(b) shows the intensity distribution for the standard FFT algorithm.The other drawing (Figures 2(c+d) are generated with the introduced chirp-Z-transformation.The difference between the both figures is the adopted boundary condition.The Figure 2(c) is calculated with an entire aperture of 8.4 x 8.4 mm and is sampled with 2750 points in both axis directions.On these terms is the resulting intensity profile of Figure 2(c) homogenous.The Figure 2(d) is developed with a smaller mask aperture shown in Figure 2(a) so that boundary effects become visible.It has an aperture size of 1 x 1 mm 2 and is sampled with 1001 x 1001 points.
Figure 2(b) are 801 sampling points used for the generation of the mask, the frequency domain and the resulting intensity profile.It is sensible to take an odd sampling rate to detect correct the zero order in the frequency spectrum.Comparing the results of the FFT (Figure 2(b)) and the chirp-

FIG. 3
FIG. 3 (a) sinusoidal amplitude grating with p = 20 µm and (b) the corresponding three-dimensional Talbot carpet over 10 periods in x and y direction.

Figure 3 (
Figure 3(b) shows the resulting intensity profile starting directly behind the mask.For a better visualisation of the 3Dintensity carpet, an intensity threshold was introduced, which suppress values below 25% of the maximum intensity.The sampling grid for x and y is 200 x 200 sampling points, and 125 sampling points have been chosen in the z-direction.The z-axis in Figure 3(b) has been scaled to the Talbot distance.

FIG. 4
FIG.43D intensity profile with a smaller observation range from ±15 µm for the x and y axis.

J
FIG. 5 Binary amplitude mask with p = 20 µm and a size of 400 µm x 160 µm