Ersoy O.K. Diffraction, Fourier optics, and imaging (Wiley, 2006)(ISBN 0471238163)(427s) PEo

368	NUMERICAL METHODS FOR RIGOROUS DIFFRACTION THEORY

20.5FINITE DIFFERENCE TIME DOMAIN METHOD

The FDTD method is an effective numerical method increasingly used in applications. In this method, the Maxwell’s equations are represented in terms of central-difference equations. The resulting equations are solved in a leapfrog manner. In other words, the electric field is solved at a given instant in time, then the magnetic field is solved at the next instant in time, and the process is repeated iteratively many times.

Consider Eqs. (3.3-14) and (3.3-15) of Chapter 3 for a source-free region (no electric or magnetic current sources) rewritten below for convenience:

	@D				@E
r H ¼			¼ e					ð3:3-14Þ
	@t				@t
	@B			¼ m			@H	ð3:3-15Þ
r E ¼
		@t					@t

Eq. (3.3-14) shows that the time derivative of the E field is proportional to the Curl of the H field. This means that the new value of the E field can be obtained from the previous value of the E field and the difference in the old value of the H field on either side of the E field point in space.

Equations (3.3-14) and (3.3-15) can be written in terms of the vector components as follows:

@Hx		1			@Ey		@Ez	r0Hx
	@t	¼	m		@z		@y		ð20:5-1Þ
@Hy		1			@Ez		@Ex	r0Hy
	@t	¼	m		@x		@z		ð20:5-2Þ
@Hz		1			@Ex		@Ey	r0Hz
	@t	¼	m		@y		@x		ð20:5-3Þ
	@Ex	1			@Hz		@Hy	sEx
	@t	¼	e		@y		@z		ð20:5-4Þ
	@Ey	1			@Hx		@Hz	sEy
	@t	¼	e		@z		@x		ð20:5-5Þ
	@Ez	1			@Hy		@Hx	sEz
	@t	¼	e		@x		@y		ð20:5-6Þ

These equations are next represented in terms of finite differences, using the Yee algorithm discussed below.

20.5.1Yee’s Algorithm

Yee’s algorithm solves Maxwell’s curl equations by using a set of finite-difference equations [Yee]. Using finite differences, each electric field component in 3-D is

FINITE DIFFERENCE TIME DOMAIN METHOD

369

Figure 20.3. Electric and magnetic field vectors in a Yee cell [Yee].

surrounded by four circulating magnetic field components, and every magnetic field component is surrounded by four circulating electric field components as shown in Figure 20.3. The electric and magnetic field components are also centered in time in a leapfrog arrangement. This means all the electric components are computed and stored for a particular time using previously stored magnetic components. Then, all the magnetic components are determined using the previous electric field components [Kuhl, Ersoy].

With respect to the Yee cell, discretization is done as follows:

½i; j; k& ¼ ði x; j y; k zÞ where x; y; and z are the space increments in the x, y, and z directions, respectively, and i, j, k are integers.

tn ¼ n t

uði x; j y; k z; n tÞ ¼ uni;j;k

Yee’s centered finite-difference expression for the first partial space derivative of u in the x-direction, evaluated at time tn ¼ n t is given by [Yee]

@u		uin	1=2;j;k uin		1=2; j;k		2
	ði x; j y; k z; n tÞ ¼	þ	x			þ Ohð xÞ	i	ð20:5-7Þ
@x

It is observed that the data a distance x=2 away from the point in question is used, similarly to the Crank–Nicholson method. The first partial derivative of u with respect to time for a particular space point is also given by

@u		uin;þj;k1=2 uin; j;k1=2		2
	ði x; j y; k z; n tÞ ¼		þ Ohð tÞ	i	ð20:5-8Þ
@t		t

2ei; j;k

370	NUMERICAL METHODS FOR RIGOROUS DIFFRACTION THEORY

Substituting in the space and time derivatives for the point ði; j; kÞ at time n in Eq. (20.5-4) yields

n 1 2

þ

n 1 2

0

Hzjin; jþ1=2;k Hzjin; j 1=2;k

1

Exji;þj;k=

Exji; j;k=

1

y

ð

20:5-9

Þ

t

¼ ei; j;k

Hyjin; j;kþ1=2 Hyjin; j;k 1=2

@

A

B

j

n

C

B

z

C

B

si; j;kEx

i; j;k C

All the terms on the right side of this equation are evaluated at time step n. All the magnetic field components needed at time n are known. Only the values of Ex up to time n 1=2 are stored. For a region without loss, this term is zero. If si;j;k is nonzero, then Ex can be estimated as [Yee]

