Numerical Technique for Electromagnetic Field Computation Including High Contrast Composite Material

Periodic structure plays important role with its sensitivity for frequency in microwave communication and optical integrated circuit system [1][2]. Also, in optical fiber communication[3], such periodic structure is used for supporting optical wave to be guided in core of the fiber. In such application, a defect in the structure generally works as a cavity or a waveguide by making use of its high selectivity. However, by lacking periodicity in the structure with a defect, general mathematical approaches have often difficulties, due to inconvenience in description of the problem. In such cases, computational simulation technique for electromagnetic wave propagation and scattering is very important and effective. By the recent development of computers, it is possible to model a large scale periodic structure with defect.


Introduction
Periodic structure plays important role with its sensitivity for frequency in microwave communication and optical integrated circuit system [1] [2].
Also, in optical fiber communication [3], such periodic structure is used for supporting optical wave to be guided in core of the fiber.In such application, a defect in the structure generally works as a cavity or a waveguide by making use of its high selectivity.However, by lacking periodicity in the structure with a defect, general mathematical approaches have often difficulties, due to inconvenience in description of the problem.In such cases, computational simulation technique for electromagnetic wave propagation and scattering is very important and effective.By the recent development of computers, it is possible to model a large scale periodic structure with defect.
In simulation of electromagnetic phenomena, finite difference time domain (FDTD) method [4] [5] is most widely used.However, in numerical analysis of the wave behavior near boundary, where the dielectric constant is quite different on each side, a special care is required.Supposing in lossless dielectric medium, it is well known that wavelength of the electromagnetic wave changes due to the dielectric constant.Because of this compression of wavelength in high dielectric constant medium, grid size of space becomes rather coarse compared with material with lower dielectric constant.Therefore, accuracy of finite difference approximation deteriorate in most computation with uniform grid size.In numerical analysis of periodic structure such as photonic crystal(PC), the dielectric constant is generally quite high compared with its background medium.
A constrained interpolation profile (CIP) method [6][7] is payed much attention because of its accurate simulation result compared with conventional FDTD method.In this paper, on scattering problem by a dielectric cylinder with high contrast with its background air, performance of CIP method is compared with analytical approximated method using Hertz potential [8] and conventional FDTD simulation [9].As a measure of accuracy of CIP and FDTD simulation, an normalized cross correlation function is defined and compared with it, by setting analytical approximated result as a reference.Consequently, results of CIP method showed better correlation than that of FDTD method.As applications of CIP method, analysis of electromagnetic field propagation in Y-shaped branching waveguide and Mach-Zehnder interferometer in two dimensional photonic crystal structure were demonstrated.Both of analysis results showed reasonable behaviour.Especially for asymmetrical Mach-Zehnder interferometer, the measurement result by microwave model and the numerical result of CIP corresponded to each other.Complicated output characteristics of asymmetric Mach-Zhender interferometer was interpreted very well by refering to the electric field profile obtained by CIP method.

Formulation of electromagnetic field propagation by CIP method
Applications of CIP method to a problem of electromagnetic wave scattering is reported in Ref. [7].Outline of the CIP method and its application to solve two dimensional Maxwell's equation is briefly explained in following subsections, based on the reference.

CIP method and advection equation
We treat wave propagation with velocity u in isotropic and uniform space.The wave form is expressed by f (x, t) as a function of space x and time t.When the function f (x, t) satisfies following equation Eq.( 1) is called an advection equation.This equation holds for arbitrary wave function f with variable x ∓ ut, including electromagnetic wave propagation.This equation means that each point on the wave form moves to +Δx or −Δx from current point with velocity u after passage of Δt.This property is important to calculate electromagnetic field both in uniform space and on boundary.

Interpolation by cubic polynomial
In CIP method, wave function is approximated by following cubic polynomial, where a i , b i , c i and d i are unknown coefficients to be determined by making use of function value at each discrete point.Suppose that data at discrete points F i (x i ) and F i (x i−1 ) are known with its derivatives dF i (x i )/dx and dF i (x i−1 )/dx.Substituting x i and x i−1 into Eq.( 2) and its derivative function, we obtain where x i − x i−1 = Δx is used.According to property of advection Eq.( 1), the latest values of where, ξ = −uΔt.By solving Eq.(3) to Eq.( 6) with renewal of coefficients a i , b i , c i and d i , we have updated value of the function f and its derivative g at each discrete point from Eq.( 7) and Eq.( 8).

