2.1

Theory of a Fresnel diffraction for N rectangular apertures

Let us consider the case of N rectangular parallel apertures (which is the basis of an amplitude diffraction grating) being illuminated by a light source with a wavelength λ at a distance p away. The observation screen is located at a distance q away. The width of each aperture is a and the opaque regions between apertures has a width b, the inter-aperture separation being (a+b). Each aperture has a height of c in the z-direction. The two edges of the central aperture of the N-aperture system (called aperture 0) are then located at y₀= - a/2 and y₀’=a/2 respectively, and the edges of the aperture for the next aperture to the right (called aperture +1) are located at y₁=a/2+b and y₁’=3a/2+b respectively. The next aperture to the right (called aperture +2) are located at y₂=3a/2+2b and y₂’=5a/2+2b. Finally the edges of the aperture +n will be located at y_n=(2n-1)a/2+nb and y_n’=(2n+1)a/2+nb. Similarly, the edges of the first left aperture (aperture -1) are located at y_-1’= -a/2-b and y_-1= -3a/2-b, and, the edges of the aperture -n will be located at y_-n’ = -(2n-1)a/2-nb and y_-n = -(2n+1)a/2-nb.

As shown in Fig. 1 (shown for the case of N=5), the coordinate systems on the aperture and on the image planes are chosen to be centered on the optical axis passing through the center of the central aperture (aperture 0) and normal to it, and are denoted by (y, z) and (Y, Z) axes, respectively. The Huygens−Fresnel principle is then invoked to calculate the total electric field at any given point of the image plane (Y,Z) by summing up all the contributions of all the elementary wavelets emitted by different area elements inside each of the rectangular apertures. The net complex electric field at the point P (located at the origin of the image plane or screen) due to waves emitted from the central aperture is given by^1,9

Figure 1.

Geometric configuration for Fresnel diffraction for N=5. a is the individual aperture width and (a+b) is the center-to-center aperture separation. The element at (y,z) is shown inside aperture 1, but it can be anywhere in any of the apertures. (Reproduced from Abedin et.al. Optik Vol.126. pp. 3743-3751, Copyright (2015) with permission from Elsevier)

where ε₀ is the electric field of the source, k is the wavenumber (=2π/λ₀) and K(θ) is the obliquity factor, which accounts for secondary Huygens wavelets emitted from area elements dS in an oblique direction. By integrating the contributions of the Huygens wavelets over the entire area of the central aperture, it can be shown^1,9 that the total electric field at P is,

where E_u is the unobstructed electric field at P (i.e. the electric field that would have existed if the entire aperture were absent), and C(u) and S(u) are the Fresnel cosine and sine integrals, being defined by,

Here w represents either of the two dimensionless variables u or v,

The variables u and v are proportional to Cartesian coordinates x and y. The limits u₀ and u₀’ in Eq. (2) are the values of u corresponding to two left and right edges of the central aperture, i.e. for y₀=-a/2 and y₀’=a/2, respectively. Similarly, the limits v’ and v are the values of the variable v corresponding to the lower and the upper edges of this aperture i.e. for z =-c and z’=+c respectively.

For the rest of these Proceedings, for simplicity, we consider plane-wave illumination of the aperture. This will effectively move the source S to infinity, which can be achieved by placing S in the focal point of a positive (convex) lens and allowing the collimated light from the lens to fall on the apertures. For the simple case of plane-wave illumination, p₀ is effectively infinity and, equations (4) can be simplified to,

Upon illumination from the source, the total electric field at the center P of the image plane consists of the contributions from all of the (2n+1) apertures. The electric field contribution from the aperture 1 is given, in analogy to Eq. (2), by

Here, the limits u₁ and u₁’ are the values of the dimensionless variable u corresponding to two edges of aperture 1, i.e. for y₁=a/2+b and y₁’=3a/2+b, respectively. Similarly, the limits v and v’ are the values of v for z=-c and z’=+c respectively, as before. The electric field contribution by the aperture n is given by

