has to be devoted to the fact that Ez, Dz is defined in grid points while Hx, Bx, Hy, By are shifted by half grid points, thus two kinds of conductivity vectors should be used: σi and σi+1/2 depending on which field is calculated. TMz equations in two dimensions can be derived similarly. One-dimensional equations can be obtained by the choice κy = 1, σy = 0.
Eqs . (2.20a), (2.20b) and (2.20c) take the form:
Dye0ez =µxxh˜x → h˜x= (µxx)−1Dhy0ez, (2.21a)
−Dxe0ez =µyy˜hy → ˜hy =−(µyy)−1Dxh0ez, (2.21b)
Dhx0˜hy−Dhy0˜hx =εzzez. (2.21c)
where the symbol Dey0 represents the spatial difference operator in the ydirection for the ez field vector using the appropriate Yee mesh. The derivation of the spatial difference operators D,i is presented next. We define the following difference operator matrices:
De(N) =
−1 1 · · · 0 0 0 −1 · · · 0 0 ... ... . .. ... ...
0 0 · · · −1 1 0 0 · · · 0 −1
N×N
(2.22)
Dh(N) =
1 0· · · 0 0
−1 1· · · 0 0 ... ... . .. ... ...
0 0· · · 1 0 0 0· · · −1 1
N×N
(2.23)
The discrete difference operator matrices in theFDFD method are then:
Dex0 = INy ⊗De(Nx)
/(k0∆x), (2.24a)
Dye0 = (De(Ny)⊗INx)/(k0∆x), (2.24b)
Dhx0 = INy ⊗Dh(Nx)
/(k0∆x), (2.24c)
Dyh0 = Dh(Ny)⊗INx
/(k0∆x), (2.24d)
where⊗denotes Kronecker matrix multiplication. Matrices with size ofN×N are denoted by straight bold symbols. For the sake of clarity we denote the scaled up NxNy ×NxNy matrices differently by italic capital letters.
Substituting Eqs. (2.21a) and (2.21b) back to Eq. (2.21c) one obtains for aTEz mode:
Dhy0(µxx)−1Dye0 +Dhx0(µyy)−1Dex0 +εzz
ez = 0, (2.25)
and similarly for TMz:
Dey0(εxx)−1Dhy0+Dex0(εyy)−1Dxh0 +µzz˜hz = 0, (2.26)
2.2.1 Numerical phase velocity and anisotropy
When Maxwell’s equations are solved by numerical simulations, the purpose is to de-termine the electromagnetic fields at grid points. These values should be the same or nearly the same as those obtained by evaluation of the continuous or analytical solu-tion at the grid points. Evidently, the sampling density influences the accuracy of the numerical simulation. There are two critical dimensions which have to be considered when the grid resolution is determined: one is the minimum structural dimension, an-other is the wavelength. The grid resolution, ∆x needs to be small enough to resolve the minimum structural dimension, the exact value depends on the actual problem. It can be optimized (selected larger) by smoothing the structure boundaries [85]. In order to quantify the grid resolution compared to the wavelength the grid sampling resolution is defined: Nλ = λ0/∆x [81], Eq. (2.28b) on page 30. A natural limit for Nλ comes from the Nyquist-Shannon sampling theorem, which states that to avoid aliasing, the smallest spatial frequency should be less than 1/2∆x (see [81] page 151):
1 λ0
< 1
2∆x → Nλ >2 (2.27)
Not only the resolution but also the way how the discretization was performed influences the accuracy of the simulation.
In Sec.2.2 Eqs. (2.22) and (2.23) show that the discretization of the spatial derivative is carried out in the FDFDmethod by first order approximation, which causes numerical dispersion and anisotropy like in the FDTD method (see Sec. 2.6 in [81]). A detailed calculation for numerical dispersion and anisotropy is presented in Appendix F for the FDFD method based on the FDTD formulation (see Sec. 2.1.2). By introducing the nu-merical wave vector ˜k, and the numerical phase velocity ˜vp =ω/k, using Eq. (6.37) from˜ Appendix F we obtain the relative numerical phase velocity along the major grid axis:
˜
vap/vp = π Nλ
sin−1 π Nλ
−1
(2.28)
and, using Eq. (6.38), the same along the grid diagonal axis (φ= 45◦):
˜
vdp/vp = π
√2Nλ
sin−1 π
√2Nλ
−1
(2.29)
The numerical anisotropy is defined by
A= ˜vap/˜vpd. (2.30)
2.2.2 Boundary Condition
At the boundaries of the computational domain we are going to applyUPMLboundary condition [84]. The formulation is presented in Sec. 2.1.4. The same quantities s can be used for FDFD as for FDTD. In two dimensions each component of s forms a matrix Nx×Ny. The matrix Sx contains the UPML boundary for xdirection. In the simulation area all elements are ones, and in the UPML area at both ends in the x direction the absorption is increased terminating the boundary end (see Eqs. (2.17) and (2.18)). The matrixSy is similar but for theydirection. The matrixSx changes only inxdirection and is y-independent, while vice versa Sy changes only in y direction and is x-independent.
Fig.2.14 shows the geometric interpretation of the above description, where red indicates higher absorption.
Using the definition ofs from Eq. (2.16) in two dimensions:
Sxx =Sy./Sx, Syy =Sx./Sy, Szz =Sx·Sy, (2.31) where./and·are pointwise division and multiplication. Eqs. (2.21a) - (2.21c) indicate that evaluation of these matrices should be done at positions corresponding to the equations:
Sxx where ˜hx is defined, Syy where ˜hy is defined, and Szz where ez is defined. For the positions of the fields Fig.2.5provides a great help, by using it Eq. (2.31) can be rewritten as:
Sxx = [Sy(j)./Sx(i)]|i,j+1/2 =Sy(j + 1/2)./Sx(i), (2.32a) Syy = [Sx(i)./Sy(j)]|i+1/2,j =Sx(i + 1/2)./Sy(j), (2.32b)
Szz =Sx(i)·Sy(j). (2.32c)
Thus these matrices should be evaluated at half grid positions too. We will merge the UPML matrices with the permittivity and permeability (see Eqs. (2.15a) and (2.15b)).
We have to rearrange these Nx×Ny matrices (µxx,µyy, εzz) into diagonalNxNy×NxNy
matrices (µxx, µyy,εzz).
µxx → µxx =µxxSxx → µxx (2.33a)
µyy → µyy =µyySyy → µyy (2.33b)
εzz → εzz =εzzSzz → εzz (2.33c)
So Eqs. (2.25) and (2.26) have still the same form but the meaning of µxx, µyy, εzz has changed.
The way how anNx×Ny matrix is scaled up into NxNy×NxNy block matrix comes from the fact that normal matrix multiplication by block matrix should give the same result on column vector ez like pointwise multiplication by an Nx ×Ny matrix on field matrix Ez:
εzz·Ez ≡εzzez. (2.34)
In practice it can be executed by taking the two-dimensionalNx×Ny matrix, rearranging the components the same way as described in Fig. 2.15, and then inserting this vector into diagonal elements of a block matrix with size NxNy ×NxNy.
2.2.3 Implementation of source
Eqs. (2.25) and (2.26) show the formAx= 0 and have only trivial solution (fields are zeros). The reason is that the source has not yet been implemented. The source can be an electric or magnetic current in Maxwell’s equations, which would be situated on the right-hand side of the matrix equation after arranging the terms. As described in Sec.
2.1.3 there are plenty of advantages using TF-SF source condition. Luckily it is possible to implement it in the FDFD method and there is a sophisticated way to work out the corrections between TF and SF regions [83], these corrections will play the role of the source.
First we define a plane wave over the whole simulation area: Fsrc, which is scaled up to column vector fsrc the same way as depicted in Fig. 2.15. For distinguishing the SF and TF areas a mask matrix Q → Q is introduced, where Q is a diagonal block matrix with the same conversion as described in the previous section. This matrix contains zeros in the TF region and ones in the SF region. So the source field can be isolated to theSF region:fSF =Qfsrcand to the total field region:fTF = (I−Q)fsrc. We have to compensate the total field source which is in the SF region and the scattered field source which is in the TF region:
1. Total field source in the SF region: QAfTF (it should be subtracted).
2. Scattered field source in the TF region:(I−Q)AfSF (it should be added).
So that theFDFD equation by using TF-SF source condition reads:
Ax−QAfTF+ (I−Q)AfSF = 0 (2.35)
where we can identify the source term:
b =QAfTF−(I−Q)AfSF =QA(I−Q)fsrc−(I−Q)AQfsrc= (QA−AQ)fsrc, (2.36) Finally the solution is the field distribution obtained by matrix inversion: x=A−1b.
Fig. 2.16. Realization of TF-SF source condition by masking procedure. fs(x, y) is a plane wave, Q(x, y) is a mask matrix containing zeros in the TF region and ones in the SF region.<[b(x, y)]is the real part of the source term defined by Eq. (2.36) and hx(x, y) is the resulting field distribution.
Fig. 2.16 depicts how the TF-SF source condition works: using the incoming plane wave fsrc, and mask Q, the source b is easily constructed. This source field yields a wave
distribution which is nonzero only in the total field region (indicated by the mask) and the field is zero in the scattered field region since there is no scattering structure in this example.
2.2.4 Realization of oblique incidence
Simulation of oblique incidence is quite easy in the FDFD method, since only one frequency component is present. We only have to compensate the phase at PBC bound-aries. Assuming a plane wave with wave number k and angle of incidence θ , we have ky = ksinθ. Let us recall the matrices in Eqs. (2.22) and (2.23) and modify them by implementing oblique incidence:
De(N) =
−1 1 · · · 0 0 0 −1· · · 0 0 ... ... . .. ... ...
0 0 · · · −1 1 e−jkyLy 0 · · · 0 −1
N×N
(2.37)
Dh(N) =
1 0· · · 0 −e−jkyLy
−1 1· · · 0 0 ... ... . .. ... ... 0 0· · · 1 0 0 0· · · −1 1
N×N
(2.38)
The tilt of the incoming wave usingTF-SF source condition should agree with θ used in ky for this boundary condition.
2.2.5 Final matrix equations in the FDFD method
To summarize all features described above here we conclude those findings and present the final version of the FDFDequations. The general form of the FDFDmatrix equation in two dimensions is
Ax=b, (2.39)
whereA is the FDFDmatrix, andxdenotes the field column vectorez forTEmode, and
˜hz for TM mode. ForTE mode this equation stands for Dhy0(µxx)−1Dye0 +Dhx0(µyy)−1Dex0 +εzz
ez =b , (2.40)
and similarly for TM mode:
Dey0(εxx)−1Dhy0+Dex0(εyy)−1Dxh0 +µzz˜hz =b , (2.41) where b = (QA−AQ)fsrc is the TF-SF source (see in Sec. 2.2.3), the matrices ε and µ include theUPMLboundary condition (see in Sec.2.2.2), and matricesDxe/h0/y0 may contain the PBC boundary condition for normal or oblique angle of incidence (see in Sec.2.2.4).