n		Exjin;þj;k1=2 þ Exjin; j;k1=2		ð20:5-10Þ
Exji; j;k	¼		2

which is the average of the known value of Ex at time n 1=2 and the unknown value at time n þ 1=2. Using Eq. (20.5-10) in Eq. (20.5-9) yields

Hz

in; jþ1=2;k Hzjin; j 1=2;k

0

j

nþ1=2

n 1=2

t

y

2

Ex i; j;k

Ex i; j;k

B

n 1 2

n 1 2

j

þ j

¼ ei; j;k

@

Ex

i;þj;k=

Ex

i; j;k=

B

þ

B

si; j;k

j

j

n 1 2

=

results in

Collecting like terms and solving for Ex i;þj;k

Hy	in; j;kþ1=2 Hyjin; j; k 1=2
j	z	1
			C
			C
			C
			A
		ð20:5-11Þ

				k t
	0	1		si; j;	1Ex
Ex nþ1=2				2ei; j;k		n 1=2
ji; j;k	B		þ	si; j;k t	C	ji; j;k
	@				A
	¼ B1				C

t

Hz

in; jþ1=2;k Hzjin; j 1=2;k

0

10

j

1

ei; j;k

y

B

CB

C

þ

si; j;k t

z

@

AB

C

@

A

þ

1

B

Hy

in; j;k

1=2

Hy

in; j;k

1=2 C

B

2ei; j;k

C

j

þ

j

ð20:5-12Þ

All the other unknown field components can be similarly evaluated as

si; j;k t

0

Hx

in; j;kþ1=2 Hxjin; j;k 1=2

1

1

t

j

z

2ei; j;k

n 1=2

n 1=2

ei; j;k

Ey i;þj;k

B

C

i;j;k

B

CB

j þ

j

C

þ 2ei; j;k

þ

2ei; j;k

0

si; j;k

t1Ey

0

si; j;k

t1

n

x

n

C

@

A

@

AB

@

A

j

¼

1

j

þ

1

B Hz i 1=2; j;k

Hz i 1=2; j;k C

ð20:5-13Þ

2mi;j;k

COMPUTER EXPERIMENTS

				k t
		0	1	si; j;		1Ez
Ez nþ1=2				2ei; j;k			n 1=2
ji; j;k	¼	@		si; j;k t	A		ji; j;k
		B1			C
		B þ		2ei; j;k	C

0		t
		ei; j;k
@		si; j;k t
þ B1	þ
B		2ei; j;k

						371
	Hy		in	1=2; j;k Hyjin		1=2; j;k
10		j	þ	x			1
B		Hxjin; jþ1=2;k Hxjin; j 1=2;kC
CB				y				C
CB								C
A@								A

ð20:5-14Þ

0

t

t

10

Ey

nþ1=2

Ey

nþ1=2

1

ri; j;k

i; j;k

1=2

ji; j;k

1=2

01

1

0

j

þ

n 1

2mi; j;k

n

mi; j;k

z

Hxji;þj;k

¼

Hxji; j;k

þ

B

C

B1

ri0; j;k t

C

B1

ri0; j;k t

Ez

nþ1=2

Ez nþ1=2

CB

ji; j

1=2;k

ji; j

1=2;k

C

B

þ

2mi; j;k

C

B

þ 2mi; j;k

CB

þ

y

C

B

C

B

CB

C

@

A

@

A@

A

ð20:5-15Þ

Hyjni;þj;1k

Hzjni;þj;1k

			r0		k t
	01			i; j;		1
				2mi; j;k			n
¼	B					CHyji; j;k
		1	ri0; j;k t
			þ 2
	B				mi; j;k	C
	B					C
	@					A
			r0		k t
	01			i; j;		1
				2mi; j;k			n
¼	B					CHzji; j;k
		1	ri0; j;k t
	B		þ			C
	B					C
	@					A

0

mi; j;k

10

j

nþ1=2

x

nþ1=2

1

þ

t

Ez

i 1=2; j;k

Ezji

1=2; j;k

B

CB

C

z

þ

@

ri0; j;k t

B

Ex nþ1=2

Ex

nþ1=2

C

A@

A

B

þ

2mi; j;k

CB

þ

1=2

C

B1

CB

ji; j;k

1=2 ji; j;k

C

ð20:5-16Þ

10

nþ1=2

nþ1=2

1

Ex

i; j

1=2;k

Exji; j 1=2;k

0

t

j

þ

y

mi; j;k

B

CB

x

C

þ

@

ri0; j;k t

B

E nþ1=2

E

nþ1=2

C

þ 2mi; j;k

A@

j

þ

j

A

B