where u_n and u_n’ are the values of u corresponding to two edges of aperture n, i.e. for y_n=(2n-1)a/2+nb and y_n’=(2n+1)a/2+nb, respectively. The electric field contributions by the apertures on the other side, (e.g., for the -1 and –n apertures denoted by E_P-1 and E_P-n respectively) are given by equations similar to Eqs. (6) and (7), where all the limits are correspondingly negative. The net complex electric field at P contributed by the all the N (=2n+1) apertures is given by the arithmetic sum of all complex amplitudes E_Pn, …. E_P1, E_P0, E_P-1, …..E_P-n; i.e. by,

Separating the cosine and sine integrals, and using the summation notation for the Fresnel sines and cosines, we can write the total complex electric field at P as,

The summation from –n to +n for both the Fresnel cosine and sine integrals (for the u variable only), in effect, carries out the summation of the complex electric field contributions from aperture –n to aperture n, a total of N=(2n+1) apertures in the system. No such summation is involved for the v variable, since the two edges (upper and lower) of the apertures are involved in this z-direction. The net intensity at P is proportional to (EE*), i.e.,

where I₀ denotes the intensity of the unobstructed wave due to E_u, i.e. I₀= E_u².

To calculate the intensity at any off-axis point P' on the image plane, one can fix the SOP line and instead of moving the point P, one can move the entire aperture system by small amounts in the yz plane, so that the relative positions of the aperture edges in the new position and P remains unchanged. For example, to find the intensity a point P' 1mm below P, one can keep the screen undisturbed and move the apertures 1mm upwards, and find the intensity, at P, in this situation instead. The point P will now see a new set of values for z, and therefore, for v in Eq.(8). In principle, the intensity at any point P' on the image plane can be found in this way by making suitable translations of the aperture in the y and z directions, and in Eq.(8), substituting the corresponding values of the edges of each aperture in the u direction u_n′..n₀′.…u_-n′ and u_n..n₀..n-_n, (which indicate the positions of the edges of the aperture in the y-direction as seen from P) and v′ and v, (which indicate the positions of the edges of the apertures in the z-direction). This means that new values of Fresnel cosine and sine integrals need to be calculated. Using this method, the entire intensity distribution in the image plane can be mapped out in principle. From Eqs. (8) or (9), it is clear that the calculation of electric field or intensity at a point P requires the evaluation of 2N Fresnel cosine and 2N Fresnel sine integrals (corresponding to the 2N edges of the N aperture system for the u(y) variable), plus two Fresnel cosine and sine integrals for the v(z) variable. The cosine integrals form the real parts of the electric field, and the sine integrals form the imaginary parts. After calculating the complex electric field, the intensity at P is obtained in Eq. (10) by multiplying the field E_N by its complex conjugate.

These equations are the basis of calculation of the complete intensity distribution of the Fresnel diffraction pattern from an N-aperture, as will be explained next.

2.2

Simulation strategy and methodology