Conversion of Maxwell's equation into advection equation
Maxwell's curl equations for electric field vector E and magnetic field vector H in lossless, isotropic and nonconductive material are given as follows, where ε is permittivity and μ is permeability.
Assuming two dimensional uniform space along with z axis (i.e.∂/∂z = 0), Maxwell equations are decomposed into two sets of polarization.We treat E-wave which includes (H x , H y , E z ) as the component.For E-wave, Maxwell's equations are reduced to followings: From Eqs.(11) to (13), we obtain where Split step procedure is used for Eq.( 14), then we have where superscripts n and n + 1 denote time step, and W * means middle value of EM components.These equations are converged into the advection equation.For plane wave with component (E z , H y ), Eq.( 14) can be written down as By using wave impedance Z = μ/ε and velocity of light c = 1/ √ εμ, above equations can be rewritten as where, Hy = ZH y .By addition and subtraction of these equations, we obtain for propagation along with x axis.In the same manner, we have advection equations for y-direction.Solving these equation with adequate Δx, Δy, and Δt, it is possible to renew electromagnetic field at each discrete point.Similarly, three dimensional analysis can be performed.

Scattering of plane wave by a dielectric rod with high contrast permittivity to background
In this section, scattered wave by a dielectric cylinder, which is obtained by CIP and FDTD method, is compared with analytical approximated solution [8].Model of a dielectric cylinder and the coordinate system is illustrated in Fig. 1.Two dimensional analysis region is 80.0 × 80.0 mm 2 and diameter of the cylinder is 20.0 mm.Plane wave of E-polarization with frequency f = 10.0 GHz is given as continuous incident wave.
In analytical approximated approach, wave function is expanded into summation of Bessel and Hankel functions as basis function.Electromagnetic field is determined so that boundary condition on surface of a dielectric cylinder is satisfied.In this demonstration, convergence of electromagnetic field is confirmed by increasing number of truncation of basis function.
Simulation for the structure was demonstrated by setting Δx = Δz = 0.5mm and Δt = 0.17ps, both for CIP and FDTD method.As examples of total electric field intensity near a cylinder, the profiles are shown in Fig. 2 and Fig. 3, where relative dielectric constant ε r of the cylinder is chosen to be 25.0, and 36.0,respectively [9].The background medium is air and its relative dielectric constant ε r = 1.0 is supposed.These high values of ε is needed for constructing photonic crystal structure, which requires high contrast with respect to the background medium.
From these field profiles, we can evaluate normalized cross correlation function defined as, where (x i , y j ) is discrete point in analyzing region, E z,Ana (x i , y j ) is E z field by analytical approximated solution, E z,{CIP|FDTD} (x i , y j ) is that by CIP and FDTD method, respectively.Using Eq.( 25), we obtain Table .1.From Table .1,values of normalized cross correlation by CIP method is always better than that by FDTD method for typical four kinds of relative dielectric constant ε r of the cylinder.

Y-shaped branch waveguide in two-dimensional photonic crystal with triangular lattice
In Fig. 4, two dimensional scale model of Y-shaped branch waveguide is illustrated.From port 1, continuous electric field with Gaussian profile along with vertical axis are given.The electromagnetic field is composed of unique electric field which is oriented to vertical direction to the paper surface and magnetic fields which are parallel to L 0 and W 0 axes.Therefore, the incident wave is E-polalization.Spot size 2w 0 of the Gaussian profile corresponds to waveguide full width √ 3P, where P = 26.5mm is lattice period.Dielectric rods have same parameters with previous numerical examples.Example of the total electric field profile with frequency f = 4.0GHz by CIP method is shown in Fig. 5, where parameters were set as follows; space discretization Δ = 0.5mm, time step Δt = 2.5 × 10 −13 sec, total cell numbers 648 × 659, cells on diameter of the cylinder 10, respectively.From the Fig. 6, it is found that incident wave is equally divided by the branch circuit.In Fig. 7, maximum electric field profile at output after filtering is shown.Two pairs of rods were inserted simultaneously in each output waveguide for filtering different frequencies.