CB

C

B1

CB

y

i 1=2; j;k

y

i

1=2; j;k

C

ð20:5-17Þ

20.6COMPUTER EXPERIMENTS

The software package XFDTD by Remcom was used in the computer experiments [Kuhl, Ersoy]. In this environment, the region of interest is a cubical mesh, where each mesh edge can be given different material properties in order to simulate a specified geometry. For each cell, the material may be a perfect conductor or free space, or may be defined in other ways. The sampling in space has sub-wavelength format, typically in the range of 1/10–1/30 of a wavelength. The region of interest is excited by either a plane wave or multiple voltage sources. The excitation may be pulsed or sinusoidal. The duration of the simulation is set by specifying the number of desired time steps.

372	NUMERICAL METHODS FOR RIGOROUS DIFFRACTION THEORY

Figure 20.4. 1-D FZP with m ¼ 3.

When the modeled region extends to infinity, absorbing boundary conditions (ABCs) are used at the boundary of the grid. This allows all outgoing waves to leave the region with negligible reflection. The region of interest can also be enclosed by a perfect electrical conductor.

Once the fields are calculated for the specified number of time steps, near zone transient and steady state fields can be visualized as color intensity images, or a field component at a specific point can be plotted versus time. When the steady-state output is desired, observing a specific point over time helps to indicate whether a steady-state has been reached.

In the computer experiments performed, a cell size of l=20 was used [Kuhl, Ersoy]. The excitation was a y-polarized sinusoidal plane wave propagating in the z-direction. All diffracting structures were made of perfect electrical conductors. All edges of the diffracting structures were parallel to the x-axis to avoid canceling the y-polarized electric field.

A 1-D FZP is the same as the FZP discussed in Section 15.10 with the x-variable dropped. The mode m is defined as the number of even or odd zones which are opaque.

An example with three opaque zones (m ¼ 3) is shown in Figure 20.4.

Using XFDTD, a 1-D FZP with a thickness of 0:1l and focal length 3l was simulated. Its output was analyzed for the first three modes. The intensity along the z-axis passing through the center of the FZP is plotted as a function of distance from the FZP in Figure 20.5. The plot shows that the intensity peaks near 3l behind the plate, and gets higher and narrower as the mode increases. The peak also gets closer to the desired focal length of 3l for higher modes.

The plot of the intensity along the y-axis at the focal line is shown in Figure 20.6. The plot shows that the spot size decreases with increasing mode and side lobe intensity is reduced for higher modes.

These experiments show that the FDTD method provides considerable freedom in generating the desired geometry and specifying material parameters. It is especially useful when the scalar diffraction theory cannot be used with sufficient accuracy due to reasons such as difficult geometries and/or diffracting apertures considerably smaller than the wavelength.

COMPUTER EXPERIMENTS	373

Figure 20.5. FDTD results for intensity along the z-axis.

Figure 20.6. FDTD results for intensity along the y-axis on the focal line.

374	NUMERICAL METHODS FOR RIGOROUS DIFFRACTION THEORY

20.7FOURIER MODAL METHODS

Fourier modal methods are among the several most versatile rigorous methods to analyze surface relief DOEs. The initial work in this area is known as rigorous coupled-wave analysis (RCWA) [Moharam, 82]. It was mostly used initially with diffraction gratings, but they were later generalized to aperiodic structures as well [Lalanne–Silberstein]. In the introductory discussion in this section, we will assume a grating structure.

The FMM is a frequency domain method involving computations of the grating modes as eigenvectors, Fourier expansions of the permittivity and of the EM fields inside the grating region. The Fourier expansion of the field inside the grating generates a system of differential equations. After finding the eigenvalues and eigenvectors of this system, the boundary conditions at the grating interfaces are matched to compute the diffraction efficiencies [Lalagne and Morris 1996].

The geometry to be used to explain the method in 2-D is shown in Figure 20.7. The periodic grating with period is along the x-axis with a permittivity function eðxÞ. The TM mode is considered such that the magnetic field is polarized in the y- direction, and there are no variations along the y-direction. The incident plane wave makes an angle y with the z-direction. In addition, the following Fourier series

definitions are made:

	eðxÞ		1				e jKkx			20:7-1
			¼ k¼X1 ek						ð		Þ
	eo
	o		1
e			¼ k¼X1 ake jKkx						ð20:7-2Þ
eðxÞ

Figure 20.7. The geometry used in the FMM method [Lalanne and Morris 1996].

FOURIER MODAL METHODS

375

where K ¼ 2p= . The field components can be written as [Lalanne and Morris 1996]