The flow chart of the algorithm is shown in fig. 2. For calculation of the electric field distribution using Eq. (10) at all the points (pixels) on the image plane, a large number of the Fresnel cosine and sine integrals will need to be evaluated quickly and efficiently. This is accomplished in the MATLAB program using the special functions mfun (‘FresnelC', i:s:f) and mfun (‘FresnelS', i:s:f). For example, the MATLAB statement R=mfun (‘FresnelC', i:s:f) generates an array R of the Fresnel cosine integrals with arguments staring from i and ending in f, with a step size s.

Figure 2.

Flow chart of the algorithm for the N-aperture problem. (Reproduced from Abedin et.al. Optik Vol.126. pp. 3743-3751, Copyright (2015) with permission from Elsevier)

As shown in fig. 3 for N=3 (three apertures), the aperture was displaced (virtually) by an amount W, and the extreme limits of the 2N edges of the displaced aperture was determined and compared to the corresponding limits for the un-displaced aperture to find the range of u and v values. In a more general N-aperture system, when the aperture is moved by W to the left, the ranges for the 2N edges of the N=(2n+1) apertures will be established as:

Figure 3.

Apparent limits of the virtually displaced aperture seen from the observation plane for N=3. (Reproduced from Abedin et.al. Optik Vol.126. pp. 3743-3751, Copyright (2015) with permission from Elsevier)

Arrays of Fresnel cosine and sine integrals need to be calculated corresponding to these input ranges by the mfun statements, with a given step size s. By using these Fresnel arrays, the electric field or intensity dependence in the y (u) direction can be calculated (first factor in Eq. 9). Following a similar procedure for the entire aperture, i.e. by moving the aperture by W downwards, the z′ and z ranges are determined: c to (W+c) for z′ and –c to (W-c) for z. Four arrays of Fresnel integrals are to be evaluated for these two ranges, giving numerical values for the calculation of the second factor in Eq. (9). Further details of the computation method and explanation of the MATLAB program is found in Ref. (8) and Ref. (5).

3.1

General

To use the program, the values of the required inputs, i.e. aperture parameters a and b, height c, image area W, step size s, the aperture-image distance q₀ (all in mm), the wavelength l (in nm), the exposure and the number N of the slits (apertures) are input through a GUI (fig. 4a). The exposure factor ex is used to adjust the apparent visual intensity in the generated image. The computer generated Fresnel image due to N-slits for wavelength λ=500 nm, q₀=400mm, image area W=5mm, aperture width a=0.1 mm, b=0.1 mm (with aperture separation 0.2 mm), c=4mm, s=0.01mm and N=31 is shown in fig. 4b. Since (a+b)=0.2mm, this system is equivalent to an amplitude grating N=31 with 5 lines/mm. The generated diffraction image resembles that of a grating in the near-field, with Fresnel-like characteristics exhibited in the top and bottom edges. Some interference effects can be seen in the left and right edges of the image.

Figure 4.

(a) Screenshot of the MATLAB GUI, showing the 9 input parameters. (b) The computer simulated Fresnel image from N=31 apertures. The input parameters correspond to those in the GUI. (Reproduced from Abedin et. al. Optik Vol.126. pp. 3743-3751, Copyright (2015) with permission from Elsevier)

To find the effect of a change of parameters on the diffracted images, we generated a series of images. In the first series, the aperture separation (a+b) was made successively smaller, while keeping the aperture width (a=0.5mm) constant. A series of diffraction images were simulated as shown in figs 5(a)-(d) for a value of N=7. In fig. 5(a), in the region between the apertures, no significant interference could be seen between the light from each of the aperture diffraction patterns. But, in fig. 5(b) and (c), some distinct interference between diffracted light from the individual apertures was observed as the separation of the slits is further reduced. As the separation was decreased to zero [fig. 5(d) for b=0], a typical single aperture Fresnel diffraction image was obtained for an aperture size 3.5mm×3mm.

Figure 5.

Computer-simulated Fresnel diffraction images for decreasing (center-to-center) aperture separation (a+b), while keeping the aperture width (a =0.5mm) constant. for N=7. Other parameters are: λ=500 nm, q0=400mm, W=8 mm, s=0.01mm and c=3mm. Fig. 5(d) corresponds to b=0 (single aperture). (Reproduced from Abedin et.al. Optik Vol.126. pp. 3743-3751, Copyright (2015) with permission from Elsevier)

In the next series, for N=7, the width a of the apertures was made successively narrower while their separation (a+b) was kept constant at 1.5 mm [fig. 6]. From figs 6(a) and (b), as the apertures become narrower, the light is diffracted over a wider region, and some interference between diffracted light can be detected. In figs. 6(c) and (d), at smaller aperture widths, strong mutual interference between diffracted light has generated fringes which look like typical Young-like fringes. Light from each aperture is seen to be diffracted over a wider portion of the image area. In figs 6(c) and (d), fringe density is the same, since this depends only on aperture separation (constant at 1.5mm).

Figure 6.

Computer-simulated Fresnel diffraction images for decreasing aperture width a, while keeping the aperture separation (a+b)=1.5mm constant for a seven-slit system. λ=500 nm, q₀=400mm, s=0.01 mm and c=3mm. (Reproduced from Abedin et.al. Optik Vol.126. pp. 3743-3751, Copyright (2015) with permission from Elsevier)

3.2

Comparison with N-slit far-field diffraction

If the aperture-screen distance q₀ is increased, it is well-known in optics that a gradual transition from the Fresnel (near-field) regime to Fraunhofer regime (far-field) should occur. As a general rule of thumb, if D is the lateral dimension of an aperture, then Fresnel diffraction occurs if D satisfies the following relation¹

Otherwise, Fraunhofer diffraction should be observed. In the next simulations, we use a practical diffraction grating with a=0.001 and b=0.003 (with aperture separation a+b=0.004 mm, i.e. 250 lines/mm). We selected the following parameters for an 11-slit system, λ=500 nm, c=2mm, image area size W variable, s=0.01 mm or 0.05 mm, as appropriate. We then increased the aperture screen distance q₀ in steps from 1 mm to 100 mm while keeping other parameters constant [fig. 7]. A transition from the Fresnel regime to Fraunhofer regime is clearly observed, being consistent with the predictions of Eq. (11). For example, for a 11-slit system with (a+b)=0.004mm, if we take the total lateral dimension of the slit system to be D=0.04 mm, then D²/λ is about 3.2 mm. We expect the diffraction to be Fresnel-like if q₀ is less than 3.2 mm. If, q₀> (D²/λ)=3.2mm, and we expect the diffraction to be Fraunhofer-like. In between these extremes, a transition from Fresnel to Fraunhofer regime should occur [as in figs 7(c) and 7(d)]. The usual mathematical analysis of an N-slit system^1,9 deals with the N-slit pattern only in the far-field Fraunhofer regime. The far-field diffraction pattern of an N-aperture system in the image plane (Y, Z) can be exactly calculated by an analytical formula^1,9

Figure 7.

Computer simulated Fresnel diffraction images for a 7-aperture system. Aperture-screen distance q₀ is increased in steps, while keeping other parameters constant, with a=0.001mm, b=0.003mm, c=2mm, λ=500 nm. (a) and (b) are in the Fresnel regime, while (e) and (f) are in the Fraunhofer regime, according to Eq. (12).

where β, γ and α are defined as,

Here, a is aperture width, (a+b) is the inter-aperture separation, c is aperture height, λ is the wavelength and R is the aperture-screen distance (assumed to be sufficiently large). As stated above, in the computer simulations, this far-field conditions should be amply satisfied in figs 7(e) and 7(f), for R (or q₀)=50 mm and R=100 mm, respectively.

In the usual analysis of far-filed Fraunhofer diffraction from an N-slit system,⁹ it is shown that the principal maxima of the diffraction pattern occurs when the [sin Nγ/sin γ]² factor have maximum values, i.e. when,

Here, θ is the diffraction angle and m is the diffraction order. For m=0, we have the zeroth-order, for m=1, the first order, and so on. In addition to these principal maxima, several much weaker secondary maxima, and several minima of zero intensity are predicted to occur with equal separation between the principal maxima. If N is the number of grating slits (apertures), then the number of secondary maxima is equal to (N-2), and the number of minima is equal to (N-1).⁹ To make comparisons, we reproduce fig. 7(f) next, with the intensity values multiplied by 10 to bring out the faint details. We observe 7 principal maxima, with m=0 (at Y=0), with m=+1 (at Y=12.5 mm), m=+2 (at Y=25.03 mm), m=+3 (at Y=37.59 mm), with m=-1 (at Y=-12.5 mm), m=-2 (at Y=-25.0 mm), and m= -3 (at Y=-37.5 mm). In addition, 9 (i.e., N-2) secondary maxima can be discerned between the principal maxima at the zeroth order and the first and the second orders, along with 10 minima of zero intensity, interspersed between them. The values of sin θ can be estimated by finding the ratios (Y/R, with R=100 mm) from the computer-simulated image (Fig. 8). These calculated values are: 0.125 (for m=+1 and m= -1), 0.2503 (for m=+2 and m= -2), and 0.375 (for m=+3 and m= -3) From Eq. (14), on the other hand, for λ =500 nm, (a+b)= 0.004 mm, and the following values of sin θ are calculated: 0.125 (for m=+1 and m= -1) 0.250 (for m=+2 and m= -2), and 0.375 (for m=+3 and m= -3). We therefore conclude: the computer-simulated results of the diffraction field, in the far-filed Fraunhofer limit, precisely agree, both qualitatively and quantitatively, with the results from the exact Fraunhofer theory, as represented by Eqs. (12) and (14). This agreement gives strong support to the validity and correctness of the present simulation methods.

Figure 8.

Computer-simulated Fresnel diffraction images from an 11-slitgrating, reproduced for clarity for a distance q₀=100 mm (Fraunhofer regime), with a =0.001mm, b=0.003mm, c=2mm and λ=500 nm. The intensity values have been multiplied by 10. The principal diffraction orders can be seen, as well as the secondary maxima and minima.