Incident Gaussian continuous wave
Length of the cavity is 2.5P and another one is 3.0P, respectively, where P is lattice period.As expected, frequencies with half-wavelength-along-waveguide which corresponding to each cavity length are filtered and obtained in each port.

Mach-Zehnder interferometer in two dimensional photonic crystal with triangular lattice
In Fig. 8, Mach-Zehnder(MZ) type interferometers are depicted.The structure is situated in two dimensional photonic crystal with triangular lattice.To investigate interference at combining point, asymmetrical structure with different arm lengths are compared with symmetrical one.At first, electric field profile in symmetrical MZ structure is shown in Fig. 9 The output characteristics of the asymmetric structure was investigated by experiment using a model in same microwave frequency range.The measurement results are shown in Fig. 10.

3.974GHz
Output electric field profile after filtering by cavities in arbitrary unit  In Fig. 10, it is found that the asymmetric MZ interferometer (MZI) shows maximum and minimum transmission at frequencies 3.816GHz and 4.196GHz, which are indicated by arrows and dots in the figure.Measurement of symmetrical MZI showed almost flat and relatively high transmission characteristics over same frequency range, comparing with result of asymmetrical one.Therefore, simulations by CIP method were demonstrated for these two frequencies and additional one frequency 4.006GHz as middle level of output.
The electric field profiles by CIP method along with MZI are shown for these three frequencies in Fig. 11.The electric field profile at output port is also indicated in Fig. 12. From these figures, it is found that electric field propagated along with two arms comes to combining point with relatively small phase difference in Fig. 11(a) and (b), while electric field become extinct for the two fields comes to combining point with out of phase in Fig. 11(c) .In Fig. 12(c), the maximum electric field is quite small compared with output with small phase difference in MZI as shown in Fig. 12(a) and (b).From Fig. 10 to 12, complicated output characteristics were qlearly interpreted by numerical results by CIP method.This is because CIP method provides more precise results of electromagnetic wave scattering compared with FDTD method, as we saw in subsection 3.1.From these results, it was shown that superiority and significance of the CIP method for designing photonic crystal structure which is composed of periodic structure with high contrast of material constant.

Conclusion
By using CIP method, numerical analysis of scattered electromagnetic field by a dielectric cylinder was demonstrated.Referring to result of analytical approximated approach, result of CIP method showed good accuracy in comparison with result by conventional FDTD method.As examples of designing of photonic crystal structure with high-contrast material profile, Y-shaped branch waveguide and Mach-Zehnder interferometer in photonic crystal structure with triangular lattice were numerically demonstrated.The CIP method showed reasonable results of branch circuit and the filtering characteristics by using cavity.Also for Mach-Zehnder interferometer, numerical results of electric field profile by CIP method implemented experimental result of the structure for typical frequencies.
For example, designing of filtering device for microwave communication or guiding device for optical fiber communication system, the CIP method is expected to show superior performance in accuracy compared with conventional FDTD method.Application of CIP method to design electromagnetic or optical signal processing devices with some defects in periodic structure can be designed by CIP method.

Figure 1 .
Figure 1.Illustration of the analysis region with coordinate system.A black circle indicates a dielectric cylinder with dielectric constant ε r .

Figure 2 .
Figure 2. Field profile of E z by (a) analytical approximated solution, (b) CIP method and (c) FDTD method, respectively, for a case that ε r of the cylinder is 25.0.

Figure 3 .
Figure 3. Field profile of E z by (a) analytical approximated solution, (b) CIP method and (c) FDTD method, respectively, for ε r = 36.0.

Table 1 .
Values of normalized cross correlation.

Figure 5 .
Figure 5. Electric field profile of E-polarised wave in Y-shaped branch waveguide.

Figure 6 .
Figure 6.Maximum electric field profile at output port.Two identical peaks locates on center of each output waveguide.0.10

Figure 7 .Figure 8 .Figure 9 .
Figure 7. Electric field profile at output after filtering.Fields with maximum amplitude in each output port were depicted.

Figure 12 .
Figure 12.Maximum electric field profile at output port for typical three frequencies.