X

		Ex ¼			m			SmðzÞe jðKmþbÞx		ð20:7-3Þ
		Ez ¼		X						ð20:7-4Þ
					m			fmðzÞe jðKmþbÞx
		Hy ¼		X						ð20:7-5Þ
					m			UmðzÞe jðKmþbÞx
where b ¼ k sin y ¼	2p	y.
	l
Maxwell’s curl equations in this case are given by
			dEz	þ			dEx		¼ jwmoHy	ð20:7-6Þ
			dx					dx
							dHy		¼ jweEx	ð20:7-7Þ
								dz
				1 dHy					¼ jwEz	ð20:7-8Þ
					e			dx

Prime and double prime quantities will be used to denote the first and second partial derivatives with respect to the z-variable. Using Eqs. (20.7-3)–(20.7-5) in Eqs. (20.7-6)–(20.7-8) results in

jðKm þ bÞfm þ Sm0 ¼ jkoUm							ð20:7-9Þ
		Um0 ¼ jko Xp				em pSp	ð20:7-10Þ
1				Xp ðpK þ bÞam pUp
		fm ¼					ð20:7-11Þ
			ko
Using Eq. (20.7-11) in Eq. (20.7-9) yields
Sm0 ¼ jkoUm þ	j	ðKm þ bÞ Xp			ðpK þ bÞam pUp		ð20:7-12Þ
	ko

In practice, the sums above are truncated with jpj M. Then, Eqs. (20.7-10) and (20.7-12) make up an eigenvalue problem of size 2ð2M þ 1Þ. However, it turns out to be more advantageous computationally to compute the second partial derivative of Um with Eq. (20.7-10) and obtain [Caylord, Moharam, 85], [Peng]

X				X	1
Um00 ¼ ko2( p	1			)( r		ðrK þ bÞap rUr) ð20:7-13Þ
	em p Up	ko	ðpK þ bÞ		ko

376

	72
	70
(%)	68
intensity	66
Transmitted	64
	62

60

580

NUMERICAL METHODS FOR RIGOROUS DIFFRACTION THEORY

69.822

69.229

50

100

150

200

250

300

350

400

Number of retained orders

Figure 20.8. Diffraction efficiency of the transmitted zeroth order beam from a metallic grating with TM polarized light [Courtesy of Lalanne and Morris, 1996].

Equation (20.7-13) can be written in matrix form as
1		U00 ¼ EðKxE 1Kx IÞU	ð20:7-14Þ
	k2
	o

where all the bold upper case letters represent matrices; I is the identity matrix, E is the permittivity harmonic coefficients matrix, and Kx is a diagonal matrix whose ith diagonal element equals ðiK þ bÞ=ko.

A more efficient form of Eq. (20.7-14) computationally is given by [Lalanne

1996]
1		U00 ¼ EðEKxE 1Kx IÞU	ð20:7-15Þ
	k2
	o

As an example, Figure 20.8 shows the diffraction efficiency of the zeroth order beam of a metallic grating with TM polarized light with respect to the number of Fourier coefficients kept [Lalagne, Morris, 1996]. The solid curve and the curve with circles were obtained with Eqs. (20.7-14) and (20.7-15), respectively.

The coverage above is introductory at best. There are currently a number of refinements of the Fourier modal method for different applications, and further refinements are expected

Appendix A

The Impulse Function

The impulse (delta or Dirac delta) function dðtÞ can be regarded as the idealization of a very narrow pulse with unit area. Consider the finite pulse shown in Figure A.1. It is defined by

	8	1		a		< t <	a
x t								ð	A:1-1	Þ
		a		2			2
	:	0	otherwise
ð Þ ¼ <

The area under the pulse is 1 and remains as 1 for all values of a. The impulse function can be defined as

dðtÞ ¼ lim xðtÞ

ðA:1-2Þ

a!0

dðtÞ satisfies

1ð

dðtÞdt ¼ 1

ðA:1-3Þ

1

dðtÞ ¼ 0 t 6¼ 0

In mathematics, dðtÞ is considered to be not an ordinary function, but a generalized function or a distribution, as discussed in the following sections. dðtÞ has the sampling (sifting) property given by

t2	t2
tð1 dðt tÞhðtÞdt ¼ tð1 dðt tÞhðtÞdt ¼ hðtÞ		ðA:1-4Þ

where h(t), a given function, is continuous at t ¼ t, and t1 < t < t2.

It is observed that if h(t) is the impulse response of a LTI system, application of dðtÞ to the system leads to the identification of h(t).

Diffraction, Fourier Optics and Imaging, by Okan K. Ersoy

Copyright # 2007 John Wiley & Sons, Inc.

377